# Converting hybrid levels to pressures

4 posts / 0 new
rentianyang@...
Converting hybrid levels to pressures

Hi,

I've obtained CESM output data in hybrid levels, and wish to convert to pressure levels. I know that NCL can do the job, but I don't really understand how to use it. Do I need to write a script in order to do the conversion or there is some options to NCL (e.g. ncl ifile ofile) such that the conversion can be done quickly? Please help me. Thank you so much.

Regards
Tianyang

mai

Here is an example of what you want. Once you understand every line, you can modify it for your case.

begin

; ..Initializations
levpc = (/300.,850./)                     ; specify pressure levels for PSI and CHI
levt  = (/850./)                          ; specify pressure levels for T
levh  = (/500./)                          ; specify pressure levels for Z3
klvpc = dimsizes(levpc)                   ; number of pressure levels for PSI and CHI
klvt  = dimsizes(levt)                    ; number of pressure levels for T
klvh  = dimsizes(levh)                    ; number of pressure levels for Z3
nyrs  = 62
nflv  = 6
locP  = "b.e11.B1850C5CN.f09_g16.011"
locF1 = ".cam.h0."
locFe = ".nc"

; ..Allocate space for needed arrays based on first history tape
loc   = locP+locF1+"2000-01"+locFe
hyam  = f->hyam                           ; hybrid coef (used below)
hybm  = f->hybm                           ;
t     = f->T                              ; temperature
dimc  = dimsizes(t)                       ; get dimension sizes for multilevel fields
ntim  = dimc(0)
nlev  = dimc(1)
nlat  = dimc(2)
nlon  = dimc(3)

; ..Get PHIS and P0 (needed for vertical interpolation to pressure surfaces)
phis  = f->PHIS(0,:,:)
intyp = 1                                 ; 1=linear, 2=log, 3=log-log
kxtrp = True                              ; True=extrapolate
p0    = f->P0                             ; units of P0 is Pa
p0mb  = p0/100.                           ; need ref sur pressure P0 in mb [mb=hPa]
delete(t)
delete(f)

; ..Step through the input files
do imo = 01,12
do iyr = 2000,2061
loc  = locP+locF1+locF2+sprinti("%0.4i",iyr)+"-"+sprinti("%0.2i",imo)+locFe
psfc = f->PS                        ; surface pressure [Pa]
t    = f->T                         ; temperature
u    = f->U                         ; zonal wind
v    = f->V                         ; meridional wind
z    = f->Z3                        ; geopotential height
psl  = f->PSL                       ; sea level pressure
ts   = f->TS                        ; surface temperature
tbot = t(:,nlev-1,:,:)              ; bottom T level [for varflg /= 0]

; ..Calculate Z3,T,PSI,CHI on desired pressure surfaces
varflg = -1                         ; = -1 for geopotential
zp = vinth2p_ecmwf(z,hyam,hybm,levh,psfc,intyp,p0mb,1,kxtrp,\
varflg,tbot,phis)
varflg = 1                          ; =  1 for temperature
tp = vinth2p_ecmwf(t,hyam,hybm,levt,psfc,intyp,p0mb,1,kxtrp,\
varflg,tbot,phis)
varflg = 0                          ; not T or Z3
up = vinth2p_ecmwf(u,hyam,hybm,levpc,psfc,intyp,p0mb,1,kxtrp,\
varflg,tbot,phis)
vp = vinth2p_ecmwf(v,hyam,hybm,levpc,psfc,intyp,p0mb,1,kxtrp,\
varflg,tbot,phis)
psip = up
chip = up
uv2sfvpg(up,vp,psip,chip)           ; U,V  ==> PSI,CHI

end do                                ; year
end do                                  ; month

print("Normal termination of exam.ncl")
end

I strongly suggest that you read the documentation for the NCL functions that you will be calling. For example:

http://www.ncl.ucar.edu/Document/Functions/Built-in/vinth2p_ecmwf.shtml

leilanepassos@...

Dear Mai,

I made the NCL script below for interpolate from hybrid sigma levels to pressure. But my intetion is to interpolate from hybrid sigma levels to 10m height. What NCL function I can to use to make this interpolation? I tried the function int2p, but didn't work.

; --------------  LOAD FUNCTIONS AND PROCEDURES ----------------

; --------------  BEGINING OF NCL SCRIPT ----------------

begin

outfile = "press_heig_teste.nc"

if (isfilepresent(outfile)) then
system("rm -rf "+outfile)          ;-- make sure that file does not exist
end if

;-- open data file

;-- get variable t and its dimensions and dimension sizes
;-- variables at midpoints, so hyai e hybi

hyai         = in_config->hyai      ; hybrid a coefficients
hybi         = in_config->hybi      ; hybrid b coefficients
PSFC      = in_ps->PS
P0mb      = 0.01*in_config->P0   ; passa Pa para milibar
PHIS       = in_phis->PHIS
PHIS       = lonFlip(PHIS)        ; convert from 0:360 to -180:180
phis_sub = PHIS(1,{-62:12},{-73:0})  ; sub-domain

U = in_u->U ;(time, lev, lat, lon)                       ;-- get variable u
V = in_v->V ;(time, lev, lat, lon)                       ;-- get variable v
Z = in_z->Z3 ;(time, lev, lat, lon)                      ;-- get variable u
T = in_t->T                          ; temperature at bottom hybrid level

tbot = T(:,0,:,:)               ; bot temp level [clarity]
nlev = dimsizes(hyai)

lev_p   = (/ 1000., 995., 990., 985., 980. /)
lev_p!0         = "lev_p"            ; variable and dimension name the same
lev_p&lev_p     = lev_p              ; create coordinate variable
lev_p@long_name = "pressure"         ; attach some attributes
lev_p@units     = "hPa"
lev_p@positive  = "down"

intyp = 1                             ; 1=linear, 2=log, 3=log-log
kxtrp = True                          ; True=extrapolate

varflg = 1                           ; temperature is variable
Tp     = vinth2p_ecmwf(T,hyai,hybi,lev_p,PSFC,intyp,P0mb, \
1,kxtrp,varflg,tbot,phis_sub)
varflg = -1                          ; geo pot hgt is variable [tbot is used]
Zp     = vinth2p_ecmwf(Z, hyai,hybi, lev_p ,PSFC, intyp, P0mb, \
1,kxtrp,varflg,tbot,phis_sub)
varflg = 0                           ; not T or Z
Up     = vinth2p_ecmwf(U,hyai,hybi,lev_p,PSFC,intyp,P0mb,\
1,kxtrp,varflg,tbot,phis_sub)

; Pressure to height
xlvl = (/10, 50, 100 /) ; height levels
xlvl@units = "m"

;  Reorder to interpolation

Uh = Up(time|:,lat|:,lon|:,lev_p|:)

linlog = -1 ; or 1 and negative for extrapolation
uHeight = int2p(lev_p \
,Uh(time|:,lat|:,lon|:,lev_p|:) \
,xlvl,linlog)

uHeight@long_name = U@long_name
uHeight@units = U@units

uHeight!0 = "time"
uHeight!1 = "lat"
uHeight!2 = "lon"
uHeight!3 = "height"
printVarSummary(uHeight)

uHeight&time = U&time
uHeight&lat  = U&lat
uHeight&lon  = U&lon
uHeight&height  = xlvl

height = xlvl
height!0 = "height"
printVarSummary(height)

; Reorder
uH = uHeight(time|:,height|:,lat|:,lon|:)
delete(uHeight)

time  =  in_u->time                    ;-- get dimension time
lat   =  in_u->lat                     ;-- get dimension lat
lon   =  in_u->lon                     ;-- get dimension lon

ntim  =  dimsizes(time)              ;-- get dimension sizes of time
nlat  =  dimsizes(lat)               ;-- get dimension sizes of lat
nlon  =  dimsizes(lon)               ;-- get dimension sizes of lon

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;-- create new netCDF file

;-- begin output file settings
setfileoption(fout,"DefineMode",True) ;-- explicitly declare file definition mode

;-- create global attributes of the file
fAtt                  =  True        ;-- assign file attributes
fAtt@title            = "NCL Efficient Approach to netCDF Creation"
fAtt@source_file      = "CCSM4.nc"
fAtt@Conventions      = "CF"
fAtt@creation_date    =  systemfunc ("date")
fAtt@history          =  "NCL script: ex_netcdf_2.ncl"
fAtt@comment          = "Convert variable t: Kelvin to Celsius"
fileattdef(fout,fAtt)                ;-- copy file attributes

;-- predefine the coordinate variables and their dimensionality
nhei  =  dimsizes(height)
dimNames = (/"time", "height", "lat", "lon"/)
dimSizes = (/ -1   , nhei, nlat,  nlon /)
dimUnlim = (/ True , False, False, False/)
filedimdef(fout,dimNames,dimSizes,dimUnlim)

;-- predefine the the dimensionality of the variables to be written out
filevardef(fout, "time" ,typeof(time),getvardims(time))
filevardef(fout, "height",typeof(height), getvardims(height))
filevardef(fout, "lat"  ,typeof(lat), getvardims(lat))
filevardef(fout, "lon"  ,typeof(lon), getvardims(lon))
filevardef(fout, "uH" ,typeof(uH),  getvardims(uH))
;  filevardef(fout, "v1000p" ,typeof(v1000p),  getvardims(v1000p))

;-- copy attributes associated with each variable to the file
filevarattdef(fout,"time" ,time)       ;-- copy time attributes
filevarattdef(fout,"height",height)       ;-- copy time attributes
filevarattdef(fout,"lat"  ,lat)        ;-- copy lat attributes
filevarattdef(fout,"lon"  ,lon)        ;-- copy lon attributes
filevarattdef(fout,"uH" ,uH)       ;-- copy u10m attributes
; filevarattdef(fout,"v1000p" ,v1000p)       ;-- copy v10m attributes

;-- explicitly exit file definition mode (not required)
setfileoption(fout,"DefineMode",False)

;-- output only the data values since the dimensionality and such have been predefined.
;-- The "(/", "/)" syntax tells NCL to only output the data values to the predefined
;-- locations on the file.
fout->time   =  (/time/)               ;-- write time to new netCDF file
fout->height =  (/height/)
fout->lat    =  (/lat/)                ;-- write lat to new netCDF file
fout->lon    =  (/lon/)                ;-- write lon to new netCDF file
fout->uH   =  (/uH/)               ;-- write variable to new netCDF file
; fout->v1000p   =  (/v1000p/)               ;-- write variable to new netCDF file

end

Thank you very much.

Leilane Passos

s.undorf@...

Hello, just to add to the reply above so it might be of use to a future user, here a combination of cdo and nco funtions also doing the task: If you have sea level pressure (PSL) and surface geopotential height (PHIS) in the same file as the coordinate, you do

ncap2 -s "hyai=hyai*100000;hyam=hyam*100000;" \$f \$f_tmp

cdo ml2pl,100000,92500,85000,70000,60000,50000,40000,30000,25000,20000,15000,10000,7000,5000,3000,2000,1000, -chname,PSL,sp -chname,PHIS,geosp \$f_tmp \$f_out

(see https://code.mpimet.mpg.de/boards/1/topics/2997 for why the first command is necessary). Having PSL and PHIS are not needed for all variables, see https://code.mpimet.mpg.de/boards/1/topics/35 . Note that this is an arbitrary range of pressure levels (the CMIP5 ones), independent of the original coordinates.