運用GFS/CWB/CAMS數值預報數進行台灣地區CMAQ模擬

(DEPRECATED !! 本文已不再更新，最新版本請至https://sinotec2.github.io/Focus-on-Air-Quality/GridModels/ForecastSystem/查詢, 2022-10-13 20:17:00)

(todo list加在此頁之末)

背景

• 以光化學網格模式進行高解析度空氣品質數值預報、對空品不良狀況的預告以及應變措施有重要的參考價值。確定性(deterministic)模式主要的特點與好處包括：
  • 基本的背景空氣品質預報
  • 區分境外與本土排放的貢獻
  • 區分背景與特定污染源的貢獻
  • 區分人為污染或天候因素
  • 緊急排放、或是假設性排放之投入與增量分析
• 背景預報的功能過去大多依靠統計模式、類型預報、類神經網路、機器學習等等模式系統，邏輯系統與因果關係模式較少，主要限制與解決對策列表如下：

高解析度空氣品質數值預報的限制與對策

整體作業流程

graph TD
    A((GFS</br>forecasting)) -- WPS --> B(met_em)
    C((CWBWRF</br>forecasting)) -- get_M-A0064.cs --> D(wrfout_d01/3)
    B --> E{mk_metoa}
    D --> E
    E --> F{real,WRF}
    F --> K{mcip}
    G((CAMS</br>forecasting)) -- grb2ic, grb2bc--> H(CMAQ ICON,BCON)
    H --> I{CMAQ}
    K --> I
    I --> J((displaying</br>on website))

下載作業

CWBWRF之下載

• CWBWRF的時間約與實際時間落差6~10小時，因此如需00Z（0800LST）起始的預報結果，則需在1800LST起動下載。（下載檔案的時間標籤可以使用[EAC4_Times]來檢視。）
• CWB沒有提供下載選項、所幸變數、層數有限，也僅提供每6小時預報，檔案容量與個數並不多。
• CWB WRF程式結果每6小時更新，分別為每天的2/8/14/20時(LST)等4次。各次預報起始時間與安排如下：

GFS預報檔下載與分析

• GFS預報時間較長(384 hr)，層數、變數項目較完整，然而NOAA也提供了部分下載的功能，可供展示使用。
• GFS預報數據是否存檔的考量
  1. 畢竟是模式計算結果，仍有一定程度的不準度。
  2. 相同解析度另有FNL再分析數據，結合觀測及預報之三維場，有更高的正確性。只是稍晚才會上架，約延遲2天。
  3. 如果GFS一定要存檔，僅需渡過前述空窗期即可刪除。或可考慮下載模式起始所使用的分析檔。(https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.YYYYMMDD/HH/atmos/gfs.tHHz.pgrb2.1p00.anl)、YMDH分別為UTC之年月日時，約35M
  4. NOAA 也只開放下載最近10天內的預報作業結果。
• 下載的作業細節可以參考GFS數據驅動WRF->GFS預報結果的下載
• 下載安排如下：

• 下載後即按照WRF的執行流程進行準備。

^* real的執行本身是平行，但在整體腳本流程中只是一個動作，不會由腳本另外開啟需等候完成的程序。

CAMS預報結果之下載

• 使用get_All.py
• 引數：
  1. dt：YYYY-MM-DD， dt=$(date -d "now -20hours" +%Y-%m-%d)
  2. hr:HH:00，hr=$(date -d "now -20hours" +%H:00)
  3. 1 ~ 3部分物質種類。詳見CAMS_ic
• 由於CAMS數據較實際時間晚12小時，再加上時差，因此差距為20小時。如要下載0時數據，只能
  1. 在前1天的20時(LST)，以預報第1天(第9筆)的結果來替代當天0時值。(get_All.py每3小時下載一筆)
  2. 在8時下載前1天12:00起始的預報，以預報第5筆)的結果來替代當天0時值。

CAMS預報數據之下載與分析

• CAMS預報數據除了寫成CMAQ邊界檔及初始場之外，也是earth展示的重點項目。
• CAMS預報的更新頻率為12小時，上架時間較模式起始時間晚20個小時。如表所示：

執行下載分析之腳本

• get_all.cs
  1. 0 ~ 120小時，每3小時下載
  2. 27項污染物質，分3批下載
  3. 下載後轉nc檔有較快的讀取速度
  4. 抽出第8筆，為當日0時之預報值作為CMAQ的初始值(注意ncks -O以覆蓋同一檔名)
  5. 執行橫向整合merge，後將做為CMAQ的初始檔。注意必須先刪除結果檔以利初始化。
  6. 執行grb2icon.py
  7. 執行grb2bcon.py

#kuang@node03 /nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022
#$ cat get_all.cs
src=/nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022
ncks=/usr/bin/ncks
cd $src
dt=$(date -d "now -1days" +%Y-%m-%d)
hr=$(date -d "now -20hours" +%H:00)
hr=00:00
for i in 1 2 3;do
  rm allEA_$i.grib
  ./get_All.py $dt $hr $i >& /dev/null
  PATH=/bin:/opt/miniconda3/envs/ncl_stable/bin ncl_convert2nc allEA_$i.grib -nc4c >& /dev/null
  if [[ -e allEA_$i.nc4 ]];then mv allEA_$i.nc4 allEA_$i.nc;fi
  $ncks -O -d forecast_time0,8 allEA_$i.nc allEA_$i.nc_0
done
if [[ -e AllEA.nc ]];then rm AllEA.nc;fi
./merge.cs
./grb2icon.py &
./grb2bcon.py &

自動定期執行

• crontab：每天0時執行，以配合CMAQ模擬之需求。

0 0 * * * /nas1/ecmwf/CAMS/CAMS_global_atmospheric_composition_forecasts/2022/get_all.cs

計算與核心之分配

• 短時間的計算需求，直接進行即，可此處特別處理需要較長時間的計算。
• 雖然大多數模式、次模式都可以平行運作，但因為各個模式使用的編譯器、mpi種類版本略有差異，因此需先個別調控(tuning)到最佳狀態，再行組合、納入到自動化作業流程之中。
• 目標及優先順序是
  1. 穩定計算達到正確的結果
  2. 使用到單機最多的計算核心
  3. 使用到跨機最多的計算核心

REAL and WRF

• 注意要正確指定執行檔所需程式庫的路徑($LD_LIBRARY_PATH)
• 執行腳本如逐日WRF及CMAQ執行腳本所示
  1. real與wrf必須循序進行
  2. mcip隨後是執行CCTM，然而因第3層mcip執行完是執行第1層的CCTM，因此可在背景獨立平行執行，並不會造成需要瓶頸。
• WRF的計算核心安排可以參考GFS數據驅動WRF->執行核心數之考量的詳細說明，此處說明bash技巧。
• mpirun選項：儲存在環境變數序列 ${MPI[$i]}之中，$MPI的定義及應用如下

MPI=( '-f machinefile -np 200' '-f machinefile -np 196' '-f machinefile -np 140')
...
  # real
  LD_LIBRARY_PATH=/nas1/WRF4.0/WRFv4.3/WRFV4/LIBRARIES/lib:/opt/intel_f/compilers_and_libraries_2020.0.166/linux/compiler/lib/intel64_lin:/opt/mpich/mpich-3.4.2-icc/lib /opt/mpich/mpich-3.4.2-icc/bin/mpirun ${MPI[$i]} /nas1/WRF4.0/WRFv4.3/WRFV4/main/real.exe >& /dev/null
  # wrf
  LD_LIBRARY_PATH=/opt/netcdf/netcdf4_gcc/lib /opt/mpich/mpich3_gcc/bin/mpirun ${MPI[$i]} /opt/WRF4/WRFv4.2/main/wrf.exe >& /dev/null

• 跨機核心數檔案：以相同名稱存放在每個Domain的目錄下，平均分配所有的負荷到2台超微工作站機器上。

kuang@master /nas1/backup/data/NOAA/NCEP/GFS/YYYY
$ for i in $(ls */machinefile);do echo $i;cat $i;done
CWBWRF_45k/machinefile
devp:100
dev2:100
SECN_9k/machinefile
devp:98
dev2:98
TWEPA_3k/machinefile
devp:70
dev2:70

• 總核心數必須等於namelist.input中nproc_x與nproc_y的乘積，也必須小或等於machinefile中2台機器核心數的總合。

2022-09-16 修正

• d02速度太慢的問題，可能出在FDDA太過強烈，系統不容易達成動力的平衡。將3層domain同時進行模擬(two way (tw) nesting)，似能解決此一問題。
• 3層套疊全部需時約1.5小時，較原先約3小時縮減一半時間。需修改：
  1. namelist.input中的max_dom，由1->3，相應網格數、起始網格位置、網格間距比值等等。
  2. 循序執行real_wrf的迴圈，只需執行一次。
  3. 新增網格定義與目錄。（./tw_CWBWRF_45k）
  4. tw 所需的核心數120，machinefile：

$ grep proc tw_CWBWRF_45k/namelist.input
 nproc_x                             = 10,
 nproc_y                             = 12,
$ cat tw_CWBWRF_45k/machinefile
devp:60
dev2:60

• 核心數大幅減少(受限於3層網格系統之公約數)，計算速率並沒有下降，似乎還有所提昇。應為網路傳輸速度造成瓶頸所致。

每日預報mcip

腳本特色

• 由於執行完CCTM後，會需要下一層的METBDY3D檔案，來產生下一層的邊界條件，因此在執行CCTM之前，將所有網格的mcip一次執行完，會是較好的作法。

mcip 執行的核心控制

• 相較WRF只檢討最高核心數，mcip對核心數的工作安排有較高的敏感性，如果不是能整除網格數的核心數，會導致計算的不穩定與錯誤。
• 加以mcip的腳本是csh指令，必須在腳本中逐一設定。
• 因mcip只是檔案轉換，並沒有實質的數值積分，所以此處僅安排以單機方式處理。

1. 如前所述，不整除的核心數會導致計算不穩定。d03因公版模式設定NROW=131，是個質數，只能屈就使用46個核心。
2. 2022-09-16 修正：d01及d02會因為核心數太多造成寫出檔案錯誤，發生在較小的檔案如GRIDBDY2D.nc或GRIDCRO2D.nc，可能是輸出不同步所造成。將d01及d02的NP分別改成16及20

CCTM

• 與mcip類似的作法，cctm是每個domain有個別的腳本，其NPCOL_NPROW設定如下：

• 跨機核心數檔案：以Domain網格間距識別名稱存放在CCTM的根目錄，平均分配所有的負荷到2台機器上。

kuang@dev2 /nas2/cmaqruns/2022fcst
$ for i in $(ls machinefile??);do echo $i;cat $i;done
machinefile03
devp:92
dev2:92
machinefile09
devp:98
dev2:98
machinefile45
devp:88
dev2:88

自動執行腳本與執行時間

下載fcst.cs腳本

• 詳見逐日WRF及CMAQ執行腳本

分段時間分析

• 總計時間約需8小時

成果展示

• 動態展示以earth套件，靜態則以wrf-python。
• 不同解析度：個別建立自己的套件，以8084/6等3個port來展示(修改最少、最快上手)。
• 需發展cmaq_json.py來解讀wrfout與CCTM_ACONC檔案。
• 自動執行排程
  1. 每3小時連結current檔案
  2. 每天處理wrfout與CCTM_ACONC預報結果

cmaq_json.py

• 詳見[cmaq_json.py說明][cmaq_json.md]
• 下載cmaq_json.py

lnk_curr.cs

• 依序進行3個解析度結果的連結
• 只連結wind及ozone
• 粒狀物另參照[2022-09-14-PM_earth]、跨日的處理詳[2022-09-15-join_nc]

y=$(date -d -8hours +%Y)
m=$(date -d -8hours +%m)
d=$(date -d -8hours +%d)
h=$(date -d -8hours +%H)
h=$(( 10#$h / 3 * 3|bc -l ))
h=$(printf "%02d" $h)
for grd in 45 09 03;do
  weather=/nas1/Data/javascripts/D3js/earthFcst${res}/public/data/weather
  path=$weather/$y/$m/$d
  for itm in wind ozone;do
    lev=-surface-level
    fn=${h}00-${itm}${lev}-fcst-${res}.json
    cr=${fn/${h}00/current}
    if [ -e $path/$fn ]; then
      ln -sf $path/$fn $weather/current/$cr
    fi
  done
done

預報結果比較

2022/09/01 1500空氣品質實測比較

• 時間：2022/09/01/06Z，模式起始時間9/1及環保署空品監測網1500(LST)AQI實測值

CMAQ預報東亞範圍地面臭氧濃度結果	環保署AQI值

東南亞強颱「軒嵐諾」污染事件預報

• 時間：2022/09/03/06Z，模式起始時間8/30(CAMS)及8/31(CMAQ)

CAMS預報東南亞地面臭氧濃度結果	CMAQ預報東亞範圍臭氧濃度

• CAMS似乎對颱風的風速與位置有所低估。
• 對大陸東南沿海城市的臭氧煙陣也有低估。
• 由於CAMS預報時間不及第5天。CMAQ可補足此一部分，結果如下：

CMAQ預報2022/09/04 06Z東亞範圍臭氧濃度

東北季風境外影響

• 2022/09/06 1700臭氧8小時值之預報與觀測

臺灣地區臭氧8小時值預報濃度	環保署AQI值

• 圖中顯示台灣本地的排放量已經比TEDS11減很多了，因為模式結果顯示有高估本地的傾向。

作業過程之問題與解決

CCTM_CONC寫出時遇到的格式問題

• 發生時機：新日或換日，寫出CCTM_CONC時，發生在DEV2結點。grid09、grid03。
• 可能原因：跨機平行運作，因網路速度耽誤寫出時間差異
• 改善方式：避免跨機運作

東亞排放量REAS高/低估

• 經確認在[reas2cmaqD2.py][reas2cmaqD2.py]沒有做網格面積的校正，修正如下。
• REAS 0.25度網格面積的計算。原有網格系統與經緯度：

  lon=list(set(lon+[float(i.split()[0]) for i in l[9:]]))
  lat=list(set(lat+[float(i.split()[1]) for i in l[9:]]))
# generate the x and y arrays for REAS datafile
# def nlat, nlon here
for ll in ['lon','lat']:
  exec(ll+'.sort()')
  exec('n'+ll+'=nn=int(('+ll+'[-1]-'+ll+'[0])/0.25)+1')
  exec(ll+'M=['+ll+'[0]+0.25*i for i in range(nn)]')
  exec(ll+'n={l:'+ll+'M.index(l) for l in '+ll+'M}')
lonm, latm = np.meshgrid(lonM, latM)

• 面積的計算
  1. 按照新網格系統的原點，將經緯度轉成直角座標系統
  2. 計算dx與dy。並假設最末行及最末列，與前一行及前一列相同。
  3. 面積即為dx與dy的乘積

pnyc = Proj(proj='lcc', datum='NAD83', lat_1=nc.P_ALP, lat_2=nc.P_BET, lat_0=nc.YCENT, lon_0=nc.XCENT, x_0=0, y_0=0.0)
x,y=pnyc(lonm,latm, inverse=False)
dx,dy=np.zeros(shape=x.shape),np.zeros(shape=x.shape)
dx[:,:-1]=x[:,1:]-x[:,:-1];dx[:,-1]=dx[:,-2]
dy[:-1,:]=y[1:,:]-y[:-1,:];dy[-1,:]=dy[-2,:]
area=dx[:]*dy[:]

• 面積的應用，將REAS排放量單位轉成每m²。因為面積是2維，排放量nc檔案是5階，因此其他階必須加上None，見矩陣階層numpy.newaxis(None)的用法。

var[:]*=1.E6/31/24/3600. #Ton/Mon to g/s
var[:,:,:,:,:]/=area[None,None,None,:,:] #change to g/s/m^2

• 面積的校正

 zz[t,:,: ] = griddata(xyc, c[:], (x1, y1), method='linear') *nc.XCELL*nc.XCELL

todo list

TEDS與grid03結果高估

• (TEDS11, 2019待更新)
• [2022-09-06-ptse_today]

1day BCON

• 1dayBCON

[run_mcip]: <> “” [1dbcon]: https://sinotec2.github.io/FAQ/2022/08/27/1dayBCON.html “逐日循序執行bcon.exe” [add0]: https://sinotec2.github.io/Focus-on-Air-Quality/GridModels/MCIP/add_firstHr/ “CMAQ Model System->Met. Chem. Interface Proc.->mcip結果初始小時值的延伸”