背景
- 此程式系列之目的在於整併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
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 |
- Carbon # is the precise number of carbon atoms for each model species.
- 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).
-
https://sinotec2.github.io/FAQ/2022/07/20/aok.html “ CAMx模擬結果之比對(aok)” ↩