Hi ALL,
I kindly request for your help if you have some knowledge about the NCL. I am plotting zenith hydrostatic delay, zenith wet delay, total zenith path delay and integrated water vapor over the entire domain from the attached script.
; NCL script to plot zenith hydrostatic delay, zenith wet delay, total zenith path delay and
; integrated water vapor over the entire domain
;############################################################################################
; Load functions and procedures
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/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;############################################################################################
begin
;############################################################################################
; first remove any existing plots
system("/bin/rm -f " + "plot_ZHD.pdf")
system("/bin/rm -f " + "plot_ZWD.pdf")
system("/bin/rm -f " + "plot_ZTD.pdf")
system("/bin/rm -f " + "plot_IWV.pdf")
system("/bin/rm -f " + "plot_terrain_ht.pdf")
;-------------------------------------------------------------------------------------------
; open workstations, output to pdf
wkszhd = gsn_open_wks ("pdf","plot_ZHD")
wkszwd = gsn_open_wks ("pdf","plot_ZWD")
wksztd = gsn_open_wks ("pdf","plot_ZTD")
wksiwv = gsn_open_wks ("pdf","plot_IWV")
wkstht = gsn_open_wks ("pdf","plot_terrain_ht")
;-------------------------------------------------------------------------------------------
; set basic resources to apply to all the plots
; change the default map
mpres = True
mpres@mpOutlineOn = True
mpres@mpGeophysicalLineColor = "Black" ; outlines of all continents, islands, and inland water bodies
mpres@mpGeophysicalLineThicknessF = 1.5
mpres@mpNationalLineColor = "Black" ; outlines of all political boundaries
mpres@mpNationalLineThicknessF = 1.5
mpres@mpOutlineBoundarySets = "National"
mpres@mpLimitMode = "LatLon"
mpres@mpMinLatF = -12.5
mpres@mpMaxLatF = 4.6
mpres@mpMinLonF = 28.1
mpres@mpMaxLonF = 41.9
; set other basic resources
res = True
res@InitTime = False
res@Footer = False
res@cnFillOn = True ; create a color fill plot
res@cnLineLabelsOn = True ; Turn on line labels
res@cnLineLabelPerimOn = False
res@tmXBLabelFont = 0.01
res@tmYRLabelFont = 0.01
res@tmXBMajorLengthF = 0.01
res@tmYRMajorLengthF = 0.01
res@lbTitleOn = True
res@lbTitleFontHeightF = 0.013
res@lbLabelOffsetF = 0.02
res@lbLabelFontHeightF = 0.01
;res@lbOrientation = "vertical"
;res@pmLabelBarSide = "right"
; overlay plot onto a map (using default plot and map setting):
pltres = True ; set plot options
mpres = True ; set map options
pltres@NoTitles = False
;##########################################################################################
; make list of the wrfout files and loop through them extracting the variables and making the plots
fList = systemfunc("ls /media/richard/cliffe0704544658/2016nwmdata/march/cliffe01/wrfprd/*")
Files = dimsizes(fList)
; start looping through the *.nc files
do File = 0, Files-1
nc_file = addfile(fList(File),"r")
time = 0 ; the netcdf files are individual for each time variable
;*******************************************************************************************
; extract the variables
; lat and lon coordinates
xlat = nc_file->XLAT(time,:,:) ; latitude
xlon = nc_file->XLONG(time,:,:) ; longitude
;-------------------------------------------------------------------------------------------
; surface pressure
psfc = nc_file->PSFC(time,:,:) ; surface pressure in Pascals at all model levels
;-------------------------------------------------------------------------------------------
; surface temperature (2m temperature)
t2 = nc_file->T2(time,:,:) ; temperature at 2m in Kelvin
;-------------------------------------------------------------------------------------------
; water vapor mixing ratio (also spec. humidity) at 2m (near surface)
q2 = nc_file->Q2(time,:,:) ; specific humidity at 2m in Kg/kg
;-------------------------------------------------------------------------------------------
; perturbation and base-state pressure
p = nc_file->P(time,:,:,:) ; perturbation pressure variable in Pascals at all model levels
pb = nc_file->PB(time,:,:,:) ; base-state pressure variable in Pascals at all model levels
;-------------------------------------------------------------------------------------------
; perturbation potential temperature
t = nc_file->T(time,:,:,:) ; perturbation potential temperature (theta-t0) in Kelvin at all model levels
;-------------------------------------------------------------------------------------------
; water vapor mixing ratio at model levels
qvapor = nc_file->QVAPOR(time,:,:,:) ;water vapor mixing ratio in kg/kg at all model levels
;-------------------------------------------------------------------------------------------
; perturbation and base-state geopotentials (staggered in Z)
ph = nc_file->PH(time,:,:,:) ; perturbation geopotential in m2/s2 at FULL levels
phb = nc_file->PHB(time,:,:,:) ;base_state geopotential in m2/s2 at FULL levels
;*******************************************************************************************
; calculations
; total pressure over the entire domain and model levels
p_total = p + pb ; total pressure over the entire domain
;-------------------------------------------------------------------------------------------
; temperature at model levels
theta = t + 300 ; potential temperature
tk = wrf_tk (p_total, theta) ; temperature in Kelvin for all model levels over the entire domain
;-------------------------------------------------------------------------------------------
; geopotential heights
g = 9.80665 ; standard acceleration due to gravity at mean sea level in ms-2
z_geop = (ph + phb)/g ; geopotential heights on FULL levels in m
dimv = dimsizes(z_geop)
z = 0.5*(z_geop(0:dimv(0)-2,:,:)+z_geop(1:dimv(0)-1,:,:)) ; geopotential heights on HALF levels in m
z!0 = "bottom_top" ; name
z!1 = "south_north" ; the
z!2 = "west_east" ; dimensions
;--------------------------------------------------------------------------------------------
; local surface gravity on the entire domain
ge = 9.780356 ; gravity at the earth's equator in ms-2
dtr = 3.14159/180 ; degrees to radians
g_sfc = ge * (1 + 0.0052885*(sin(xlat*dtr))^2 - (0.0000059*(sin(2*xlat*dtr)^2))) ; gravity on the surface of entire domain
;--------------------------------------------------------------------------------------------
; gravity at model levels over the entire domain
n_level = 44 ; levels 0 - 43 (indices,unstaggered)
x_dim = 193
y_dim = 241
xlat_new = new((/x_dim,y_dim/),float) ; create new variable with two dimensions
xlat_new!0 = "south_north" ; assign
xlat_new!1 = "west_east" ; the dimensions
copy_VarAtts(xlat,xlat_new) ; copy variable attributes from xlat to xlat_new
do i = 0, x_dim-1
do j = 0, y_dim-1
xlat_new(i,j) = (/xlat(i,j)/) ; copy the data only (without dimensions) from xlat to xlat_new
end do
end do
; replicate the xlat_new variable to all the model levels
xlat_level = new((/n_level,x_dim,y_dim/),float)
xlat_level!0 = "bottom_top" ; assign
xlat_level!1 = "south_north" ; the
xlat_level!2 = "west_east" ; dimensions
do i = 0, n_level-1
xlat_level(i,:,:) = (/xlat_new/) ; copy the data in xlat_new to xlat_level
end do
g_level = ge * (1 + 0.0052885*(sin(xlat_level*dtr))^2 - (0.0000059*(sin(2*xlat_level*dtr)^2))) - (0.003086*z)/1000 ; in ms-2
g_level!0 = "bottom_top" ; name
g_level!1 = "south_north" ; the
g_level!2 = "west_east" ; dimensions
;--------------------------------------------------------------------------------------------
; pressure difference at half pressure levels
delta_p = new((/n_level,x_dim,y_dim/),float)
do i = 0, n_level-2
delta_p(i,:,:) = p_total(i,:,:) - p_total(i+1,:,:) ; n_level-2 ignores the last value(unpaired), difference with 0, pressure difference in Pa
end do
delta_p!0 = "bottom_top" ; name
delta_p!1 = "south_north" ; the
delta_p!2 = "west_east" ; dimensions
; define constants
k1 = 77.6 ; refraction constant in K/hPa ; refraction constants after Boudouris, Newell and Baird
k2 = 70.4 ; refraction constant in K/hPa ;refs. R\FCeger.pdf, Refraction constants - best average.pdf
k3 = 373900 ; in K\B2/hPa
Mw = 18.0152*10^-3 ; molar mass of wet air in kg/mol (kg/kmol converted to kg/mol)
Md = 28.9644*10^-3 ; molar mass of dry air in kg/mol (kg/kmol converted to kg/mol)
e = Mw/Md
R = 8.31434 ; universal gas constant in J mol-1 K-1,
Rd = 287.054 ; gas constant of dry air in J kg-1 K-1
; level variables, reference to the write-up clarifies this step
pg_levels = new((/n_level,x_dim,y_dim/),float) ; necessary for ZHD in model levels ; pressure and gravity variables
pgt_levels = new((/n_level,x_dim,y_dim/),float) ; necessary for ZWD in model levels ; pressure, gravity, and temperature variables
pgr_levels = new((/n_level,x_dim,y_dim/),float) ; necessary for IWV in model levels ; pressure, gravity, and water vapor mixing ratio variables
do i = 0, n_level-1
pg_levels(i,:,:) = delta_p(i,:,:)*0.01/g_level(i,:,:)
pgt_levels(i,:,:) = (1/g_level(i,:,:)) * (k2 - k1*e + (k3/tk(i,:,:))) * (qvapor(i,:,:)/(qvapor(i,:,:)+1)) * delta_p(i,:,:)*0.01
pgr_levels(i,:,:) = (1/g_level(i,:,:)) * (qvapor(i,:,:)/(qvapor(i,:,:)+1)) * delta_p(i,:,:)
end do
; name the dimensions
pg_levels!0 = "bottom_top"
pg_levels!1 = "south_north"
pg_levels!2 = "west_east"
pgt_levels!0 = "bottom_top"
pgt_levels!1 = "south_north"
pgt_levels!2 = "west_east"
pgr_levels!0 = "bottom_top"
pgr_levels!1 = "south_north"
pgr_levels!2 = "west_east"
; sigma, replicated at all model levels
pg_levels_sum = new((/n_level,x_dim,y_dim/),float)
pgt_levels_sum = new((/n_level,x_dim,y_dim/),float)
pgr_levels_sum = new((/n_level,x_dim,y_dim/),float)
do j = 0, x_dim-1
do k = 0,y_dim-1
pg_levels_sum(:,j,k) = sum((pg_levels(:,j,k)))
pgt_levels_sum(:,j,k) = sum((pgt_levels(:,j,k)))
pgr_levels_sum(:,j,k) = sum((pgr_levels(:,j,k)))
end do
end do
; ; name the dimensions
pg_levels_sum!0 = "bottom_top"
pg_levels_sum!1 = "south_north"
pg_levels_sum!2 = "west_east"
pgt_levels_sum!0 = "bottom_top"
pgt_levels_sum!1 = "south_north"
pgt_levels_sum!2 = "west_east"
pgr_levels_sum!0 = "bottom_top"
pgr_levels_sum!1 = "south_north"
pgr_levels_sum!2 = "west_east"
; 2D sigma obtained by taking on the data on a single level (same in all levels), in this case, level index 0.
pg_sum = (/pg_levels_sum(0,0:193,0:241)/) ;copy only the data
pgt_sum = (/pgt_levels_sum(0,0:193,0:241)/) ;copy only the data
pgr_sum = (/pgr_levels_sum(0,0:193,0:241)/) ;copy only the data
;--------------------------------------------------------------------------------------------
; calculate zenith hydrostatic delay in m
ZHD_lower = 10^-6 * k1 * Rd * (psfc - p_total(0,:,:))*0.01/g_sfc ; below model orography
ZHD_levels = 10^-6 * k1 * Rd * pg_sum ; within model orography
ZHD_upper = 10^-6 * k1 * Rd * p_total(43,:,:)*0.01/g_level(43,:,:) ; above model orography
ZHD_NWP = ZHD_lower + ZHD_levels + ZHD_upper ; total ZHD
ZHD_NWP@units = "m" ; assign units to ZHD
;--------------------------------------------------------------------------------------------
; calculate zenith wet delay in m
ZWD_lower = 10^-6 * (R/Mw) * (1/g_sfc) * (k2 - k1*e + (k3/t2)) * (q2/(q2+1)) * (psfc - p_total(0,:,:))*0.01
ZWD_levels = 10^-6 * (R/Mw) * pgt_sum
ZWD_upper = 10^-6 * (R/Mw) * (1/g_level(43,:,:)) * (k2 - k1 + (k3/tk(43,:,:))) * (qvapor(43,:,:)/(qvapor(43,:,:)+1)) * p_total(43,:,:)*0.01
ZWD_NWP = ZWD_lower + ZWD_levels + ZWD_upper
ZWD_NWP@units = "m" ; assign units to ZWD
;--------------------------------------------------------------------------------------------
; calculate total zenith path delay in m
ZTD_NWP = ZHD_NWP + ZWD_NWP
ZTD_NWP@units = "m" ; assign units to ZTD
;--------------------------------------------------------------------------------------------
; calculate integrated water vapor in km/m2
IWV_lower = (1/g_sfc) * (q2/(q2+1)) * (psfc - p_total(0,:,:))
IWV_levels = pgr_sum
IWV_upper = (1/g_level(43,:,:)) * (qvapor(43,:,:)/(qvapor(43,:,:)+1)) * p_total(43,:,:)
IWV_NWP = IWV_lower + IWV_levels + IWV_upper
IWV_NWP@units = "m" ; assign units to IWV
;############################################################################################
; set plotting options for ZHD and make plots
res@FieldTitle = "Zenith Hydrostatic Delay" ; overwrite the field title
contour_zhd = wrf_contour(nc_file,wkszhd,ZHD_NWP,res) ; create ZHD contour plot
plot_zhd = wrf_map_overlays(nc_file,wkszhd,(/contour_zhd/),pltres,mpres) ; plot the data over a map background
;--------------------------------------------------------------------------------------------
; set plotting options for ZWD and make plots
res@FieldTitle = "Zenith Wet Delay" ; overwrite the field title
contour_zwd = wrf_contour(nc_file,wkszwd,ZWD_NWP,res) ; create ZWD contour plot
plot_zwd = wrf_map_overlays(nc_file,wkszwd,(/contour_zwd/),pltres,mpres) ; plot the data over a map background
;--------------------------------------------------------------------------------------------
; set plotting options for ZTD and make plots
res@FieldTitle = "Zenith Total Delay" ; overwrite the field title
contour_ztd = wrf_contour(nc_file,wksztd,ZTD_NWP,res) ; create ZTD contour plot
plot_ztd = wrf_map_overlays(nc_file,wksztd,(/contour_ztd/),pltres,mpres) ; plot the data over a map background
;--------------------------------------------------------------------------------------------
; set plotting options for IWV and make plots
res@FieldTitle = "Integrated Water Vapor" ; overwrite the field title
contour_iwv = wrf_contour(nc_file,wksiwv,IWV_NWP,res) ; create IWV contour plot
plot_iwv = wrf_map_overlays(nc_file,wksiwv,(/contour_iwv/),pltres,mpres) ; plot the data over a map background
;*******************************************************************************************
end do ; done with all files
;*******************************************************************************************
; plot terrain height over the entire domain using the geopotential height variable at full-levels obtained
; from the last netcdf file opened through the loop above
; full-level geopotential corresponding to level index 0 is equivalent to the terrain height
terrain_height = (/z_geop(0,0:190,0:240)/) ;copy only the data
terrain_height@units = "m" ; assign units to terrain height
; set plotting options for the terrain height and make a plot
res@FieldTitle = "Terrain Height" ; overwrite the field title
contour_tht = wrf_contour(nc_file,wkstht,terrain_height,res) ; create terrain height contour plot
plot_tht = wrf_map_overlays(nc_file,wkstht,(/contour_tht/),pltres,mpres) ; plot the data over a map background
;##########################################################################################
end ; exit program
But when i run it, it fails to plot. Here is the message received from the terminal;
Copyright (C) 1995-2015 - All Rights Reserved
University Corporation for Atmospheric Research
NCAR Command Language Version 6.3.0
The use of this software is governed by a License Agreement.
See http://www.ncl.ucar.edu/ for more details.
fatal:Subscript out of range, error in subscript #1
fatal:An error occurred reading xlat
fatal:["Execute.c":8575]:Execute: Error occurred at or near line 212 in file 7_tropodelay_iwv_and_terrainht_plots.ncl
I kindly request for help since i am stuck a this stage. I also want to have these plots seasonally like a plot DJF, MAM, JJA, SON
I will be grateful to have a solution.
Thank you very much.
NCL script to plot zenith hydrostatic delay, zenith wet delay, total zenith path delay and ; integrated water vapor over
Who is online
Users browsing this forum: No registered users and 1 guest