compare wrf output with MEERA

Hi,

I'm trying to plot the wrf output and MERRA data to compare the surface temperature. could pleas help me with the script. below are the details of files and script that i'm using. i've highlight the line with red color where i've some doubt.

wrfoutput =wrfout_d02_2011-01-01_00:00:00
Dir =/Desktop/wrfout-2011

MEERA =MERRA300.prod.assim.inst3_3d_asm_Cp.20100101.SUB.nc
Dir=/Desktop/MEERA/2011

sample usage:
./ncl month=1 domain=2 year=2011 'WRFDataDir='/Desktop/wrfout-2011/"' 'RNLDataDir="/Desktop/MEERA/2011/"' ts_compare.ncl

NCL Script

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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/popRemap.ncl"
begin

if (.not. isvar("year")) then ; is year on command line?
      year = 2011
  end if

if (.not. isvar("month")) then ; is month on command line?
      month = 1
  end if

if (.not. isvar("domain")) then ; is domain on command line?
      domain = 2
end if

;Dir where wrfout is stored
if (.not. isvar("WRFDataDir")) then ; is DataDir on the command line ?
WRFDataDir = "/Desktop/wrfout-2011/"
end if

;Dir where reanalysis data is stored
if (.not. isvar("RNLDataDir")) then ; is DataDir on the command line ?
RNLDataDir = "/Desktop/MEERA/2011/"
end if

monstring = (/ "Zero", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" /)

; Plot monthly surface temperature

if (month .le. 9) then
    all_files1 = systemfunc ("ls " + WRFDataDir + "wrfout_d0" + domain + "_" + year + "-0" + month + "*")
    all_files4 = systemfunc ("ls " + RNLDataDir + "MERRA300.prod.assim.inst3_3d_asm_Cp." + year + "0" + month + "*.nc")
else
    all_files1 = systemfunc ("ls " + WRFDataDir + "wrfout_d0" + domain + "_" + year + "-" + month + "*")
    all_files4 = systemfunc ("ls " + RNLDataDir + "MERRA300.prod.assim.inst3_3d_asm_Cp." + year + month + "*.nc")
end if

l = dimsizes(all_files1)
do i = 0, l-1
   all_files1(i) = all_files1(i)+".nc"
end do

f1 = addfiles (all_files1, "r")

vardims = getfilevardimsizes(f1[0],"TSK")
sn = vardims(1)
we = vardims(2)

latwrf = f1[0]->XLAT
maxlat = max(latwrf)
minlat = min(latwrf)

lonwrf = f1[0]->XLONG
maxlon = max(lonwrf)
minlon = min(lonwrf)
;print(maxlat)
;print(minlat)
;print(maxlon)
;print(minlon)

ListSetType (f1, "cat")
p1 = dim_avg_n_Wrap(f1[:]->TSK,0)

nlat = fspan(minlat, maxlat, sn)
nlon = fspan(minlon, maxlon, we)
nlat@units = "degrees_north"
nlon@units = "degrees_east"

p1&south_north = nlat
p1&west_east = nlon

f4 = addfiles (all_files4, "r")
ListSetType (f4, "cat")
l = dimsizes(all_files4)
p41 = new((/l,181,360/),"float")
do i=0,l-1
 p41(i,:,:) = f4[i]->TMP_3_SFC_10(::-1,:)
end do
p41!0 = "time"
p42 = dim_avg_n_Wrap(p41,0)
p4 = p42({minlat:maxlat},{minlon:maxlon})

dims_p4 = dimsizes(p4)
dims_p4_lat = dims_p4(0)
dims_p4_lon = dims_p4(1)

Lat = fspan(minlat, maxlat, dims_p4_lat)
Lon = fspan(minlon, maxlon, dims_p4_lon)
Lat@units = "degrees_north"
Lon@units = "degrees_east"

v1 = area_hi2lores_Wrap (p1&west_east,p1&south_north, p1 , True, 1, Lon, Lat, False)
v4 = p4

v1!0 = "lat"
v1!1 = "lon"

v1&lat = Lat
v1&lon = Lon

v4!0 = "lat"
v4!1 = "lon"
v4&lat = Lat
v4&lon = Lon

wks = gsn_open_wks("pdf",monstring(month))
gsn_define_colormap(wks,"rainbow")
plot = new(2,graphic)
plot2 = new(1,graphic)

res = True ; plot mods desired
res@cnFillOn = True ; turn on color
res@cnLinesOn = False ; turn off contour lines
res@gsnSpreadColors = True ; use full colormap
res@lbLabelBarOn = False ; turn off individual label bars
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
res@gsnAddCyclic = False ; data already has cyclic point
res@cnLineLabelsOn = False

;res@cnLevelSelectionMode = "ManualLevels"; manual set levels so lb consistent
;res@cnMinLevelValF = 82000 ; min level
;res@cnMaxLevelValF = 102000 ; max level
;res@cnLevelSpacingF = 100 ; contour interval
res@gsnLeftString = "Surface temperature of " + monstring(month)
res@gsnRightString = "K"

res@mpMinLatF = minlat
res@mpMaxLatF = maxlat
res@mpMinLonF = minlon
res@mpMaxLonF = maxlon

res@tiMainString = "cu_physics = Betts-Miller"
plot(0) = gsn_csm_contour_map_ce(wks,v1,res)
res@tiMainString = "Reanalysis"
plot(1) = gsn_csm_contour_map_ce(wks,v4,res)

res_U = True
res_U@gsnFrame = False
res_U@gsnPanelTop = 0.90
res_U@gsnPanelBottom = 0.50
res_U@gsnPanelLabelBar = True
res_U@lbLabelAutoStride = True
res_U@txString = "Surface temperature of " + monstring(month) + " " + year

gsn_panel(wks,plot,(/1,2/),res_U) ; now draw as one plot

diff1 = v1-v4
diff1!0 = "lat"
diff1!1 = "lon"
diff1&lat = Lat
diff1&lon = Lon

gsn_define_colormap(wks,"hotcold_18lev")
res@cnLevelSelectionMode = "ManualLevels"; manual set levels so lb consistent
res@cnMinLevelValF = -2
res@cnMaxLevelValF = 2
res@cnLevelSpacingF = .25

res@tiMainString = "(cu_physics = Betts-Miller) - (Reanalysis)"
plot2(0) = gsn_csm_contour_map_ce(wks,diff1,res)

res_L = True
res_L@gsnFrame = False
res_L@gsnPanelTop = 0.40
res_L@gsnPanelLabelBar = True
res_L@lbLabelAutoStride = True
res_L@txString = "Difference plot"

gsn_panel(wks,plot2,(/1,1/),res_L) ; now draw as one plot
frame(wks)
end

tks & rgds
lily
                                               

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Mar 20 22:18:44 2013