CAMx模擬結果之壓縮

 

背景

  • 此程式系列之目的在於整併CAMx的模擬結果,產生測站量測之空品項目,以利後續分析及比對(aok,see also[[2022-07-20-aok]]1)。整併項目包括:
    • NMHC:濃度乘上碳數之sumproduct
    • PM:PM2.5及PM10
  • 除了fortran程式之外,同樣功能也可以在python中實現

shkavgcb6.f90程式設計

一般氣狀物及粒狀物

  • 使用二維疏鬆矩陣ridx做為對照表
    • 第1維是欲輸出之污染物項目。其順序為A10
    • 第2維是CAMx.in中所列之模擬輸出項目,順序列於最右邊駐解
    • 對應到VOC之ridx全為0,另行計算。
    real ridx(nout,nlst)
...
!                O3  NO2  SO2  VOC PM25 PM10 PNO3,PSO4
    data A10/'O3','NO2','SO2','VOC','PM25','PM10','PNO3','PSO4','CO'/
    data ridx /   0.,  0.,  0.,  0.,  0.,  0.,  0.,   0.,   0., &  !NO
                  0.,  1.,  0.,  0.,  0.,  0.,  0.,   0.,   0., &  !NO2
                  1.,  0.,  0.,  0.,  0.,  0.,  0.,   0.,   0., &  !O3
                  0.,  0.,  1.,  0.,  0.,  0.,  0.,   0.,   0., &  !SO2
                  0.,  0.,  0.,  0.,  0.,  0.,  0.,   0.,   0., &  !NH3
                  0.,  0.,  0.,  0.,0.58,  1.,  1.,   0.,   0., &  !PNO3Tsai & Cheng
...
                  0.,  0.,  0.,  0.,  1.,  1.,  0.,   0.,   0., &  !SOPB
                  0.,  0.,  0.,  0.,  1.,  1.,  0.,   0.,   0./    !SOA7
  • 氣狀物之單位轉換:ppm -> ppb
    • 對象:’O3’,’NO2’,’SO2’
    • CO仍為ppm
ridx(1:3,:)=ridx(1:3,:)*1000.
  • ridx之應用
 !non-voc
           sumc=0
           do lnv=1,nspec_nvoc
             l=SEQ(lnv)
             sumc=sumc+rval(ii,kk,l,it)*ridx(ll,iin(l))
           enddo
           oconc(ii,kk,ll,it)=sumc

VOC’s之計算

  • 碳數參考CAMx手冊表格(如)
  • VOC’s項目之化學物質名稱(lsvname),及其碳數rvidx
    • 與反應機制有關。當化學機制選項不同、程式更新時,此2序列需重新檢查並予以更新。
    data lsvname/'HCO3', 'HO2', 'MEO2', 'ROR', 'CO', 'CH4', 'FACD', 'FORM', 'KET', 'MEOH', &
                 'MEPX', 'PAR', 'ECH4', 'XPAR', 'C2O3', 'AACD', 'ALD2', 'ALDX', 'ETH', 'ETHA',&
                 'ETHY', 'ETOH', 'GLY', 'GLYD', 'PACD', 'PAN', 'CXO3', 'ACET', 'MGLY', 'OLE', &
                 'PANX', 'PRPA', 'XPRP', 'OPO3', 'IOLE','ISPD', 'NTR', 'OPAN', 'OPEN', 'EPX2',&
                 'ISO2', 'EPOX', 'HPLD', 'INTR', 'ISOP', 'ISPX', 'XOPN','BZO2', 'BENZ','CRO',&
                 'TO2', 'CAT1', 'CRES', 'CRON', 'TOL', 'CRNO', 'CRN2', 'CRPX', 'CAO2','XLO2',&
                  'XYL','TERP'/
    data rvidx / 1, 1, 1, 1, 0, 0, 1, 1, 1, 1,&
                 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,&
                 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,&
                 3, 3, 3, 4, 4, 4, 4, 4, 4, 5,&
                 5, 5, 5, 5, 5, 5, 5, 6, 6, 7,&
                 7, 7, 7, 7, 7, 7, 7, 7, 7, 8,&
                 8, 10/
  • rvidx之應用
!voc
          sumc=0
          do lnv=1,nspec_voc
            l=SEQV(lnv)
            sumc=sumc+rval(ii,kk,l,it)*rvidx(iin(l))
          enddo
          oconc(ii,kk,4,it)=sumc

shkavg4.4.f90

  • 這個版本的VOC’s項目(nvlst)較少,只有36項,適用較早期的cb5反應機制。
  • 輸出項目也少1項CO。這是因為早期反應機制中沒有輸出這個項目。
$ diff shkavg4.4.f90 shkavgcb6.f90
9,10c9,10
<    integer,parameter::nlst=32,nvlst=36  
<    integer,parameter::nout=8  
---
>    integer,parameter::nlst=32,nvlst=62 
>    integer,parameter::nout=9  

shkavgPar.f90

  • 這個最早版本與shkavg4.4.f90的差異在於部分VOC’s的名稱在反應機制中略有改變。
  • 影響所及也包括個別物質的碳數
diff shkavg4.4.f90 shkavgPar.f90
130,131c131,132
<    lsvname(17)  = 'ISPD'
<    lsvname(18)  = 'ISPX' !
---
>    lsvname(17)  = 'ISP'  !ISP = ISOP
>    lsvname(18)  = 'ISPD' !
145c146
<    lsvname(32)  = 'TO2' !TO2 =  Peroxy radical from OH addition to TOL 7 173.1
>    lsvname(32)  = 'TOLA' !TOLA = TOL or ARO1
149c150
<    lsvname(36)  = 'XLO2' !Peroxy radical from OH addition to XYL 8 187.1
---
>    lsvname(36)  = 'XYLA' !XYLA = XYL or ARO2

shk3.py

  • 這支程式會讀取測站上的測值來修正CAMx的瞬間檔,做為以觀測濃度重啟模式之用。

程式IO

  • 引數:年月日時(YYMMDDHH共8碼)
  • 粒狀物名稱:從'/nas1/camxruns/2013_enKF/inputs/part.txt'中讀取
    • 此檔案由反應機制文字檔中拮取
  • 測站觀測值(同一網格測站取平均值):'/home/backup/data/epa/pys/2013ByMonth/'+YMDH[:4]+'IJ4.csv'
  • CAMx模擬結果瞬時檔(YYMMDDHH.finst)
    • 會複製成輸出檔(YYMMDDHH R.finst)
    • 檔案格式:uamiv

分配邏輯

  • 氣狀物['CO', 'NO2', 'O3', 'SO2']
    • 測站網格所在位置的濃度,直接換成測站(平均)值
  • NMHC
    • 先計算模擬之NMHC(ppbc)
    • 按照個別物種所佔ppbc比例整比分配觀測值
    • 除回碳數成為ppm濃度
  • 粒狀物
    • 每個成分的PM2.5部分有其比例。乘上比例加總成PM2.5模擬值,修正值為觀測值按模擬值正比分配。
    • PM10為所有粒狀物的總合。計算粗粒部分之觀測與模擬值。將粗粒分配給'CPRM CCRS PH2O'.split()這3項

shk_camx.cs

  • 這個版本以pncgen進行變數之提取(沒有加總)
  • 如果是nc格式,則以ncks進行提取,選項略有差異
kuang@master ~/bin
$ cat shk_camx.cs
if [ $HOSTNAME == '114-32-164-198.HINET-IP.hinet.net' ];then NCO='/opt/anaconda3/bin'
elif [ $HOSTNAME == 'master' ];then NCO='/cluster/netcdf/bin'
elif [ $HOSTNAME == 'node01' ];then NCO='/cluster/netcdf/bin'
elif [ $HOSTNAME == 'centos8' ];then NCO='/opt/anaconda3/envs/py37/bin'
elif [ $HOSTNAME == 'node03' ];then NCO='/opt/miniconda3/bin'
else NCO='/usr/bin'
fi
NCKS='pncgen -f uamiv'
NCRCAT=${NCO}/ncrcat

VAR='TFLAG,CO,NO2,SO2,O3,PAR,PSO4,PNO3,PNH4,FPRM,FCRS'

$NCKS -O --slice LAY,0,0 -v $VAR $1 $2
#$NCKS -O -v $VAR -d LAY,0,0 $1 $2

shk.cs

  • 原則同上,程式自行判斷是CMAQ或者是CAMx產出之結果。

cs腳本

kuang@master ~/bin
$ cat shk.cs
nc=$1
if [ $HOSTNAME == '114-32-164-198.HINET-IP.hinet.net' ];then NCO='/opt/anaconda3/bin'
elif [ $HOSTNAME == 'master' ];then NCO='/cluster/netcdf/bin'
elif [ $HOSTNAME == 'node01' ];then NCO='/cluster/netcdf/bin'
elif [ $HOSTNAME == 'node02' ];then NCO='/cluster/netcdf/bin'
elif [ $HOSTNAME == 'centos8' ];then NCO='/opt/anaconda3/envs/py37/bin'
elif [ $HOSTNAME == 'node03' ];then NCO='/opt/miniconda3/bin'
else NCO='/usr/bin'
fi
NCKS=${NCO}/ncks
NCRCAT=${NCO}/ncrcat
NCDUMP=${NCO}/ncdump
v=$($NCDUMP -h $nc|grep PM25|wc -l)
a=$($NCDUMP -h $nc|grep AVERAGE|wc -l)
if [ $v != 0 ];then
  if  [ $a == 0 ];then
    #cmaq combined files
    VAR='TFLAG,CO,NO2,SO2,O3,PM25_NO3,PM25_SO4,PM25_TOT,PM10,VOC'
  else
  #version 700 nc directly from CAMx
    VAR='TFLAG,CO,NO2,SO2,O3,PNO3,PSO4,PNH4,PAR,CCRS,FCRS,CPRM,FPRM'
  fi
#camx>400, nc generated by pncgen from uamiv file
else
  VAR='TFLAG,CO,NO2,SO2,O3,PNO3,PSO4,PNH4,PAR,CCRS,FCRS,CPRM,FPRM'
fi
if ! [ -e $2 ];then
  touch $2
  res=${VAR//[^,]}
  nvars=$(echo ${#res})
  $NCKS -O -v $VAR -d LAY,0,0 -d VAR,1,$nvars $1 $2
  ncatted -a NVARS,global,o,i,$nvars $2
fi
  • $res$VAR中的所有逗號(從開使到逗號中的內容置換成null,[...]則會重覆執行)
  • ${#res}則為$res的長度

do_shk

  • 這項應用是將全年的combine結果一次提取
kuang@master /nas1/cmaqruns/2019base/Annual
$ cat ./do_shk
#note the ncks or ncrcat can not run in parallel, will conflick in memory
ls -r /nas1/cmaqruns/2016base/data/out*/POST/COMBINE_ACONC*_sCh*10.nc > fnames.txt
sort fnames.txt>a;mv a fnames.txt
for i in $(cat fnames.txt);do ymd=$(echo $i|cut -d_ -f11);shk.cs $i $ymd.nc;done &

shk_Days_DM.cs

  • 這支程式是專為CMAQ模擬結果所用的。由於combine結果是逐日分檔,而且含有許多不需要的污染項目,因此如果需要一段期間、特定項目的濃度檔,就必須一一執行ncks提取。再以ncrcat予以串連。
  • 程式有3個引數,分別為起、迄的月日,以及domain id
kuang@master ~/bin
$ cat shk_Days_DM.cs
#$1:first mmdd
#$2:last  mmdd
#$3:domain d01~d04

if [ $HOSTNAME == '114-32-164-198.HINET-IP.hinet.net' ];then NCO='/opt/anaconda3/bin'
elif [ $HOSTNAME == 'master' ];then NCO='/cluster/netcdf/bin'
elif [ $HOSTNAME == 'centos8' ];then NCO='/opt/anaconda3/envs/py37/bin'
elif [ $HOSTNAME == 'node03' ];then NCO='/opt/miniconda3/bin'
else NCO='/usr/bin'
fi
NCKS=${NCO}/ncks
NCRCAT=${NCO}/ncrcat

VAR='TFLAG,CO,NO2,SO2,O3,PM25_NO3,PM25_SO4,PM25_TOT,PM10,VOC'

CASE=10
ROOT=$3
test $ROOT == 'd04' && GRID='TWN_3X3'
test $ROOT == 'd01' && GRID='EAsia_81K'
test $ROOT == 'd02' && GRID='sChina_27k'

yymm=$(echo $PWD|cut -d'_' -f5|cut -d'/' -f1)
bdate=`echo $(ls -rt COMBINE_ACONC*${GRID}_${CASE}.nc|head -n1)|cut -d'_' -f7`
bjul=$(date -d "$bdate" +%Y%j)
fjul=$(date -d "2016$1" +%Y%j)
ljul=$(date -d "2016$2" +%Y%j)
bi=$(( $fjul - $bjul + 1))
ei=$(( $ljul - $bjul + 1))
for ((i=$bi;i<=$ei;i+=1));do
  j=$(( 10#$i - 1 ))
  r=$( echo $j/4+5|bc )
  cdate=$(date -d "$bdate +${j} day" +%Y%m%d)
  file=COMBINE_ACONC_v53_gcc_${yymm}_run${r}_${cdate}_${GRID}_${CASE}.nc
  echo $file
  $NCKS -O -v $VAR -d LAY,0,0 $file ${cdate}_$ROOT.nc
done
$NCRCAT -O 20*_$ROOT.nc $ROOT.nc

CB6化學物質名稱碳數及分子量

Table 5-2. CAMx species names and descriptions common to all Carbon Bond Mechanisms.

Model Species Description Carbon #1 Mol. Wt.2
BZO2 Peroxy radical from OH addition to benzene 6 159.1
C2O3 Acetylperoxy radical 2 75.0
CRO Alkoxy radical from cresol 7 107.1
CXO3 C3 and higher acylperoxy radicals 3 89.0
EPX2 Peroxy radical from EPOX reaction with OH 5 149.1
HCO3 Adduct from HO2 plus formaldehyde 1 63.0
HO2 Hydroperoxy radical 1 28.0
ISO2 Peroxy radical from OH addition to isoprene 5 117.1
MEO2 Methylperoxy radical 1 47.0
NO3 Nitrate radical 0 62.0
O Oxygen atom in the O3 (P) electronic state 0 16.0
O1D Oxygen atom in the O1 (D) electronic state 0 16.0
OH Hydroxyl radical 0 17.0
OPO3 Peroxyacyl radical from OPEN 4 115.0
RO2 Operator to approximate total peroxy radical concentration 0 87.1
ROR Secondary alkoxy radical 1 71.1
TO2 Peroxy radical from OH addition to TOL 7 173.1
XLO2 Peroxy radical from OH addition to XYL 8 187.1
XO2 NO to NO2 conversion from alkylperoxy (RO2) radical 0 87.1
XO2H NO to NO2 conversion (XO2) accompanied by HO2 production 0 87.1
XO2N NO to organic nitrate conversion from alkylperoxy (RO2) adical 0 87.1
AACD Acetic acid 2 60.0
ACET Acetone 3 58.1
ALD2 Acetaldehyde 2 44.0
ALDX Propionaldehyde and higher aldehydes 2 58.1
BENZ Benzene 6 78.1
CAT1 Methyl-catechols 7 124.1
CO Carbon monoxide 1 28.0
CH4 Methane 1 16.0
CRES Cresols 7 108.1
CRON Nitro-cresols 7 153.1
EPOX Epoxide formed from ISPX reaction with OH 5 118.1
ETH Ethene 2 28.0
ETHA Ethane 2 30.1
ETHY Ethyne 2 26.0
ETOH Ethanol 2 46.1
FACD Formic acid 1 46.0
FORM Formaldehyde 1 30.0
GLY Glyoxal 2 58.0
GLYD Glycolaldehyde 2 60.0
H2O2 Hydrogen peroxide 0 34.0
HNO3 Nitric acid 0 63.0
HONO Nitrous acid 0 47.0
HPLD hydroperoxyaldehyde 5 116.1
INTR Organic nitrates from ISO2 reaction with NO 5 147.1
IOLE Internal olefin carbon bond (R-C=C-R) 4 56.1
ISOP Isoprene 5 68.1
ISPD Isoprene product (lumped methacrolein, methyl vinyl ketone, etc.) 4 70.1
ISPX Hydroperoxides from ISO2 reaction with HO2 5 118.1
KET Ketone carbon bond (C=O) 1 72.1
MEOH Methanol 1 32.0
MEPX Methylhydroperoxide 1 48.0
MGLY Methylglyoxal 3 72.0
N2O5 Dinitrogen pentoxide 0 108.0
NO Nitric oxide 0 30.0
NO2 Nitrogen dioxide 0 46.0
NTR Organic nitrates 4 119.1
O3 Ozone 0 48.0
OLE Terminal olefin carbon bond (R-C=C) 3 42.1
OPAN Peroxyacyl nitrate (PAN compound) from OPO3 4 161.0
OPEN Aromatic ring opening product (unsaturated dicarbonyl) 4 84.0
PACD Peroxyacetic and higher peroxycarboxylic acids 2 76.0
PAN Peroxyacetyl Nitrate 2 121.0
PANX C3 and higher peroxyacyl nitrate 3 135.0
PAR Paraffin carbon bond (C-C) 1 72.1
PNA Peroxynitric acid 0 79.0
PRPA Propane 3 44.1
ROOH Higher organic peroxide 0 90.1
SO2 Sulfur dioxide 0 64.0
SULF Sulfuric acid (gaseous) 0 98.0
TERP Monoterpenes 10 136.2
TOL Toluene and other monoalkyl aromatics 7 92.1
XOPN Aromatic ring opening product (unsaturated dicarbonyl) 5 98.1
XYL Xylene and other polyalkyl aromatics 8 106.2
NTR1 Simple organic nitrates 0 119.1
NTR2 Multi-functional organic nitrates 0 135.1
ECH4 Emitted methane (to enable tracking separate from CH4) 1 16.0
XPRP Operator for organic nitrates from PRPA 3 89.1
XPAR Operator for organic nitrates from PAR 1 117.1
CRNO Nitro-cresol oxy radical 7 152.1
CRN2 Nitro-cresol peroxy radical 7 168.1
CRPX Nitro-cresol hydroperoxide 7 169.1
CAO2 Ring-opening product from methyl catechol 7 173.1
  1. Carbon # is the precise number of carbon atoms for each model species.
  2. Mol. Wt. is a representative molecular weight, intended only for estimating molecular diffusivity, e.g. in dry deposition calculations. Diffusivity requires a different interpretation (complete molecules) than Carbon Bond chemistry (chemical groups).
  1. https://sinotec2.github.io/FAQ/2022/07/20/aok.html “ CAMx模擬結果之比對(aok)”