NCL Home >
Documentation >
Functions >
CESM
pres_hybrid_jra55
Calculates the "full" hybrid pressure levels for models using the Simmons/Burridge formulations (eg: Japanese ReAnalysis).
Available in version 6.4.0 and later.
Prototype
function pres_hybrid_jra55 ( ps : numeric, hyai [*] : numeric, hybi [*] : numeric ) return_val : numeric
Arguments
psAn array of at least 2 dimensions equal to surface pressure (Pa) data. The rightmost dimensions must be latitude and longitude.
hyaiA one-dimensional array containing the hybrid A interface coefficients. The levels are expected to be ordered from bottom to top.
hybiA one-dimensional array containing the hybrid B interface coefficients. The levels are expected to be ordered from bottom to top.
Return value
If ps is two-dimensional [e.g. (nlat,mlon)] or three-dimensional [e.g. (ntim,nlat,mlon)] then the return array will have an additional level dimension of size sixty (60,nlat,mlon) or (ntim,60,nlat,mlon). The returned type will be double if ps is double, float otherwise.
Description
This function is being tested. Contact ncl-talk@ucar.edu if you wish to try it.
Calculates full pressure levels using the formula found on page 17 of the "JRA-55 Product User's Handbook". Also, see the Simmons and Burridge (1981) reference below,
Nominally, the formula is:
p(k) = exp( (1/dp(k))*(p(k-0.5)*log(p(k-0.5)) - p(k+0.5)*log(p(k+0.5))-1.0))
where the p(k-05) and p(k+0.5) represent the interface coefficients. The levels are expected to be ordered from bottom to top.
The original reference for the above equation is:
Simmons. A.J., and D. M. Burridge:
An Energy and Angular-Momentum Conserving Vertical Finite-Difference Scheme
and Hybrid Vertical Coordinates
MON WEATHER REV , vol. 109, no. 4, 1981
http://dx.doi.org/10.1175/1520-0493(1981)109<0758:AEAAMC>2.0.CO;2
Examples
For JRA55, hyai and hybi are one-dimensional arrays ordered bottom-to-top and ps is a two- or three-dimensional array of size (nlat,mlon) or (ntim,nlat,mlon) in units of pascals. The full pressure levels pm) will be returned as a three-dimensional array of size (60,nlat,nlon) or a four-dimensional array of size (ntim,60,nlat,nlon).
Example 1 The simplest case is illustrated.
nlat = 1 ; function expects (..,nlat,mlon) structure mlon = 1 ps = new((/nlat,mlon/), "double") ps = 100000. ; ps(1,1): surface pressure [Pa] dirhy = "./" filhy = "JRA55.hybrid_levels.nc" ; contains 61 hybrid interface coef finhy = addfile(dirhy+filhy,"r") hyai = finhy->hyai hybi = finhy->hybi p = pres_hybrid_jra55(ps,hyai,hybi) ; [60] x [1] x [1] ; (klev,nlat,mlon) ; verification when ps=100000 lev = finhy->lev diff = p(:,0,0)-lev print(lev+" "+p(:,0,0)+" "+diff)The edited output looks like
lev p diff
(0) 99849.96244364697 99849.96244364697 0
(1) 99549.96233052592 99549.96233052592 0
(2) 99149.89494018865 99149.89494018865 0
(3) 98549.7928287735 98549.7928287735 0
(4) 97699.57352197972 97699.57352197972 0
(5) 96599.37887679219 96599.37887679219 0
(6) 95299.14304700246 95299.14304700246 0
(7) 93698.55920278581 93698.55920278979 3.987224772572517e-09
(8) 91798.18441228721 91798.18441228606 -1.149601303040981e-09
(9) 89697.75169188292 89697.75169188037 -2.546585164964199e-09
(10) 87347.01856643189 87347.01856643189 0
(11) 84696.14302236082 84696.14302236518 4.365574568510056e-09
(12) 81795.41531387174 81795.41531386811 -3.623426891863346e-09
(13) 78694.5780828073 78694.5780828073 0
(14) 75443.98546020423 75443.98546020543 1.193257048726082e-09
(15) 72042.91489346673 72042.91489346455 -2.182787284255028e-09
(16) 68441.66534734554 68441.66534734869 3.157765604555607e-09
(17) 64741.1889174707 64741.18891747024 -4.583853296935558e-10
(18) 60990.13453779992 60990.13453779992 0
(19) 57189.47883614307 57189.47883614226 -8.076312951743603e-10
(20) 53388.72974167685 53388.72974167523 -1.615262590348721e-09
(21) 49587.86576445695 49587.86576445792 9.749783203005791e-10
(22) 45837.55467777611 45837.55467777644 3.274180926382542e-10
(23) 42187.19874312158 42187.19874312113 -4.511093720793724e-10
(24) 38636.78800522001 38636.78800522021 2.037268131971359e-10
(25) 35186.30936535118 35186.30936535118 0
(26) 31886.61757076395 31886.61757076395 0
(27) 28736.06368117695 28736.06368117695 0
(28) 25736.38222305299 25736.38222305299 0
(29) 22885.72350887325 22885.72350887309 -1.637090463191271e-10
(30) 20186.04356927362 20186.04356927405 4.292814992368221e-10
(31) 17686.42714647552 17686.42714647502 -5.020410753786564e-10
(32) 15386.89025254726 15386.8902525474 1.364242052659392e-10
(33) 13287.45328330309 13287.45328330316 7.09405867382884e-11
(34) 11388.14185796683 11388.14185796685 2.000888343900442e-11
(35) 9688.987180267748 9688.987180267834 8.549250196665525e-11
(36) 8164.265567858072 8164.265567858072 0
(37) 6763.767298472668 6763.767298472656 -1.182343112304807e-11
(38) 5515.002875033884 5515.002875033893 9.094947017729282e-12
(39) 4466.576216266381 4466.576216266381 0
(40) 3608.149325172584 3608.14932517259 6.366462912410498e-12
(41) 2914.501376923828 2914.501376923823 -5.002220859751105e-12
(42) 2352.981404124263 2352.981404124263 0
(43) 1898.89882407333 1898.89882407333 0
(44) 1532.036534874198 1532.036534874198 0
(45) 1235.129618271835 1235.129618271835 0
(46) 997.1208628290149 997.1208628290149 0
(47) 804.9498872016828 804.9498872016828 0
(48) 649.2593496916936 649.2593496916936 0
(49) 524.0019775689651 524.0019775689651 0
(50) 422.1641249658988 422.1641249658988 0
(51) 338.288125071407 338.288125071407 0
(52) 268.363593649281 268.363593649281 0
(53) 208.8961054950331 208.8961054950331 0
(54) 158.4429547310325 158.4429547310325 0
(55) 115.9526640886324 115.9526640886324 0
(56) 80.46871791830324 80.46871791830324 0
(57) 51.45078427193928 51.45078427193928 0
(58) 28.97808727180794 28.97808727180794 0
(59) 10 10 0
Example 2 Read from JRA55 GRIB file and interpolate to user specified isobaric levels. Also, create a netCDF file.
;*****************************************************************
; Two local utility routines: created for convenience
;*****************************************************************
undef ("changeDimNames")
procedure changeDimNames(v3:numeric, v4:numeric)
; convenience: change fron NCL generated names to more 'common' (friendly) names
begin
v3!0 = "time"
v3!1 = "lat"
v3!2 = "lon"
v4!0 = "time"
v4!1 = "lev"
v4!2 = "lat"
v4!3 = "lon"
end
;--
undef ("deleteSelectedAtts")
procedure deleteSelectedAtts(v4:numeric)
; esthetics: these pertain more to the original GRIB file; not to pressure data
begin
delete([/ v4@sub_center, v4@level_indicator, v4@parameter_table_version \
, v4@parameter_table_version, v4@parameter_number, v4@gds_grid_type /])
end
; ---------------------------> MAIN DRIVER <----------------------
;*****************************************************************
; User specified levels
;*****************************************************************
plvl = (/10.,20.,30.,50.,70.,100.,150.,200.,250.,300.,400.,500.,700.,850.,925.,1000./)*100
plvl = plvl(::-1) ; reorder ... make same vertical order as JRA55
plvl!0 = "plvl"
plvl@units = "Pa"
plvl@long_name = "pressure level"
plvl&plvl = plvl
printVarSummary(plvl)
nplvl = dimsizes(plvl)
netCDF = True
vari = "TMP_GDS4_HYBL"
varo = "TMP"
linlog = -1 ; type vertical interpolation (see: int2p_n)
dirx = "./"
filx = systemfunc("cd "+dirx+" ; ls fcst_mdl*")
nfilx = dimsizes(filx)
print(filx)
dirps = "./"
filps = systemfunc("cd "+dirps+" ; ls jra55PS6hr*")
nfilps = dimsizes(filx)
print(filps)
;*****************************************************************
; end user input
;*****************************************************************
; input error checks
;*****************************************************************
if (nfilx.ne.nfilps) then
print("number of files do not match: nfilx="+nfilx+" nfilps="+nfilps)
exit
end if
;*****************************************************************
; Read hybrid information; external file containg hyai and hybi coefficients
;*****************************************************************
dirhy = "./" ; read file with interface hybrid info
filhy = "JRA55.hybrid_levels.nc"
finhy = addfile(dirhy+filhy,"r")
hyai = finhy->hyai
hybi = finhy->hybi
P0 = finhy->P0
ilev = finhy->ilev
lev = finhy->lev
;*****************************************************************
; loop over all GRIB files
;*****************************************************************
fext = ".grb" ; all the input are GRIB
do nf=0,nfilx-1
loopTime = get_cpu_time()
print("-------------------------------------------------")
print("nf="+nf+" dirx+filx(nf)+fext="+dirx+filx(nf)+fext)
finx = addfile(dirx+filx(nf)+fext,"r")
;;x := finx->$vari$(0:4,:,:,:) ; x( 5, lv_HYBL1, g4_lat_2, g4_lon_3 ) ; testing
x := finx->$vari$ ; x( initial_time0_hours, lv_HYBL1, g4_lat_2, g4_lon_3 )
printVarSummary(x)
finps = addfile(dirps+filps(nf)+fext,"r")
;;ps := finps->PRES_GDS4_SFC(0:4,:,:) ; ( 5 , g4_lat_2, g4_lon_3 )
ps := finps->PRES_GDS4_SFC ; x( initial_time0_hours, g4_lat_2, g4_lon_3 )
printVarSummary(ps)
;*****************************************************************
; error checking
;*****************************************************************
dimx = dimsizes(x)
rankx = dimsizes(dimx)
if (rankx.ne.4)then
print("FATAL: rank of variable must be 4: rankx="+rankx)
exit
end if
dimps = dimsizes(ps)
rankps= dimsizes(dimps)
if (rankps.ne.3)then
print("FATAL: rank of PS must be 3: rankps="+rankps)
exit
end if
ntimx = dimx(0)
ntimps= dimps(0)
if (ntimx.ne.ntimps) then
print("FATAL: number of time steps do not match: ntimx="+ntimx+" ntimps="+ntimps)
exit
end if
if (.not.all(x&initial_time0_hours .eq. ps&initial_time0_hours)) then
print("FATAL: time mismatch between variable and sfc pressure")
exit
end if
;*****************************************************************
; Miscellaneous
;*****************************************************************
changeDimNames(ps, x) ; convenience ... user friendly names
xMin = min(x) ; source min
xMax = max(x) ; max
ntim = dimx(0)
klev = dimx(1)
nlat = dimx(2)
mlon = dimx(3)
klevi= klev+1 ; number of interface levels
;*****************************************************************
; derive interface pressure levels at each grid point
;*****************************************************************
p = pres_hybrid_jra55(ps,hyai,hybi) ; (ntim,60,nlat,mlon)
;*****************************************************************
; Interpolate from hybrid pressure levels to isobaric levels ('plvl' )
;*****************************************************************
x_p := int2p_n_Wrap(p,x,plvl,linlog,1) ; (ntim,nplvl,nlat,mlon)
deleteSelectedAtts(x_p)
x_p@info = "interpolated to pressure levels via NCL script"
x_p = where(x_p.lt.xMin .or. x_p.gt.xMax, x_p@_FillValue, x_p) ; extrapolated values
delete(p)
;*****************************************************************
; Write netCDF
;*****************************************************************
if (netCDF) then
diro = "./"
filo = "jra55_2_so."+filx(nf)+".nc" ; arbitrary
ptho = diro+filo
system("/bin/rm -f "+ ptho)
ncdf = addfile(ptho,"c")
ncdf->$varo$ = x_p
end if
print("LOOP time: " + (get_cpu_time() - loopTime))
end do