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" begin a = addfile("Southwest_ATM.1991070100.time_brief.nc","r") xlat = a->xlat xlon = a->xlon q = a->qv ; kg per kg ; qv(time, kz, iy, jx) printVarSummary(q) ; printMinMax(q, True) ps = a->ps ; hPa ; ps(time, iy, jx) printVarSummary(ps) printMinMax(ps, True) ptop = a->ptop ; hPa ; (1) ... scalar print(ptop) sigma = a->sigma ; sigma(kz) print(sigma) ; ------------------------------------------------- ; create "interface" sigma levels to mimic the ; CAM models interface hybrid levels ksig = dimsizes(sigma) ; 'super' sigma .. include top and bottom SIGMA = new(ksig+2, typeof(sigma),"No_FillValue") SIGMA(0) = 0.0 ; anything less than sigma(0) SIGMA(1:ksig) = (/ sigma /) SIGMA(ksig+1) = 1.0 print(SIGMA) ; interface sigma levels sigmai = (SIGMA(1:) + SIGMA(0:ksig))*0.5 sigmai(0) = SIGMA(0) ; force top for full thickness sigmai(ksig) = 1.0 ; force bot print(sigmai) ; ------------------------------------------------- ; use the sigmai in the dpres_hybrid_ccm function ; to get layer thicknesses A = sigmai A = 0.0 P0 = 0.0 dp_sigma = dpres_hybrid_ccm(ps,P0,A,sigmai) dp_sigma = dp_sigma*100. dp_sigma@units = "Pa" dp_sigma!0 = ps!0 dp_sigma!1 = "plevel" dp_sigma!2 = ps!1 dp_sigma!3 = ps!2 printVarSummary(dp_sigma) printMinMax(dp_sigma, True) pw = prcwater_dp(q(time|:,iy|:,jx|:,kz|:) \ ,dp_sigma(time|:,iy|:,jx|:,plevel|:)) copy_VarCoords(ps,pw) pw@long_name = "total column precipitable water" pw@units = "kg/m2" pw@coordinates = "xlat xlon" pw@grid_mapping = "rcm_map" printVarSummary(pw) printMinMax(pw, True) ; --------------------------------------------------- g = 9.81 PW = dim_sum_n(q*dp_sigma, 1)/g ; v5.1.1 onward copy_VarCoords(ps, PW) PW@long_name = "total column precipitable water" PW@units = "kg/m2" PW@coordinates = "xlat xlon" PW@grid_mapping = "rcm_map" printVarSummary(PW) printMinMax(PW, True) exit out=addfile("pw.nc","c") out->pw=pw out->xlat=xlat out->xlon=xlon end