Link Search Menu Expand Document
Focus on Air Quality

MOZARD/WACCM模式輸出轉成CMAQ初始條件_水平內插與污染項目對照

Table of contents

背景

  • 全球模式模擬結果要使用在地區的空品模擬，需要經過空間、時間的內插、以及空品項目的對照等作業。此處進行水平內插與污染物的對照。
  • 前置作業包括
  • 由於CMAQ濃度檔案可以使用combine予以整併，(或直接)就可以用VERDI或進行檢視。
  • 後續作業
    • 邊界濃度的製作，包括時間的內插、邊界框濃度的解析等，將在bcon程式內進行。
    • 直接以初始濃度引用：在run_cctm.csh內指定即可
  • 有鑒於全球空品模擬結果越來越多，且Mozart模式不再維護，有更多直接轉換之程式與作業方式(如ecmwf之EAC4)，可以簡化程序。

程式說明

執行方式

  • (eg.)執行2018年4月1日0時之內插：python moz2cmaqHd1J.py 1809100

引數

  • 年代+Julian day+小時共7碼

I/O檔案

同步執行

  • 範例同時執行程式由2018/4/1~4/30日(Juian date=091~120)之濃度內插
    • 使用Julian day$j和小時$h雙層迴圈
    • sub=$1 $2 $3 $4 $5 $6 $7 $8 $9 ${10} ${11} ${12} ${13} ${14} ${15} ${16} ${17} ${18} ${19} ${20} &
for j in {091..120};do 
  for h in 00 06 12 18;do 
    sub python moz2cmaqHd1J.py 18$j$h
  done
done

執行成果的整合(ncrcat)

  • 逐6小時的結果，需要整合成批次時間範圍內的濃度檔序列，此處使用ncrcat程式
  • 濃度檔序列可以連結到一新的目錄，然後用ncrcat將其整合，或使用迴圈累加檔名字串，在將字串內的檔案整合。範例如下：
    • 將全月檔案連結到目錄、再以ncrcat整合
for i in {12..12};do
  cdate=2019${i}
  pdate=`date -d "${cdate}01 -1  days" +%Y%j`
  ndate=`date -d "${cdate}01 +1 month" +%Y%j`
  mkdir -p $cdate
  cd $cdate
  for ((d=$pdate;d<=$ndate;d+=1));do
    for h in 00 06 12 18;do
      ff=../ICON_${d}${h}.d1
      if ! [ -e $ff ];then continue;fi
      ln -sf $ff .
    done
  done
  cd ..
  ncrcat -O ${cdate}/* ICON_d1_${cdate}.nc &
done
  • 按照執行批次的日期範圍，進行字串累加， 再以ncrcat整合
for mo in {01..12};do
begd=$(date -ud "2019${mo}15 - 1 month" +%Y%m%d)
for r in {5..12};do
    d0=$(( ($r-1)*4 ))
    d=d0
    II=ICON_
  for ((d=d0;d<=$(( $d0+5));d+=1));do
    j=$(date -ud "$begd + $d days" +%Y%j)
    if [ $d == $d0 ]; then
      files=$II${j}*.d1
    else
      files=${files}' '$II${j}*.d1
    fi
  done
  sub ncrcat -O ${files} ICON_d1_${y}${mo}_run${r}.nc
done
done

程式分段說明

  • 調用模組
import numpy as np
import netCDF4
import os,sys, datetime, json
from pyproj import Proj
  • 使用pyproj的Proj來進行LCC轉經緯度。此處採經緯度網格的內插(較為方便)。
Latitude_Pole, Longitude_Pole = 23.61000, 120.990
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)
  • 讀取引數與時間
if (len(sys.argv) != 2):
  print ('usage: '+sys.argv[0]+' YYJJJHH(eg. 1636500)')
yrjulhh=sys.argv[1]
yrmn=datetime.datetime.strptime(yrjulhh[:-2], '%y%j').strftime('%y%m')
  • 讀取mozart等全球模式輸出結果的經緯度及時間標籤字串
    • 時間標籤以字串形式，方便與引數進行比對，找到指定時間(tM)
#read the mozart model results
fname='moz_41_20'+yrmn+'.nc'
if not os.path.exists(fname):sys.exit(fname+' not exist')
ncM = netCDF4.Dataset(fname,'r')
v4M=list(filter(lambda x:ncM.variables[x].ndim==4, [i for i in ncM.variables]))
ntM,nlayM,nrowM,ncolM=(ncM.variables[v4M[0]].shape[i] for i in range(4))
lonM=[ncM.XORIG+ncM.XCELL*i for i in range(ncolM)]
latM=[ncM.YORIG+ncM.YCELL*i for i in range(nrowM)]
tflagM=[str(i*100+j//10000)[2:] for i,j in zip(ncM.variables['TFLAG'][:,0,0],ncM.variables['TFLAG'][:,0,1])]
tM=tflagM.index(yrjulhh)

#load mz2cmaq map
with open('cb6_new.json','r') as f:
  mz2cm=json.load(f)
with open('cb6_newNum.json','r') as f:
  mz2cmNum=json.load(f)
with open('lay2VGLEVLS.json','r') as f:
  lay2VGLEVLS=json.load(f)
lay2VGLEVLS.update({'40':40})
  • 讀取CMAQ濃度檔案模版
    • 讀取座標與時間標籤
#read the template(s)
fname='ICON_20'+yrjulhh+'.d1'
nc = netCDF4.Dataset(fname,'r')
v4=list(filter(lambda x:nc.variables[x].ndim==4, [i for i in nc.variables]))
nt,nlay,nrow,ncol=(nc.variables[v4[0]].shape[i] for i in range(4))
X=[nc.XORIG+nc.XCELL*i for i in range(ncol)]
Y=[nc.YORIG+nc.YCELL*i for i in range(nrow)]
lon, lat= pnyc(X,Y, inverse=True)
lon_ss=np.searchsorted(lonM,lon)
lat_ss=np.searchsorted(latM,lat)
tflag=nc.variables['TFLAG'][:,0,:]
nc.close()
  • 去除對照不到的項目
#drop keys which values not in new CMAQ spec_list
mc=list(mz2cm)
for i in mc:
  if i not in v4M:
    mz2cm.pop(i)
    continue
  if mz2cm[i] not in v4:
    mz2cm.pop(i)
  • 儲存全球模式濃度矩陣
    • 進行高度對照
    • 將濃度單位由Volume Mixing Ratio轉到CMAQ的PPM
#save the matrix
v4M=list(mz2cm)
A5=np.zeros(shape=(len(v4M),nlayM,nrowM,ncolM))
for ix in range(len(v4M)):
  for k in range(nlayM):
    A5[ix,k,:,:]=ncM.variables[v4M[ix]][tM,lay2VGLEVLS[str(k)],:,:]*1000.*1000. #Volume Mixing Ratio to PPM
ncM.close()
  • 清空模版的濃度值
    • 避免netCDF4自動遮蔽
    • 因為模版只有單一時間點，沒有時間的維度
#perform the horizontal interpolation and write results
fname='ICON_20'+yrjulhh+'.d1'
nc = netCDF4.Dataset(fname,'r+')
for x in v4M:
  nc.variables[mz2cm[x]][0,:,:,:]=0
  • 空間內插與污染項目對照
    • 使用簡單的線性內插
    • 由於不同物質可能都會貢獻到同一碳鍵項目，因此需要累加
for jcrs in range(nrow):
  jmz=lat_ss[jcrs]
  for icrs in range(ncol):
    imz=lon_ss[icrs]
    rx=(lon[icrs]-lonM[imz-1])/(lonM[imz]-lonM[imz-1])
    ry=(lat[jcrs]-latM[jmz-1])/(latM[jmz]-latM[jmz-1])
    A2x=A5[:,:,jmz,imz]*rx+A5[:,:,jmz,imz-1]*(1-rx)
    A2y=A5[:,:,jmz,imz]*ry+A5[:,:,jmz-1,imz]*(1-ry)
    A2=(A2x+A2y)/2.
    for x in v4M:
      nc.variables[mz2cm[x]][0,:,jcrs,icrs]+=A2[v4M.index(x),:-1]*mz2cmNum[x]
nc.close()

程式下載

Reference

  • lizadams, Visualization Environment for Rich Data Interpretation (VERDI): User’s Manual, github, August 03, 2021
  • sinotec2, NC矩陣遮罩之檢查與修改, FAQ,Dec 10 2021