環境工程篇

AERMOD - AERMET 氣象處理模組

AERMET 是 AERMOD 裡面很重要的一環!因為氣象參數的影響對於空品模式是相當重要的!

最近更新時間: 2019/07/28

影響擴散模式最重要的一環,當然就是氣象場資料。AERMOD 的氣象場資料,主要是透過相關附屬模組,也就是 AERMET 的協助來處理。 如前所述,AERMOD 運作時,需要一組地面風場資料,以及一組探空資料。然後將這兩個資料丟進 AERMET 整理成為 AERMOD 可以讀取的數據, 才能夠開始進行模式的運作。就讓我們來理解一下如何操作 AERMET 吧!

AERMET 的處理流程概述 (ps1)

AERMET 對於 AERMOD 來說,是很重要的!因為氣象條件不準確的話,後續的所有模擬與分析都是垃圾!

  • 前言

AERMET 模式是 AERMOD 的附屬模組之一,主要的目的在提供 AERMOD 運作時,所需要的 (1)地面風場資料與 (2)探空資料等。 但是 AERMET 一開始是針對美國的氣象場設計而來的,因此,讀取的資料格式主要是美國國家氣象服務 (National Weather Service, NWS) 所觀測的氣象資料格式。 NWS 主要提供兩種資料,一個是每天兩次的探空資料 (soundings),一個是地面的觀測站資料。在取得上述的資料後,AERMET 可以進行如下的動作:

  • 首先,將 NWS 提供的氣象資料讀入,並進行氣象資料的品質確認 (QA/QC) 任務,確認數據的可靠度
  • 再將上述的資料以小時逐時記載的方式,儲存到單一檔案中
  • 由上述的檔案中,讀入所需要的氣象場資料,開始整合所需要的地面風場資料,以及計算出模擬區域的 PBL 資料,當然, 這些探空資料是多層次的。

最終 AERMET 在完成上述的動作之後,會產生兩個檔案:

  • 地面小時風場資料與 PBL 的特徵參數;
  • 多層次的風速、風向資料
  • AERMET 概觀

US EPA 指定了 AERMOD 作為法規標準的擴散模式,而 AERMOD 擴散模式主要在模擬近場的污染濃度擴散模擬,所謂的近場大致上範圍就是 50 公里內。 為了要模擬擴散行為,AERMET 就被開發用來作為擴散模擬的氣象場準備資料。AERMET 可以評估計算大氣邊界層 (ABL) 或者是所謂的行星邊界層 (PBL) 的特徵參數, 據以提供給 AERMOD 進行模擬。

AERMET 主要運作三種資料格式:

  1. 由 NWS 或 FAA (Federal Aviation Administration) 所提供的地面測站的逐時觀測資料;
  2. 由 NWS 每天兩次施放監測到的探空氣球的探空資料 (主要計算 PBL)
  3. 由類似 WRF 的氣象資料,並經過 MMIF 模組轉換之後,取得的大氣模式輸出的資料 (並不是觀測數據喔!)

基本上的運作主要分為兩種,一種是完全由觀測值取得的地面與高空資料,就是上面的 1, 2 項資料。如果你手邊有大氣模式的輸出檔,例如 WRF 的全年輸出資料, 那就可以透過 AERMET 的外掛,亦即是 MMIF 這個模組,直接從已經網格化的氣象場當中取出所需要的數據即可!因此,上述的第 3 項資料,其實就已經包含了全部所需要的地面與探空資料了。 而且數據可能比觀測值還要更可靠!

因為 WRF 所輸出的氣象場資料非常龐大,除了原本就有在玩氣象模式的朋友之外,大部分的朋友大概很難取這些資料。所以,US EPA 所提供的使用手冊, 也只是針對上述的 1, 2 項資料來進行說明。整體的 AERMET 流程如下圖所示:

AERMET 的整體處理流程
AERMET 的整體處理流程

因為 NWS 提供的數據資料相當龐大,因此地面觀測資料與探空資料,就得要針對不同地點來擷取 (extracted) 出來!你可以看到上圖中, 從左而右的最上方兩條線路,第一條就是取出地面風場資料 (surface observations),第二條就是探空資料 (soundings)。第一條跟第二條中間的一個小區塊, 名稱為『 1-min ASOS 』是比較新的數據,是一分鐘一筆的數據資料,並不是小時資料!因為台灣似乎沒有這方面的數據,所以,這一條 (1-min ASOS) 就可以先忽略。 至於最後一條『 Raw site-specific data 』則是針對測站的相關設定,則可以直接指定到設定檔當中,也不用事先讀出。

對於台灣這種非 NWS 資料格式的系統來說,比較討厭的,就是得要事先將我們自己的觀測資料格式轉成 AERMET 可以處理的格式!也就是說, 要先將格式轉為 NWS 的格式之後,後續就能直接使用 AERMET 的相關模式功能。否則,連開始都無法開始啊!

如上圖所示,基本上, AERMET 共需要進行三個階段 (stage),第一階段就是在進行地面測站及探空資料的擷取,一直到 QA 之後,都是屬於第一階段的工作範圍。 第二階段則是在整理 24 小時的資料,將資料彙整成為單一輸出檔,如上圖中右側數來的第二個欄位:『 Merged Data: 24-h blocks 』就是這個階段的資料。 第三階段就是將剛剛的單一氣象檔案轉成地面與 PBL 參數,還有垂直剖面檔案兩者。這三個階段必需要分別進行,因為每個階段都會產生下個階段的輸入資料! 所以,當然需要依序進行三次才行!

step 1: 中央氣象局的資料轉換

AERMET 可以讀入的氣象場參數資料主要格式為 TD-3280 與 CD-144 等,不過,咱們中央氣象局的氣象場資料,並不是這兩種格式喔! 你可以自行到『大氣水文研究資料庫』(ps2)去查閱相關的氣象局下載資料。因此,如果需要讓 AERMET 可以讀入氣象場資料, 就得要將中央氣象局的資料進行轉換。不過,要轉換也有點麻煩~因為 AERMET 有規定自己的讀取格式!相關格式可以參考 AERMET 使用者手冊的附錄 C (ps1),裡面就有談到相關的讀取格式。

不過,光是從大氣水文研究資料庫下載資料,可能就花費許多時間了,更何況還得要參考上述的文件去轉換成為 AERMET 可以認識的格式! 鳥哥考慮到大家的困難~就先根據上述的文件,編譯好了一組可以在 Linux 系統上面運作的 cwb2aermet-static.exe 檔案。這個檔案只能在 Linux 底下運作! 不過運作之前還是得要注意輸入資料的格式問題!基本上,就是直接從中央氣象局的資料去分析就是了。

  • 步驟一:先從大氣水文資料庫下載中央氣象局資料與探空資料

一般研究人員可以從大氣水文資料庫所提供的界面去申請下載,當然,您得要有所謂研究人員的資格,否則,可能就得要發文去給相關的單位申請。 在 2000 年以前,我們想要做這種研究,是得樣向中央氣象局申請購買,然後對方會寄光碟片來,一打開...裡面都是圖片檔~然後一張一張列印出來, 再重複輸入到類似 Excel 的檔案裡面後,才能做後續的分析。現在已經都有電子檔,格式也都調整成固定的模式,已經改善非常非常多了。

上述的網站裡面順利的登入之後,點選『資料庫』的『基本查詢』之後,點選中央氣象局的『局屬逐時』、『探空』這兩個格式就好, 往下拉動後,到達時間設定的部份。這次我們選擇 2016 年的數據來解釋的話,那就勾選 2016 年的 12 個月的資料。 之後就按下『送出查詢』,你就會看到類似如下的畫面:

大氣水文資料庫的中央氣象局資料查詢結果
大氣水文資料庫的中央氣象局資料查詢結果

如上圖,按下『全選』然後再按下右上角的『產生壓縮檔』之後,過不多久,就可以下載壓縮檔了!這是因為資料量非常龐大,沒有壓縮的話,很耗頻寬啊! 等到你下載之後,在將該檔案解開,裡面就有 2016 年全年的地面測站資料與探空資料。至於檔名的制定方式,通常是這樣的:

  • 地面資料: 201601_cwb_hr.txt (2016年1月份的逐時資料)
  • 探空資料: 201601_upair.txt (2016年1月份的探空資料)

最後,將下載的檔案透過 zip 解壓縮之後,就可以得到正確的檔案了。在 Linux 系統底下,將 24 個檔案通通放置到同一個目錄去,假設為 cwb 目錄好了, 那麼就會有點像這樣:

# 1. 假設使用者名稱為 user,預計操作 AERMET 目錄在家目錄的 aermet 底下:
[user@localhost ~]$ mkdir aermet
[user@localhost ~]$ cd aermet
[user@localhost aermet]$ mkdir cwb
[user@localhost aermet]$ cd cwb

# 2. 假設你透過各種方式,已經將名為 2018091217311031398.zip 放置到本目錄中,並將他解壓縮
[user@localhost cwb]$ ll
-rw-rw-r-- 1 user user 55361553  7月 24 16:46 2018091217311031398.zip  <==必須存在這個檔案

[user@localhost cwb]$ unzip 2018091217311031398.zip  <==進行解壓縮!
Archive:  2018091217311031398.zip
  inflating: 201601_upair.txt
  inflating: 201602_upair.txt
  ...

[user@localhost cwb]$ ll
-rw-rw-r-- 1 user user  7867084  7月 24 16:48 201601_cwb_hr.txt
-rw-rw-r-- 1 user user   685822  7月 24 16:48 201601_upair.txt
-rw-rw-r-- 1 user user  7359676  7月 24 16:48 201602_cwb_hr.txt
-rw-rw-r-- 1 user user   700462  7月 24 16:48 201602_upair.txt
-rw-rw-r-- 1 user user  7867084  7月 24 16:48 201603_cwb_hr.txt
-rw-rw-r-- 1 user user   731511  7月 24 16:48 201603_upair.txt
-rw-rw-r-- 1 user user  7613380  7月 24 16:48 201604_cwb_hr.txt
-rw-rw-r-- 1 user user   811665  7月 24 16:48 201604_upair.txt
-rw-rw-r-- 1 user user  7867084  7月 24 16:48 201605_cwb_hr.txt
-rw-rw-r-- 1 user user   915487  7月 24 16:48 201605_upair.txt
-rw-rw-r-- 1 user user  7466068  7月 24 16:48 201606_cwb_hr.txt
-rw-rw-r-- 1 user user   816301  7月 24 16:48 201606_upair.txt
-rw-rw-r-- 1 user user  7613380  7月 24 16:48 201607_cwb_hr.txt
-rw-rw-r-- 1 user user   854182  7月 24 16:48 201607_upair.txt
-rw-rw-r-- 1 user user  7613380  7月 24 16:48 201608_cwb_hr.txt
-rw-rw-r-- 1 user user   844361  7月 24 16:48 201608_upair.txt
-rw-rw-r-- 1 user user  7367860  7月 24 16:48 201609_cwb_hr.txt
-rw-rw-r-- 1 user user   980635  7月 24 16:48 201609_upair.txt
-rw-rw-r-- 1 user user  7613380  7月 24 16:48 201610_cwb_hr.txt
-rw-rw-r-- 1 user user  1029069  7月 24 16:48 201610_upair.txt
-rw-rw-r-- 1 user user  7367860  7月 24 16:48 201611_cwb_hr.txt
-rw-rw-r-- 1 user user   839664  7月 24 16:48 201611_upair.txt
-rw-rw-r-- 1 user user  7613380  7月 24 16:48 201612_cwb_hr.txt
-rw-rw-r-- 1 user user   851132  7月 24 16:48 201612_upair.txt
-rw-rw-r-- 1 user user 55361553  7月 24 16:46 2018091217311031398.zip

[user@localhost cwb]$ rm 2018091217311031398.zip

這樣就完成下載檔案在 Linux 的解壓縮工作了!如果你不熟 Linux 的話,那麼,就先不要玩模式了!大部分的模式操作都是在 Linux 系統底下運作的! 如果不熟悉 Linux ,會很難入門啊!所以,基本功 (蹲馬步) 是需要的!請在本站慢慢的先學習 Linux 基礎吧!

  • 步驟二:根據模擬區域的環境,選擇正確的測站點,並分別擷取地面、探空資料

因為 AERMET 只會分析一個地面測站以及一個高空測站的資料,但是中央氣象局的資料,通常是全台的資料!所以,你得要先知道,你要取得的是哪一個縣市的地面測站! 你可以參考底下的連結:

舉例來說,我們可能想要知道高雄市附近的擴散模擬,因此,選擇了高雄測站的『 467440 』這一個站址,然後,探空資料只有兩個, 一個是台北 (466920) 一個是花蓮 (466990),畢竟我們是在西岸,所以,探空資料也只能選擇 466920 這一個測站位址。知道了測站資訊之後, 接下來就可以透過底下的方式來擷取所需要的探空與地面測站資料了!

  • 範例年份: 2016
  • 範例測站: 467440
    • 經度: 120.3157
    • 緯度: 22.5660
    • 海拔高度: 2.3
  • 範例探空: 466920
  • 範例地面觀測檔案: cwb_2016_surface
  • 範例探空資料檔案: cwb_2016_upair
# 先進入到剛剛建立的 aermet 目錄中,並擷取地面觀測資料
[user@localhost ~]$ cd ~/aermet
[user@localhost aermet]$ grep -h '^467440' ./cwb/2016*cwb*.txt | grep -v '\/\/\/' > cwb_2016_surface
[user@localhost aermet]$ grep -h '^466920' ./cwb/2016*up*.txt  | grep -v '\/\/\/' > cwb_2016_upair
[user@localhost aermet]$ ll
drwxrwxr-x 2 user user    4096  7月 24 16:49 cwb
-rw-rw-r-- 1 user user 2995344  7月 24 17:07 cwb_2016_surface  <==地面測站的觀測資料
-rw-rw-r-- 1 user user 5035855  7月 24 17:07 cwb_2016_upair    <==探空資料的觀測

有興趣的朋友可以使用 vim 去看一下每一個檔案的內容,否則就直接略過~這樣就已經處理好等等我們需要的中央氣象局的資料了!

  • 步驟三:從鳥站下載需要的轉檔程式,開始進行轉檔

接下來就是轉檔!透過鳥站下載程式後,直接將剛剛取得的資料轉檔即可。轉檔程式在這裡

[user@localhost ~]$ cd ~/aermet
[user@localhost aermet]$ wget https://linux.vbird.org/enve/images/cwb2aermet-static.exe
[user@localhost aermet]$ chmod a+x cwb2aermet-static.exe
[user@localhost aermet]$ ./cwb2aermet-static.exe cwb_2016_surface cwb_2016_upair SF_2016_467440.EXT UP_2016_466920.FSL 22.5660 120.3157 2.3
 地面測站檔名:cwb_2016_surface
 探空測站檔名:cwb_2016_upair
 完成探空資料
 開始處理地面資料的讀寫
 完成地面資料的處理

[user@localhost aermet]$ ll
drwxrwxr-x 2 user user    4096  7月 24 16:49 cwb
-rw-rw-r-- 1 user user 2995344  7月 24 17:34 cwb_2016_surface
-rw-rw-r-- 1 user user 5020300  7月 24 17:34 cwb_2016_upair
-rwxrwxr-x 1 user user 1116552  7月 24 11:54 cwb2aermet-static.exe
-rw-rw-r-- 1 user user 1352736  7月 24 17:36 SF_2016_467440.EXT     <==用於 AERMET 的地面測站資料
-rw-rw-r-- 1 user user 1020756  7月 24 17:36 UP_2016_466920.FSL     <==用於 AERMET 的探空資料

轉檔程式的語法也挺單純的:

[user@localhost ~]$ ./cwb2aermet-static.exe 地表觀測值 探空觀測值 地表輸出檔 探空輸出檔 測站緯度 測站精度 測站高度

這樣就準備好了 AERMET 要讀取的大氣資料囉!

進行 AERMET 的程式編譯

在 Linux 伺服器上面運作模式有許多的好處,除了你能夠讓 aermet 的資料與所有同事共享之外 (只要有帳號), 也可以從任何地方取得同一份模式運作完的資料!此外,大家統一使用 Server 的資源,也可以不用擔心會影響到你的 windows 工作用機器! 你也無須擔心,如果你的 windows 死掉,或者是需要系統升級,會導致磁碟內的模式資料遺失!最後, Linux 系統的穩定性與資源分配都比 windows 來的好!因此,要跑模式,首選就是 Linux。

可惜的是, US EPA 編譯完成的 AERMET 僅針對 Windows 系統,所以該程式碼無法在 Linux 上面運作。 為此,我們得要從 US EPA 下載原始碼 (純文字的程式碼檔案),並且重新編譯成為我們自己 Linux 可以運作的模樣, 這樣就能在我們的 Linux 上面運作了。整體流程也不難,鳥哥是以 CentOS 7 為範本進行編譯, 你可以依據你的系統來進行重新改寫的動作!

  • 步驟一:以管理員身份,安裝好所需要的編譯器以及函式庫

首先,你必須要取得管理員的身份才行!如果不清楚如何處理,請向你的 Linux server 管理員詢問。 假設你的身份已經是管理員了,接下來可以這樣進行任務:

# 直接使用 CentOS 的 yum 功能,直接安裝所需要的軟體即可
[root@localhost ~]# yum -y groupinstall "Development Tools"
[root@localhost ~]# yum -y install libgfortran-static and glibc-static wget
[root@localhost ~]# gfortran -v
使用內建 specs。
COLLECT_GCC=gfortran
COLLECT_LTO_WRAPPER=/usr/libexec/gcc/x86_64-redhat-linux/4.8.5/lto-wrapper
目的:x86_64-redhat-linux
配置為:../configure --prefix=/usr --mandir=/usr/share/man .....
執行緒模型:posix
gcc version 4.8.5 20150623 (Red Hat 4.8.5-36) (GCC)

這樣就搞定了我們編譯過程中所需要的各項參數了!

  • 步驟二:檔案下載與重新編輯編譯腳本

你可以自己前往 US EPA 的官網下載所需要的軟體,也可以透過底下的方式,直接使用你的 Linux server 線上直接下載原始碼。 鳥哥這裡使用的原始碼為 v18081 這個版本~直接這樣做即可:

# 
[user@localhost ~]$ cd ~/aermet
[user@localhost aermet]$ mkdir build
[user@localhost aermet]$ cd build
[user@localhost build]$ wget https://www3.epa.gov/ttn/scram/7thconf/aermod/aermet_source.zip
[user@localhost build]$ ll
-rw-rw-r--. 1 user user 475062  4月 27  2018 aermet_source.zip

[user@localhost build]$ unzip aermet_source.zip
[user@localhost build]$ vim compile.sh
#!/bin/bash

COMPILE_FLAGS=" -fbounds-check -Wuninitialized -Ofast -static -march=native -ffast-math -funroll-loops"
LINK_FLAGS=" -static -Ofast -march=native -ffast-math -funroll-loops"
sources="mod_AsosCommDates.for AERMET.FOR AERSURF.FOR AERSURF2.FOR ASOSREC.FOR AUDIT.FOR
	AUTCHK.FOR AVGCRD.FOR BANNER.FOR BULKRI.FOR CALMS.FOR CBLHT.FOR
	CHRCRD.FOR CHRCRD2.FOR CHROND.FOR CLHT.FOR CLMCRD.FOR CLOUDS.FOR
	COMPDT.FOR CUBIC.FOR CVG.FOR D028LV.FOR D144HD.FOR D144LV.FOR
	D3280H.FOR D3280L.FOR D6201H.FOR D6201L.FOR DATCRD.FOR DATER.FOR
	DEF256.FOR DEFINE.FOR DOCLDS.FOR DTCRD.FOR EQ_CCVR.FOR ERRHDL.FOR
	FDKEY.FOR FDPATH.FOR FETCH.FOR FLIWK1.FOR FLIWK2.FOR FLOPEN.FOR
	FLOS.FOR FLSDG.FOR FLSFC.FOR FLWRK1.FOR FLWRK2.FOR FMTCRD.FOR
	FNDCOMDT.FOR GEO.FOR GET620.FOR GETASOS.FOR GETCCVR.FOR GETFIL.FOR
	GETFLD.FOR GETFSL.FOR GETSFC.FOR GETTEMP.FOR GETWRD.FOR GMTLST.FOR
	GREG.FOR HDPROC.FOR HEADER.FOR HEAT.FOR HGTCRD.FOR HR0024.FOR
	HTCALC.FOR HTKEY.FOR HUMID.FOR HUSWX.FOR ICHRND.FOR INCRAD.FOR
	INTEQA.FOR INTHF.FOR ISHWX.FOR JBCARD.FOR LATLON.FOR LOCCRD.FOR
	LWRUPR.FOR MANDEL.FOR MDCARD.FOR MERGE.FOR MIDNITE.FOR MODEL.FOR
	MPCARD.FOR MPFIN.FOR MPHEAD.FOR MPMET.FOR MPOUT.FOR MPPBL.FOR
	MPPROC.FOR MPTEST.FOR MRCARD.FOR MRHDR.FOR MRPATH.FOR NETRAD.FOR
	NR_ANG.FOR NWSHGT.FOR OAUDIT.FOR OSCARD.FOR OSCHK.FOR OSDTCD.FOR
	OSDUMP.FOR OSFILL.FOR OSFILL2.FOR OSHRAV.FOR OSNEXT.FOR OSPATH.FOR
	OSPRNT.FOR OSQACK.FOR OSQAST.FOR OSRANGE.FOR OSREAD.FOR OSSMRY.FOR
	OSSUMS.FOR OSSWAP.FOR OSTEST.FOR OSTRA.FOR OSWRTE.FOR OTHHDR.FOR
	P2MSUB.FOR PRESET.FOR PTAREA.FOR PTGRAD.FOR RDHUSW.FOR RDISHD.FOR
	RDLREC.FOR RDSAMS.FOR READRL.FOR REALQA.FOR RHOCAL.FOR RNGCRD.FOR
	SAMWX.FOR SAUDIT.FOR SBLHT.FOR SCNGEN.FOR SECCRD.FOR SECCRD2.FOR
	SETHUS.FOR SETSAM.FOR SETUP.FOR SFCARD.FOR SFCCH.FOR SFCCH2.FOR
	SFCCRD.FOR SFCCRD2.FOR SFCHK.FOR SFCWXX.FOR SFEXST.FOR SFEXT.FOR
	SFPATH.FOR SFQASM.FOR SFQAST.FOR SFQATM.FOR SFTRA.FOR SMTHZI.FOR
	STONUM.FOR SUBST.FOR SUMHF.FOR SUMRY1.FOR SUMRY2.FOR SUNDAT.FOR
	TDPEST.FOR TEST.FOR THRESH1MIN.FOR UACARD.FOR UACHK.FOR UAEXST.FOR
	UAEXT.FOR UAMOVE.FOR UAPATH.FOR UAQASM.FOR UAQAST.FOR UATRA.FOR
	UAUDIT.FOR UAWNDW.FOR UCALCO.FOR UCALST.FOR VALCRD.FOR VARCRD.FOR
	VRCARD.FOR WRTCRD.FOR YR2TOYR4.FOR YR4TOYR2.FOR XDTCRD.FOR XTNDUA.FOR"

for file in ${sources}
do
        gfortran -m64 -c ${COMPILE_FLAGS} ${file}
done

gfortran -m64 -o aermet.exe ${LINK_FLAGS}  mod_AsosCommDates.o AERMET.o AERSURF.o  AERSURF2.o ASOSREC.o AUDIT.o  AUTCHK.o  AVGCRD.o  \
        BANNER.o  BULKRI.o  CALMS.o  CBLHT.o  CHRCRD.o  CHRCRD2.o  CHROND.o  CLHT.o  CLMCRD.o  CLOUDS.o  \
        COMPDT.o  CUBIC.o  CVG.o  D028LV.o  D144HD.o  D144LV.o  D3280H.o  D3280L.o  D6201H.o  D6201L.o   \
        DATCRD.o  DATER.o  DEF256.o  DEFINE.o  DOCLDS.o  DTCRD.o  EQ_CCVR.o  ERRHDL.o  FDKEY.o FDPATH.o  \
        FETCH.o  FLIWK1.o  FLIWK2.o  FLOPEN.o  FLOS.o  FLSDG.o  FLSFC.o  FLWRK1.o  FLWRK2.o  FNDCOMDT.o  \
        FMTCRD.o  GEO.o  GET620.o  GETASOS.o  GETCCVR.o  GETFIL.o  GETFLD.o  GETFSL.o  GETSFC.o  GETTEMP.o  GETWRD.o  GMTLST.o  GREG.o  HDPROC.o  HEADER.o  HEAT.o \
        HGTCRD.o  HR0024.o  HTCALC.o  HTKEY.o  HUMID.o  HUSWX.o  ICHRND.o  INCRAD.o  INTEQA.o  INTHF.o  ISHWX.o  JBCARD.o \
        LATLON.o  LOCCRD.o  LWRUPR.o  MANDEL.o  MDCARD.o  MERGE.o  MIDNITE.o  MODEL.o  MPCARD.o  MPFIN.o  MPHEAD.o  MPMET.o \
        MPOUT.o  MPPBL.o  MPPROC.o  MPTEST.o  MRCARD.o  MRHDR.o  MRPATH.o  NETRAD.o  NR_ANG.o  NWSHGT.o  OAUDIT.o  OSCARD.o \
        OSCHK.o  OSDTCD.o  OSDUMP.o  OSFILL.o  OSFILL2.o  OSHRAV.o  OSNEXT.o  OSPATH.o  OSPRNT.o  OSQACK.o  OSQAST.o  OSRANGE.o  OSREAD.o  OSSMRY.o \
        OSSUMS.o  OSSWAP.o OSTEST.o  OSTRA.o  OSWRTE.o  OTHHDR.o  P2MSUB.o  PRESET.o  PTAREA.o  PTGRAD.o  RDHUSW.o  RDISHD.o  RDLREC.o \
        RDSAMS.o  READRL.o  REALQA.o  RHOCAL.o  RNGCRD.o  SAMWX.o  SAUDIT.o  SBLHT.o  SCNGEN.o  SECCRD.o  SECCRD2.o  SETHUS.o \
        SETSAM.o  SETUP.o  SFCARD.o  SFCCH.o  SFCCH2.o  SFCCRD.o  SFCCRD2.o  SFCHK.o  SFCWXX.o  SFEXST.o  SFEXT.o  SFPATH.o  \
        SFQASM.o  SFQAST.o  SFQATM.o  SFTRA.o  SMTHZI.o  STONUM.o  SUBST.o  SUMHF.o  SUMRY1.o  SUMRY2.o  SUNDAT.o  TDPEST.o \
        TEST.o  THRESH1MIN.o UACARD.o  UACHK.o  UAEXST.o  UAEXT.o  UAMOVE.o  UAPATH.o  UAQASM.o  UAQAST.o  UATRA.o  UAUDIT.o  UAWNDW.o  \
        UCALCO.o  UCALST.o  VALCRD.o  VARCRD.o  VRCARD.o  WRTCRD.o  YR2TOYR4.o  YR4TOYR2.o  XDTCRD.o  XTNDUA.o

[user@localhost build]$ sh compile.sh
[user@localhost build]$ ll aermet.exe
-rwxrwxr-x. 1 user user 2900600  7月 25 00:16 aermet.exe

[user@localhost build]$ cd ~/aermet/
[user@localhost aermet]$ ln -s build/aermet.exe .
[user@localhost aermet]$ ll
lrwxrwxrwx. 1 user user      16  7月 25 00:21 aermet.exe -> build/aermet.exe   <==這就是重點!
drwxrwxr-x. 2 user user   12288  7月 25 00:16 build
drwxrwxr-x. 2 user user    4096  7月 24 16:49 cwb
-rw-rw-r--. 1 user user 2995344  7月 24 23:57 cwb_2016_surface
-rw-rw-r--. 1 user user 5020300  7月 24 23:57 cwb_2016_upair
-rwxrwxr-x. 1 user user 1116552  7月 24 11:54 cwb2aermet-static.exe
-rw-rw-r--. 1 user user 1352736  7月 24 23:57 SF_2016_467440.EXT
-rw-rw-r--. 1 user user 1020756  7月 24 23:57 UP_2016_466920.FSL

這樣就編譯好了 AERMET 主程式碼!接下來的三個步驟,全部都需要用到這個傢伙喔!

step 2: 處理 AERMET 的 stage1 流程

如同上述的資料,要完成 AERMET 的處理,需要有三個基本的階段 (stage),第一個階段是解析出我們需要的數據, 第二個階段則是整合地面資料與探空資料成為唯一的一個氣象檔,第三個階段則是處理出地面資料+探空資料兩個 AERMOD ready 的檔案。 全部的三個階段都是使用剛剛我們編譯成功的 aermet.exe 來執行的。那 aermet 怎麼知道要處理那一個階段呢?這就需要我們撰寫設定檔來處置了。

至於設定檔的檔名,通常副檔名會指定為 .inp 就是了,而設定檔內容大致分為 6 個部份,這 6 部份可稱為設定關鍵字或路徑 (pathway):

  • JOB:說明整體運作的資訊
  • UPPERAIR:處理探空資料 (sounding data) 的解析與數據保證工作
  • SURFACE:處理地面觀測站逐時氣象資料
  • ONSITE:使用者提供的測站當地大氣資訊
  • MERGE:整合大氣資料
  • METPREP:評估給 AERMOD 使用的邊界層參數

假設以我們這次模擬的 2016 年全年高雄地區的 AERMET stage1 的設定為例,它的設定內容,有點類似這樣:

JOB
    MESSAGES       MSG.STAGE1
    REPORT         RPT.STAGE1

UPPERAIR
    DATA           UP_2016_466920.FSL FSL
    EXTRACT        UP_2016_466920.EXT
    QAOUT          UP_2016_466920.QA
    XDATES         2016/01/01 TO 2016/12/31
    LOCATION       466920 25.04N 121.51E -8
    AUDIT          UAPR UAHT UATT UATD UAWD UAWS

** 這底下是關於地面測站,其中 SF_2016_467440.EXT 就是直接從氣象局數據轉換得來的檔案
SURFACE
    EXTRACT        SF_2016_467440.EXT
    QAOUT          SF_2016_467440.QA
    XDATES         2016/01/01 TO 2016/12/31
    LOCATION       467440 22.566N 120.316E 0 2.3

基本語法說明如下:

  • 每一個路徑在同一個設定檔中,基本上只能出現一次,且該路徑全部的細部設定,都需要連續寫在同一個區塊中。以上述的內容為例, 若以 JOB 來看,它裡面有兩個設定,分別是 MESSAGES 以及 REPORT,那麼這兩個設定就一定要在該 JOB 底下!不能寫在其他地方,是這樣的意思。
  • 每一行設定值,最大只能到 132 個字元,這 132 個字元包含前面的空白字元喔!
  • 可以擁有空白行,為了讓整個設定檔更容易瞭解。
  • 在每一行的第 1, 2 個字元,如果出現兩個星號 (**) 時,那一行就是註解。因此,如果你針對該行有什麼特別的說明,擔心自己忘記時, 就可以使用這樣的方式處理。但請注意,不是每一行出現的第 1, 2 個字串喔,是每一行第 1, 2 個字元,就是行首的意思! 例如上述範例中的說明喔!
  • 英文字可以是大寫也能是小寫,因為預設是給 Windows 系統運作的。不過,因為我們將之改版成為 Linux 版本,Linux 大小寫是不同的。 因此,在 AERMET 的說明文件中,它建議檔名最好都使用大寫!這是因為 AERMET 會主動的將設定檔內的文字都自動更改為大寫, 所以,如果檔名使用小寫英文的話,AERMET 會找不到檔名!所以,都用大寫吧!
  • 在每一個 pathway (如 JOB, UPPERAIR) 裡面,有很多的設定值 (keyword),某些設定值是主要的,一定要存在於 pathway 裡面, 某些則是參數類型,可有可無。某些設定是可以重複出現的!底下我們在 stage3 時,就會看到重複出現的一些設定值喔!

再回頭來看一下,我們現在有的檔案有兩個,其中探空資料是尚未經過 QA 的類型,但是地面資料則已經是經過 QA 之後的類型了。 也就是說:

  • UP_2016_466920.FSL :內容是尚未經過 QA 的探空資料型態,格式為 FSL 格式;
  • SF_2016_467440.EXT :內容是已經經過 QA 的地面觀測值,所以直接用在 AERMET 即可。

現在,讓我們來建立這個 stage1.inp 的輸入檔案:

[user@localhost ~]$ cd ~/aermet
[user@localhost aermet]$ vim stage1.inp
** 這個部份在說明整體工作的意思,訊息輸出在 MSG.STAGE1 的內容,報告則在 RPT.STAGE1 囉!
JOB
    MESSAGES       MSG.STAGE1
    REPORT         RPT.STAGE1

** 這個部份在講解地面觀測值!其中 SF_2016_467440.EXT 是鳥哥程式轉換過來的,AERMET 自動生成 SF_2016_467440.QA
SURFACE
    EXTRACT        SF_2016_467440.EXT
    QAOUT          SF_2016_467440.QA
    XDATES         2016/01/01 TO 2016/12/31
    LOCATION       467440 22.566N 120.316E 0 2.3
    AUDIT          PRCP SLVP PRES TMPD DPTP RHUM WDIR WSPD

** 這個部份在講探空資料,其中 UP_2016_466920.FSL 是經過鳥哥程式轉換過來的,AERMET 會主動生成 UP_2016_466920.EXT 與 UP_2016_466920.QA
UPPERAIR
    DATA           UP_2016_466920.FSL FSL
    EXTRACT        UP_2016_466920.EXT
    QAOUT          UP_2016_466920.QA
    XDATES         2016/01/01 TO 2016/12/31
    LOCATION       466920 25.04N 121.51E -8
    AUDIT          UAPR UAHT UATT UATD UAWD UAWS

上述的檔案中,除了 JOB 的內容大致上就是轉檔過程中的錯誤訊息告知 (例如原始資料是否有缺值、是否有其他數值不合理而改變等等), 在地面觀測值與探空資料的說明如下:

  • DATA:指的就是 AERMET 所需要的觀測值未格式化之前的數據,因為我們的探空資料只能做成符合 AERMET 的 INPUT 資料, 因此,只有在 UPPERAIR 的路徑裡面擁有這個資訊而已,而且,主要的格式為 FSL 喔!
  • EXTRACT:解開我們需要的資料後,轉出的檔案。因為地面觀測值 (SURFACE) 的資料已經是從中央氣象局轉換過來的, 所以,地面測站的資料就是從中央氣象局過來的資料,無須重新解開。至於探空資料則需要重新解開取得所需要的內容。 因此,在 UPPERAIR 的部份,這個檔名得要重新另存一個新檔,通常副檔名就會是 .EXT (EXTRACT) !
  • QAOUT:進行數據品質保證工作的過程檔,也是 AERMET 自己生成的!通常副檔名會取 .QA
  • XDATES:就是取得測站的資料的年月日!通常都是整年的數據!因此,就給開頭 (2016/01/01) 到結尾 (2016/12/31) 的數據即可! 中間使用 TO 來連結~
  • LOCATION:就是測站點的位置!所以,請繼續前往底下的網站,找出來探空資料以及地面測站的經緯度座標。LOCATION 的內容比較需要注意:
    • 氣象局測站資料:https://e-service.cwb.gov.tw/wdps/obs/state.htm
    • 第一個數值就是測站代碼,測站代碼請由上面的連結找出來!測站代碼不能寫錯,寫錯的話,AERMET 會找不到數據喔!
    • 第二個數值是緯度,需要填寫北緯或南緯。台灣在北緯,所以才會是 22.566N 或 25.04N 之類的。
    • 第三個數值是精度,需要填寫東經會西經。台灣在東經,所以才會是 120.316E 或 121.51E。
    • 第四個數值是時間轉換。地面觀測資料原本就是台灣地區本地時間,所以無須轉換。但是探空資料全球都統一使用格林威治時間 (GMT or UTC), 因此需要轉換時間成為台灣時間。不過由於 AERMET 是美國環保署發展的,為了方便他們自己,所以時區全部為正值。 台灣時間則需要填寫為負值 (-8) 才行!這個第一次接觸 AERMET 的朋友,可能會完全搞混啦!要注意!要注意!
    • 第五個數值則是地表高度,這個數值只對地面測站有意義,探空資料的這個數值會被忽略!
  • AUDIT:剛剛上面不是講到 QA (數據的品質保證) 嘛?那是針對哪些資料進行呢?資料的簡寫可以參考 AERMET 使用者手冊 B-1, B-2 的表格。 我們這裡僅列出我們有興趣的部份資料檢測而已:
    • PRCP:降水量,數值為: mm * 1000
    • SLVP:水平面氣壓,單位為 hpa*10
    • PRES:測站點的氣壓,單位 hpa*10
    • TMPD:乾球溫度 °C * 10
    • DPTP:露點溫度 °C * 10
    • RHUM:相對濕度,單位 %
    • WDIR:風向,360 度角,北風為 0 度,東風 90 度這樣。
    • WSPD:風速,單位 m/s * 10
    • UAPR:大氣壓力, hpa * 10
    • UAHT:高度, m
    • UATT:乾球溫度 °C * 10
    • UATD:露點溫度 °C * 10
    • UAWD:風向,360 度角
    • UAWS:風速 m/s*10

這個設定檔設定妥當之後,接下來就是需要執行!執行過程也是很簡單!下達這樣的指令即可:

[user@localhost aermet]$ ./aermet.exe stage1.inp
   AERMET Version  18081

   Pre-processing the Setup Information

   Processing the Setup Information

              ********************************************************
              ***        AERMET Setup Finished Successfully        ***
              ********************************************************

+  Stage 1: Extracting upper air data for month/day/year 01/01/16
+  Stage 1: Extracting upper air data for month/day/year 01/01/16
.....
+  Stage 1: Extracting upper air data for month/day/year 12/31/16
+  Stage 1: Extracting upper air data for month/day/year 12/31/16

+  Stage 1: QA'ing upper air data for month/day/year 01/01/16
+  Stage 1: QA'ing upper air data for month/day/year 01/01/16
....
+  Stage 1: QA'ing upper air data for month/day/year 12/31/16
+  Stage 1: QA'ing upper air data for month/day/year 12/31/16

+  Stage 1: QA'ing surface data for month/day/year 01/01/16
+  Stage 1: QA'ing surface data for month/day/year 01/02/16
....
+  Stage 1: QA'ing surface data for month/day/year 12/31/16
   Processing completed; writing summary files


              ********************************************************
              ***   AERMET Data Processing Finished Successfully   ***
              ********************************************************

   The Summary Report Generated by AERMET Is In:
     RPT.STAGE1
# 對!這樣就跑完了!過程會有相當相當多的資訊,每天會跑出一個資訊!所以要注意查看。

[user@localhost aermet]$ ll
-rw-rw-r--. 1 user user   10346  7月 26 15:03 MSG.STAGE1           <==主要的訊息輸出檔
-rw-rw-r--. 1 user user   13091  7月 26 15:03 RPT.STAGE1           <==成果報告檔
-rw-rw-r--. 1 user user 1352736  7月 24 23:57 SF_2016_467440.EXT   <==中央氣象局地面觀測轉出檔
-rw-rw-r--. 1 user user 1353656  7月 26 15:03 SF_2016_467440.QA    <==AERMET QA 後的檔案
-rw-rw-r--. 1 user user     987  7月 26 14:51 stage1.inp           <==我們剛剛上頭撰寫的參數設定檔
-rw-rw-r--. 1 user user  767656  7月 26 15:03 UP_2016_466920.EXT   <==經由 FSL 轉檔過來的資料
-rw-rw-r--. 1 user user 1020756  7月 24 23:57 UP_2016_466920.FSL   <==中央氣象局取得的探空資料
-rw-rw-r--. 1 user user  767620  7月 26 15:03 UP_2016_466920.QA    <==AERMET QA 後的檔案

有興趣的話,我個人建議你一定要去 MSG.STAGE1 看看,因為所有的可能錯誤訊息都在裡面!比較多的資料都是『缺值』! 例如中央氣象局的探空資料,可能因為某些緣故有缺乏~缺乏資料時,AERMET 會略過啦!

step 3: 處理 AERMET 的 stage2 流程

Stage2 的流程很簡單,就是要將剛剛建立好的經過 QA 的檔案整合成為一個氣象場檔案資料!連設定檔的內容都很單純!

[user@localhost aermet]$ vim stage2.inp
JOB
    MESSAGES       MSG.STAGE2
    REPORT         RPT.STAGE2

UPPERAIR
    QAOUT          UP_2016_466920.QA  <==由 stage1 產生的 QA 檔案

SURFACE
    QAOUT          SF_2016_467440.QA  <==由 stage1 產生的 QA 檔案

MERGE
    OUTPUT         MERGE.STAGE2       <==預計要建立的唯一氣象場檔案
    XDATES         2016/01/01 TO 2016/12/31

比較重要的是 MERGE 這個 pathway,這個路徑主要在處理 stage2 的整合工作!最重要的就是輸出 (OUTPUT) 這個設定值就是了! 所以整體設定很單純!接下來,你只要這樣執行即可:

[user@localhost aermet]$ ./aermet.exe stage2.inp
   AERMET Version  18081

   Pre-processing the Setup Information

   Processing the Setup Information

              ********************************************************
              ***        AERMET Setup Finished Successfully        ***
              ********************************************************

+  Stage 2: Merging month/day/year 01/01/16
+  Stage 2: Merging month/day/year 01/02/16
.....
+  Stage 2: Merging month/day/year 12/30/16
+  Stage 2: Merging month/day/year 12/31/16
   Processing completed; writing summary files


              ********************************************************
              ***   AERMET Data Processing Finished Successfully   ***
              ********************************************************

   The Summary Report Generated by AERMET Is In:
     RPT.STAGE2

[user@localhost aermet]$ ll
-rw-rw-r--. 1 user user 4072136  7月 26 15:14 MERGE.STAGE2
-rw-rw-r--. 1 user user     543  7月 26 15:14 MSG.STAGE2
-rw-rw-r--. 1 user user   19539  7月 26 15:14 RPT.STAGE2
-rw-rw-r--. 1 user user     240  7月 26 15:14 stage2.inp

新產生的資料也不過就是上述的檔案而已。比較重要的是第一個 MERGE.STAGE2 檔名!這個是等等 stage3 會用到的檔案喔!

step 4: 處理 AERMET 的 stage3 流程

AERMET 的第三階段,就是要產生 AERMOD 能夠讀取的氣象資料了!而當中最重要的,大概就是 PBL 的相關資訊!其中就以模擬區域的地表特徵最重要! 我們在前一章節有提到流體力學的特性,在表面的時候,可能會有流速為 0 的情境發生!但是這個流體的情況與表面粗糙度有關! 因此,在設計模擬區域的特徵時,就得要提供這方面的模擬表面資料才行。

基本上,模擬範圍內的地表面,不可能都是一模一樣的!因此,US EPA 有另外提供一個名為 AERSURFACE 的模組來處理! 不過,AERSURFACE 還是需要其他地形特徵的輸入參數才可以運作,因此,這裡我們還是將模擬區域假設為單一樣貌的地表特徵! 這樣比較好處理。

台灣地區的土地利用型態研究,透過環保署的專案計畫結果,可以使用底下的資料來設定喔!

環保署模式技術支援計畫提供的參考值
環保署模式技術支援計畫提供的參考值 (ps3, page 4-209)

如果以我們這個案例當中想要模擬的高雄市區來說,應該是可以選擇 7 號都市區的設定值。那麼就來瞧一瞧 stage3.inp 的內容應該要有哪些設定吧!

[user@localhost aermod]$ vim stage3.inp
JOB
    MESSAGES       MSG.STAGE3
    REPORT         RPT.STAGE3

METPREP
    DATA           MERGE.STAGE2
    XDATES         2016/01/01 TO 2016/12/31
    MODEL          AERMOD
    METHOD         WIND_DIR RANDOM
    METHOD         REFLEVEL  SUBNWS
    METHOD         STABLEBL BULKRN
    METHOD         ASOS_ADJ NO_ADJ
    METHOD         UASELECT SUNRISE
    NWS_HGT        WIND      14

    OUTPUT         aermet_2016_467440.SFC
    PROFILE        aermet_2016_466920.PFL

    FREQ_SECT      ANNUAL 1
    SECTOR         1    0   360
    SITE_CHAR      1   1  0.16  0.8  1.2

METPREP 這個 pathway 主要是在進行大氣資訊的評估計算,裡面的設定值大概有這些:

  • DATA:就是來自 stage2 整合的那個唯一氣象資料檔案!
  • XDATES:就是 2016 一整年
  • MODEL:主要的輸出資料是給 AERMOD 用的!這一般也是預設值
  • METHOD:主要的大氣特徵評估方式,若需要詳細的資訊,則請參考 AERMET 手冊
  • NWS_HGT:主要就是測站風速採樣點的高度,一般採樣口高度大概在 10m 左右!
  • OUTPUT:就是地面測站的觀測資料以及相關 PBL 的資訊。
  • PROFILE:就是高空資料的風速風向等大氣資訊
  • FREQ_SECT:就是地表資料的設定值,後續的 ANNUAL 指全年一樣,且設定號碼為 1 號
  • SECTOR:模擬區域的分區,針對上述的 1 號區域設計,範圍為 0~360 度,代表全區!
  • SITE_CHAR:針對第 1 號的第 1 個分區,設計的三個數值。數值請參考上表!

設計妥當之後,自然就是開始跑囉!執行過程也是很簡單:

[user@localhost aermet]$ ./aermet.exe stage3.inp
   AERMET Version  18081

   Pre-processing the Setup Information

   Processing the Setup Information

              ********************************************************
              ***        AERMET Setup Finished Successfully        ***
              ********************************************************

+  Stage 3: Now processing month/day/year 01/01/16
+  Stage 3: Now processing month/day/year 01/02/16
....
+  Stage 3: Now processing month/day/year 12/30/16
+  Stage 3: Now processing month/day/year 12/31/16
   Processing completed; writing summary files


              ********************************************************
              ***   AERMET Data Processing Finished Successfully   ***
              ********************************************************

   The Summary Report Generated by AERMET Is In:
     RPT.STAGE3

[user@localhost aermet]$ ll
-rw-rw-r--. 1 user user  579744  7月 27 00:48 AERMET_2016_466920.PFL  <==這個是高空氣象場
-rw-rw-r--. 1 user user 1554894  7月 27 00:48 AERMET_2016_467440.SFC  <==這個是地面測站資料喔!
-rw-rw-r--. 1 user user   60385  7月 27 00:48 MSG.STAGE3
-rw-rw-r--. 1 user user   23565  7月 27 00:48 RPT.STAGE3
-rw-rw-r--. 1 user user     572  7月 27 00:22 stage3.inp

如上所示,跑出來的最終兩個檔案,就是我們需要的 AERMOD 的輸入參數了!就這麼簡單喔!

風花圖的製作

完成了 AERMET 的輸出之後,再來就是最好拿個分析軟體,探索一下風向變化是否符合實際的環境!畢竟,許多轉換軟體, 或者是自己撰寫的軟體,可能會有部份錯誤的資料轉換。鳥哥寫的第一版轉檔資料,也是出現過數據不符合的情況 (這是因為 AERMET 是給美國本土的數據處理,跟我們認知的中央氣象局的數具有點不同!),跑完 AERMOD 害我嚇一跳!所以, 還是分析一下比較好。最簡單的分析,就是探索風向的風花圖!

繪製風花圖最簡單的軟體就是由 Lakes Environmental Software 所開發的 WRPLOT View 這個軟體! 相關的連結位置在底下的超連結上面,請自行前往下載正確的檔案來安裝處理。不過, WRPLOT View 是 windows 的軟體, 這意味著,您得要將剛剛在 Linux 上面產生的檔案下載到你的 windows 系統上面才行。

軟體下載完畢之後,先前往該上述連結去按下『 Registration 』,並進入到註冊的頁面,開始進行註冊! 註冊完畢之後,Lakes 公司會提供一組驗證碼給妳,該驗證碼可以讓你的軟體操作使用一年 (詳細情況,請參考 Lakes 提供給您的認證信的說明)。之後就開始在你的 windows 上面點擊兩下剛剛下載的軟體, 安裝妥當之後,直接點選『 WRPLOT View 』,該軟體會提醒『Do you want to register WRPLOT Viewer now?』, 點選『Yes』之後,輸入剛剛收到的信件內容的啟動碼,這樣就可以開始使用該軟體了。

如果經過一年之後,還想使用這個軟體呢?那就請再次的前往上述的網站註冊 (register),然後再次的啟動該軟體, 就又可以再用一年喔!軟體啟動之後會有點像這樣:

WRPLOT viewer 軟體操作與執行
WRPLOT viewer 軟體操作與執行

如上,選擇『 Add File 』的按鈕,並將剛剛我們建立好的 aermet 屬於 SFC 的檔名載入即可!載入之後, 就會出現如下的圖示!注意看到測站的代碼以及日期有沒有問題!同時,預設的情況下,風花圖會是全年的平均! 如果你只想要看某些日期 (不同的季節),就得要在底下的 Date Range 項目底下去指定了!

WRPLOT viewer 軟體操作與執行
WRPLOT viewer 軟體操作與執行

如果一切都 OK 沒啥大問題,你就可以點選『 Wind Rose 』的按鈕,即可看到如下的風花圖的圖示囉!相當有趣! 如果還想要知道其他的頻率相關問題,就直接點選其他的頁面查閱即可!

WRPLOT viewer 軟體操作與執行
WRPLOT viewer 軟體操作與執行

這個軟體非常簡單!希望大家用的開心!

參考資料

其他連結
環境工程模式篇
鳥園討論區
鳥哥舊站

今日 人數統計
昨日 人數統計
本月 人數統計
上月 人數統計