分享一下去年底寫的一隻程式,雖然是去年寫好的,但是因為工作所需今年才上線,就在這邊分享一下了

源起:會做這工作主要是因為工作所需,這邊的博士後把去年的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檔的時候,為了怕矩陣位置被覆蓋或者記憶體重複使用才這樣做的吧,所以如果只是要重新做一些初始資料的應用的話,建議直接用取代的方式比較簡單。

arrow
arrow
    全站熱搜

    廖董不懂 發表在 痞客邦 留言(1) 人氣()