Tm plot and extract

All other topics about postprocessing model data (GrADS and other software), about other numerical weather prediction software (including WRF-NMM and WRF-ARW discussion unrelated to UEMS/WRF EMS), and general meteorology talk go in this forum.
Post Reply
cliffe
Posts: 17
Joined: Fri Sep 08, 2017 9:44 am

Tm plot and extract

Post by cliffe » Wed Dec 20, 2017 8:12 pm

Hi all,
I kindly request for your help. I am trying to plot and extract weighted mean temp (Tm) and surface temp (Ts) from my files extracted using uem. Below is the script i am trying to run:

;
; NCL script to plot weighted mean temperature (num_int and bevis) over the entire domain
; It also extracts surface temperature and together with Tm_num_int and Tm_bevis saves results into a data file (for comparison)
; Note that the script uses the beforeDA output files only

;############################################################################################

; 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_Tm1_num_int.pdf")
system("/bin/rm -f " + "plot_Tm2_bevis.pdf")
system("/bin/rm -f " + "Tm_and_Ts.dat")

;-------------------------------------------------------------------------------------------

; open workstations, output to pdf

wkstm1 = gsn_open_wks ("pdf","plot_Tm1_num_int")
wkstm2 = gsn_open_wks ("pdf","plot_Tm2_bevis")

;-------------------------------------------------------------------------------------------

; 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 = 45.857
mpres@mpMaxLatF = 54.5656
mpres@mpMinLonF = 5.36002
mpres@mpMaxLonF = 24.378

; 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

;##########################################################################################

; create list of the wrfout beforeDA netcdf files

fList = systemfunc("ls /usr/nwp/wrfprd/*")
Files = dimsizes(fList)

; start looping through the *.nc files
do iFile = 0, Files-1
nc_file = addfile(fList(iFile),"r")
time = 0

;*******************************************************************************************

; 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

pgr1_levels = new((/n_level,x_dim,y_dim/),float) ; pressure, gravity, and water vapor mixing ratio variables
pgr2_levels = new((/n_level,x_dim,y_dim/),float) ; pressure, gravity, and water vapor mixing ratio variables

do i = 0, n_level-1
pgr1_levels(i,:,:) = (1/g_level(i,:,:)) * (1/tk(i,:,:)) * (qvapor(i,:,:)/(qvapor(i,:,:)+1)) * (delta_p(i,:,:)*0.01)
pgr2_levels(i,:,:) = (1/g_level(i,:,:)) * (1/(tk(i,:,:)*tk(i,:,:))) * (qvapor(i,:,:)/(qvapor(i,:,:)+1)) * (delta_p(i,:,:)*0.01)

end do

; name the dimensions

pgr1_levels!0 = "bottom_top"
pgr1_levels!1 = "south_north"
pgr1_levels!2 = "west_east"

; sigma, replicated at all model levels
pgr1_levels_sum = new((/n_level,x_dim,y_dim/),float)
do j = 0, x_dim-1
do k = 0,y_dim-1
pgr1_levels_sum(:,j,k) = sum((pgr1_levels(:,j,k)))
end do
end do

; name the dimensions

pgr1_levels_sum!0 = "bottom_top"
pgr1_levels_sum!1 = "south_north"
pgr1_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.

pgr1_sum = (/pgr1_levels_sum(0,0:156,0:206)/) ;copy only the data


pgr2_levels!0 = "bottom_top"
pgr2_levels!1 = "south_north"
pgr2_levels!2 = "west_east"

; sigma, replicated at all model levels
pgr2_levels_sum = new((/n_level,x_dim,y_dim/),float)
do j = 0, x_dim-1
do k = 0,y_dim-1
pgr2_levels_sum(:,j,k) = sum((pgr2_levels(:,j,k)))
end do
end do

; name the dimensions

pgr2_levels_sum!0 = "bottom_top"
pgr2_levels_sum!1 = "south_north"
pgr2_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.

pgr2_sum = (/pgr2_levels_sum(0,0:156,0:206)/) ;copy only the data

;--------------------------------------------------------------------------------------------

; calculate mean temperature in Kelvin (using numerical integration)

Tm_lower1 = ((1/g_sfc) * (1/t2) * (q2/(q2+1)) * ((psfc - p_total(0,:,:))*0.01))
Tm_lower2 = ((1/g_sfc) * (1/(t2*t2)) * (q2/(q2+1)) * ((psfc - p_total(0,:,:))*0.01))

Tm_levels1 = pgr1_sum
Tm_levels2 = pgr2_sum

Tm_upper1 = ((1/g_level(43,:,:)) * (1/tk(43,:,:)) * (qvapor(43,:,:)/(qvapor(43,:,:)+1)) * (p_total(43,:,:)*0.01))
Tm_upper2 =((1/g_level(43,:,:)) * (1/(tk(43,:,:)*tk(43,:,:))) * (qvapor(43,:,:)/(qvapor(43,:,:)+1)) * (p_total(43,:,:)*0.01))

Tm_num_int = (Tm_lower1 + Tm_levels1 + Tm_upper1)/(Tm_lower2 + Tm_levels2 + Tm_upper2)

Tm_num_int = Tm_num_int - 273.15 ; convert Kelvin to deg Celsius

;------------------------------------------------------------------------------------------
; set plotting options and make plots

res@FieldTitle = "Mean temperature (num_int) [deg C]" ; overwrite the field title

contour_Tm = wrf_contour(nc_file,wkstm1,Tm_num_int,res) ; create Tm contour plot

plot_Tm = wrf_map_overlays(nc_file,wkstm1,(/contour_Tm/),pltres,mpres) ; plot the data over a map background

delete([/contour_Tm,plot_Tm/])

;------------------------------------------------------------------------------------------

; calculate mean temperature in Kelvin (using Bevis et al. formular). Note that the surface temp is required to be in Kelvin

Tm_bevis = 70.2 + (0.72 * t2) ; in units of Kelvin


Tm_bevis = Tm_bevis - 273.15 ; convert deg Kelvin to deg Celsius

;------------------------------------------------------------------------------------------

; plotting options

res@FieldTitle = "Mean temperature (Bevis Tm-Ts relation)[deg C]" ; overwrite the field title

contour_Tm = wrf_contour(nc_file,wkstm2,Tm_bevis,res) ; create contour plot

plot_Tm = wrf_map_overlays(nc_file,wkstm2,(/contour_Tm/),pltres,mpres) ; plot the data over a map background

delete([/contour_Tm,plot_Tm/])

;------------------------------------------------------------------------------------------

;*******************************************************************************************

t2 = t2 - 273.15 ; convert deg Kelvin to deg Celsius

results = [/Tm_num_int,Tm_bevis,t2/]

write_table("Tm_and_Ts.dat","a",results,"%4.2f %4.2f %4.2f")


;*******************************************************************************************

end do ; done with all files

;#######################################################################################################

end ; exit program




But when i run it, it fails to plot and extract the Ts. 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 200 in file 11_Tm_plot.ncl

I kindly request for help to solve this problem. In my thinking, i believe the cause of the problem is from these two statements from the script;

x_dim = 193
y_dim = 241

I also request to know what do these mean and where to get these two values i.e. the x_dim and y_dim.
I kindly request for help since i am stuck a this stage.
I will be grateful to have a solution as soon as possible.
Thank you very much.
Richard Cliffe

Post Reply

Who is online

Users browsing this forum: Baidu [Spider] and 13 guests