;***********Load the ncl internal libraries*********; 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/shea_util.ncl" begin ;get the data era_all_pas = addfile("/media/SAMSUNG/REGCM_SST/EIN15/2015/ERAINT_0.75/MSLP/pas_eraint_0.75_daily.nc","r") era_all_ta = addfile("/media/SAMSUNG/REGCM_SST/EIN15/2015/ERAINT_0.75/TA/ta_eraint_0.75_daily.nc","r") era_all_vwnd = addfile("/media/SAMSUNG/REGCM_SST/EIN15/2015/ERAINT_0.75/VERTWND/vertwnd_eraint_0.75_daily.nc","r") era_all_qas = addfile("/media/SAMSUNG/REGCM_SST/EIN15/2015/ERAINT_0.75/QAS/qas_eraint_0.75_daily.nc","r") ;which region latsouthraw = -67 latnorthraw = 32 lonwestraw = -160 loneastraw = 10 ;select pressure levels lev1 = 10 lev2 = 36 ;select dates seltime0 = 20 seltime1 = 21;22 mar seltime2 = 22;23 mar seltime3 = 23;24 mar seltime4 = 24 seltime5 = 25 seltime6 = 26 ;pressure levels levera = era_all_ta->level(lev1:lev2) levera_pte = era_all_ta->level(lev1:lev2) levera_pte = levera_pte*100 time = era_all_ta->time latd = era_all_ta->latitude({latitude|latsouthraw:latnorthraw}) ;extract surface pressure and apply lonFlip era_var_ps = era_all_pas->sp era_flip_ps = lonFlip(era_var_ps) era_ps = short2flt(era_flip_ps(time|:,{latitude|latsouthraw:latnorthraw},{longitude|lonwestraw:loneastraw})) era_ps@_FillValue = -32767 lond = era_ps&longitude fint = era_ps*0.01 fint!0 = "time" fint!1 = "latitude" fint!2 = "longitude" fint&time = time fint&latitude = latd fint&longitude = lond ;extract air temperature and apply lonFlip era_var_ta = era_all_ta->t era_flip_ta = lonFlip(era_var_ta) era_ta = short2flt(era_flip_ta(time|:,level|lev1:lev2,{latitude|latsouthraw:latnorthraw},{longitude|lonwestraw:loneastraw})) era_ta@_FillValue = -32767 fint_ta = era_ta ;extract omega and apply lonFlip era_var_vwnd = era_all_vwnd->w era_flip_vwnd = lonFlip(era_var_vwnd) era_vwnd = short2flt(era_flip_vwnd(time|:,level|lev1:lev2,{latitude|latsouthraw:latnorthraw},{longitude|lonwestraw:loneastraw})) era_vwnd@_FillValue = -32767 fint_vwnd = era_vwnd P = conform(fint_vwnd,levera,1) ;omega to vertical velocity fint_vwnd_ms = omega_to_w(fint_vwnd,P,fint_ta) ;extract potential temperature and apply lonFlip era_var_pottemp = era_all_ta->t era_flip_pottemp = lonFlip(era_var_pottemp) pottemp_all = short2flt(era_flip_pottemp(time|:,level|lev1:lev2,{latitude|latsouthraw:latnorthraw},{longitude|lonwestraw:loneastraw})) pottemp_all@_FillValue = -32767 ;extract spec.hum. and apply lonFlip era_var_pr = era_all_qas->q era_flip_pr = lonFlip(era_var_pr) era_precw = short2flt(era_flip_pr(time|:,level|lev1:lev2,{latitude|latsouthraw:latnorthraw},{longitude|lonwestraw:loneastraw})) era_precw@_FillValue = -32767 fint_pr = era_precw fint_pr!0 = "time" fint_pr!1 = "level" fint_pr!2 = "latitude" fint_pr!3 = "longitude" fint_pr&time = time fint_pr&level = levera fint_pr&latitude = latd fint_pr&longitude = lond ;theta-e calculation using wrf_eth P1 = conform(fint_pr,levera_pte,1) pottemp_all_e = wrf_eth(fint_pr,pottemp_all,P1) pottemp_all_e@_FillValue = pottemp_all@_FillValue ;cross section coordinates latant1 = -23.65;-23.65 lonant1 = -71;-70.4 ;attributes and dimensions of theta-e array pottemp_all_e!0 = "time" pottemp_all_e!1 = "level" pottemp_all_e!2 = "latitude" pottemp_all_e!3 = "longitude" pottemp_all_e&time = pottemp_all&time pottemp_all_e&level = pottemp_all&level pottemp_all_e&latitude = pottemp_all&latitude pottemp_all_e&longitude = pottemp_all&longitude ;masking based on pressure level do i=0,dimsizes(levera)-1 pottemp_all(time|:,level|i,latitude|:,longitude|:) = fint_ta(time|:,level|i,latitude|:,longitude|:)*tofloat((1000/levera(i))^0.286) end do do tt = 0,dimsizes(time)-1 do gg = 0,dimsizes(levera)-1 print("Now working on "+levera(gg)+"hPa") pottemp_all(time|tt,level|gg,latitude|:,longitude|:) = where(fint(time|tt,latitude|:,longitude|:).lt.levera(gg),pottemp_all@_FillValue,pottemp_all(time|tt,level|gg,latitude|:,longitude|:)) end do end do do tt = 0,dimsizes(time)-1 do gg = 0,dimsizes(levera)-1 print("Now working on "+levera(gg)+"hPa") fint_vwnd_ms(time|tt,level|gg,latitude|:,longitude|:) = where(fint(time|tt,latitude|:,longitude|:).lt.levera(gg),fint_vwnd_ms@_FillValue,fint_vwnd_ms(time|tt,level|gg,latitude|:,longitude|:)) end do end do do tt = 0,dimsizes(time)-1 do gg = 0,dimsizes(levera)-1 print("Now working on "+levera(gg)+"hPa") pottemp_all_e(time|tt,level|gg,latitude|:,longitude|:) = where(fint(time|tt,latitude|:,longitude|:).lt.levera(gg),pottemp_all_e@_FillValue,pottemp_all_e(time|tt,level|gg,latitude|:,longitude|:)) end do end do ;longitude pressure pt= pottemp_all(time|seltime4,level|:,{latitude|-23.65},{longitude|-90:-60}) vwnd_ms = fint_vwnd_ms(time|seltime4,level|:,{latitude|-23.65},{longitude|-90:-60}) ;time pressure pt_time= pottemp_all(time|14:30,level|:,{latitude|-23.65},{longitude|-70.4}) vwnd_time = fint_vwnd_ms(time|14:30,level|:,{latitude|-23.65},{longitude|-70.4}) print (pt_time&time) dims_sta = getvardims(pt_time) pt_new = pt_time($dims_sta(1)$|:,$dims_sta(0)$|:) vwnd_new = vwnd_time($dims_sta(1)$|:,$dims_sta(0)$|:) pottemp_all_e!0 = "time" pottemp_all_e!1 = "level" pottemp_all_e!2 = "latitude" pottemp_all_e!3 = "longitude" pottemp_all_e&time = pottemp_all&time pottemp_all_e&level = pottemp_all&level pottemp_all_e&latitude = pottemp_all&latitude pottemp_all_e&longitude = pottemp_all&longitude precw = pottemp_all_e(time|seltime4,level|:,{latitude|-23.65},{longitude|-90:-60}) precw_time = pottemp_all_e(time|14:30,level|:,{latitude|-23.65},{longitude|-70.4}) precw_new = precw_time($dims_sta(1)$|:,$dims_sta(0)$|:) wks=gsn_open_wks("ps","ERAINT_POT_TEMP_VERTWIND_POTE") gsn_define_colormap (wks,"BlueDarkRed18") res = True res = True ; plot mods desired res@cnFillOn = True ; color fill res@cnLinesOn = False res@cnLineLabelsOn = False ; no contour labels res@cnInfoLabelOn = False ; no contour info label res@tiYAxisString = "" res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = -16;0.0010;-2 ; set min contour level res@cnMaxLevelValF = 16;0.0080; 2 ; set max contour level res@cnLevelSpacingF = 2;0.0005; 0.2 ; set contour interval res@gsnSpreadColors = True res@lbLabelBarOn = False res@tiMainString = "" res@gsnCenterString = "" res@gsnCenterStringFontHeightF = 0.02 res@tmXBMode = "Explicit" res@tmXBValues = (/-90,-85,-80,-75,-70,-65,-60/) res@tmXBLabels = (/"90W","85W","80W","75W","70W","65W","60W"/) res@gsnLeftString = "" res@gsnRightString = "" res@tmXTOn = False res@vpXF = 0.13 ; change x-scale res@vpWidthF = 0.75 ; change height and width res@vpHeightF = 0.45 res@cnSmoothingOn = True res@cnMissingValFillColor = "burlywood4" res@gsnDraw = False ; do not draw the plot res@gsnFrame = False ; do not advance the frame sres = True ; set up a second resource list sres@gsnDraw = False ; do not draw the plot sres@gsnFrame = False ; do not advance the frame sres@cnLevelSelectionMode = "ManualLevels" ; manually set the contour levels with the following 3 resources sres@cnMinLevelValF = 300. ; set the minimum contour level sres@cnMaxLevelValF = 400. ; set the maximum contour level sres@cnLevelSpacingF = 5. ; set the interval between contours sres@cnLineDashSegLenF = 0.2 ; assist in controlling concentration sres@cnLineLabelInterval = 3 ; default = 2 sres@cnLineLabelPlacementMode = "constant" ; choose constant label method sres@cnLineLabelFontHeightF = 0.014 sres@cnLineLabelsOn = True ; no contour labels sres@cnInfoLabelOn = False ; no contour info label sres@gsnLeftString = "" sres@gsnRightString = "" sres@tmYRMode = "Automatic" sres@cnSmoothingOn = True sres@cnFillMode = "RasterFill" ;sres@cnLineThicknessF = 2.5 res@tiYAxisString = "Pressure (hPa)" res@tiYAxisFontHeightF = 0.02 res@tiYAxisOffsetYF = -0.004 tpwres = True ; set up a second resource list tpwres@gsnDraw = False ; do not draw the plot tpwres@gsnFrame = False ; do not advance the frame tpwres@cnLevelSelectionMode = "ManualLevels" ; manually set the contour levels with the following 3 resources tpwres@cnMinLevelValF = 300. ; set the minimum contour level tpwres@cnMaxLevelValF = 400. ; set the maximum contour level tpwres@cnLevelSpacingF = 5. ; set the interval between contours tpwres@cnLineDashSegLenF = 0.24 ; assist in controlling concentration tpwres@cnLineLabelInterval = 2 ; default = 2 tpwres@cnLineLabelPlacementMode = "constant" ; choose constant label method tpwres@cnLineLabelFontHeightF = 0.014 tpwres@cnLineLabelsOn = True ; no contour labels tpwres@cnInfoLabelOn = False ; no contour info label tpwres@gsnLeftString = "" tpwres@gsnRightString = "" tpwres@tmYRMode = "Automatic" tpwres@cnSmoothingOn = True tpwres@cnFillMode = "RasterFill" tpwres@cnLineLabelFontColor= "Gray45" tpwres@cnLineColor = "Gray45" tpwres@cnLineDashPattern = 2 tpwres@cnLineThicknessF = 3 plot = new(2,graphic) res@gsnCenterString = "2015-03-25, 23.65S" plotA = gsn_csm_pres_hgt(wks,vwnd_ms,res ) ; place holder ;plotB = gsn_csm_pres_hgt(wks,pt,sres ) ; place holder ;overlay (plotA,plotB) plotK = gsn_csm_pres_hgt(wks,precw,sres) overlay (plotA,plotK) plot(0) = plotA res@tmXBMode = "Explicit" delete(res@tmXBValues) delete(res@tmXBLabels) ;res@tmXBValues =pt_time&time res@tmXBValues = (/1009842,1009890,1009938,1009986,1010034,1010082,1010130,1010178,1010226/) res@tmXBLabels = (/"15 Mar","17 Mar","19 Mar","21 Mar","23 Mar","25 Mar","27 Mar","29 Mar","31 Mar"/) res@gsnCenterString = "Antofagasta" plotC = gsn_csm_pres_hgt(wks,vwnd_new,res ) ; place holder ;plotD = gsn_csm_pres_hgt(wks,pt_new,sres ) ; place holder ;overlay (plotC,plotD) plotM = gsn_csm_pres_hgt(wks,precw_new,sres) overlay (plotC,plotM) plot(1) = plotC resP = True resP@gsnPanelLabelBar = True ;resP@gsnFrame = False resP@gsnMaximize = True resP@lbOrientation = "Vertical" resP@lbLabelAutoStride = True ;resP@gsnPanelBottom = 0.1 ;resP@gsnPanelLeft = 0.1 resP@pmLabelBarWidthF = 0.045 resP@pmLabelBarHeightF = 0.55 ;resP@lbLabelAlignment = "InteriorEdges" resP@lbTitleOn = True ; turn on title resP@lbTitleString = "cm/s" ; title string (Pa*100/100) resP@lbTitlePosition = "Bottom" ; title position resP@lbTitleFontHeightF= .014 ; make title smaller resP@lbTitleDirection = "Across" ; title direction resP@lbLabelStride = 2 resP@pmLabelBarParallelPosF = -0.01 resP@lbLabelFontHeightF = 0.014 resP@pmLabelBarOrthogonalPosF = 0.01 ;resP@gsnPanelRowSpec = True ;resP@gsnPanelCenter = False ;resP@gsnPanelYWhiteSpacePercent = 0.3 resP@txString = "" ;resP@gsnPanelFigureStrings= (/"a)","b)","c)","d)","e)"/) ; add strings to panel gsn_panel(wks,plot,(/2,1/),resP) ;draw (plot) ;frame(wks) end