Dr. Kuang's Utilities

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

 
  • Jul 06, 2022
  • views
>

背景

目標

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

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

  1. 經緯度系統之網格 → 正交、直角座標系統
  2. 網格範圍解析度 → d1、d2、等指定範圍解析度
  3. 污染物質名稱(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
  1 REAS,MOZ_/cluster/src/CAMx/mozart2camx_v3.2.1/src/G2Lconv_CB6r4_CF_from_WACCM_18mar18.EXT,wt
  2 Acetylene,C2H2,26
  3 BC_,bc_a1,4030.
  4 Benzene,BENZENE,78
  5 Butanes,BIGALK,58
  6 CH4,CH4,16
  7 CO2,CO2,44
  8 CO_,CO,28
  9 Ethane,C2H6,30
 10 Ethene,C2H4,28
 11 Ethylene,C2H4,28
 12 Formaldehyde,CH2O,30
 13 Halocarbons,None,1
 14 Internal_Alkenes,BIGENE,28
 15 Ketones,MEK,72.11
 16 N2O,N2O,44
 17 NH3,NH3,17
 18 NMHC,C3H6,58.08
 19 NOX,NO2,46
 20 OC_,pom_a1,8867.
 21 Other_Aldehyde,CH3CHO,44.05
 22 Other_Alkanes,BIGALK,48
 23 Other_Aromatics,TOL,92.14
 24 Others_emissions,BIGENE,84.1608
 25 PM10_,dst_a1,10.0
 26 PM2.5,dst_a2,10.0
 27 Pentanes,BIGALK,72.15
 28 Propane,C3H8,44.1
 29 Propene,C3H6,42.1
 30 SO2,SO2,64
 31 Terminal_Alkenes,BIGENE,28
 32 Toluene,TOLUENE,92.14
 33 Xylenes,XYLENES,106.16

解壓與轉換

解壓

  1. 解開REAS提供的各個gz檔
  2. 確定所要的年代是否出現 確定解開之文字檔格式是否正常(如下，分別為經緯度及每個月的排放量(噸/月)
    • 注意：檔案內版本與年代沒有更新到。
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][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：
    • 每個等經緯度網格範圍的面積(約略值,KM2)(line 88~102,)
    • 由於mozart2camx 轉檔過程會考慮網格差異所造成的質量守恆，因此以intensive 型態表示，將經緯度網格排放量除以面積，未來視網格的面積為何，再乘回來即可。(line 195, 201)
    88  #generate the area(km^2) of each grid cells, 0.25 squre = 500~800Km^2 
    89  Latitude_Pole, Longitude_Pole = 23.61000, 120.9900 
    90  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) 
    91  lonG, latG = np.meshgrid(x0, y0) 
    92  x,y=pnyc(lonG, latG, inverse=False) 
    93  x,y=x/1000.,y/1000. 
    94  area=np.zeros(shape=(ny,nx)) 
    95  for j in range(ny-1): 
    96    for i in range(nx-1): 
    97      dx=(x[j+1,i+1]+x[j,i+1]-x[j+1,i]-x[j,i])/2. 
    98      dy=(y[j+1,i+1]+y[j+1,i]-y[j,i+1]-y[j,i])/2. 
    99      area[j,i]=dx*dy 
   100    area[j,nx-1]=area[j,nx-2] 
   101  for i in range(nx): 
   102    area[ny-1,i]=area[ny-2,i]

applying…

   190      try: 
   191        f.createVariable(io,dtmap[str(a.variables[i].dtype)],a.variables[i].dimensions) 
   192      except: #already create, 
   193        zs=np.array(f.variables[specM][:,:,:,:]) 
   194        for m in range(nm): 
   195          f.variables[specM][m]=zs[m]+z[m,0,:,:]/area[:,:]/unit_REAS[spec]/HRS[m] 
   196 
   197      else:       #new to f 
   198        f.variables[specM].units='gmoleAir/hr/km^2 or AirDens*mg/hr/km^2'#like MOZART in Volume Mix-ratio(to PPM) 
   199        f.variables[specM].long_name=specM+' monthly emission from REAS' 
   200        for m in range(nm): 
   201          z2d=z[m,0,:,:]/area[:,:]/unit_REAS[spec]/HRS[m]  #in unit of mole for splitting and summing 
   202          f.variables[specM][m]=z2d

單位轉換

  • 由原本的Ton/mon/0.25deg2 ->g(mole)/hr/km2
  • gas：公式 VMR = 28.9644 / mw * 1e6 * MMR 可以參考ecmwf網站。但mozart2camx 程式中只有乘上1e6 Particle：粒狀物ug/M3=MMR* 1e9 * (g/M3_air) 。mozart2camx 程式中還有分子量轉換，乘上不同係數，詳aero_fac.txt
    35  facG=10**6 /1E6 
    36  unit_REAS={i:facG*mw for i,mw in zip(df.REAS,df.wt)} 
    37  with open('aero_fac.txt','r') as f: 
    38    facP={i.split()[0]:float(i.split()[1])/1E6 for i in f}
    39  REAS_part=['BC_','OC_','PM2.5','PM10_'] 
    40  MOZT_part=['CB1','OC1','SA1','DUST1'] 
    41  R2M={i:j for i,j in zip(REAS_part,MOZT_part)}
    42  for i in REAS_part: 
    43    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)
   164    for fname in fnames: 
   165      if 'POWER_PLANTS_POINT' in fname:continue 
   166      if cate not in fname: continue 
   167  #read the txt file 
   168      with open(fname) as text_file: 
   169        d=[line.strip('\n').split() for line in text_file] 
   170      print (fname+' '+cate) 
   171      f1=int(d[0][0]) 
   172      spec=d[1][0] 
   173      if spec=='Total':spec='NMHC' 
   174      if spec=='Total_NMV':spec='NMHC' 
   175      if spec=='Others':spec=d[1][0]+'_'+d[1][1] 
   176      specM=specREAS_MOZ[spec] 
   177      if specM=='None':continue 
   178 
   179      z=np.zeros(shape=(nm,1,ny,nx)) 
   180      for l in range(f1,len(d)): 
   181        xx=min(max(float(d[l][0]),xmin),xmax) 
   182        yy=min(max(float(d[l][1]),ymin),ymax) 
   183        i=x0.index(xx) 
   184        j=y0.index(yy) 
   185        for m in range(12): 
   186          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][NCF2IOAPI]。此舉將mozart等一般的nc檔案(前述的模版frame.nc)，轉成USEPA的標準格式(model3/IOAPI檔案，可以用verdi開啟)。其後再進行切割或類別間的合併。

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

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

  1. EXE ：為執行檔，可以由Environ公司下載、編譯、使用。
  2. OUTFILE3D ：為結果檔案，在隔壁join_spec目錄下。
kuang@114-32-164-198 /Users/TEDS/REAS3.1/origins
$ cat -n nc2m3.cs
     1  export EXECUTION_ID=mz2camx.job
     2  export PROMPTFLAG=N
     3  export IOAPI_ISPH=20
     4
     5  EXE=/Users/camxruns/src/mozart2camx_v3.2.1/ncf2ioapi_mozart/NCF2IOAPI
     6  for i in $(ls [ADEFIMOPRSW]*.nc|cut -d'.' -f1) ;do
     7    export INFILE=$i.nc
     8    export OUTFILE3D=../join_spec/$i".m3.nc"
     9    echo $i
    10    $EXE|tail -n5
    11  done

來源數據檢核

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

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

主控腳本add.cs

  • mz2camx.job ：
    • 主要呼叫的腳本，輸入年代月份($YYMM )、巢狀層數($D )以及m3.nc檔案名稱($nc )(line 13)
  • addavrg2 ：
    • uamiv格式檔案相加之小程式
    • 分別將點、線、面源的污染源予以加總
    • 點源（地面部分）處理詳參[[2022-07-06-REASPtRd]]
  • multavrg ：
    • uamiv格式乘上scalar($mult)的小程式，會覆蓋原有檔案
    • scalar包括該層網格的的面積(單位為KM2)
    • 結果為gmole/月
  • [add_ttime.f][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
     1  #performed the nz2camx.job
     2  for YYMM in 150{1..9} 151{0..2};do
     3  #YYMM=$1
     4  j0=`date -v+0m -j -f "%Y%m%d"  "20${YYMM}01" +%j`
     5  j1=`date -v+1m -j -f "%Y%m%d"  "20${YYMM}01" +%j`
     6  mnday= $(( 10#$j1 - 10#$j0 ))
     7  if [ $mnday -lt 0 ];then
     8    j2=`date -j -f "%Y%m%d"  "20${YYMM}31" +%j`
     9    mnday=$(( $mnday + 10#$j2 ))
    10  fi
    11  for D in d1 d2;do
    12    for nc in $(ls *m3.nc);do
    13      mz2camx.job $YYMM $D $nc
    14    done
    15
    16    #sum_up the REAS emission from different categories into area,line,avi, and ind
    17    ln $YYMM$D.AVIATION $YYMM$D.avi
    18
    19    cp $YYMM$D.DOMESTIC $YYMM$D.area
    20    for nc in ENTERIC_FERMENTATION EXTRACTION FERTILIZER INTNNV \
    21      MANURE_MANAGEMENT SOIL_INDIRECT WASTE SOIL_DIRECT SOIL RICE_CULTIVATION MISC;do
    22      if [ -e $YYMM$D.$nc ];then addavrg2 $YYMM$D.$nc $YYMM$D.area tmp;mv tmp $YYMM$D.area;fi
    23    done
    24
    25    cp $YYMM$D.INDUSTRY $YYMM$D.ind
    26    for nc in FUGITIVE_COAL FUGITIVE_GAS FUGITIVE_OIL POWER_PLANTS_NON-POINT_JPN \
    27      POWER_PLANTS_NON-POINT SOLVENTS;do
    28      if [ -e $YYMM$D.$nc ];then addavrg2 $YYMM$D.$nc $YYMM$D.ind tmp;mv tmp $YYMM$D.ind;fi
    29    done
    30
    31    cp $YYMM$D.OTHER_TRANSPORT $YYMM$D.line;addavrg2 $YYMM$D.ROAD_TRANSPORT $YYMM$D.line tmp;mv tmp $YYMM$D.line
    32
    33    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
    35    km2=$(( 81 * 81 ))
    36    test $D == 'd2' && km2=$(( 27 * 27 ))
    37    mult=`echo $km2` # mult=`echo $km2/$mnday/24 |bc -l`
    38    for nc in area avi ind line tot;do multavrg $YYMM$D.$nc $mult;mv $YYMM$D.$nc$mult $YYMM$D.$nc
    39    done
    40    done
    41    
    42    #cbin all month
    43    for D in d1 d2;do
    44      for nc in line area avi ind tot;do
    45        cbin_all "15*$D.$nc" $D.$nc
    46        pncgen --format=uamiv -O --out-format=uamiv --slice=LAY,0 -a NAME,global,o,c,'EMISSIONS ' $D.$nc ${D}L.$nc>&/dev/null
    47        for m in 0{1..9} 1{0..2};do add_ttime ${D}L.$nc $m;done
    48      done
    49    done

核心腳本mz2camx.job

  1. EXE ：此腳本為c-shell指令，重點在執行mozart2camx程式(line 6)，注意與前述REAS物質名稱對照關係，必須是同一個設定方式。
  2. MET ：主要用來抓垂直網格的設定。由於mozart系統的時間都是UTC，因此MET也必須是UTC(time zone grid 設為空白) (line 16)
  3. 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
     1  #!/bin/csh -f
     2  #$1->yymm
     3  #$2->m3.nc
     4  setenv PROMPTFLAG N
     5  setenv IOAPI_ISPH 20
     6  set EXE = /Users/camxruns/src/mozart2camx_v3.2.1/src/mozart2camx_CB6r4_CF__WACCM
     7
     8  set YYMM = $1
     9  set d = $2
    10  set OU = `echo $3|cut -d'.' -f1`
    11  set MM = `date -j -f "%Y%m%d"  "20${YYMM}01" +%m`
    12  set YY = `echo $YYMM|cut -c 1-2`
    13  @ MM = $MM - 1
    14  @ ND = $MM * 30
    15  set MET = /Users/WRF4.1/201601/wrfout/1601${d}
    16
    17  # DEFINE OUTPUT FILE NAMES
    18  setenv EXECUTION_ID mz2camx.job
    19  foreach i ( 0 )#1 2 3 )
    20  foreach j ( 1 )#1 2 3 4 5 6 7 8 9 )
    21  if ( $i == '3' && $j > '1' ) goto BYPASS
    22  set k = $i$j
    23  if ( $k == '00' ) goto BYPASS
    24  set DATE = `date -v+${ND}d -j -f "%Y%m%d"  "20150101" +%Y%m%d`
    25  # DEFINE INPUT MOZART FILES
    26  # IF MORE THAN 1 MOZART FILE IS NEEDED, ADD setenv INFILE2
    27  set NINFILE = 1
    28
    29  foreach t ( 00 )#06 12 18 )
    30  setenv OUTFILEIC  $YYMM$d"."$OU
    31  setenv OUTFILEBC  tmp
    32  foreach INFILE ($3)# m3/$2 m3/$3 m3/$4 m3/$5 m3/$6 )
    33  setenv INFILE1 $INFILE
    34  set fs = `ls -l "$OUTFILEIC"|awk '{print $5}'`
    35  if ( $fs > '10000000' ) goto BYPASS2
    36  echo $OUTFILEIC
    37  rm -f $OUTFILEBC $OUTFILEIC
    38
    39  echo $DATE$t
    40  set YYYYMMDD = $DATE
    41  $EXE << IEOF
    42  CAMx5,CAMx6,CMAQ   |CAMx 6
    43  ProcessDateYYYYMMDD|$YYYYMMDD
    44  Output BC file?    |.false.
    45  Output IC file?    |.true.
    46  If IC, starting hr |$t
    47  Output TC file?    |.false.
    48  Max num MZRT files |$NINFILE
    49  CAMx 3D met file   |$MET.3d
    50  CAMx 2D met file   |$MET.2d
    51  IEOF
    52  mv OUTFILEIC $OUTFILEIC
    53  echo $INFILE1
    54  BYPASS2:
    55  end
    56  end
    57  BYPASS:
    58  end
    59  end
    60  exit 0

粒狀物的單位轉換

  1. 氣狀物在mozart2過程中會乘上28/Mw*106將mass mixing ration 轉換成PPM，因此：
    • 在txt2sum.py中事先除以facG轉成莫耳數及單位面積。
    • 在add.cs進行面積(網格面積)的轉換
  2. 粒狀物，在mozart2camx過程中會乘上109、空氣密度與分子量比例，將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
1853 !           CONVERT TO PPM OR UG/M3 
1854             IF ( MZRT_CMAQ_MAP(VAR) .LE. N_CMAQ_GAS_SPC ) THEN 
1855             ! convert gas species into ppmV 
1856               concgrd_buf = concgrd_buf * 10**6 
1857               tconcgrd_buf = tconcgrd_buf * 10**6 
1858             ELSE 
1859             ! convert aerosol species into micrograms/m**3 
1860               concgrd_buf = concgrd_buf * rho_grd_buf * 10**9 * 
1861      &                MWvar(MZRT_CMAQ_MAP(VAR))/MWair 
1862               tconcgrd_buf = tconcgrd_buf * trho_grd_buf * 10**9 * 
1863      &                MWvar(MZRT_CMAQ_MAP(VAR))/MWair 
1864               print *,'kuang',MOZART_SPECIES(VAR), 
1865      &         rho_grd_buf(1,1,1)* 10**9 * MWvar(MZRT_CMAQ_MAP(VAR))/MWair, 
1866      &         trho_grd_buf(1,1) * 10**9 * MWvar(MZRT_CMAQ_MAP(VAR))/MWair 
1867             ENDIF

kuang@114-32-164-198 /Users/TEDS/REAS3.1/origins
$ 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
  • [add_ttime.f][add_ttime.f]
    • 從月均值延長成24小時之逐時值(CAMx可以接受典型日變化之日排放檔)

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

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

氣溫反比變化

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

  1. 全年逐時T2擷取(T2.nc )
    • 以ncks由逐日檔案擷取T2逐時值
    • 逐日檔案合併(ncrcat)成全年逐時T2
  2. 全年月均值之氣溫(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)
  1. CMAQ排放檔案模板之準備
    • 修改d1.tot等檔案的TFLAG值，camx2cmaq不會認TSTEP=10000以外的數值（此處為30日=720）
    • 執行camx2cmaqd1.job進行12筆排放量之轉換
    • 將此12筆檔案複製(ncrcat)成至少744筆之模版
  2. 計算網格月均值之回歸公式（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][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][uamiv]時間標籤
    • 由[pncgen][pncgen]轉成nc檔之後，還需要再加上必須的全域屬性，因此放棄不持續發展，改由REAS2hr_TimVar.py為之。

add_ncatt.cs

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

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

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

  

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/
  • here：REASv3.1排放檔案之處理 *
  • relevent：REASv3.1點源之讀取、格式轉換與合併
  • parent：Dr. Kuang’s Evernotes_Grid Models
Search