分享一下去年底寫的一隻程式，雖然是去年寫好的，但是因為工作所需今年才上線，就在這邊分享一下了
源起：會做這工作主要是因為工作所需，這邊的博士後把去年的TEDS分析好了，接下來就是要更新我們那用了十年的TEDS排放資料了，但是基本上程式資料都是人家給的，根本沒有自己做的，而且裡面的物種也是不完整，這位博士後很認真的幫我們分析資料，並加入了亞洲的排放，點線面都完成了，把汙染物從26種提升到50種，所以排放資料有了，接下來就是預報了。所以我的工作就是把他處理完的資料(TEDS)跟自然排放源(MEGAEN)加在一起，寫成一組新的netcdf檔案，提供給模式(CMAQ、WRFCHEM)使用，廢話不多說先來開始說程式吧。
其實一開始我並不想用NCL來寫的，我一開始用fortran去寫nc檔案，但是寫了一兩天就發生問題了，寫入nc矩陣很簡單，也一下就寫好了，但是寫矩陣說明跟global_attribute的時候格式都不能被模式接受，有很多小地方的錯誤，後來我就放棄了，使用NCL是因為我之前做氣象模擬的時候，有利用NCL去改變初始資料，測試溫度的敏感度，所以想說有了基礎，那就用這個好了，所以就去查了一下官網的說明https://www.ncl.ucar.edu/Applications/netcdf4.shtml，就從基礎開始了。
以下是我寫的source code
;=========================================================
;;;; read TEDS and MEGAN file to combine a NewFile ;;;;;
;========================================================
TEDS="2015335_D1.nc" !!整理好的 TEDS file
MEGAN="MEGAN_emission_SAPRC99_2015364_D1.nc" !!用WRF當天氣條件每日跑出來的自然排放
TEDS_MEGAN = addfile("Emission_15364_d1.nc","c") !!要產生的TEDS+MEGAN file
NTSTEP=73 三天預報
NDATE=2 NDATE指d d+1 d+2三天
NLAY=14 垂直14層
NVAR=50 50組變數
NROW=67 y軸網格
NCOL=85 x軸網格
dim_names = (/ "TSTEP","DATE-TIME", "LAY", "VAR", "ROW","COL" /) dimension title
dim_sizes = (/ NTSTEP, NDATE ,NLAY, NVAR, NROW,NCOL /) dimension number
dimUnlim = (/ True, False , False, False, False, False /) 固定dimension unmber,False固定，true不固定
filedimdef(TEDS_MEGAN, dim_names, dim_sizes, dimUnlim ) file dimension def(dimensions:
TSTEP = UNLIMITED ; // (25 currently)
DATE-TIME = 2 ;
LAY = 14 ;
VAR = 50 ;
ROW = 67 ;
COL = 85 ;)
ST_DATE=2015364
print(ST_DATE)
T_d04 = addfile(TEDS,"rw") 讀取TEDS file
M_d04 = addfile(MEGAN,"rw") 讀取Megan file
;==========================================================
; add species a new Time flag in the New file (TFLAG)
;==========================================================
T_TFLAG = T_d04->TFLAG 把TEDS的排放時間輸出
M_TFLAG = M_d04->TFLAG 把Megan的排放時間輸出
TOTAL_TFLAG = new((/73,50,2/),integer) 創立新的Teds+Megan矩陣
do TFLAG=0 ,23 ======================
do var=0 , 49 把
TOTAL_TFLAG(TFLAG,var,0)=M_TFLAG(TFLAG,0,0) 一
TOTAL_TFLAG(TFLAG,var,1)=T_TFLAG(TFLAG,var,1) 天
end do 的
end do 矩
do TFLAG= 24, 47 陣
do var=0 , 49 變
TOTAL_TFLAG(TFLAG,var,0)=M_TFLAG(TFLAG,0,0) 成
TOTAL_TFLAG(TFLAG,var,1)=T_TFLAG(TFLAG-24,var,1) 三
end do 天
end do
do TFLAG=48 ,72
do var=0 , 49
TOTAL_TFLAG(TFLAG,var,0)=M_TFLAG(TFLAG,0,0)
TOTAL_TFLAG(TFLAG,var,1)=T_TFLAG(TFLAG-48,var,1)
end do
end do =====================
varatt = getvaratts(T_TFLAG) 寫入變數說明，我是吃入TEDS TFLAG的說明
TOTAL_TFLAG@$varatt(2)$=T_TFLAG@$varatt(2)$
TOTAL_TFLAG@$varatt(1)$=T_TFLAG@$varatt(1)$
TOTAL_TFLAG@$varatt(0)$=T_TFLAG@$varatt(0)$
TEDS_MEGAN->TFLAG = TOTAL_TFLAG =====================
delete(T_TFLAG) 關閉矩陣減少記憶體耗損
delete(M_TFLAG)
delete(TOTAL_TFLAG)
print(systemfunc("/bin/date "))
print("TOTAL_TFLAG")
;==========================================================
; add species a new content in the New file (NO2)
;==========================================================
T_NO2 = T_d04->NO2 (同上變數，只是矩陣不同，因為純時間變數跟物理變數的矩陣不同)
M_NO2 = M_d04->NO2
TOTAL_NO2 = new((/73,14,67,85/),float)
TOTAL_NO2 = 0.0
do i = 0,23
TOTAL_NO2(i,:,:,:)=T_NO2(i,:,:,:)
end do
do i = 24,47
TOTAL_NO2(i,:,:,:)=T_NO2((i-24),:,:,:)
end do
do i = 48,72
TOTAL_NO2(i,:,:,:)=T_NO2((i-48),:,:,:)
end do
TOTAL_NO2(:,0,:,:) = M_NO2(:,0,:,:)+TOTAL_NO2(:,0,:,:)
TOTAL_NO2!0 = "TSTEP"
TOTAL_NO2!1 = "LAY"
TOTAL_NO2!2 = "ROW"
TOTAL_NO2!3 = "COL"
varatt = getvaratts(T_NO2)
TOTAL_NO2@$varatt(2)$=T_NO2@$varatt(2)$
TOTAL_NO2@$varatt(1)$=T_NO2@$varatt(1)$
TOTAL_NO2@$varatt(0)$=T_NO2@$varatt(0)$
TEDS_MEGAN->NO2 = TOTAL_NO2
delete(T_NO2)
delete(M_NO2)
delete(TOTAL_NO2)
print(systemfunc("/bin/date "))
print("TOTAL_NO2")
;==============================================
看你要增加幾個變數就自己增加吧，最後也是最重要的就是global_attribute
我是把 TEDS整段global_attribute複製過來貼上去，再去改一些細部的說明
========================
---------- global_attribute--------------
=========================
global_attnames = getvaratts(T_d04) ;-- retrieve the global attributes from input file
printVarSummary(global_attnames)
print(global_attnames)
do i=0,dimsizes(global_attnames)-1
print(T_d04@$global_attnames(i)$) ;-- print global attribute contents
TEDS_MEGAN@$global_attnames(i)$ = T_d04@$global_attnames(i)$ ;-- write global attributes to new file
end do
fAtt = True
fAtt@SDATE =ST_DATE
fAtt@FILEDESC = "this is a emission file for CMAQ by carls "
fAtt@WDATE =ST_DATE
fAtt@CDATE =ST_DATE
fAtt@$"VAR-LIST"$ = "NO2 NO HONO CO SO2 SULF HCHO MEOH MEK PROD2 CCO_OH CCHO RCO_OH ACET PHEN HCOOH RCHO GLY MGLY BACL CRES BALD MACR MVK IPROD ETHENE ISOPRENE TRP1 ALK1 ALK2 ALK3 ALK4 ALK5 ARO1 ARO2 BENZENE OLE1 OLE2 SESQ NH3 POA PSO4 PNO3 PEC PMFINE PMC HAC GLYC C2H2 PRPE "
fileattdef( TEDS_MEGAN, fAtt )
;================================================
然後運氣很好的，就可以每天自動產出emssion給CMAQ WRFCHEM吃了。
感想：NCL其實還滿好用的，只是輸出NC檔的時候非常的笨，因為他不會直接加進去後面的矩陣，當你要在寫入B的時候，他要先打開A，再把B放進去，所以寫入C的時候，就要先把A&B打開，再把C放進去，如此一來變數越多其實越花時間，不過根據側面了解ncgen這隻執行檔在做nc檔的時候也是一樣的方式，所以可能不是ncl問題，而是寫入nc檔的時候，為了怕矩陣位置被覆蓋或者記憶體重複使用才這樣做的吧，所以如果只是要重新做一些初始資料的應用的話，建議直接用取代的方式比較簡單。
留言列表