Dr. Kuang's Utilities

CAMS高度之內插方式

 
  • Sep 17, 2022
  • views
>

背景

  • 此處詳細處理CAMS之IFS模式與CMAQ之高度對照及內插議題。

pk+1/2 = Ak+1/2 + Bk+1/2ps

pk = 1/2 (pk-1/2 + pk+1/2)

  • 式中的A/B值，如下表所示。
n a [Pa] b ph [hPa] pf [hPa] Geopotential Altitude [m] Geometric Altitude [m] Temperature [K] Density [kg/m^3]
137 0.000000 1.000000 1013.2500 1012.0494 10.00 10.00 288.09 1.223803
136 0.000000 0.997630 1010.8487 1009.5363 30.96 30.96 287.95 1.221341
135 3.757813 0.995003 1008.2239 1006.7900 53.92 53.92 287.80 1.218650
134 22.835938 0.991984 1005.3562 1003.7906 79.04 79.04 287.64 1.215710
133 62.781250 0.988500 1002.2250 1000.5165 106.54 106.54 287.46 1.212498
2 3.102241 0.000000 0.0310 0.0255 73721.58 74584.91 209.21 0.000042
1 2.000365 0.000000 0.0200 0.0100 79301.79 80301.65 198.05 0.000018
0 0.000000 0.000000 0.000000 - - - - -
  • 由於CAMS並不提供各層壓力的下載、而按照上述計算式各層的壓力會隨時間、位置而異，如需精確校正壓力，必須計算各個時間框架的3維壓力。
  • 壓力的內插，可以使用wrf-python的函數（詳範例）。
  • WRF 24層對應到ifs之層數（約略對照地理高度求得）
    • kk= ['137', '135', '133', '129', '125', '122', '120', '117', '114', '112', '110', '107', '105', '101', '96', '92', '87', '83', '78', '73', '67', '61', '56', '51']
  • 影響所及包括[[2022-08-16-CAMS_bc]]、[[2022-08-16-CAMS_ic]]

程式設計

from pandas import *
import netCDF4
from wrf import (getvar, interplevel)

df=read_csv('/nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022/heights.csv')
cname={'A':'a [Pa]','B':'b'}
for s in 'AB':
  for i in '01':
    exec(s+i+'=np.array([df.loc[df.n==int(k)-'+i+',"'+cname[s]+'"].values[0] for k in kk])')

root='/nas2/cmaqruns/2022fcst/grid45/wrfout/wrfout_d01_'
fnames=[root+str(i) for i in range(6)]
nc = netCDF4.Dataset(fnames[0], 'r')
nt,nlay,nrow,ncol=(nc.dimensions[i].size for i in ['Time','bottom_top','south_north','west_east'])
slp=np.zeros(shape=(nt*5,nrow,ncol))
pres=np.zeros(shape=(nt*5,nlay,nrow,ncol))
dens=np.zeros(shape=(nt*5,nlay,nrow,ncol))
p=np.zeros(shape=(nt*5+1,24,nrow,ncol))

t0=0
for fname in fnames:
  nc = netCDF4.Dataset(fname, 'r')
  nt=nc.dimensions['Time'].size
  for i in range(nt):
    slp[t0+i,:,:]=getvar(nc, "slp",timeidx=i)
    pres[t0+i,:,:,:]=getvar(nc, "pressure",timeidx=i)
  nc.close()
  t0+=nt
fname='/nas2/cmaqruns/2022fcst/grid45/mcip/DENS/METCRO3D.'+bdate.strftime('%Y%m%d')
nc = netCDF4.Dataset(fname,'r')
nt1,nlay1,nrow1,ncol1=nc.variables['DENS'].shape
dens1=nc.variables['DENS'][bt,:,:,:]
nc.close()
dens[:,:,:,:]=np.mean(dens1,axis=(0,2,3))[None,:,None,None]
dens[:,:,1:-1,1:-1]=dens1[:,:,:,:]  
p[:,:,:,:] =(A0[None,:,None,None]/100+B0[None,:,None,None]*slp[:,None,:,:])/2.
p[:,:,:,:]+=(A1[None,:,None,None]/100+B1[None,:,None,None]*slp[:,None,:,:])/2.
dens1 = interplevel(dens, pres, p)[::3,:,2:-2,2:-2]
  • getvar:[[2022-08-11-wrf_pythonTAB]]

[2022-08-11-wrf_pythonTAB]: https://sinotec2.github.io/FAQ/2022/08/11/wrf_pythonTAB.html “wrf-python [//end]: # “Autogenerated link references”

Search