CAMS高度之內插方式

Table of contents

背景
程式設計

背景

此處詳細處理CAMS之IFS模式與CMAQ之高度對照及內插議題。
- CMAQ並不定義高度，而是沿用WRF之設定，在real/wrf所讀取的namelist.input檔中（eta_levels）
- CAMS空品預報濃度提供2種高度下載方式：25層定壓層、與模式內設層數（137層）。其定義在IFS Documentation

p_k+1/2 = A_k+1/2 + B_k+1/2p_s

p_k = 1/2 (p_k-1/2 + p_k+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']

程式設計

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]