逐3小時啟動CMAQ之嘗試

背景

東亞或中國空氣品質的模擬，ecmwf的[CAMS預報][CAMS]已經有非常優異的表現。此處仍然重複進行CMAQ模擬，目標主要是在建立台灣地區高解析度模擬所需的初始與邊界條件(IC/BC)。
- 但似乎5天的模擬使得濃度表現越來越遠離[CAMS的預報][CAMS]，並且
- 在台灣地區的初始化與前一天的模擬結果，發生[突兀的變化]join_nc，表示CAMS和CMAQ系統之間存在不小的差異。
縮短重新模擬的時距或許是一個不錯的策略。最極端的作法，就是將所有[CAMS預報][CAMS]時間框架結果（逐3小時、共120小時、41筆。），都拿來進行初始化，直接將結果以線性內插連結，形成符合CMAQ化學與傳輸機制、同時也符合[CAMS預報][CAMS]的濃度場，再以此作為下層的IC/BC。
此一做法的可以達到的好處：
1. 降低境外排放量不確定性對CMAQ系統造成的影響。境外排放量確實是一個好議題。
2. 提高模擬結果與[CAMS預報場][CAMS]的相似性，包括空間分布與數量級。
3. 確保有足夠的空氣品質項目、化學平衡、符合該時間框架的傳輸機制。
4. BCON必須有逐時數據、直接線性內插可減省不少CMAQ的積分時間。
當然也可能會引進不少的誤差
1. 雖然解析度不高，只模擬導入IC的單一小時，將會失去CMAQ模式的整體特色。（特別是夜間不利擴散在局部網格累積的較高濃度分布）
2. 時間線性內插的部份，不必然達成化學與傳輸機制的平衡
3. 即使以沙盒方式、另外進行系統的建置與模擬，整體作業仍然具有[原本流程][fcst]的複雜性。

前處理

與[預報系統][fcst]共用氣象等檔案
1. 氣象部份：沿用WRF及mcip系統的結果([GFS數據自動下載轉換][GFS]、[GFS數據驅動WRF][gfs2wrf]及[mcip][mcip])
2. 排放部份：與[預報系統][fcst]相同
3. 邊界條件：也是由[grb2bcon.py][grb2bcon]所產生。
4. 下載CAMS預報結果。詳見
  - earth套件貼上CAMS臭氧濃度、
  - [CAMS預報數據寫成CMAQ初始檔][CAMS_ic]及
  - [CAMS預報結果之下載][get_all]。
每時間框架處理成CMAQ的ICON檔案
- 刪除非屬CAMS預報相關的空品項目，只剩50項，以縮小檔案容量(模版檔案templateCWBWRF_45k.ncV50K24)。
- 因CAMS會更新預報內容，轉檔作業也需配合更新。

CAMS檔案轉寫成ICON檔

由每日執行的[grb2icon.py][CAMS_ic]進行改寫成grb2iconV50.py
1. 檔案名稱中含有數據的日期時間，由grib檔中的forecast_time0及initial_time等2項資訊計算而得。
2. 套用不同模版檔案

$kuang@node03 /nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022
$diff grb2iconV50.py grb2icon.py
44,45c44,46
< fcst_hr=float(np.array(nc['forecast_time0'])[0])
< bdate=datetime.datetime.strptime(nc.variables[V[3][0]].initial_time,'%m/%d/%Y (%H:%M)')+datetime.timedelta(hours=fcst_hr)
---
> td=datetime.datetime.today()
> bdate=datetime.datetime(td.year,td.month,td.day)
77,79c78,79
< path='/nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022/'
< fname=path+'templateCWBWRF_45k.ncV50K24'
< fnameO='/nas2/cmaqruns/2022fcst/grid45/icon/ICON_'+bdate.strftime('%Y%m%d%H')+'_CWBWRF_45k'
---
> fname='/nas1/cmaqruns/2022fcst/data/icon/ICON_template_CWBWRF_45k'
> fnameO=fname.replace('template','yesterday')

執行迴圈腳本(grb2iconV50.cs)

沿用get_all.cs的邏輯，應用在每一時間框架
1. 因所有項目檔案龐大，如要減少對ecmwf網站的存取次數，貌似以空品項目橫向切割較為合理(未比較按預報時間切割方案所需處理時間，至少減省merge.cs的時間)。
2. python所需執行迴圈的時間仍是工作瓶頸，放在背景個別執行可以提高整體批次的效率。
3. CAMS預報會較實際時間晚20小時(12小時模式時間、8小時的時差)將近1天，因此41個時間框架還是不能滿足所有未來5天CMAQ預報所需。

#kuang@node03 /nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022
#$ cat grb2iconV50.cs
src=/nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022
ncks=/usr/bin/ncks
cd $src
for t in {0..40};do
for i in 1 2 3;do
  $ncks -O -d forecast_time0,$t allEA_$i.nc allEA_$i.nc_0
done
if [[ -e AllEA.nc ]];then rm AllEA.nc;fi
./merge.cs
./grb2iconV50.py >& /dev/null &
sleep 10
done

CCTM執行腳本

必須執行CMAQ的理由是因為BCON需要完整的空品項目，不是ICON的50項可以滿足。
除了項目之外，也必須保持化學場的平衡及穩定、同時可以符合WRF氣象場在各個時間框架當下造成的擴散與傳輸作用。

profile.config模版差異

因為只為了達成平衡，CCTM執行時間($runlen)取1小時。
檔名標籤($startdate)：必須加上小時
起始日期($START_DATE)：隨著時間框架而移動

kuang@DEVP /nas2/cmaqruns/2022fcst
$ diff project.config_loop project.config_loopic
5,6c5,8
<  set startdate = YYYYJJJ      # YYYYDDD
<  set runlen    = 1210000      # HHH0000
---
>  set    hh = HH
>  setenv sdate YYYYMMDD
>  setenv startdate ${sdate}${hh} # YYYYMMDDHH
>  setenv runlen    10000               # HHH0000
13c15
<  set cmaqbcdate = ${startdate}
---
>  set cmaqbcdate = today
19c21
<  set START_DATE   = "mcip_start"
---
>  set START_DATE   = "BEGD"

逐3小時執行迴圈

配合前述profile.config模版進行更新。
1. 10碼的時間+日期標籤$startdate=${sdate}${hh}
2. $BEGD需隨著時間框架移動
CCTM跨日執行似乎會清空之前的檔案，因此還是將小時平均結果(ACONC、APMDIAG)備份到其他目錄(../daily2)較為穩當

#kuang@DEVP /nas2/cmaqruns/2022fcst
#$ cat doic.cs
today=$(date -d -0day +%Y%m%d)
BEGD=$(date -d "$today -0days" +%Y-%m-%d)
mcip_start=$BEGD
mcip_end=$(date -d ${BEGD}+4days +%Y-%m-%d)
dir1=/nas2/cmaqruns/2022fcst/grid45/cctm.ic/daily
dir2=/nas2/cmaqruns/2022fcst/grid45/cctm.ic/daily2
mkdir -p $dir2
for id in {0..4};do
  for HH in {00..21..3};do
    BEGD=$(date -d "$today +${id}days" +%Y-%m-%d)
    YYYYJJJ=$(date -d ${BEGD} +%Y%j)
    YYYYMMDD=$(date -d ${BEGD} +%Y%m%d)
    if ! [ -e /nas2/cmaqruns/2022fcst/grid45/icon/ICON_${YYYYMMDD}${HH}_CWBWRF_45k ];then continue;fi
    cp project.config_loopic project.configIC
    for cmd in 's/YYYYJJJ/'$YYYYJJJ'/g' \
           's/YYYYMMDD/'$YYYYMMDD'/g' \
           's/HH/'$HH'/g' \
           's/BEGD/'$BEGD'/g' \
           's/mcip_start/'$mcip_start'/g' \
           's/mcip_end/'$mcip_end'/g';do
      sed -ie $cmd project.configIC
    done
    csh ic.csh
    for sp in CONC PMDIAG;do
      aconc=$dir1/CCTM_A${sp}_v532_intel_CWBWRF_45k_${YYYYMMDD}${HH}.nc
      if [ -e $aconc ];then mv $aconc $dir2;fi
    done
  done
done

run.cctm.45.csh腳本差異

呼叫字尾有IC版本的profile.config及cctm主程式腳本，儲存到cctm.ic目錄。沙盒作法。

#kuang@DEVP /nas2/cmaqruns/2022fcst
#$ diff ic.csh run.cctm.45.csh
7,8c7,8
< set sfile       = ./project.configIC
< set sourcefile  = ./cctm.source.v5.3.1.ae7_IC
---
> set sfile       = ./project.config
> set sourcefile  = ./cctm.source.v5.3.1.ae7
30c30
<  set myjob   = ic                   # which cases
---
>  set myjob   = fcst                   # which cases

cctm.source.v5.3.1.ae7 差異

主程式腳本有4項修改
1. 執行長度：呼叫前述執行時間($runlen)1小時
2. 模擬結果的時間標籤($CTM_APPL)：增加小時標籤，以避免被覆蓋。
3. IC檔名：增加時間標籤，以配合grb2iconV50.py轉檔結果，BC檔名，套用${cmaqbcdate}，讓其值為today,並未改變實質作法。
4. 換日關上$NEW_START：每一次都是True，CCTM才會讀取ICON

kuang@DEVP /nas2/cmaqruns/2022fcst
$ diff cctm.source.v5.3.1.ae7 cctm.source.v5.3.1.ae7_IC
60,61c60,62
<  set STTIME     = 000000           #> beginning GMT time (HHMMSS)
<  set NSTEPS     = 240000           #> time duration (HHMMSS) for this run
---
>  set HH         = `echo $startdate|cut -c9-`
>  set STTIME     = ${HH}0000        #> beginning GMT time (HHMMSS)
>  set NSTEPS     = ${runlen}        #> time duration (HHMMSS) for this run
218d218
<
253c253
<   setenv CTM_APPL ${RUNID}_${YYYYMMDD}
---
>   setenv CTM_APPL ${RUNID}_${YYYYMMDD}${HH}
264,265c264,265
<   set icname = ICON_yesterday_${GRID_NAME}
<   set bcname = BCON_today_${GRID_NAME}
---
>   set icname = ICON_${cmaqicdate}_${GRID_NAME}
>   set bcname = BCON_${cmaqbcdate}_${GRID_NAME}
621c621
<   setenv NEW_START false
---
> #  setenv NEW_START false

後處理

逐3小時結果串成逐日檔

CCTM的BCON必須是逐時的時間框架，因此需進行時間內插。
逐3小時的框架在換日前的2個小時會遇到困難，解決方法是延長時間維度一個小時，儲存隔日第一小時的結果，來進行最後2個小時的內插。
1. 每天讀取9筆時間框架(當天8筆與隔天第一筆)來進行內插
2. 最後在倒入新檔案時，要注意不要倒最後一筆。

#kuang@dev2 /nas2/cmaqruns/2022fcst/grid45/cctm.ic/daily
#$ cat join_icon.py

import numpy as np
import netCDF4
import sys, os, subprocess
import datetime
from dtconvertor import dt2jul, jul2dt

bdate=datetime.datetime.strptime(sys.argv[1],"%Y-%m-%d")
YMDs=[bdate+datetime.timedelta(days=i) for i in range(5)]
path='/nas2/cmaqruns/2022fcst/grid45/cctm.ic/daily/'
tmps={i:path+'CCTM_A'+i+'_temp.nc' for i in ['PMDIAG','CONC']}
for sp in list(tmps)[1:]:
  fcst=path+'../daily2/CCTM_A'+sp+'_v532_intel_CWBWRF_45k_YYYYMMDDHH.nc'
  for day in range(1,5):
    YMDsd=YMDs[day].strftime('%Y%m%d')
    fnameO=path+'CCTM_A'+sp+'_v532_intel_CWBWRF_45k_'+YMDsd+'.nc'
#   os.system('cp '+tmps[sp]+' '+fnameO)
    nc1 = netCDF4.Dataset(fnameO, 'r+')
    V=[list(filter(lambda x:nc1.variables[x].ndim==j, [i for i in nc1.variables])) for j in [1,2,3,4]]
    nt,nlay,nrow,ncol=(nc1.variables[V[3][0]].shape[i] for i in range(4))
    nv=len(V[3])
    if day==1:var=np.zeros(shape=(nv,nt+1,nlay,nrow,ncol))
    for h in range(0,25,3):
      sdate=YMDs[day]+datetime.timedelta(hours=h)
      YMDsd,HH=sdate.strftime('%Y%m%d'),sdate.strftime('%H')
      fname=fcst.replace('YYYYMMDD',YMDsd).replace('HH',HH)
      if not os.path.isfile(fname):continue
      nc = netCDF4.Dataset(fname, 'r')
      for v in range(nv):
        var[v,h,:,:,:]=nc.variables[V[3][v]][0,:,:,:]
      if h==0:
        nc1.SDATE=nc.SDATE
        nc1.variables['TFLAG'][:,:,0]=nc.SDATE
      nc.close()
    for t in range(0,24,3):
      var[:,t+1,:,:,:]=2./3.*var[:,t,:,:,:]+1./3.*var[:,t+3,:,:,:]
      var[:,t+2,:,:,:]=1./3.*var[:,t,:,:,:]+2./3.*var[:,t+3,:,:,:]
    for v in range(nv):
      nc1.variables[V[3][v]][:,:,:,:]=var[v,:nt,:,:,:]
    nc1.close()

因為ACONC逐日檔非常大(14G)，必須先在程式外複製(5個檔同步複製約需10分鐘)。
本程式記憶體需求龐大，需以超微工作站執行，DEC會卡住。

combine差異

kuang@DEVP /nas2/cmaqruns/2022fcst
$ diff combine.sh combine.shIC
5c5
< ymd=$(echo $nc|cut -d'_' -f7|cut -d'.' -f1)
---
> ymd=$(echo $nc|cut -d'_' -f2|cut -d'.' -f1)
11c11
< export cctmout=${BASE}/grid$ii/cctm.fcst/daily
---
> export cctmout=${BASE}/grid$ii/cctm.ic/daily

結果討論

比較傳統積分結果與該小時CCTM初始化後的CAMS預報差異。此處比較PM₁以顯示空品項目、WRF氣象場的綜合效果。


(a)東亞範圍地面CMAQ PM₁之模擬結果(CCTM 積分48小時)


(b)東亞範圍地面CAMS PM₁之模擬結果(CCTM初始化後積分1小時)

(a)積分結果的濃度值較(b)CAMS值低1/2，CAMS值高到不可接受的程度。檢討可能原因：
1. 因grb2icon系列程式是引用一組平均粒徑分布比例來分配粒狀物的濃度，此系統性的不一致尚待解決。尤其在海面上的粒狀物，如果套用高度平均的分配比例，將會造成細粒部份嚴重高估。
2. CAMS粒狀物單位(Kg/Kg Mixing Ratio)轉換需要空氣的密度，而密度會隨著大氣壓力(高度)、溫度而變化(n/V=P/RT)。grb2icon系列程式為求降低記憶體需求同時因應不同步執行進度、將4階維度的密度對時間軸予以平均（預報範圍內的5天平均值）。此外，CAMS與CMAQ之間高度的對照，也可能是誤差的來源之一。
3. 東西側邊界，圖(a)因入流邊界濃度較低。
4. 日本南方奄美大島的颱風範圍與太平洋中硫磺島，圖(a)因雲雨現象而較低。
陸地上、中國東半壁圖(a)較圖(b)有較高的解析度，這是CCTM模式的效果。雖然CAMS解析度有0.4度，但圖(b)似乎無法模擬出城市及交通排放在垂直擴散不良情況(00Z)下的效果。