;*************************************************************** ; Load Libaries ;*************************************************************** 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/wrf/WRF_contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW_add.ncl" ;*************************************************************** ; Begin user input: Change directories, filenames, labels ;*************************************************************** version = "Default" forecast_days = "12" dates = "9-2-2006 0z to 9-14-2006 0z" diri = "~/myWorkDir/Current_work/12-Day_forecast_Default_folder/" fili = "WRF-"+version+"_"+forecast_days+"day-forecast_rain.nc" filA = "Original_R2_12_day_forecast_plevel_surface_and_rain_diagnostics.nc" ; End user input ;*************************************************************** ; Format of output file ;*************************************************************** ;type = "x11" type = "pdf" ;type = "eps" ;type = "ps" ;*************************************************************** ; WRF model = addfile(diri+fili,"r") ;******************************** ; WRF Precipitation WRF_prate_day = model->rain WRF_lat = model->south_north WRF_lon = model->west_east WRF_prate_day = lonFlip(WRF_prate_day) WRF_lat=WRF_lat(::-1) printVarSummary(WRF_prate_day) ;******************************** ; NCPR2 obs = addfile(diri+filA,"r") R2_prate = obs->prate ;*************************************************************** ; Adjust R2 Data. ;*************************************************************** ; Raw data is in units of kg/(m^2*s) but the data is supposed to be a 6-hourly mean ; We want units in mm/day. Conversion of units requires: ; 1) kg/m^2 = 1mm ; 2) kg/(m^2*s) = 1mm/s ; 3) kg/(m^2*s) * 3600s = 3600kg/(m^2*hr) = mm/hr R2_prate_hr = ((R2_prate * 3600) *6) copy_VarMeta(R2_prate, R2_prate_hr) R2_prate_hr =R2_prate_hr(:,::-1,:) R2_prate_hr@units = "mm/hr" ;*************************************************************** ; Create Daily precipitation from NCPR2 file ;*************************************************************** p = R2_prate_hr dimp = dimsizes( p ) ntim = dimp(0) nlat = dimp(1) mlon = dimp(2) nhr = 4 ; 6 hourly files, 4 x day nhrdim = (ntim-1)/nhr ; number of nhr-hour segments; disregard last ts ptot = new((/nhrdim,nlat,mlon/), typeof(p), getFillValue(p)) ntStrt = 0 ntLast = nhr-1 do nt=0,ntim-1,nhr ptot(nt/nhr,:,:) = dim_sum_n( p(ntStrt:ntLast,:,:) , 0) ;SUM of 4xdaily timesteps, mm/day ntStrt = ntStrt + nhr ntLast = ntLast + nhr end do R2_prate_day = ptot copy_VarMeta(R2_prate_hr,R2_prate_day) ; copy Metadata R2_prate_day@units = "mm/day" ; Correct units printVarSummary(R2_prate_day) ; Check variable delete(p) delete(R2_prate) delete(R2_prate_hr) ;*************************************************************** ; INTERPOLATION SECTION ;*************************************************************** ; Interpolate WRF data to the NCPR2 lat-lon grid, 2.5x2.5 degrees, using areahi2lores interpolation ; We must interpolate the WRF from a 0.5x0.5 degree grid to 2.5x2.5 degrees ; Note: We could use the lat and lons from the NCPR2 file but this method should be more robust. ;*************************************************************** latGRID = fspan(-90, 90, 73) latGRID!0 = "latGRID" latGRID@units = "degrees_north" latGRID&latGRID = latGRID NLAT = dimsizes(latGRID) lonGRID = fspan( 0, 357.5, 144) lonGRID!0 = "lonGRID" lonGRID@units = "degrees_east" lonGRID&lonGRID = lonGRID MLON = dimsizes(lonGRID) ;;;Now that we have the latitudes and longitudes of the new grid, interpolate the WRF variable opt = True ;opt@critpc = 50 ; default is 100% newWRF_prate_day = area_hi2lores_Wrap (WRF_lon,WRF_lat, WRF_prate_day, False, 1, lonGRID, latGRID, opt) newWRF_prate_day!0 = "time" newWRF_prate_day!1 = "lat" newWRF_prate_day&lat = latGRID newWRF_prate_day!2 = "lon" newWRF_prate_day&lon = lonGRID newWRF_prate_day@lat = "degrees_north" newWRF_prate_day@lon = "degrees_east" ;**************************** newR2_prate_day = area_hi2lores_Wrap (R2_prate_day&lon_98,R2_prate_day&lat_98, R2_prate_day, True, 1, lonGRID, latGRID, opt) newR2_prate_day!0 = "time" newR2_prate_day!1 = "lat" newR2_prate_day&lat = latGRID newR2_prate_day!2 = "lon" newR2_prate_day&lon = lonGRID newR2_prate_day@lat = "degrees_north" newR2_prate_day@lon = "degrees_east" ;****************************** ; Clipping ;****************************** newR2_prate_day = where(newR2_prate_day .lt.0.0, -9999, newR2_prate_day ) newWRF_prate_day = where(newWRF_prate_day .lt.0.0, -9999,newWRF_prate_day) printMinMax(newR2_prate_day,True) printMinMax(newWRF_prate_day,True) ;******************************* ;Longitude FLIP ;******************************* newR2_prate_day = lonFlip(newR2_prate_day) newWRF_prate_day = lonFlip(newWRF_prate_day) ;*************************************************************** ; Adjust WRF to fit to NCPR2 data ; Fit each other ;*************************************************************** newR2_prate = newR2_prate_day(1::,{5:15},{-20:30}) newWRF_prate = newWRF_prate_day(:,{5:15},{-20:30}) printVarSummary(newR2_prate) printVarSummary(newWRF_prate) delete(newR2_prate_day) delete(newWRF_prate_day) ;*************************************************************** ; Hov ;*************************************************************** R2_HOV = dim_avg_n_Wrap(newR2_prate,1) wrfHOV = dim_avg_n_Wrap(newWRF_prate,1) delete(newR2_prate) delete(newWRF_prate) diff = wrfHOV - R2_HOV copy_VarMeta(wrfHOV,diff) diff@description= "Difference: WRF minus TRMM" printVarSummary(R2_HOV) printMinMax(R2_HOV,True) printVarSummary(wrfHOV) printMinMax(wrfHOV,True) printVarSummary(diff) printMinMax(diff,True) ;*************************************************************** ; create plot ;*************************************************************** type@wkColorModel = "cmyk" ; make output cmyk wks = gsn_open_wks(type,"wrf_R2_rainHov_Panel_abc-2") ; open a ps file ;;gsn_define_colormap(wks,"BlAqGrYeOrRe") ; choose colormap gsn_merge_colormaps(wks,"BlAqGrYeOrRe","BlWhRe") res = True ; plot mods desired plot = new(3,graphic) res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@gsnMaximize = True ; biggest Plot Possible ;************************************************ ; Use These next set of commands if using color res@cnFillOn = True ; turn on color fill res@gsnSpreadColors = True ; use full range of colors res@gsnSpreadColorStart = 0 res@gsnSpreadColorEnd = 101 res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off line labels ; CMYK colors ;************************************************ res@lbOrientation = "vertical" ; vertical labelbar res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res@cnMinLevelValF = 0. res@cnMaxLevelValF = 30. res@cnLevelSpacingF = 2. ;res@tiMainString = "Daily Accumulated Precipitation, "+forecast_days+" days: ~C~ "+dates res@gsnMajorLonSpacing = 5. res@tmXBTickSpacingF = 5 ;res@tiXAxisString = "Longitude" res@tmXBLabelsOn = False res@tmXTOn = False ; top tickmarks res@tmYROn = False res@tiYAxisString = "Elapsed Time: SEP 2 to 14, 2006" n=(dimsizes(R2_HOV(:,0))) values = ispan(1,n-1,1) labels = ispan(2,n,1) res@tmYLMode = "Explicit" res@tmYLValues = values res@tmYLLabels = labels res@tmYLTickSpacingF = 1 res@gsnLeftString = "Lats: 5N to 15N" res@gsnCenterString = "NCPR2" res@gsnRightString = "Units: mm/day" res@pmTickMarkDisplayMode = "Always" plot(0) = gsn_csm_hov(wks,R2_HOV,res) res@tiYAxisString = "Elapsed Time: SEP 2 to 14, 2006" ;res@gsnLeftString = "Lats: 5N to 15N" res@gsnRightString = "Units: mm/day" res@gsnCenterString = "WRF" plot(1)=gsn_csm_hov(wks,wrfHOV,res) ;********************************* ; Experimental: Difference Plot ;********************************* res@gsnSpreadColorStart = 131 res@gsnSpreadColorEnd = 171 res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res@cnMinLevelValF = -30. res@cnMaxLevelValF = 30. res@cnLevelSpacingF = 5. ;res@tiMainString = "Hovmoller: Difference of Daily Precipitation ~C~ "+forecast_days+" days: "+dates delete(res@tmXBLabelsOn) res@tiXAxisString = "Longitude" res@gsnCenterString = "WRF - NCPR2" plot(2) = gsn_csm_hov(wks,diff,res) ;*************************************************************** ; create panel ;*************************************************************** resP = True ; modify the panel plot ;resP@ganPanelTop = .95 ;resP@gsnPanelBottom = .35 resP@txString = "Precipitation: NCPR2 and WRF "+forecast_days+"-Day Forecast";;; ~C~ "+dates ;resP@lbLabelFontHeightF = 0.007 ; make labels smaller resP@gsnMaximize = True resP@gsnPanelBottom = 0.05 ;resP@gsnPaperOrientation = "Portrait" ; force portrait resP@gsnPanelFigureStrings = (/"a)","b)","c)"/) ; add strings to panel gsn_panel(wks,plot,(/3,1/),resP) ; now draw as one plot