BCON南東北西4面2維濃度檔之轉接程式

背景

雖說模擬範圍外BCON檔案的濃度值，來自於CAMS空品預報的結果，然也需有顯示軟體，來證實濃度檔之間的轉換不會發生問題。
BCON邊界檔案的順序是南->東->北->西，逆時針方向排序的濃度檔，目前並沒有專用的顯示軟體，主要的困難：
1. BCON將4個側面邊界的水平方向、以一條逆時針方向的軌跡線來整合，方向與一般定義不同，詳EAC4檔案轉成4階邊界檔案。
2. 垂直方向的範圍與水平範圍不成比例，網格解析度也差異很大。
3. 垂直網格的數目遠遠不如水平方向，即使以網格為單位，也無法均衡展示。
4. 就檔案的格式而言，BCON的全域屬性FTYPE=2,ICON則為1，需另存新檔。
CAMx時代曾以某個單位高度(50m)作為單一垂直間距，重新將非等間距的垂直向轉成等間距系統
- 最後結果再以SURFER軟體一併放大垂直高度，以與水平網格匹配
VERDI
- 雖然可以畫3維濃度的垂直剖面。但是BCON檔案不是一個3維濃度檔。
- 並沒有某一維度放大的功能。需以外部程式先予以處理。
結果檔案之模版由程式內製、或由外部ncks產生，何者為佳?
- 雖然BCON檔案也將對應之濃度檔案網格數、格距等訊息儲存在其全域屬性之內，然而由程式內產生新的模版，不單需要產生濃度矩陣、時間標籤矩陣，還有許多必要的全域屬性，並不會比較輕鬆。
- 相反的，自一既有相同模擬範圍的ICON檔案來裁剪出需要的模版，會容易很多。(詳下說明)

ICON模版之製作

水平網格數大於垂直數情況

ICON_today_CWBWRF_45K範圍的水平網格約有200，為一portrait形式之範圍，可按需要進行裁剪成模版。
BCON檔案垂直網格僅24格，因此轉成ICON檔的y軸後，還需予以增加。
模版x方向設定(icon.NCOLS)
- 南面及北面邊界：原東西向(bcon.NCOLS)
- 東面及西面邊界：原南北向(bcon.NROWS)
模版y方向設定(icon.NROWS)為原垂直網格數(bcon.NLAYS)的3倍 icon.NLAYS=1;icon.NROWS=72
ncks指令
- icon=../icon/ICON_today_CWBWRF_45k
- ncks -O -d LAY,0 -d ROW,0,71 -d COL,1,$icon.NCOLS $icon $tname
結果模版檔名($tname)之規則
- tpl={0:'SN',1:'WE',2:'SN',3:'WE'}
- path='/nas1/cmaqruns/2022fcst/data/bcon'
- tname=path+'/template'+tpl[i]+'_CWBWRF_45k.nc', for i in range(2)
- 只需做SN及WE等2種類的模版即可。

水平網格數不足的情況

儘量避免延長一個受限制之維度
- 過程會很複雜，需要：1.將其排在第1順位、2.打開限制、3.執行ncrcat或python延長維度、4.恢復限制、5.恢復順位
- 詳參加長一個limited維度
以原始3維的ICON為底版進行旋轉與裁剪
使用4個NCO指令：
- ncrename -d：維度之更名
- ncks -d：切割特定維度
- ncpdq -a：維度之排序
- ncatted -a：屬性之更新

cp 20221212/00/2022121200.ic template.nc #as a basic template
ncrename -O -d COL,A template.nc a #COL is spared
ncrename -O -d ROW,COL a b # turn the coord. sys. anti-clockwise 90deg. ROW become COL
ncrename -O -d A,ROW b a # ROW become 72 layers
ncks -O -d LAY,0 -d ROW,1,72 a b
ncpdq -O -a TSTEP,LAY,ROW,COL b templateWE_CWBWRF_45k.nc #re-sequence
ncatted -a NROWS,global,o,i,72 templateWE_CWBWRF_45k.nc # 2d-western and eastern BC templates
python /nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022/bcon2icon.py BCON_20221204_CWBWRF_45k

程式設計

IO

引數：BCON檔案
模版：詳前述說明，共2個檔案
結果：m3nc格式之濃度檔，共南、東、北、西等4個檔案

不同範圍解析度的應用

從輸入檔名中解析其解析度

...
res={'CWBWRF_45k', 'SECN_9k', 'TWEPA_3k'}
#read a BC file as rate base
fname=sys.argv[1]
ipas=0
for r in res:
  if r in fname:
    ext=r
    ipas=1
if ipas==0:sys.exit('not right resolution extension')

引用不同解析度模版
輸出檔名中有2個WE：TWEPA，因此增加前一碼以提高取代的精準度。

for i in range(4):
  tname=path+'/template'+tpl[i]+'_'+ext+'.nc'
  fnameO=tname.replace('template','today').replace('y'+tpl[i],'y'+out[i])

BCON水平網格順序之轉換

BCON將4面邊界的水平方向、以一條逆時針方向的軌跡線來整合4個側面，因此在北側及西側，BCON指標順序與傳統的軸線方向是相反的。
水平網格數nbnd=(nc.NROWS+nc.NCOLS)*2+4
反向順序放在等式的左邊或右邊
- 左邊、新矩陣、範圍為0 ～ $NCOLS（南面）或${NCOLS}-1 ~ 0 （北面、翻轉）
  - np.array的訖點指標無法表示0以下的位置（-1不是0以下，而是最末項）
  - 當反轉順序時無法指定到指標0，頂多到1（python的索引規則）。
- 右邊、舊矩陣BCON
  - 北側與西側需翻轉，雖然有點複雜，卻不是不能做。

水平網格對照之指標

各側邊界在BCON之起迄指標(ibnd)與方向(drn)
1. 南側邊界：左下角(西南) → 右下(東南)，起訖指標：(1,nc.NCOLS+1)，順向
2. 西側邊界：右下角(東南) → 右上(東北)，起訖指標：(nc.NCOLS+2,nc.NCOLS+nc.NROWS+2)，順向
3. 北側邊界：左上角(西北) → 右上(東北)，起訖指標：(nc.NCOLS*2+nc.NROWS+2,nc.NCOLS+nc.NROWS+2)、方向相反。
4. 西側邊界：左下角(西南) → 左上(西北)，起訖指標：(nc.NCOLS*2+nc.NROWS*2+3,nc.NCOLS*2+nc.NROWS+3)、方向相反。
等式右邊、BCON檔案指標之通式，i=0,1,2,3
- [ibnd[i][0]:ibnd[i][1]:drn[i]]

drn={0:1,1:1,2:-1,3:-1}
ibnd={0:(1,nc.NCOLS+1),
1:(nc.NCOLS+2,nc.NCOLS+nc.NROWS+2),
2:(nc.NCOLS*2+nc.NROWS+2,nc.NCOLS+nc.NROWS+2),
3:(nc.NCOLS*2+nc.NROWS*2+3,nc.NCOLS*2+nc.NROWS+3)}

等式左邊、ICON檔案之起訖指標。因涉指標0，因此只能自0開始，只能有一個方向。

i1s={0:0,1:0,2:0,3:0}
i2s={0:nc.NCOLS,1:nc.NROWS,2:nc.NCOLS,3:nc.NROWS}

新舊矩陣之對照

nc1[v][:,0,1::3,i1s[i]:i2s[i]]=nc[v][:,:,ibnd[i][0]:ibnd[i][1]:drn[i]]

垂直方向的轉換

24層轉72層。
- 72層每3層的中間、對照到BCON的24層。
- 新檔案的y軸[1::3]對照到既有BCON檔案z軸[:]
底層及最上層，直接令為BCON的底層及最上層
中間以外的2層：線性加權漸變

    nc1[v][:,0,0,:] =nc1[v][:,0,1,:]
    nc1[v][:,0,-1,:]=nc1[v][:,0,-2,:]
    nc1[v][:,0,2:-1:3,:]=(nc1[v][:,0,1:-2:3,:]*2+nc1[v][:,0,4::3,:]  )/3
    nc1[v][:,0,3:-1:3,:]=(nc1[v][:,0,1:-2:3,:]  +nc1[v][:,0,4::3,:]*2)/3

北方與西方水平方向之轉置

因BCON檔案逆時針旋轉軌跡的特性，在北方、西方的順序與傳統座標軸的方向是相反的，必須進行轉置(np.flip)
- 北方邊界(i=2)的順序是由東向西
- 西方邊界(i=3)的順序是自北向南

    if i>=2:
      nc1[v][:]=np.flip(nc1[v][:],axis=3)

結果討論

南北側邊界

一般污染物濃度在對流層底部擴散，但SO2有高空之濃度。由於南側邊界已經接近赤道，可能發生全球性傳輸的現象。
北側邊界值得關注的項目是CAMS模擬砂塵的效果。由圖來看，西北方的濃度並不低，且已擴散到高空。由圖中也可以發現局部地區因雲雨現象而出現較低濃度範圍，顯見CAMS也具有模擬移除現象的能力。


2022/10/27~31南側邊界SO₂平均濃度之分布


2022/10/27~31北側邊界粗顆粒物平均濃度之分布

東西側邊界

東側太平洋以海鹽濃度最具特徵，如圖顯示東北方地面濃度較高，應為該處海面風速較大所致。濃度向上遞減，因污染源在海面上。
西側邊界濃度較低，西南方為印度洋範圍，西側邊界北半部為陸域，海鹽濃度趨近於零。


2022/10/27~31東側邊界海鹽顆粒平均濃度之分布


2022/10/27~31西側邊界海鹽顆粒物平均濃度之分布

程式碼下載

Download: BCON檔案4面2維濃度檔之轉接程式bcon2icon.py