大氣化學(xué)模式的數值計算方法

環(huán)境科學(xué)研究第13卷第1期esearch of environmental sciencesVol.13No.12000大氣化學(xué)模式的數值計算方法劉峻峰,李金龍,白郁華北京大學(xué)技術(shù)物理系北京100871)摘要∶大氣化學(xué)模式計算的一個(gè)重要研究熱點(diǎn)是如何快速、精確地對相互關(guān)聯(lián)的非線(xiàn)性常微分方程組(OEs進(jìn)行數值積分求解。由于大氣中督物種壽命長(cháng)短相差很大,導致ωDEs具有了很強的剛性因而在數值求算過(guò)程中存在著(zhù)計算精度和計算效率的平衡問(wèn)題。筆者介紹了當前幾種常見(jiàn)的大氣化學(xué)模式的計算方法包括 GEAR QSSA( MQSSA)Y&BMY&B),RADM, LSODE SMVGEAR小參數法 Gong and Cho及 TWOSTEP等并結合當前本領(lǐng)域內的部分研究成果對一些重要的方法進(jìn)行了比較和分析。關(guān)鍵詞:大氣化學(xué)模式;剛性常微分方程組;數值積分中圖分類(lèi)號:X11文獻標識碼:A文章編號:1001-692×200001-0044-06The Numerical Solvers for Atmospheric Chemistry ModelsLIU Jun-feng, LI Jin-long, BAI Yu-huaDepartment of Technical Physics Peking University Beijing 100871)Abstract: The numerical integration of nonlinear systems of ordinary differential equations( ODEs )is a major problem to atmosphechemistry models. Because of the great disparities in the rates of the atmospheric chemistry reactions it is importantd a properway to accurately and efficiently solve the stiff systems of ODEs. In this paper several solvers including GEAR QSSA( MQSSA)Y&E MY&B) RADM SODE SMVGEAR SPM ong and Cho and TWOSTEP aare presented and compared the based on current research achievements in this fieldKey words: atmospheric chemistry model stiff ODE numerical integration自70年代以來(lái),模式研究已成為大氣環(huán)境化學(xué)領(lǐng)域中或(1)代表了一組數量較大、非線(xiàn)性且互相關(guān)聯(lián)的常微不可缺少的組成部分。一—個(gè)性能優(yōu)異的大氣化學(xué)模式應具分方程細ODEs)故無(wú)法求出其解析解,只能在給定每種化備2個(gè)必要條件⑩合理的模式機理。此機理能較準確地反合物初值后用各物種濃度對時(shí)間數值積分的方法,求岀各映實(shí)際大氣中物種數量及彼此間的化學(xué)反應同時(shí)給出較準物種在特定時(shí)刻的濃度。確的動(dòng)力學(xué)數據。②高效精確的計算方法。該方法在滿(mǎn)足最早的方法是對具有初值的式(1)在固定時(shí)間步長(cháng)的條定誤差范圍的前提下應盡可能快速地計算出最終結果件下進(jìn)行 Euler積分。但由于各反應物種在大氣中的壽命長(cháng)根據研究的對象及目的不同,人們設計出了不同的機短不一,物種間的反應速率相差高達幾個(gè)數量級,造成了理。常見(jiàn)的幾種大氣化學(xué)模式機理及計算方法見(jiàn)表1。ODEs的剛性。為了保證求解方程組的準確性和穩定性由表1可以看出,般機理包含的化合物種類(lèi)由幾十種須將時(shí)間步長(cháng)定得很小造成計算效率極低。因此在數值積到上百種反應方程式由幾百個(gè)到上千個(gè)。由這些化合物及分過(guò)程中存在著(zhù)積分時(shí)間步長(cháng)同各物種的壽命是否相匹配反應方程式根據質(zhì)量作用定律可以寫(xiě)出每種化合物的反的問(wèn)題。若積分時(shí)間步長(cháng)取的較大,對于壽命短的化合物會(huì )應速率方程產(chǎn)生較大誤差時(shí)間步長(cháng)取得過(guò)小,又會(huì )使計算過(guò)慢造成dh=F(C)=P(Ct)-L(CC(1)對計算資源的浪費中國煤化工中新數值求解OEs的方法式中:P(Cd)第i種化合物的產(chǎn)生速率HCNMHG應用,但是周看人們研究物L(fēng)(C,)第i種化合物的去除速率。質(zhì)量模式中所取的網(wǎng)絡(luò )數成干上萬(wàn),一般又要模擬較長(cháng)收稿日期:1999-09-14段時(shí)間的反應因而對計算速度的要求也在提高,即需要作者簡(jiǎn)介劉峻峰(1974-)男北京市人碩士研究生種快速且較為精確的計算方法。GEAR法逐漸成為人們對新研究的快速方法精確度的檢驗手段。1977年Y&B法161第1期劉峻峰等大氣化學(xué)模式的數值計算方法45和1978年QSsA"擬穩態(tài)假設)法的提出標志著(zhù)快速計算較18-201,并提出了一批新的或是改進(jìn)的方法,如法的誕生。此后直到現在又有許多科學(xué)家對大氣化學(xué)模式 LSODEL211,RADM21, Gong andTWOSTEP 241中剛性ODEs數值解法的效率和精度進(jìn)行了研究和比 SMVGEAR25] MQSSA和MY&B201表1幾種常見(jiàn)大氣化學(xué)模式及其計算方法Table 1 Atmospheric chemistry models and their numerical solvers機理縮寫(xiě)名參考物種數反應方程數計算法RADM2-IFUStockwell et al.[(1990)RADM2-FZKStockwell et al. [( 1990)RADMRADM2-KFAStockwell et al. [I(1990)Implicit EulaEuro-RADMStockwell and Kley?( 1994)CBM-ⅣLOTOGery et al. [3( 1989 )Builties+( 1992)63CBM-ⅣCBM-N-TNOGery, et al. 3(1989 Sterr'( 1994)GearC4.1Gery et al. [3(1989EMEPEMEPSimpson et al. [7( 1993), Simpson&( 1995140MCMPilling et al. [(1998)>2500>7000ADOM-ⅡlADOMLurmann et al. 10( 1986)47114Strand and HoyI( 1994)Andersson-Skold 12(1995RuhnkeRuhnke 13( 1995)145342資料主要來(lái)源 Kuhn et al.41998筆者在前面簡(jiǎn)要介紹了GEAR, LSODE, SMVGEARhF(Cm)Cm-1為前一步的C值,CQSSA, MQSSA,Y&B,MY&B,小參數法, RADM Gong and為本步的C值Cho及 TWOSTEP等計算方法。在下面結合各方面的研究對這些方法進(jìn)行了分類(lèi)并對重要方法的精確性和計算效率n——計算步數做了比較和分析每步校正次數1方法介紹積分開(kāi)始時(shí),Z0=(CC′o),C0為初值,C'o=dC/d1.1GEAR法和改進(jìn)GEAR法(C)第一步步長(cháng)取值較小P取1以后根據計算的誤GEAR法及與GEAR法類(lèi)似的 LSODE, LSODES,差大小可按一定的方法自動(dòng)調整h和P。SMVGEAR等代表了一類(lèi)較精確的數值計算方1.1.2 LSODE與 LSODES1.1.1GEAR法LSODE( Livermore solver for ordinary differential equaGEAR數值計算法151具有適用范圍廣、可以自動(dòng)調節ons521是建立在GEAR法基礎上,被廣為應用計算ODEs時(shí)間步長(cháng)與誤差,以及計算結果準確度高等特點(diǎn),是一種的較精確方法。此方法允許使用者根據ODEs的剛性程度為理想的計算方法。但其最大的缺點(diǎn)是計算步驟繁瑣要頻選擇Adam多步計算公式(對于非剛性體系)或BDF公式繁計算 Jacobian行列式,當計算的物種比較多時(shí)會(huì )耗費大 backward difference multistep formulas對于剛性體系)大量的機時(shí)和內存。GEAR法的數值計算公式如下氣化學(xué)模式的ODEs體系剛性較大故多用BDF法。BD)F的a.預測公式計算公式為(2)C)+△tbf(Cb.校正公式式中:C"+1—當前t"1的濃度Zm+1=Zm+1-[(l1-lnh)(xZm+1)](3)中國煤化工C.最終公式(4)CNMHG邇以刊酊炬陣結構稀少的不同來(lái)源(是式中zn—(k)×(P+1矩陣來(lái)自于用戶(hù)直接輸入,還是來(lái)自于調用 Jacobian子程序或是A-(P+1)×(P+1)矩陣l≤P≤6多次調用方程評估子程序)分別采用不同方法簡(jiǎn)化計算步Ⅰ——單位矩陣驟提高計算效率26。(xm})k×1矩陣,G(Zm+1)=AF(Cm1.1. 3 SMVGEAR46環(huán)境科學(xué)研究第13卷SMVGEAR( Sparse-Matrix Vectorized Gear Code25是另(△tm≤Mt≤△tmx)(9)種GEAR式的解ODEs的方法。它與 LSODE有一定的相式中:最大允許△t的改變量取值為1%~100%。似性但也有其特點(diǎn)1.3Y&B和MY&B其最大的特點(diǎn)是利用了矢量式計算機在計算內循環(huán)時(shí)1.3.1Y&B可以顯著(zhù)提高浮點(diǎn)運算效率這一現象設計計算代碼的。矢Y&B數值計算法又稱(chēng) Hybrid法l6,由 Young和Bori量化是指計算機在計算內部DO…LOOP循環(huán)時(shí),可以提速于1977年提出其目的是建立一種先預測再迭代校正的積到其最大速度的90%以上。當然若循環(huán)中包含一些特殊分方法。此法后被 McRae odman等采用并改進(jìn)。 Yamani作如遞歸、矩陣運算時(shí))矢量化也會(huì )減慢因而此法也減1o等將其應用于 CALGRID大尺度模式計算中271少了不必要的矩陣操作以提高速度Y&B法開(kāi)始時(shí)按反應物的特征時(shí)間與積分時(shí)間步長(cháng)的對于化學(xué)模式計算中的ODEs可以認為是二維的即關(guān)系將反應物分成2類(lèi):①剛性物質(zhì)(τ<△t),其在大氣中維是化學(xué)物種數,一維是網(wǎng)格數。 Jacobson等研究認為油于的壽命較短非剛性物質(zhì)(r;>△t)在大氣中的壽命較長(cháng)般情況下網(wǎng)格數會(huì )遠遠大于物種數故以網(wǎng)格數矢量化的公式如下速度比以化學(xué)物種數矢量化的速度快。計算每步的階和時(shí)a.剛性(x;≥△tn)間步長(cháng)的方法同于 LSODE。減少矩陣運算的措施是引入預測值JSPARSE法。此法先對物種進(jìn)行排序,將出現頻率較少、和其他物種關(guān)系簡(jiǎn)單的物種排在矩陣的頂部而出現頻率C:= Ci+(P(10)多、關(guān)系復雜的物種排在矩陣的底部這樣處理可以最大限修正值度地提高矩陣和零相乘的幾率洏而后 JSPARSE將所有零和矩陣相乘的項刪除。這是其相對于 LSODE法的主要差別。(11)當然有一點(diǎn)SMVGEAR運行在非矢量化的計算b.非剛性(r<△tn)機上時(shí)和 LSODE法在性能上是相差不大的預測值1.2QSsA和 MQSSAC(2r"-4tn)+2P?r"△1.2.1 QSSAQSSA數值計算法17又稱(chēng)擬穩態(tài)近似法是 Hesstwedu修正值等于1978年提出的。其主要思想是按照積分時(shí)間步長(cháng)△t△t(P"+P;r"+r;y2與壽命τ的比將計算公式分成3種r+r;-△tna.Δt/τ>10,此時(shí)化合物的壽命很短按穩態(tài)假設處(13)理即將P;和L;看成是常數1.3.2MY&BPMathur20等人指出倘若對n的物種應用穩態(tài)假(6)設則會(huì )大大提高計算效率。他們提出了將反應物分成3類(lèi)b.0.01<△t/r<10,此時(shí)化合物的壽命與所取步長(cháng)的方法即快速物質(zhì)τ;<0.2△n中速物質(zhì)0.2△tn≤r;接近5△tn和慢速物質(zhì)r2>5△tn。對于中速和慢速物質(zhì)應用公式(7)(10~13)快速物質(zhì)則應用公式(6) Mathur等同樣將自動(dòng)c,△t/x<0.01此時(shí)化合物的壽命較長(cháng)時(shí)間步長(cháng)控制機制引入到了Y&B法見(jiàn)式(9+1=CLn+△(P-LC)(8)14小參數法(SPM)公式(7是建立在P和L;是常數的基礎上的公式6廂(8)小參數法由我國韓天敏提出281適用于中等剛性的是為了增加計算速度。當△t/r>10時(shí)用公式(6算所引ODEs。計算公式為起的誤差小于公式(7)所引起誤差的4.5×10-3%而公式M= p[ef(tn+I Ch+1)+Ch+I-C] (14)(8X簡(jiǎn)單 Euler積分公式)雖引入了誤差項(△t/r;)/2但引+1=Cn+(1-a)xn+ax+1(15)起的誤差小于公式(7)所引起誤差的5.0×103%。式中1.2. 2 MQSSAp=hh+∈)Mh2等引入矯正因子P和L:P=-2,中國煤化工L;=L;。每次計算出的結果用上式校正后將P和L;代入1.5CNMHG公式(6~8)中繼續迭代RADM法22是針對RADM( Regional acid depositon mod此外 Mathur等還引入了時(shí)間步長(cháng)控制機制:el)模式提出。其對18種自由基作了穩態(tài)假設,與QSSA法類(lèi)似對時(shí)間步長(cháng)也規定了上下界。計算公式如下a.OH自由基濃度由解下式方程得到第1期劉峻峰等大氣化學(xué)模式的數值計算方法47al oh]+ b[ OH]+c=16)各不相同但最終目的都是為了提高計算速度和精確程度對于其他自由基應用穩態(tài)假設見(jiàn)公式13然而就象測不準原理中描述的時(shí)間和空間不可能同時(shí)測準C,H,O,與HO計算公式一樣剛性ODEs的數值解法也存在著(zhù)計算速度和計算精度(L:)(17)間難以同時(shí)滿(mǎn)足的矛盾。計算速度的提高是建立在犧牲計算精度的基礎上的而高精度的計算方法必然伴隨著(zhù)占有大d.其他物種量的計算機資源。在模式的實(shí)際應用過(guò)程中有時(shí)對計算速C(1-Ltn2)+P△t度要求很高而計算精度只要滿(mǎn)足一定的誤差范圍即可故L"△tn/2(18)采用一定的方法簡(jiǎn)化計算過(guò)程可以大大提高計算速度1.6 Gong and Cho筆者介紹的方法基本上可以分成3類(lèi):第1類(lèi)是GEARng and Cho法23將時(shí)間步長(cháng)定為0.5h根據物質(zhì)的類(lèi)包括 GEAR LSODE LSODES) SMVGEAR等方法其基壽命將其分為2類(lèi):①r≥104s稱(chēng)為慢物質(zhì)熜τ<104s稱(chēng)本思路來(lái)源于GEAR法;第2類(lèi)是穩態(tài)假設類(lèi),此類(lèi)方法的為快物質(zhì)。計算公式為基本思路是將物種按其壽命大間步長(cháng)之比分類(lèi),對不C=(C1C2)同類(lèi)的物種分別采用不同的計算公式,包括QSSA(MQSSA),Y&BMY&B) RADM Gong and Cho等方法第3類(lèi)是(19)d迭代法包括小參數法和 TWOSTEP法慢物質(zhì)在策(n+1)s濃度用式(20)算在各種數值計算方法中GEAR法是一種較為精確的計Cr= Ci +f(C )At(20)算方法也是研究其他快速方法的基礎。在研究模式機理階對于快物質(zhì)則用牛頓迭代公式計算出n-(n+1)s的濃度段對計箅結果的精確程度要求較高,般使用該方法這點(diǎn)增量δC可從表1看出。但是該法的計算速度較慢不適于大尺度的C2+-1+C2+f(Cm1+1,C2+-1)氣質(zhì)模式計算。 LSODE和 SMVGEAR兩者為了提高計算速(21)度都減少了計算 Jacobian矩陣的量,即只有在收斂失敗的情況下才計算。 SMVGEAR則又利用了矢量計算機的特點(diǎn)來(lái)其中C2=C21+8C2而送代計算出的C2又可以利用加快計算速度?？傮w而言第1類(lèi)的3種方法基本上可以滿(mǎn)式(22增加C1的精度足較高的計算精度,其時(shí)間步長(cháng)可根據計算過(guò)程動(dòng)態(tài)改變。C+=C+(C1)(22)穩態(tài)假設法代表了一類(lèi)高速計算方法。由于大氣體系中物1.7 TWOSTEP種的壽命相差很大是造成ODEs剛性的來(lái)源故其核心內容WOSTEP法2428是專(zhuān)為解大氣化學(xué)中的ODEs設計是對壽命很小和很大的物種采用較為簡(jiǎn)單的公式而對壽命的。其用 Gauss-Seidel迭代代替了牛頓迭代公式如下適中的物種進(jìn)行復雜計算這樣可以提高計算速度。但是此C+l= Y+ yf( tn+I Cn+I)(23)類(lèi)方法的精度依賴(lài)于時(shí)間步長(cháng)取值大小,當Δt的取值較大式中+5=時(shí)計算精度較差而Δt取值較小雖然提高了精度但計算c+2速度較低。QSSA法和Y&B法從形式上可以看到是比較相似的,尤其經(jīng)過(guò)改進(jìn)后的 MQSSA和MY&B法可近似認為是一種方法。( ong and Chou法介于 LSODE法和Y&B法之間其對壽命較長(cháng)的物種采用一般的積分方法,而對壽命較利用 Gauss-Seidel迭代式(23)可寫(xiě)成短物質(zhì)的計算方法類(lèi)似于 LSODE法。RADM法建立在C=C)Y+mxtt11、(24)RADM模式的基礎上,它雖然針對特定模式提出但也有其普遍的適用意義特別是它提出了18種自由基的計算順序TWOSTEP的迭代計算公式為眾所周知在數值計算過(guò)程中不同的計算順序會(huì )產(chǎn)生不同的結果,該算法給出的合理計算順序對后人的研究會(huì )帶來(lái)一定的幫助。第3類(lèi)方法主要是用迭代法來(lái)求解,其中(k=12(25) TWOSTEP法假定了本步所用濃度是通過(guò)前2次的濃度計式中:m—物種數。算得汁汗田工內竿剛性的體系當然式(25)中物種的計算順序對Fk的計算結果有影中國煤化工因而有必要對其進(jìn)行誤響但試驗證明采用式25廂(26)對結果的影響是不大的。差和CNMHG交方法是計算各物種的均C=F(C1C2…Ck=12,Ck…,,Cm)(26)方根(RMS)對于第k種物質(zhì)其RMS定義為2各算法間的比較RISk(27)以上所列方法雖然在基本原理、公式結構及條件參數上式中C",C"—分別代表精確值和計算值環(huán)境科學(xué)研究第13卷總步數,m=12,,M[5] Stern R M. Entwichlung und Anwendung eines ells mit ver在滿(mǎn)足了一定的誤差范圍前提下看其對相同計算機所schiedenen chemischen Mechanismen[ A ] In: Meteorologische用時(shí)間來(lái)比較速度的高低。Say19等對 LSODE Gong&Abhandlungen Vol. 8. Institut Fur Meteorologie der Freien UrCho,Y& B QSSA4種方法作了比較。其模式選用了CBMersitat Berlin[ C ] Verlag von Dietrich Reimer, Berlin, ISSNⅣ機理(包含33種化合物和81個(gè)反應),對10幾種物質(zhì)采0342-4324,1994用不同方法時(shí)的RMS與CPU耗時(shí)情況進(jìn)行了比較?？蓺w[6] Milford J B Gao D Russell A G et al. Use of sensitivity analysis納如下:①隨著(zhù)計算誤差的降低,各種方法的CPU耗時(shí)也在to compare chemical mechanisms for air-quality modeling J ]. Environ Sci Technol 1992 26: 1179-1 189成倍增加。② LSODE法的RMS隨著(zhù)CPU耗時(shí)的增加降低[ 7] Simpson D Andersson-Skold Y Jenkin M E. Updating the chemi-最快而QSSA法降低最慢。③當對CPU耗時(shí)要求較少時(shí)4cal scheme for the EMEP MSC-w oxidant model: current status種算法的RMS近似QSA法的性能較高。 ODMAN等研[R]. EMEP MSC-W Note 2/93, Norwegian Meteorological Insti究認為使用Y&B和改進(jìn)后的QSA法可以使計算精度得到tute Oslo ,Norway 1993保證的前提下大幅提高計算速度。他們同樣采用CBM-Ⅳ[8] Simpson D. Biogenic emissions in Europe2: implication for ozone機理(一些條件不同于上者),以Y&B和改進(jìn)的QSSA法分control strategie J ] Journal of Geophysical Rescarch 1995 .100別與 LSODE法作了比較結果發(fā)現在與 LSODE的誤差控22891-22906在1%之內時(shí),Y&B和QSSA法分別是 LSODE法的40與20[9] Pilling M Saunders s,Jenkin M,etal. Tropospheric Chemistry低計算過(guò)程中對QSSA法使用了質(zhì)量守恒的控制技術(shù),此Modelling[Eb/oLI.Downloadedfromhttp:/www.chem技術(shù)占時(shí)不大此外 MATHUR等與 JACOBSON等分別leeds. ac, uk/Atmospheric/MCM/ main. html. 1998對 MQSSA MY& B RADM及 SMVGEAR LSODE等方法作[10] Lurmann F W, Lloyd A C ,Atkinson R, A chemical mechanism了比較結論較為相似。for use in long-range transport/ acid deposition computer modeling結論[J]. J Geophys Res 1986 91: 10905--1093611] Strand A Hov O. A two-dimensional global study of tropospheric通過(guò)對以上各種方法的介紹和討論可以看出大氣化學(xué)ozone production[ J]. j Geophys Res 1994 99 22877-2289模式中剛性ODFs的計算問(wèn)題的解決方法是多種多樣的但[12] Andersson-Skold Y, Updationg the chemical scheme for the IVL大多都圍繞在計算精度和計算時(shí)間的取舍上建立的可分為tochemical trajectory model[ R]. IVL B-1151, IVL,Gote-2類(lèi),一類(lèi)在可以接受的誤差范圍內盡可能地提高計算速borg sweden 1995另一類(lèi)則追求在可以忍受的計算速度下提高計算精度。顯[13] Ruhnke r. Ein Verfabren zur Analyse von chemischen Reak然隨著(zhù)三維模式和大尺度模式的廣泛應用,研究前者的人占tionssystemen und zur Erstellung von Chemienodulen fur atmo-多數spharische Modelld Ph-D]. Thesis University Essen germany在各種方法中 GEAR LSODE QSSA和Y&B法一直是人們研究的熱點(diǎn)。其中GEAR和 LSODE等方法代表了一類(lèi)[14] Kuhn M Builtjes P J H Poppe D et al. Intercomparison for the精度較高的算法,是人們研究其他方法時(shí)精確度的參照法gas-phase chemistry in several chemistry and transport models而 LSODE因其在計算速度上的優(yōu)勢已被廣泛應用。QSSA[J]. Atmos Environ 1998 32 593--709和Y&B法作為快速計算方法其計箅精度與系統穩定性能15 Gear C W. Numerical Initial Value Problems in Ordinary Differ-也在不斷提高,些大型模式計算就采用了該方法。其他方ential Equation R ] Prentiedhill Englewood Cliffs , J, 197116] Young T R Boris J P. A numerical technique for solving stiff or法的建立原則也很明顯,一種是尋找新的原理,另一種是對dinary differential equations associated with the chemical kineticEAR與QSSA等已有方法進(jìn)行綜合改進(jìn)其中有許多地方blems[J]. J Phys Chem, 1977, 81: 2424值得后人借鑒。2427參考文獻[17] Hesstvedt E Hov O,Isaksen I S A. Quasi-steady-state approxi[ 1] Stockwell WR, Middleton P Chang J S et al. The second generamations in air pollution modeling: Comparison of two numericaltion regional acid deposition model chemical mechanism for regionalschemes for oxidant predictio[ J ] Int J Chem Kinet, 1978, 10air quality modeling J ]. J Geophys Res 1990 95: 16343-16367[2 Stockwell W R, Kley D. The Euro-RADM Mechanism: A Gas. 18 J Odman M T, Kumar N Russell A G. A comparison of fast chemPhase Chemical Mechanism for European Air Quality Studie[ C]eling[j]. Atmos Environ 1992Forschungszentrum Julich gmbH KFA) Julich Germany 1994中國煤化工[ 3] Gery M w Whitten G Z Killus J P et al. A photochemical kinet- [19]CN MH Gmparison of numerical methodsics mechanism for urban and regional scale computer modelinfor the integration of kinetic equations in atmospheric chemistry[J]. I Geophys Res19899412925-12956and transport models[ J ]. Atmos Environ 1995,29: 2585[ 4] Builtjes P J. The LOTOS- Long tem ozone simulation project summary report[ R ] IMW-TNO Report 92/240, TNO, Delft, The 20 Mathur R Young JO Schere K L , et al. A comparison of niNetherlands 1992merical techniques for solution of atmospheric equation[ J ].At第1期劉峻峰等大氣化學(xué)模式的數值計算方法49mos Environ 1998 32: 1535--1553torized gear code for atmospheric models[ J ]. Atmos Environ[21] Hindmarsh A C. LSODE and LSODI, Two new initial value or199428273-284dinary differential equation solver J]. ACM Newsl 1980, 15: 10 [26] Hindmarsh A C. ODEPACK: a systematized collection of ode22] Chang I S Brosr R A Isaksen I S A, et al. A three-dimensionalNorth-Holland Amsterdam 1983 55-64acid deposition model: physical concepts anion[J]. [27 Yamartino RJ Scire J S Carmichael GR et al. The CALGRIDGeophys Res98792:14681-14700rid model: I. Model formulation[J[ 23] Gong W Cho HR. A numerical scheme for the integration of the中國煤化工1512gas-phase chemical rate equations in three-dimensional atmo- [28 1CNMHG境化學(xué)M]北京高等教育spheric model![ J ] Atmos Environ 1993 27A 2147--2160出版社1990.156~158[24] VerwerJ(. Gauss-Seidel iteration for stiff ODEs from chemical[29]關(guān)治陸金甫.數值分析基砒[M].北京高等教育出版社kinetic J ]. SIAM J Sci Comput 1994, 15: 1243-12501998.445-45125] Jacobson M Z, Turco R P. SMVGEAR: A sparse-matrix ve