; ; This script plots MOC ; ; Li Qing, 20121121 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/popRemap.ncl" begin ; define parameters type = "eps" ; plot type ; type = "x11" type@wkOrientation = "portrait" ; define file names diri = "/home/liqing/data/CTLB18/" fili = "ctl.B1850.pop.avgh.0901.1000.nc" f = addfile(diri+fili,"r") comp = (/"eulerian","bolus","submeso","total"/) plot1 = new(2, graphic) plot2 = new(2, graphic) plot3 = new(2, graphic) plot12 = new(2, graphic) plot22 = new(2, graphic) plot32 = new(2, graphic) tmax = (/36,12,8,36/) dt = (/4,2,1,4/) zl = (/4.,0,0,4./) do i = 0,3 filo = "moc3_global_"+comp(i) filo2 = "moc3_atlantic_"+comp(i) filo3 = "moc3_pacind_"+comp(i) ; read variables if (i .eq. 3) then mocg = f->MOC(0,0,0,:,:) + f->MOC(0,0,1,:,:) + f->MOC(0,0,2,:,:) moca = f->MOC(0,1,0,:,:) + f->MOC(0,1,1,:,:) + f->MOC(0,1,2,:,:) else mocg = (/f->MOC(0,0,i,:,:)/) moca = (/f->MOC(0,1,i,:,:)/) end if mocpi = mocg - moca lat = f->lat_aux_grid lev = f->moc_z / 100 nlat = dimsizes(lat) nlev = dimsizes(lev) mocg = where(abs(mocg) .gt. 1.e-6, mocg, mocg@_FillValue) moca = where(abs(moca) .gt. 1.e-6, moca, moca@_FillValue) mocpi = where(abs(mocpi) .gt. 1.e-6, mocpi, mocpi@_FillValue) vmask = new((/nlev,nlat/), float) do j = 0,nlev-1 vmask(j,:) = lat end do moca = where(vmask .gt. -35, moca, moca@_FillValue) mocpi = where(vmask .gt. -35, mocpi, mocpi@_FillValue) mocg!0 = "lev" mocg&lev = lev mocg!1 = "lat" mocg&lat = lat moca!0 = "lev" moca&lev = lev moca!1 = "lat" moca&lat = lat mocpi!0 = "lev" mocpi&lev = lev mocpi!1 = "lat" mocpi&lat = lat ; get zonal mean temperature gridin = "gx1v6" gridout = "1x1d" method = "aave" ddate = "130115" tmp0 = f->TEMP(0,:,:,:) tmp = tmp0 amask = where(ismissing(tmp),0.,1.) temp0 = PopLatLon(tmp0, gridin, gridout, method, "da", ddate) mask_rg = PopLatLon(amask, gridin, gridout, method, "da", ddate) mask_rg = where(mask_rg .eq. 0, mask_rg@_FillValue, mask_rg) gtemp = temp0 gtemp = (/temp0/mask_rg/) zmgtp = dim_avg_n(gtemp,2) ndim = dimsizes(f->TLAT) lev2 = f->z_t / 100. nlev2 = dimsizes(lev2) ; print(ndim) rgmsk = new((/nlev2,ndim(0),ndim(1)/),integer) do kk = 0,nlev2-1 rgmsk(kk,:,:) = f->REGION_MASK end do tmp = where(rgmsk .eq. 6 .or. rgmsk .eq. 8 .or. rgmsk .eq. 9 \ .or. rgmsk .eq. 10,tmp0,tmp0@_FillValue) amask = where(ismissing(tmp),0.,1.) temp0 = PopLatLon(tmp, gridin, gridout, method, "da", ddate) mask_rg = PopLatLon(amask, gridin, gridout, method, "da", ddate) mask_rg = where(mask_rg .eq. 0, mask_rg@_FillValue, mask_rg) atemp = temp0 atemp = (/temp0/mask_rg/) zmatp = dim_avg_n(atemp,2) tmp = where(rgmsk .eq. 2 .or. rgmsk .eq. 3,tmp0,tmp0@_FillValue) amask = where(ismissing(tmp),0.,1.) temp0 = PopLatLon(tmp, gridin, gridout, method, "da", ddate) mask_rg = PopLatLon(amask, gridin, gridout, method, "da", ddate) mask_rg = where(mask_rg .eq. 0, mask_rg@_FillValue, mask_rg) pitemp = temp0 pitemp = (/temp0/mask_rg/) zmpitp = dim_avg_n(pitemp,2) zmgtp!0 = "lev" zmgtp&lev = lev2 zmgtp!1 = "lat" zmgtp&lat = gtemp&lat zmatp!0 = "lev" zmatp&lev = lev2 zmatp!1 = "lat" zmatp&lat = atemp&lat zmpitp!0 = "lev" zmpitp&lev = lev2 zmpitp!1 = "lat" zmpitp&lat = pitemp&lat ; open workstation fontht = 0.018 wks = gsn_open_wks(type, filo) wks2 = gsn_open_wks(type, filo2) wks3 = gsn_open_wks(type, filo3) gsn_define_colormap(wks, "BlueWhiteOrangeRed") ;BlueYellowRed") ;WhBlGrYeRe") ;BlGrYeOrReVi200") gsn_define_colormap(wks2, "BlueWhiteOrangeRed") ;BlueYellowRed") ;WhBlGrYeRe") ;BlGrYeOrReVi200") gsn_define_colormap(wks3, "BlueWhiteOrangeRed") ;BlueYellowRed") ;WhBlGrYeRe") ;BlGrYeOrReVi200") res = True res@gsnDraw = False res@gsnFrame = False res@trXMaxF = 80 res@trXMinF = -80 res@tiYAxisString = "Depth (m)" res@tiXAxisString = "Latitude" res@tiXAxisFontHeightF = fontht res@tiYAxisFontHeightF = fontht res@tmXBLabelFontHeightF = fontht res@tmYLLabelFontHeightF = fontht res@gsnRightString = "" res@tmBorderThicknessF = 2.5 res@tmXBMajorThicknessF = 2.5 res@tmXBMinorThicknessF = 1.5 res@tmXTMajorThicknessF = 2.5 res@tmXTMinorThicknessF = 1.5 res@tmYLMajorThicknessF = 2.5 res@tmYLMinorThicknessF = 1.5 res@tmYRMajorThicknessF = 2.5 res@tmYRMinorThicknessF = 1.5 ; line labels res@cnInfoLabelOn = False res@cnLineLabelsOn = False res@cnLineLabelPlacementMode = "Computed" res@cnLineLabelInterval = 1 res@cnLineColor = "gray" res@cnLineLabelBackgroundColor = 0 res@cnLineLabelFont = "helvetica" res@cnLineLabelFontHeightF = 0.012 res@cnLineThicknessF = 2.5 res@gsnContourZeroLineThicknessF = zl(i) ; doubles thickness of zero contour res@gsnContourNegLineDashPattern = 5 ; sets negative contours to dash pattern 1 res@gsnYAxisIrregular2Linear = True res@gsnXAxisIrregular2Linear = True res@tmXBLabelsOn = True res@trYReverse = True ; reverse the Y-axis res@cnFillOn = True res@gsnSpreadColors = True res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -1. * tmax(i) res@cnMaxLevelValF = tmax(i) res@cnLevelSpacingF = dt(i) res@lbLabelBarOn = True res@lbOrientation = "vertical" ; vertical label bar res@lbLeftMarginF = 0.0 res@lbLabelAutoStride = True res@lbLabelFontHeightF = 0.015 res@pmLabelBarParallelPosF = 0.8 res@pmLabelBarHeightF = 0.4 ; contour1 res1 = res ; res1@gsnLeftString = "Global_"+comp(i) res1@gsnLeftString = "" res1@tmXBLabelsOn = False res1@tiXAxisOn = False res1@tiYAxisOn = False res1@lbLabelBarOn = False res1@vpWidthF = 0.6 res1@vpHeightF = 0.15 res1@trYMaxF = 500 res1@trYMinF = 5 res2 = res res2@gsnLeftString = "" res2@vpWidthF = 0.6 res2@vpHeightF = 0.25 res2@tiYAxisOffsetYF = 0.1 res2@trYMaxF = 5000 res2@trYMinF = 500 res3 = res1 res3@gsnLeftString = "" res3@cnFillOn = False res3@cnMinLevelValF = -2 res3@cnMaxLevelValF = 36 res3@cnLevelSpacingF = 2 res3@cnLineColor = "black" ; res3@cnLineThicknessF = 1.5 res3@cnLineLabelsOn = True res4 = res3 res4@tmXBLabelsOn = False res4@vpWidthF = 0.6 res4@vpHeightF = 0.25 res4@trYMaxF = 5000 res4@trYMinF = 500 ; plot mocg = smth9(mocg,0.50,0.25,False) ; smooth mocg = smth9(mocg,0.50,0.25,False) ; smooth moca = smth9(moca,0.50,0.25,False) moca = smth9(moca,0.50,0.25,False) mocpi = smth9(mocpi,0.50,0.25,False) mocpi = smth9(mocpi,0.50,0.25,False) ; reference line pres = True pres@gsLineColor = "black" pres@gsLineThicknessF = 2.5 pres@gsLineDashPattern = 0 pres@tfPolyDrawOrder = "Postdraw" refx1 = new(34,float) refy1 = new(34,float) refx2 = new(27,float) refy2 = new(27,float) refx1 = 0. refy1 = lev(0:33) refx2 = 0. refy2 = lev(33:59) ; attach res attres11 = True attres11@gsnAttachPlotsXAxis = True attres12 = False plot1(0) = gsn_csm_contour(wks,mocg(lev|0:33,lat|:),res1) plot1(1) = gsn_csm_contour(wks,mocg(lev|33:59,lat|:),res2) ref1 = gsn_add_polyline(wks,plot1(0),refx1,refy1,pres) ref2 = gsn_add_polyline(wks,plot1(1),refx2,refy2,pres) plot12(0) = gsn_csm_contour(wks,zmgtp(lev|0:33,lat|:),res3) plot12(1) = gsn_csm_contour(wks,zmgtp(lev|33:59,lat|:),res4) newplot = gsn_attach_plots(plot1(0),plot1(1),attres11,attres12) newplot12 = gsn_attach_plots(plot12(0),plot12(1),attres11,attres12) overlay(plot1(0),plot12(0)) draw(plot1(0)) ; draw(plot12(0)) frame(wks) ; res1@gsnLeftString = "Atlantic_"+comp(i) plot2(0) = gsn_csm_contour(wks2,moca(lev|0:33,lat|:),res1) plot2(1) = gsn_csm_contour(wks2,moca(lev|33:59,lat|:),res2) ref3 = gsn_add_polyline(wks2,plot2(0),refx1,refy1,pres) ref4 = gsn_add_polyline(wks2,plot2(1),refx2,refy2,pres) plot22(0) = gsn_csm_contour(wks2,zmatp(lev|0:33,lat|:),res3) plot22(1) = gsn_csm_contour(wks2,zmatp(lev|33:59,lat|:),res4) newplot2 = gsn_attach_plots(plot2(0),plot2(1),attres11,attres12) newplot22 = gsn_attach_plots(plot22(0),plot22(1),attres11,attres12) overlay(plot2(0),plot22(0)) draw(plot2(0)) frame(wks2) ; res1@gsnLeftString = "PacInd_"+comp(i) plot3(0) = gsn_csm_contour(wks3,mocpi(lev|0:33,lat|:),res1) plot3(1) = gsn_csm_contour(wks3,mocpi(lev|33:59,lat|:),res2) ref5 = gsn_add_polyline(wks3,plot3(0),refx1,refy1,pres) ref6 = gsn_add_polyline(wks3,plot3(1),refx2,refy2,pres) plot32(0) = gsn_csm_contour(wks3,zmpitp(lev|0:33,lat|:),res3) plot32(1) = gsn_csm_contour(wks3,zmpitp(lev|33:59,lat|:),res4) newplot3 = gsn_attach_plots(plot3(0),plot3(1),attres11,attres12) newplot32 = gsn_attach_plots(plot32(0),plot32(1),attres11,attres12) overlay(plot3(0),plot32(0)) draw(plot3(0)) frame(wks3) end do end