Hi,
Thanks for all. I’m interested to Ri, so for this I made some assumption to
calculating it. The code is working but the Ri contours are not appeared in
the figure, may be my code is not correct. Please could anyone help me to
figure out the problem? My code and the result are attached.
thanks
wei
= = = = = = =
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"
begin
a = addfile("wrfout_d03_2000-03-20.nc","r")
T = a->T
T = T + 300. ;theta K
u = wrf_user_getvar(a,"ua",-1)
v = wrf_user_getvar(a,"va",-1)
z = wrf_user_getvar(a, "z",-1)
rh = wrf_user_getvar(a,"rh",-1)
U = (u*u+v*v)^(0.5) ; m/sec
copy_VarCoords(u, U )
printVarSummary(U)
printMinMax(U, True)
;
dTdz = center_finite_diff_n(T,z,False,0,1)
eps = 1e-02 ; eliminate near zero values
dTdz@_FillValue = -999.
dTdz = where(dTdz.lt.eps, dTdz@_FillValue, dTdz)
copy_VarCoords(T, dTdz )
printVarSummary(dTdz)
printMinMax(dTdz, True)
;
dUdz = center_finite_diff_n(U,z,False,0,1)
eps = 1e-02 ; eliminate near zero values
dUdz@_FillValue = -999.
dUdz = where(dUdz.lt.eps, dUdz@_FillValue, dUdz)
copy_VarCoords(u, dUdz )
printVarSummary(dUdz)
printMinMax(dUdz, True)
;
nsqr = ((9.8/T)*(dTdz)); Brunt vaisala frequency
copy_VarCoords(T, nsqr)
printVarSummary(nsqr)
printMinMax(nsqr, True)
RI = nsqr/(dUdz)^2 ; to calculate Richardson Number(RI)
copy_VarCoords(u, RI )
printVarSummary(RI)
printMinMax(RI, True)
;***********************************
wks = gsn_open_wks("pdf","ri")
gsn_define_colormap(wks,"gui_default")
res = True ; Set up some basic plot resources
res@Footer = False
res@NoHeaderFooter = True
pltres = True
;*******************************
times = wrf_user_list_times (a); get times in the file
ntimes = dimsizes(times) ; number of times in the file
FirstTime = res
mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
nd = dimsizes(mdims)
;do it = 10,10;ntimes-1
;***************************************
if ( FirstTime ) then ; get height info for labels
zmin = 0.
zmax = max(z)/1000.
nz = floattoint(zmax/2 + 1)
FirstTime = False
end if
;**************************************
ip = 1 ; Just do the one (constant y coord) plot
opts = True ; setting start and end times
plane = new(4,float)
if(ip .eq. 1) then
plane = (/ 0,84, 200,84 /) ; start x;y & end x;y point
end if
;
rh_plane = wrf_user_intrp3d(rh,z,"v",plane,0.,opts)
ri_plane = wrf_user_intrp3d(RI,z,"v",plane,0.,opts)
printVarSummary(rh_plane)
dim = dimsizes(rh_plane) ; Find the data span - for use in labels
zspan = dim(0)
; Options for XY Plots
opts_xy = res
opts_xy@tiYAxisString = "Height (km)"
opts_xy@AspectRatio = 0.75
opts_xy@cnMissingValPerimOn = True
opts_xy@cnMissingValFillColor = 0
opts_xy@cnMissingValFillPattern = 11
opts_xy@tmYLMode = "Explicit"
opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks
opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels
opts_xy@tiXAxisFontHeightF = 0.020
opts_xy@tiYAxisFontHeightF = 0.020
opts_xy@tmXBMajorLengthF = 0.02
opts_xy@tmYLMajorLengthF = 0.02
opts_xy@tmYLLabelFontHeightF = 0.015
opts_xy@PlotOrientation = ri_plane@Orientation
; Plotting options for RH
opts_rh = opts_xy
opts_rh@pmLabelBarOrthogonalPosF = -0.07
opts_rh@ContourParameters = (/ 10., 90., 10. /)
opts_rh@cnFillOn = True
; Plotting options for ri
opts_ri = opts_xy
opts_ri@cnInfoLabelOrthogonalPosF = 0.00
opts_ri@ContourParameters = (/10. /)
; Get the contour info for the rh and ri
contour_ri = wrf_contour(a,wks,ri_plane(16,:,:),opts_ri)
contour_rh = wrf_contour(a,wks,rh_plane(16,:,:),opts_rh)
; MAKE PLOTS
plot = wrf_overlays(a,wks,(/contour_rh,contour_ri/),pltres)
; Delete options and fields, so we don't have carry over
delete(opts_ri)
delete(opts_rh)
delete(ri_plane)
delete(rh_plane)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2.1.8 : Tue Feb 15 2011 - 09:43:19 MST