REASv3.1地面排放檔案之處理

Table of contents

REASv3.1地面排放檔案之處理

背景

目標

全球或東亞地區之空氣污染排放量來源有很多，此處以日本國立環境研究所2019/12/4公開的REAS3.1 (Regional Emission inventory in ASia)為對象，由於早先此一數據庫乃提供區域酸雨研究所用，近年來雖然酸雨議題不再那麼嚴重，該研究所仍然持續發展更新，以應用在地區空氣污染的研究。

本項作業的目標就是將其0.25度解析度之數據庫，轉成網格模式(CAMx、CMAQ等)各層範圍的排放源，重要轉換包括：

經緯度系統之網格 → 正交、直角座標系統
網格範圍解析度 → d1、d2、等指定範圍解析度
污染物質名稱(NMV、PM) → 按反應機制之物質名稱

架構

REAS資料的架構是以污染物_年代來區分壓縮檔，共有12種污染項目(BC、CO2、CO、CH4、N2O、NH3、NMV、NOX、OC、PM10、PM2.5及SO2)， 2010~15共6個年代。

每個壓縮檔內有該項污染物的類別，除了NH3有10個類別、其餘污染項目(以BC為例)有8個來源類別(含點源及加總)。

#kuang@114-32-164-198 /Users/TEDS/REAS3.1/origins
$ ls NH3/2015|cut -c14-
DOMESTIC_2015_0.25x0.25
FERTILIZER_2015_0.25x0.25
INDUSTRY_2015_0.25x0.25
MANURE_MANAGEMENT_2015_0.25x0.25
MISC_2015_0.25x0.25
OTHER_TRANSPORT_2015_0.25x0.25
POWER_PLANTS_NON-POINT_2015_0.25x0.25
POWER_PLANTS_POINT_2015
ROAD_TRANSPORT_2015_0.25x0.25
TOTAL_2015_0.25x0.25
$ ls BC_/2015|cut -c14-
AVIATION_2015_0.25x0.25
DOMESTIC_2015_0.25x0.25
INDUSTRY_2015_0.25x0.25
OTHER_TRANSPORT_2015_0.25x0.25
POWER_PLANTS_NON-POINT_2015_0.25x0.25
POWER_PLANTS_POINT_2015
ROAD_TRANSPORT_2015_0.25x0.25
TOTAL_2015_0.25x0.25

類別架構與其在排放模式中所應用的時間分率特性有關。
粒狀物
- 粒狀物計有BC、OC、PM10及PM2.5共4項。由於粒狀物會按照模擬結果修正排放量，應注意其類別與空間、時間變化。
非甲烷揮發性有機物
- NMH有20項成份(在年代之下有20項成分項目，其下才是來源類別 )，
REAS物質名稱對照
- 有機物20項、加上前述12項無機物(含粒狀物)，合計共32項，
- 由於本次作業是要套用Environ公司的mozart2camx程式，因此需要建立REAS物種名稱與Mozart模式的對照如表：
- (雖然功能是mozart2camx，但事實上ucar不只有mozart模式，還有CAM-chem、WACCM等作業化之全球大氣成份模式模擬，
- 因此要按照實際、選擇正確的對照關係版本，後續轉換時才能進入程式計算。)

kuang@114-32-164-198 /Users/TEDS/REAS3.1/origins
$ cat -n REAS2MOZ.csv
REAS,MOZ_/cluster/src/CAMx/mozart2camx_v3.2.1/src/G2Lconv_CB6r4_CF_from_WACCM_18mar18.EXT,wt
Acetylene,C2H2,26
BC_,bc_a1,4030.
Benzene,BENZENE,78
Butanes,BIGALK,58
CH4,CH4,16
CO2,CO2,44
CO_,CO,28
Ethane,C2H6,30
Ethene,C2H4,28
Ethylene,C2H4,28
Formaldehyde,CH2O,30
Halocarbons,None,1
Internal_Alkenes,BIGENE,28
Ketones,MEK,72.11
N2O,N2O,44
NH3,NH3,17
NMHC,C3H6,58.08
NOX,NO2,46
OC_,pom_a1,8867.
Other_Aldehyde,CH3CHO,44.05
Other_Alkanes,BIGALK,48
Other_Aromatics,TOL,92.14
Others_emissions,BIGENE,84.1608
PM10_,dst_a1,10.0
PM2.5,dst_a2,10.0
Pentanes,BIGALK,72.15
Propane,C3H8,44.1
Propene,C3H6,42.1
SO2,SO2,64
Terminal_Alkenes,BIGENE,28
Toluene,TOLUENE,92.14
Xylenes,XYLENES,106.16

解壓與轉換

解壓

解開REAS提供的各個gz檔
確定所要的年代是否出現
確定解開之文字檔格式是否正常(如下，分別為經緯度及每個月的排放量(噸/月)
注意：檔案內版本與年代沒有更新到。

kuang@114-32-164-198 /Users/TEDS/REAS3.1/origins
$ for i in $(ls gz/*2015*gz);do tar xvfz $i;done
$ head -n20 BC_/2015/REASv3.1_BC__AVIATION_2015_0.25x0.25
10
BC_ emissions on 0.25 degree by 0.25 degree grid
REASv2.1_BC__AVIATION_2008_0.25x0.25
BC_[t/mon],2008,monthly,0.25 degree by 0.25 degree
Aviation (Asian region)
min : 0.1322E-07 max : 0.9134E+01 sum : 0.2213E+04
Format:2F8.2,12E14.7(longitude, latitude, monthly emission value)
* Longitude and latitude are at lower left (southwest) corner of grid cell
Contact:
  tohara@nies.go.jp; kurokawa@acap.asia
   91.50   80.00 0.8274797E-04 0.7740939E-04 0.8274797E-04 0.8007868E-04 0.8274797E-04 0.8007868E-04 0.8274797E-04 0.8274797E-04 0.8007868E-04 0.8274797E-04 0.8007868E-04 0.8274797E-04
   91.75   80.00 0.1498731E-03 0.1402039E-03 0.1498731E-03 0.1450385E-03 0.1498731E-03 0.1450385E-03 0.1498731E-03 0.1498731E-03 0.1450385E-03 0.1498731E-03 0.1450385E-03 0.1498731E-03
  148.00   80.00 0.1219596E-03 0.1140912E-03 0.1219596E-03 0.1180254E-03 0.1219596E-03 0.1180254E-03 0.1219596E-03 0.1219596E-03 0.1180254E-03 0.1219596E-03 0.1180254E-03 0.1219596E-03
  148.25   80.00 0.1219596E-03 0.1140912E-03 0.1219596E-03 0.1180254E-03 0.1219596E-03 0.1180254E-03 0.1219596E-03 0.1219596E-03 0.1180254E-03 0.1219596E-03 0.1180254E-03 0.1219596E-03
  148.50   80.00 0.3895595E-03 0.3644266E-03 0.3895595E-03 0.3769930E-03 0.3895595E-03 0.3769930E-03 0.3895595E-03 0.3895595E-03 0.3769930E-03 0.3895595E-03 0.3769930E-03 0.3895595E-03
  148.75   80.00 0.3895595E-03 0.3644266E-03 0.3895595E-03 0.3769930E-03 0.3895595E-03 0.3769930E-03 0.3895595E-03 0.3895595E-03 0.3769930E-03 0.3895595E-03 0.3769930E-03 0.3895595E-03
  149.00   80.00 0.4312594E-03 0.4034362E-03 0.4312594E-03 0.4173478E-03 0.4312594E-03 0.4173478E-03 0.4312594E-03 0.4312594E-03 0.4173478E-03 0.4312594E-03 0.4173478E-03 0.4312594E-03
  149.25   80.00 0.9504238E-03 0.8891062E-03 0.9504238E-03 0.9197650E-03 0.9504238E-03 0.9197650E-03 0.9504238E-03 0.9504238E-03 0.9197650E-03 0.9504238E-03 0.9197650E-03 0.9504238E-03
  149.50   80.00 0.7156649E-03 0.6694930E-03 0.7156649E-03 0.6925789E-03 0.7156649E-03 0.6925789E-03 0.7156649E-03 0.7156649E-03 0.6925789E-03 0.7156649E-03 0.6925789E-03 0.7156649E-03
  149.75   80.00 0.2447359E-03 0.2289465E-03 0.2447359E-03 0.2368412E-03 0.2447359E-03 0.2368412E-03 0.2447359E-03 0.2447359E-03 0.2368412E-03 0.2447359E-03 0.2368412E-03 0.2447359E-03

將文字檔讀成nc檔(txt2sum.py)

由於未來將使用mozart2camx程式做座標系統、網格系統及檔案格式轉檔(nc/m3ioapi→uamiv)，因此有必要先將文字檔轉成nc檔，將逐月橫向的文字檔，轉成nc檔案的時間軸
應用程式為tx22sum.py，程式碼如下，需要輸入檔為：
- all_file.nam ：為所有待處理之文字檔檔名(line22~24)。讀入後依序處理。此檔之產生只須將ls結果redirect到該檔即可，範例如下。 REAS2MOZ.csv ：污染物質項目對照表，內容如上所述。讀入後以dict方式儲存備用(line25~33)
- cate.csv ：類別之列表，歷來REAS共22個類別(line34~35)，雖然REAS3.1只有12個類別(不含TOTAL與點源)，此處還是讀入此檔，如果該類別沒有檔案，程式會跳過。
- coord.txt ：為各點經緯度座標數值，由各類別網格之聯集求解，程式會自行產生方便閱讀檢查。(line34~74)
- frame.nc ：為mozart模式nc/ioapi檔範例模版
area：
- 每個等經緯度網格範圍的面積(約略值,KM²)(line 88~102,)
- 由於mozart2camx 轉檔過程會考慮網格差異所造成的質量守恆，因此以intensive 型態表示，將經緯度網格排放量除以面積，未來視網格的面積為何，再乘回來即可。(line 195, 201)

#generate the area(km^2) of each grid cells, 0.25 squre = 500~800Km^2 
Latitude_Pole, Longitude_Pole = 23.61000, 120.9900 
pnyc = Proj(proj='lcc', datum='NAD83', lat_1=10, lat_2=40, lat_0=Latitude_Pole, lon_0=Longitude_Pole, x_0=0, y_0=0.0) 
lonG, latG = np.meshgrid(x0, y0) 
x,y=pnyc(lonG, latG, inverse=False) 
x,y=x/1000.,y/1000. 
area=np.zeros(shape=(ny,nx)) 
for j in range(ny-1): 
  for i in range(nx-1): 
    dx=(x[j+1,i+1]+x[j,i+1]-x[j+1,i]-x[j,i])/2. 
    dy=(y[j+1,i+1]+y[j+1,i]-y[j,i+1]-y[j,i])/2. 
    area[j,i]=dx*dy 
  area[j,nx-1]=area[j,nx-2] 
for i in range(nx): 
  area[ny-1,i]=area[ny-2,i]

applying…

    try: 
      f.createVariable(io,dtmap[str(a.variables[i].dtype)],a.variables[i].dimensions) 
    except: #already create, 
      zs=np.array(f.variables[specM][:,:,:,:]) 
      for m in range(nm): 
        f.variables[specM][m]=zs[m]+z[m,0,:,:]/area[:,:]/unit_REAS[spec]/HRS[m] 

    else:       #new to f 
      f.variables[specM].units='gmoleAir/hr/km^2 or AirDens*mg/hr/km^2'#like MOZART in Volume Mix-ratio(to PPM) 
      f.variables[specM].long_name=specM+' monthly emission from REAS' 
      for m in range(nm): 
        z2d=z[m,0,:,:]/area[:,:]/unit_REAS[spec]/HRS[m]  #in unit of mole for splitting and summing 
        f.variables[specM][m]=z2d

單位轉換

由原本的Ton/mon/0.25deg² ->g(mole)/hr/km²
- gas：公式 VMR = 28.9644 / mw * 1e6 * MMR 可以參考ecmwf網站。但mozart2camx 程式中只有乘上1e6
- Particle：粒狀物ug/M3=MMR* 1e9 * (g/M3_air) 。mozart2camx 程式中還有分子量轉換，乘上不同係數，詳aero_fac.txt

facG=10**6 /1E6 
unit_REAS={i:facG*mw for i,mw in zip(df.REAS,df.wt)} 
with open('aero_fac.txt','r') as f: 
  facP={i.split()[0]:float(i.split()[1])/1E6 for i in f}
REAS_part=['BC_','OC_','PM2.5','PM10_'] 
MOZT_part=['CB1','OC1','SA1','DUST1'] 
R2M={i:j for i,j in zip(REAS_part,MOZT_part)}
for i in REAS_part: 
  unit_REAS.update({i:facP[R2M[i]]})

applying…

195:        f.variables[specM][m]=zs[m]+z[m,0,:,:]/area[:,:]/unit_REAS[spec]/HRS[m] 
201:        z2d=z[m,0,:,:]/area[:,:]/unit_REAS[spec]/HRS[m]  #in unit of mole for splitting and summing

讀進REAS txt檔的內容(line 164~169)
結果存成矩陣z(179~186)

  for fname in fnames: 
    if 'POWER_PLANTS_POINT' in fname:continue 
    if cate not in fname: continue 
#read the txt file 
    with open(fname) as text_file: 
      d=[line.strip('\n').split() for line in text_file] 
    print (fname+' '+cate) 
    f1=int(d[0][0]) 
    spec=d[1][0] 
    if spec=='Total':spec='NMHC' 
    if spec=='Total_NMV':spec='NMHC' 
    if spec=='Others':spec=d[1][0]+'_'+d[1][1] 
    specM=specREAS_MOZ[spec] 
    if specM=='None':continue 

    z=np.zeros(shape=(nm,1,ny,nx)) 
    for l in range(f1,len(d)): 
      xx=min(max(float(d[l][0]),xmin),xmax) 
      yy=min(max(float(d[l][1]),ymin),ymax) 
      i=x0.index(xx) 
      j=y0.index(yy) 
      for m in range(12): 
        z[m,0,j,i]=float(d[l][m+2])

mk_rec_dmn ：將時間作為可增加之維度

204    os.system(path[hname]+'/ncks -O --mk_rec_dmn time '+cate+".nc tmp;mv tmp "+cate+".nc")

excution of txt2sum.py

kuang@114-32-164-198 /Users/TEDS/REAS3.1/origins
$ python txt2sum.py 
reading nx,ny,x0,y0... 
BC_/2015/REASv3.1_BC__AVIATION_2015_0.25x0.25 AVIATION 
SO2/2015/REASv3.1_SO2_AVIATION_2015_0.25x0.25 AVIATION 
CO_/2015/REASv3.1_CO__AVIATION_2015_0.25x0.25 AVIATION 
OC_/2015/REASv3.1_OC__AVIATION_2015_0.25x0.25 AVIATION 
PM10_/2015/REASv3.1_PM10__AVIATION_2015_0.25x0.25 AVIATION 
NOX/2015/REASv3.1_NOX_AVIATION_2015_0.25x0.25 AVIATION 
CO2/2015/REASv3.1_CO2_AVIATION_2015_0.25x0.25 AVIATION
...
NMV/2015/20/REASv3.1_NMV_20_WASTE_2015_0.25x0.25 WASTE
(base) 
kuang@114-32-164-198 /Users/TEDS/REAS3.1/origins 
$ lst
-rw-r--r--   1 kuang  AQM    39M Dec 23 11:11 AVIATION.nc 
-rw-r--r--   1 kuang  AQM    51M Dec 23 11:11 DOMESTIC.nc 
-rw-r--r--   1 kuang  AQM    33M Dec 23 11:12 EXTRACTION.nc 
-rw-r--r--   1 kuang  AQM   6.2M Dec 23 11:12 FERTILIZER.nc 
-rw-r--r--   1 kuang  AQM    51M Dec 23 11:12 INDUSTRY.nc 
-rw-r--r--   1 kuang  AQM   6.2M Dec 23 11:12 MANURE_MANAGEMENT.nc 
-rw-r--r--   1 kuang  AQM   6.2M Dec 23 11:12 MISC.nc 
-rw-r--r--   1 kuang  AQM    51M Dec 23 11:12 OTHER_TRANSPORT.nc 
-rw-r--r--   1 kuang  AQM    51M Dec 23 11:12 POWER_PLANTS_NON-POINT.nc 
-rw-r--r--   1 kuang  AQM    51M Dec 23 11:12 ROAD_TRANSPORT.nc 
-rw-r--r--   1 kuang  AQM    33M Dec 23 11:12 SOLVENTS.nc 
-rw-r--r--   1 kuang  AQM    33M Dec 23 11:12 WASTE.nc

txt2sum.py輸入檔案的內容

all_file.nam
cate.csv
coord.txt
aero_fac.txt

$ head all_file.nam
BC_/2015/REASv3.1_BC__DOMESTIC_2015_0.25x0.25
BC_/2015/REASv3.1_BC__INDUSTRY_2015_0.25x0.25
...
CO2/2015/REASv3.1_CO2_POWER_PLANTS_NON-POINT_2015_0.25x0.25
$ head cate.csv
cate
AVIATION
DOMESTIC
ENTERIC_FERMENTATION
EXTRACTION
FERTILIZER
FUGITIVE_COAL
FUGITIVE_GAS
FUGITIVE_OIL
INDUSTRY
$ wc cate.csv
      23      23     289 cate.csv
$ cat coord.txt
241 185
91.0 91.25 91.5 ... 150.75 151.0
0.0 0.25 0.5 ... 45.75 46.0 (base)
$ cat -n aero_fac.txt 
     1   SO4                3.86916019E+09 
     2   NH4NO3             2.49883264E+09 
     3   NH4                725467584. 
     4   SOA                7.25467546E+09 
     5   SOA                7.25467546E+09 
     6   OC1                8.86682624E+09 
     7   OC2                8.86682624E+09 
     8   CB1                4.03037542E+09 
     9   CB2                4.03037542E+09 
    10   DUST1              4.03037542E+09 
    11   DUST2              4.03037542E+09 
    12   DUST3              4.03037542E+09 
    13   DUST4              4.03037542E+09 
    14   SA1                926986304. 
    15   SA2                926986304. 
    16   SA1                1.41063130E+09 
    17   SA2                1.41063130E+09

轉換為m3.nc檔案

程式又稱為NCF2IOAPI。此舉將mozart等一般的nc檔案(前述的模版frame.nc)，轉成USEPA的標準格式(model3/IOAPI檔案，可以用verdi開啟)。其後再進行切割或類別間的合併。

NCF2IOAPI不會改變檔案的格點數、網格系統，主要改變的是一般性的變數的名稱，如原點是nc.XORIG、間距是nc.XCELL、時間標籤是TFLAG等，IOAPI有自己成套的convention。

Environ目前尚無刪減或合併此步驟的計畫打算。

EXE ：為執行檔，可以由Environ公司下載、編譯、使用。
OUTFILE3D ：為結果檔案，在隔壁join_spec目錄下。

kuang@114-32-164-198 /Users/TEDS/REAS3.1/origins
$ cat -n nc2m3.cs
export EXECUTION_ID=mz2camx.job
export PROMPTFLAG=N
export IOAPI_ISPH=20
     4
EXE=/Users/camxruns/src/mozart2camx_v3.2.1/ncf2ioapi_mozart/NCF2IOAPI
for i in $(ls [ADEFIMOPRSW]*.nc|cut -d'.' -f1) ;do
  export INFILE=$i.nc
  export OUTFILE3D=../join_spec/$i".m3.nc"
  echo $i
  $EXE|tail -n5
done

來源數據檢核

1月份AVIATION的SO2排放量，網格數、範圍仍然為REAS原始狀態，尚未切割或整併。

直角座標排放檔之切割合併

主控腳本add.cs

mz2camx.job ：
- 主要呼叫的腳本，輸入年代月份($YYMM )、巢狀層數($D )以及m3.nc檔案名稱($nc )(line 13)
addavrg2 ：
- uamiv格式檔案相加之小程式
- 分別將點、線、面源的污染源予以加總
multavrg ：
- uamiv格式乘上scalar($mult)的小程式，會覆蓋原有檔案
- scalar包括該層網格的的面積(單位為KM²)
- 結果為gmole/月
add_ttime.f加入時間標籤(僅24小時，可供非日期鑑別、典型日變化之CAMx模擬)
重複執行add.cs即可完成全年轉檔
注意：
- 目錄下如已經有成果檔案，程式將會跳過不覆蓋
- 乘上乘數後原檔會被覆蓋掉，因此只能從新來過才能檢查確認。
- 乘數最後結果，單位為gmole/month(CAMx格式)
- macOS的數字檢查很嚴格，數字和文字不能混用，因此不會計算{01..09}，會自動略去前面的0。

kuang@114-32-164-198 /Users/TEDS/REAS3.1/join_spec
$ cat -n add.cs
#performed the nz2camx.job
for YYMM in 150{1..9} 151{0..2};do
#YYMM=$1
j0=`date -v+0m -j -f "%Y%m%d"  "20${YYMM}01" +%j`
j1=`date -v+1m -j -f "%Y%m%d"  "20${YYMM}01" +%j`
mnday= $(( 10#$j1 - 10#$j0 ))
if [ $mnday -lt 0 ];then
  j2=`date -j -f "%Y%m%d"  "20${YYMM}31" +%j`
  mnday=$(( $mnday + 10#$j2 ))
fi
for D in d1 d2;do
  for nc in $(ls *m3.nc);do
    mz2camx.job $YYMM $D $nc
  done
    15
  #sum_up the REAS emission from different categories into area,line,avi, and ind
  ln $YYMM$D.AVIATION $YYMM$D.avi
    18
  cp $YYMM$D.DOMESTIC $YYMM$D.area
  for nc in ENTERIC_FERMENTATION EXTRACTION FERTILIZER INTNNV \
    MANURE_MANAGEMENT SOIL_INDIRECT WASTE SOIL_DIRECT SOIL RICE_CULTIVATION MISC;do
    if [ -e $YYMM$D.$nc ];then addavrg2 $YYMM$D.$nc $YYMM$D.area tmp;mv tmp $YYMM$D.area;fi
  done
    24
  cp $YYMM$D.INDUSTRY $YYMM$D.ind
  for nc in FUGITIVE_COAL FUGITIVE_GAS FUGITIVE_OIL POWER_PLANTS_NON-POINT_JPN \
    POWER_PLANTS_NON-POINT SOLVENTS;do
    if [ -e $YYMM$D.$nc ];then addavrg2 $YYMM$D.$nc $YYMM$D.ind tmp;mv tmp $YYMM$D.ind;fi
  done
    30
  cp $YYMM$D.OTHER_TRANSPORT $YYMM$D.line;addavrg2 $YYMM$D.ROAD_TRANSPORT $YYMM$D.line tmp;mv tmp $YYMM$D.line
    32
  cp $YYMM$D.line $YYMM$D.tot;for nc in area avi ind;do addavrg2 $YYMM$D.$nc $YYMM$D.tot tmp;mv tmp $YYMM$D.tot;done
    34
  km2=$(( 81 * 81 ))
  test $D == 'd2' && km2=$(( 27 * 27 ))
  mult=`echo $km2` # mult=`echo $km2/$mnday/24 |bc -l`
  for nc in area avi ind line tot;do multavrg $YYMM$D.$nc $mult;mv $YYMM$D.$nc$mult $YYMM$D.$nc
  done
  done
  
  #cbin all month
  for D in d1 d2;do
    for nc in line area avi ind tot;do
      cbin_all "15*$D.$nc" $D.$nc
      pncgen --format=uamiv -O --out-format=uamiv --slice=LAY,0 -a NAME,global,o,c,'EMISSIONS ' $D.$nc ${D}L.$nc>&/dev/null
      for m in 0{1..9} 1{0..2};do add_ttime ${D}L.$nc $m;done
    done
  done

核心腳本mz2camx.job

EXE ：此腳本為c-shell指令，重點在執行mozart2camx程式(line 6)，注意與前述REAS物質名稱對照關係，必須是同一個設定方式。
MET ：主要用來抓垂直網格的設定。由於mozart系統的時間都是UTC，因此MET也必須是UTC(time zone grid 設為空白) (line 16)
DATE :雖然是每隔1個月會增加1個數字，但是m3.nc是嚴格的等間距檔案，此處設為30日(720小時)，反應在${ND}
1. macOS/centos的date略有差異，以下以macOS為主，(內為centos版本)

    13  @ MM = $MM - 1 
    14  @ ND = $MM * 30
...
    24  set DATE = `date -v+${ND}d -j -f "%Y%m%d"  "20150101" +%Y%m%d`
        (set DATE = `date -d "20150101 +${ND}day" +%Y%m%d`)
...   
    40  set YYYYMMDD = $DATE
...
    43  ProcessDateYYYYMMDD|$YYYYMMDD

kuang@114-32-164-198 /Users/TEDS/REAS3.1/join_spec
$ cat -n mz2camx.job
#!/bin/csh -f
#$1->yymm
#$2->m3.nc
setenv PROMPTFLAG N
setenv IOAPI_ISPH 20
set EXE = /Users/camxruns/src/mozart2camx_v3.2.1/src/mozart2camx_CB6r4_CF__WACCM
     7
set YYMM = $1
set d = $2
set OU = `echo $3|cut -d'.' -f1`
set MM = `date -j -f "%Y%m%d"  "20${YYMM}01" +%m`
set YY = `echo $YYMM|cut -c 1-2`
@ MM = $MM - 1
@ ND = $MM * 30
set MET = /Users/WRF4.1/201601/wrfout/1601${d}
    16
# DEFINE OUTPUT FILE NAMES
setenv EXECUTION_ID mz2camx.job
foreach i ( 0 )#1 2 3 )
foreach j ( 1 )#1 2 3 4 5 6 7 8 9 )
if ( $i == '3' && $j > '1' ) goto BYPASS
set k = $i$j
if ( $k == '00' ) goto BYPASS
set DATE = `date -v+${ND}d -j -f "%Y%m%d"  "20150101" +%Y%m%d`
# DEFINE INPUT MOZART FILES
# IF MORE THAN 1 MOZART FILE IS NEEDED, ADD setenv INFILE2
set NINFILE = 1
    28
foreach t ( 00 )#06 12 18 )
setenv OUTFILEIC  $YYMM$d"."$OU
setenv OUTFILEBC  tmp
foreach INFILE ($3)# m3/$2 m3/$3 m3/$4 m3/$5 m3/$6 )
setenv INFILE1 $INFILE
set fs = `ls -l "$OUTFILEIC"|awk '{print $5}'`
if ( $fs > '10000000' ) goto BYPASS2
echo $OUTFILEIC
rm -f $OUTFILEBC $OUTFILEIC
    38
echo $DATE$t
set YYYYMMDD = $DATE
$EXE << IEOF
CAMx5,CAMx6,CMAQ   |CAMx 6
ProcessDateYYYYMMDD|$YYYYMMDD
Output BC file?    |.false.
Output IC file?    |.true.
If IC, starting hr |$t
Output TC file?    |.false.
Max num MZRT files |$NINFILE
CAMx 3D met file   |$MET.3d
CAMx 2D met file   |$MET.2d
IEOF
mv OUTFILEIC $OUTFILEIC
echo $INFILE1
BYPASS2:
end
end
BYPASS:
end
end
exit 0

粒狀物的單位轉換

氣狀物在mozart2過程中會乘上28/Mw*10⁶將mass mixing ration 轉換成PPM，因此：
- 在txt2sum.py中事先除以facG轉成莫耳數及單位面積。
- 在add.cs進行面積(網格面積)的轉換
粒狀物，在mozart2camx過程中會乘上10⁹、空氣密度與分子量比例，將mass mixing ration 轉換成ug/M3
- 各項粒狀物之乘數如下表
- 其中OC、CB、及DUST項目，存在REAS排放清單中，因此必須先將其除回來，以避免多乘。 - 效果等同排放量除以分子量，因此OC之除數為8867、CB、及DUST除數為4030 - 其餘程序與氣狀物相同

!kuang@114-32-164-198 /Users/camxruns/src/mozart2camx_v3.2.1/src
!$ vi mozart2camx.f
!           CONVERT TO PPM OR UG/M3 
           IF ( MZRT_CMAQ_MAP(VAR) .LE. N_CMAQ_GAS_SPC ) THEN 
           ! convert gas species into ppmV 
             concgrd_buf = concgrd_buf * 10**6 
             tconcgrd_buf = tconcgrd_buf * 10**6 
           ELSE 
           ! convert aerosol species into micrograms/m**3 
             concgrd_buf = concgrd_buf * rho_grd_buf * 10**9 * 
    &                MWvar(MZRT_CMAQ_MAP(VAR))/MWair 
             tconcgrd_buf = tconcgrd_buf * trho_grd_buf * 10**9 * 
    &                MWvar(MZRT_CMAQ_MAP(VAR))/MWair 
             print *,'kuang',MOZART_SPECIES(VAR), 
    &         rho_grd_buf(1,1,1)* 10**9 * MWvar(MZRT_CMAQ_MAP(VAR))/MWair, 
    &         trho_grd_buf(1,1) * 10**9 * MWvar(MZRT_CMAQ_MAP(VAR))/MWair 
           ENDIF

kuang@114-32-164-198 /Users/TEDS/REAS3.1/origins
$ cat -n aero_fac.txt
 SO4                3.86916019E+09
 NH4NO3             2.49883264E+09
 NH4                725467584.
 SOA                7.25467546E+09
 SOA                7.25467546E+09
 OC1                8.86682624E+09
 OC2                8.86682624E+09
 CB1                4.03037542E+09
 CB2                4.03037542E+09
 DUST1              4.03037542E+09
 DUST2              4.03037542E+09
 DUST3              4.03037542E+09
 DUST4              4.03037542E+09
 SA1                926986304.
 SA2                926986304.
 SA1                1.41063130E+09
 SA2                1.41063130E+09

add_ttime.f
- 從月均值延長成24小時之逐時值(CAMx可以接受典型日變化之日排放檔)

面源時間變化(mkMon.py’s)

前述面源排放量為月均值，如何處理成模式所需的逐時值有2種做法，此處以溫度為時間變化之指標，溫度低者排放量大。另，以常數方式進行展開，詳下述。

氣溫反比變化

家戶加熱、電廠、這些與冬季排放有關，在此以WRF計算網格月平均氣溫倒數為X陣列，以污染排放為Y陣列，進行線性回歸，藉此估算逐時的排放量。架構步驟

全年逐時T2擷取(T2.nc )
- 以ncks由逐日檔案擷取T2逐時值
- 逐日檔案合併(ncrcat)成全年逐時T2
全年月均值之氣溫(wrfout_d01_2016_monthly_mean)
- 在每個月的目錄下另創wrfout目錄，在OS計算
- 使用ln_wrfout.cs腳本在2016MM/wrfout目錄下做好連結。
- 3/5/7/9/11月初缺值，須另外從前月run12最後1～2天連結過來。
- 以ncks擷取WRF地面氣溫(T2)
- 逐日檔案合併成月檔案
- 計算月均值(tmNC) - 合併12個月均值，併入前述某排放量檔案之中，儲存備用。
- 直接在python程式內計算(用for…if…條件篩選)

   beg=datetime.datetime(2015,12,31)
   wrf_time=[beg+datetime.timedelta(days=t/24.) for t in range(ntw)]
   T2Mi=np.zeros(shape=(12,nroww,ncolw))
   for m in range(12):
     mdate=np.array([t for t in range(ntw) if wrf_time[t].month==m+1])
     T2Mi[m,:,:]=np.mean(T2[mdate,:,:],axis=0)

CMAQ排放檔案模板之準備
- 修改d1.tot等檔案的TFLAG值，camx2cmaq不會認TSTEP=10000以外的數值（此處為30日=720）
- 執行camx2cmaqd1.job進行12筆排放量之轉換
- 將此12筆檔案複製(ncrcat)成至少744筆之模版
計算網格月均值之回歸公式（ems = a/T2 + b )
- 填入TFLAGs
- 填入逐時之估計值

月均溫與逐時氣溫檔案之預備

全年月均值之氣溫

for i in {01..12};do cd 2016$i/wrfout;for j in $(ls wrfout_d01_2016-${i}-??_00\:00\:00);do ncks -O -d bottom_top,0,0 $j $j.k1;done;ncrcat -O wrfout_d01_2016-${i}*.k1 wrfout_d01_2016-${i};cd ../../;done
for i in {01..12};do python ~/bin/tmNC 2016$i/wrfout/wrfout_d01_2016-$i;done
ncrcat 2016??/wrfout/*T wrfout_d01_2016_monthly_mean

WRF地面逐時氣溫(T2)之擷取

for nc in $(ls */*/*k1);do ncks -v Times,T2 ${nc}_T2 ;done
ncrcat */*/*k1_T2 T2.nc

mkMon2.py(詳codelisting)

模板為fortBE.213.teds10.base00.nc(camx2cmaq結果)
mkMon3.py(詳github)
模版為REAS2hr結果，數據為add.cs結果，格式為uamiv，
面源時間變化(REAS2hr’s)
本項作業同樣是給予面源檔案有時間變化，但著重在小時變化，而非前述著眼在日均溫度之月變化。(經證明對臭氧有更高的敏感性)
由於REAS -> MOZART -> CAMx -> CMAQ 之作業順序(相較前述REAS -> MOZART -> CAMx -> CMAQ ->CAMx)有較高的可重製特性，因此在第一次CAMx階段即給予日變化，會較方便。

uamiv格式版本(REAS2hr.py)

輸入檔案argment: 月份(2 digit)uamiv 各類別各月份d1、d2檔案(前述mozart2camx轉檔與合併結果)fortBE.?13樣版
輸出檔案，在emis/$CATE目錄下，生成fortBE.DDD_XXXXX.baseMM檔案。DDD=113~213、XXXXX=REAS3/STEAM、MM=01~12。由前述步驟所產生。
時間變化：
- 線源：將仿照fortBE.?13L之變化。
- 船舶：按照STEAM有逐日變化
- 其他：無變化
此版本直接處理uamiv格式檔案。
- 沿用原uamiv時間標籤
- 由pncgen轉成nc檔之後，還需要再加上必須的全域屬性，因此放棄不持續發展，改由REAS2hr_TimVar.py為之。

add_ncatt.cs

- REAS2hr.py完成後所生檔案為uamiv格式，再利用[pncgen][pncgen]將uamiv檔案成為nc檔之後，仍須加入必要之全域屬性，方能進行空品模擬。

先產生nc檔案再轉成uamiv格式(REAS2hr_TimVar.py

本程式與前述REAS2hr.py功能相同，然以nc檔版模為主體，而非以uamiv模版，
- 重新編寫時間標籤以符合CAMx及CMAQ模式所需
- 如此可保有完整屬性訊息，(未來)不必再轉nc檔
如要輸出nc檔案，在程式最後直接另存檔名即可，不必再進行pncgen。

面源/工業源之時間變化

由於面源/工業源之硫氧化物與細懸浮微粒排放對模擬結果具有敏感性，因此，在此以逐時氣溫做為校正。

船舶d4部分程式(STEAMd4ToHr.py)

d4的排放量主要是TEDS，船舶部分因STEAM除了近海之外，還包括公海部分，因此在邊界附近會有較高的連續性。此程式只有進行ship目錄下之新增，並沒有做area~ind等目錄。

kuang@master /nas1/camxruns/2016_v7/emis
$ diff REAS2hr.py STEAMd4ToHr.py
6,7c6,7
< for d in ['1', '2']:
<   for cate in ['area', 'avi', 'ind','line', 'ship']:
---
> for d in ['4']:
>   for cate in ['ship']:

執行批次

kuang@master /nas1/camxruns/2016_v7/emis
$ cat prep_REAS.cs
for m in {01..12};do for d in 1 2;do for i in line area avi ind;do cp fortBE.${d}13L_base$m $i/fortBE.${d}13_REAS3.base$m;done;done;done
for m in {01..12};do for d in 1 2 4;do for i in ship;do cp fortBE.${d}13L_base$m $i/fortBE.${d}13_STEAM.base$m;done;done;done
for i in {01..12};do sub python REAS2hr_TimVar.py $i >& prep_REAS.log;done
for i in {01..12};do sub python STEAMd4ToHr.py $i >& prep_STEAM.log;done

成果檢視

空間趨勢

除了tot外，尚有面源、線源、工業源、以及航空等分類。BENZ似較多來自線源，SO2則來自重工業與海上運輸。

船舶排放(http://www.evernote.com/l/AH36Psxby9dKcJ9kdRyX6azfz6m_b4Q-O3M) vs

area+avi+ind+line sources vs REAS2

d2範圍SO2海上比例更大，NH3則以陸地為主，台灣中南部較大，排放量解析度有限。

時間變化

氨具有溫度趨勢，夏季明顯高過冬季。CO則具有冬季取暖排放較大。NO2、CCRS與交通有關，變化不似CO。

Code Listing

Download: REAS文字檔讀取程式txt2sum.py

Reference

Resources

亞洲區域排放清冊 https://www.nies.go.jp/REAS/
mozart2camx http://www.camx.com/getmedia/028fbd41-c80c-49fd-aa05-8abb3fae261f/mozart2camx-26feb19_1.tgz
CAMx(UAM)的檔案格式, Yungchuan Kuang edited this page on 12 Jul 2016 · 2 revision, shttps://github.com/sinotec2/camxruns/wiki/CAMx(UAM)的檔案格式
VERDI:
- verdi usage https://www.airqualitymodeling.org/index.php/VERDI_1.5_User_Manual#3.1_Installation_Instructions_for_Linux_and_Mac
- VERDI使用說明 : http://www.evernote.com/l/AH3leuVQTuBEF7Vrs0D1C8Q-Iff5CpHl7eU
- VERDI使用說明 v2： http://www.evernote.com/l/AH2gBcV7qsJFr4E7jSW5x4A0FCoW4QW7otE Linear Regression in Python https://realpython.com/linear-regression-in-python/