;************************************************* ; PLOT PRECIPITATION (DAILY EVOLUTION) ;;ERA-INT dateoffile = 198710 months = 10 years = 1987 day_start = 15 day_end = 18 ; INPUTFILE and SETTINGS PLOT ; directory input dir_in = "/gpfs_ci/data/projects/tchristoudias/ei/" ; 12 h acc dir_in2 = "/gpfs_ci/data/adevries/APHRODITE/" ; directory output(plot) dir_out = "/gpfs_ci/data/adevries/NCL/PLOTS/ECMWF/ERA-Int/PANEL/"+dateoffile+"/" ;************************************************ 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/contrib/calendar_decode2.ncl" ;************************************************ begin ;************************************************ ; read in netCDF file ;************************************************ ;;ERA-INT a = addfile (dir_in2+"APHRO_ME_025deg_V1101."+years+".nc","r") b = addfile (dir_in+"ei_fc_totprec_gg128_step12_"+dateoffile+".grib","r") ; 12 h acc ; APHRO prec_aphro = a->precip ; read precipitation lon3 = a->longitude ; read longitude lat3 = a->latitude ; read latitude prec_aphro @_FillValue = -99.9 ; set missing value prec_aphro @long_name = "daily precipitation (Aphrodite)" ;name of simulation in title of plot prec_aphro @units = "[mm/day]" time_double = a->time time_3 = doubletointeger(time_double) prec_daily_aphro = new((/((day_end-day_start)+1),dimsizes(lat3),dimsizes(lon3)/),float) prec_daily_aphro = 0 ; ERA-INT prec_step12 = b->TP_GDS4_SFC(:,{1:61},{-1:91}) ; read convective precipitation time = b->initial_time0_hours ; read in time date = calendar_decode2(time(:),0) lon = b->g4_lon_2({-1:91}) lat = b->g4_lat_1({1:61}) yyyy = date(:,0) mm = date(:,1) dd = date(:,2) prec_daily_era = new((/((day_end-day_start)+1),dimsizes(lat),dimsizes(lon)/),float) copy_VarMeta(prec_step12(0:(day_end-day_start),:,:),prec_daily_era) prec_daily_era = 0 do d=day_start,day_end j = d-day_start print(d) print(j) ;*********************************************************************************************************************************************** ; PRECIP CALCULATION Aphrodite ;*********************************************************************************************************************************************** ddd = day_of_year(years,months,d) prec_daily_aphro(j,:,:) = prec_aphro(ddd-1,:,:) ;************************************************* ; CALCULATE PRECIPITATION DAILY ;************************************************* prec_daily_era(j,:,:) = (prec_step12(((d-1)*2),:,:)+prec_step12((((d-1)*2)+1),:,:))*1000 ; change unities (from m (12h acc). to mm (daily acc)) end do ;*********************************************************************************************************************************************** ; END CALCULATION PRECIP CALCULATION ;*********************************************************************************************************************************************** ;************************************************ ; metadata, long_name and unities ;************************************************ prec_daily_era @long_name = "tot. precip. (acc fcst 0-12 & 12-24 hrs)" ;name of simulation in title of plot prec_daily_era @units = "[mm/day]" ;************************************************* ; PLOT SETTINGS ;************************************************* res = True res@gsnAddCyclic = False res@mpGeophysicalLineThicknessF = 1.0 res@cnInfoLabelOn = False ; Turn off info label res@gsnLeftString = "" res@gsnRightString = "" res@mpFillOn = False res@mpMinLonF = 20 ; specify the plot domain res@mpMaxLonF = 50 res@mpMinLatF = 15 res@mpMaxLatF = 40 res@mpOutlineOn = True ; turn the map outline on res@gsnDraw = False ; do not draw the plot res@gsnFrame = False ; do not advance the frame res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels res@cnLevels = fspan(1,11,11) ; set the contour levels for PREC res@lbLabelBarOn = False ; turn off individual cb's res@cnLineLabelsOn = False ; do not use line labels res@cnFillOn = True ; color fill res@cnLinesOn = False ; do not draw contour lines res@gsnSpreadColors = True ; automatically choose the fill colors res@cnMissingValFillColor = "black" ;************************************************* ; Open pdf ;************************************************* ;ERA-int wks = gsn_open_wks("pdf",dir_out+"PRECIP_ERA-Int_Aphro_daily_CT2_ncl-talk_"+yyyy(day_start)+"-"+mm(day_start)+"-"+day_start+"_to_"+day_end) ; open a pdf file gsn_define_colormap(wks,"precip_11lev") ;************************************************ ; create plots ;************************************************ plot_precip_era_aphro = new(((day_end-day_start)+1)*2,graphic) ;************************************************ ; PLOTTING ;************************************************ do i=0,(day_end-day_start) d = (day_start-1)+i print("i = "+i) print("d = "+d) print("day = "+dd(d*2) ) ;; make shaded plot of PRECIP_TOT res@gsnCenterStringOrthogonalPosF = 0.02 res@gsnCenterStringFontHeightF = 0.025 res@gsnCenterString = yyyy(d*2)+"-"+mm(d*2)+"-"+dd(d*2) ; set the main title if (i.le.2) res@tmXBLabelsOn = False ; do not draw bottom labels res@tmXBOn = False ; no bottom tickmarks else res@tmXBLabelsOn = True ; do not draw bottom labels res@tmXBOn = True ; no bottom tickmarks end if res@tmXTLabelsOn = False ; do not draw bottom labels res@tmXTOn = False ; no bottom tickmarks res@tmYRLabelsOn = False ; no right labels res@tmYROn = False ; no right tickmarks res@tmYLLabelsOn = True ; no right labels res@tmYLOn = True ; no right tickmarks plot_precip_era_aphro(i*2) = gsn_csm_contour_map(wks,prec_daily_era(i,:,:),res) res@tmYLLabelsOn = False ; no right labels res@tmYLOn = False ; no right tickmarks plot_precip_era_aphro((i*2)+1) = gsn_csm_contour_map(wks,prec_daily_aphro(i,:,:),res) delete(d) end do ; loop day_start to day_end pnlres = True pnlres@gsnMaximize = True ; Maximize in frame pnlres@gsnPanelYWhiteSpacePercent = 3 pnlres@gsnPanelLabelBar = True ; add common colorbar pnlres@lbTitleFontHeightF= .015 ; make title smaller pnlres@lbLabelFontHeightF = 0.01 ; make labels smaller ;;pnlres@lbTopMarginF = 2 pnlres@gsnPanelTop = 0.97 pnlres@gsnPanelBottom = 0.01 pnlres@pmLabelBarOrthogonalPosF = -0.01 ; move whole thing down txres = True txres@txFontHeightF = 0.015 gsn_text_ndc(wks," ERA-Interim Aphrodite ",0.5,0.99,txres) gsn_panel(wks,plot_precip_era_aphro,(/4,2/),pnlres) ;************************************************ draw(plot_precip_era_aphro) delete(plot_precip_era_aphro) delete(res) delete(wks) delete(yyyy) delete(mm) delete(dd) delete(months) delete(years) delete(time) delete(date) delete(time_3) delete(time_double) delete(lon) delete(lat) delete(lon3) delete(lat3) delete(dir_in) delete(dir_out) delete(b) delete(a) delete(prec_step12) delete(prec_aphro) delete(day_start) delete(day_end) delete(dateoffile) delete(prec_daily_era) delete(prec_daily_aphro) end