從空品檔案切割邊界濃度 BNDEXTR

背景

• 空品模式需要的邊界濃度，通常是更大範圍(全球模式如CAM-chem([2022-06-27-CAMx_BC])、上層網格、粗網格)的模擬、或再分析結果，因此會需要從空氣品質檔案中切割出邊界上的濃度，除了座標系統的對照、切割之外，還需要進行內插。
• CMAQ系統有bcon可以進行切割、也有python的版本。CAMx則需要BNDEXTR。執行則需要有c-shell 的腳本。
• PseudoNetCDF也有提供CAMx邊界濃度檔案格式(pncgen -f lateral_boundary)的讀取及輸出，因此也可以在python平台來處理。
  • ([2022-06-26-pncgen])

BNDEXTR.f程式

程式下載

• BNDEXTR可以從CAMx官網下載。

編譯

• 使用makefile
• 3個編譯器可以選擇：pgf90/ifort/gfortran
• 因以IO為主，計算並沒有平行化。

bndex.job範例

輸入

• standard input：共12項設定參數或檔名
• 檔案($INPF)：.avrg.grd02檔案（uamiv格式，需要完整所有高度的數據、涵蓋所有模擬時間）

呼叫程式

• pick讀取初始日期。現在這支程式的功能完全可以被pncdump取代。
  • -H：header only
  • -f lateral_boundary：指定檔案格式，詳見pncgen/pncdump 所有可接受的格式-CAMx、pncdump及lateral_boundary格式

pncdump -H -f lateral_boundary base.grd02.1801.bc |grep SDATE
                :SDATE = 2017365 ;
$ pncdump -H -f lateral_boundary base.grd02.1801.bc |grep STIME
                :STIME = 200000 ;               

結果檔案

• base.grd02.$dy.bc (bndary格式)
• base.grd02.$dy.ic (uamiv格式)

#kuang@114-32-164-198 /Users/camxruns/2016_v7/ICBC/bndextr_d2
#$ cat bndex-d2.job 
export EXE="/cluster/src/CAMx/bndextr/src/bndextr"
export INP="/nas1/camxruns/2016_v7/ICBC/bndextr_d2"
#for INPF in $(ls $INP/*.avrgIP|grep 16);do
for INPF in $(ls $INP/*.avrg.grd02|grep 1611);do
export dy=`echo $INPF|cut -d'/' -f7|cut -c 1-4`
export jul=`pick $INPF|grep jules|awk '{print $5}'`
echo $dy $jul
#YYNN JJJ
rm base.grd01.$dy.bc
rm base.grd01.$dy.ic
rm camx.nest.$dy.diag1
rm camx.nest.$dy.diag2

$EXE << EOF
Input average file |$INPF
Output BC file     |base.grd02.$dy.bc
Make IC file?      |T
Hour/date for IC   |20  $jul
Output IC file     |base.grd02.$dy.ic
Output projection  |LAMBERT
UTM zone           |0
Center/Pole lat/lon|120.9900,23.61000
True lats          |10.000,40.000,
Grid definition    |-124.500, -205.500, 3., 3.,83,137,15
diagnostic file #1 |camx.nest.$dy.diag1
diagnostic file #2 |camx.nest.$dy.diag2
EOF

done

設定說明

BNDEXTR.py程式

• python可以用pnc模組直接開啟並讀取CAMx的粗網格濃度檔(uamiv)、開啟邊界濃度檔案(lateral_boundary格式bc檔)、寫入數據、直接儲存。但是在python程式內延長檔案的時間軸、速度會很慢。如果時間不是很長差異不會太大，有需要的讀者可以參考BNDEXTR_pnc.py版本。
• 以下還是以nc檔案格式來處理，比較有效率一些。
• pnc及netCDF4程式庫的差異說明如下

整體流程邏輯

• 執行流程
  1. 以pncgen將lateral_boundary格式檔案轉成nc模版備用(可延長時間軸)
  2. 執行python讀取粗網格濃度檔、進行切割、截取及時間、空間的內插
  3. 按照模版輸出成nc檔案
  4. (在程式內)使用pncgen將nc檔轉成lateral_boundary格式檔案

IO檔案

• PseudoNetCDF模組會將lateral_boundary格式檔案的污染物濃度變數名稱，加上東西南北方向標籤。因此變數變成原來的4倍。
  • SO2 -> SOUTH_SO2
  • NOSPEC=60 -> NVARS = 240
• 按實際大小給定矩陣的格數(注意維度的順序)
  • 東西向：WEST_NO(TSTEP, ROW, LAY)
  • 南北向：SOUTH_NO(TSTEP, COL, LAY)
• 模版的準備：pncgen可以進行切割，但不能將時間軸變成可延長，需轉成nc檔另以ncks來處理。

pncgen -s TSTEP,0 -f lateral_boundary base.grd01.base.bc bnd.nc
ncks -O --mk_rec_dmn TSTEP bnd.nc a;mv a bnd.nc

• 引數與輸入檔：CAMx粗網格模擬結果
• 執行結果：細網格邊界濃度檔。nc格式及lateral_boundary格式結果

程式設計

• 4面邊界相關變數的簡化。此處使用dict來取代if…else的程式設計，項目包括
  1. 變數名稱：sides
  2. 粗網格相對細網格角落位置的起迄點：mni/mnj/mxi/mxj。維度長度1會造成矩陣對照的困難，使用flatten().reshape()來去掉。使用連乘指令np.prod。
  3. 粗網格軸線上的座標值：X
  4. 細網格軸線上的座標值：xy1d1

sides={0:'WEST',1:'EAST',2:'SOUTH',3:'NORTH'}
mni={0:minx,  1:maxx,  2:minx,  3:minx}
mxi={0:minx+1,1:maxx+1,2:maxx,  3:maxx}
mnj={0:miny,  1:miny,  2:miny,  3:maxy}
mxj={0:maxy,  1:maxy,  2:miny+1,3:maxy+1}
X={0:np.array([y1d0[i] for i in range(miny,maxy)]), 2:np.array([x1d0[i] for i in range(minx,maxx)])}
X.update({1:X[0],3:X[2]})
xy1d1={0:y1d1,1:y1d1,2:x1d1,3:x1d1}

• 內插與填入新檔案
  • 內插方式：使用Scipy.interpolate.CubicSpline一次內插整個矩陣。
    • 充分發揮Scipy平行計算的功能。
    • 直接使用粗細網格軸線上的座標值進行內插。因此不論粗、細網格間距的比例為何，都可以直接計算。
  • 新檔案的維度順序是[時間、XY、高度]，與一般濃度檔案習慣不同，以np.transpose()來轉置。

#cubicspline interpolation of coarse grid along x and y directions(axis=2)
for face in range(4):
  for v in V0[3]:
    var=nc0.variables[v][:,:,mnj[face]:mxj[face],mni[face]:mxi[face]]
    var=var.flatten().reshape(nt0,nlay0,np.prod(var.shape)//nt0//nlay0) #reshape for dropping dimension of 1
    cs=CubicSpline(X[face],var,axis=2)
    nc1[sides[face]+'_'+v][:,:,:]=np.transpose(cs(xy1d1[face][:]),axes=(0,2,1))#shape in [NSTEPS,NCOLS/NROWS,NLAYS]

• 最後程式內執行ncgen轉檔

pncg=subprocess.check_output('which pncgen',shell=True).decode('utf8').strip('\n')
if len(pncg)>0:
  os.system(pncg+' -f netcdf --out-format=lateral_boundary '+pathO+' '+pathO.replace('.nc','.bc'))
else:
  sys.exit('pncgen not found')

pnc及netCDF4程式庫的差異