Re: compare wrf output with MEERA

Hi,

You didn't say what your doubt was. Was there an error?

By the way, if you want to append the same thing to the end of an array of NCL strings, you don't need a do loop. You can do this with one line:

all_files1 = all_files1+".nc"

--Mary

On Mar 20, 2013, at 10:18 PM, Neerju Shrestha wrote:

> 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

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Mar 22 15:07:25 2013