NCL script to plot zenith hydrostatic delay, zenith wet delay, total zenith path delay and ; integrated water vapor over

[cliffe]

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.