環(huán)境科學(xué)研究第13卷第1期esearch of environmental sciencesVol.13No.12000大氣化學(xué)模式的數(shù)值計算方法劉峻峰,李金龍,白郁華北京大學(xué)技術(shù)物理系北京100871)摘要∶大氣化學(xué)模式計算的一個重要研究熱點是如何快速、精確地對相互關(guān)聯(lián)的非線性常微分方程組(OEs進(jìn)行數(shù)值積分求解。由于大氣中督物種壽命長短相差很大,導(dǎo)致ωDEs具有了很強(qiáng)的剛性因而在數(shù)值求算過程中存在著計算精度和計算效率的平衡問題。筆者介紹了當(dāng)前幾種常見的大氣化學(xué)模式的計算方法包括 GEAR QSSA( MQSSA)Y&BMY&B),RADM, LSODE SMVGEAR小參數(shù)法 Gong and Cho及 TWOSTEP等并結(jié)合當(dāng)前本領(lǐng)域內(nèi)的部分研究成果對一些重要的方法進(jìn)行了比較和分析。關(guān)鍵詞:大氣化學(xué)模式;剛性常微分方程組;數(shù)值積分中圖分類號:X11文獻(xiàn)標(biāo)識碼: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年代以來,模式研究已成為大氣環(huán)境化學(xué)領(lǐng)域中或(1)代表了一組數(shù)量較大、非線性且互相關(guān)聯(lián)的常微不可缺少的組成部分。一—個性能優(yōu)異的大氣化學(xué)模式應(yīng)具分方程細(xì)ODEs)故無法求出其解析解,只能在給定每種化備2個必要條件⑩合理的模式機(jī)理。此機(jī)理能較準(zhǔn)確地反合物初值后用各物種濃度對時間數(shù)值積分的方法,求岀各映實際大氣中物種數(shù)量及彼此間的化學(xué)反應(yīng)同時給出較準(zhǔn)物種在特定時刻的濃度。確的動力學(xué)數(shù)據(jù)。②高效精確的計算方法。該方法在滿足最早的方法是對具有初值的式(1)在固定時間步長的條定誤差范圍的前提下應(yīng)盡可能快速地計算出最終結(jié)果件下進(jìn)行 Euler積分。但由于各反應(yīng)物種在大氣中的壽命長根據(jù)研究的對象及目的不同,人們設(shè)計出了不同的機(jī)短不一,物種間的反應(yīng)速率相差高達(dá)幾個數(shù)量級,造成了理。常見的幾種大氣化學(xué)模式機(jī)理及計算方法見表1。ODEs的剛性。為了保證求解方程組的準(zhǔn)確性和穩(wěn)定性由表1可以看出,般機(jī)理包含的化合物種類由幾十種須將時間步長定得很小造成計算效率極低。因此在數(shù)值積到上百種反應(yīng)方程式由幾百個到上千個。由這些化合物及分過程中存在著積分時間步長同各物種的壽命是否相匹配反應(yīng)方程式根據(jù)質(zhì)量作用定律可以寫出每種化合物的反的問題。若積分時間步長取的較大,對于壽命短的化合物會應(yīng)速率方程產(chǎn)生較大誤差時間步長取得過小,又會使計算過慢造成dh=F(C)=P(Ct)-L(CC(1)對計算資源的浪費(fèi)中國煤化工中新數(shù)值求解OEs的方法式中:P(Cd)第i種化合物的產(chǎn)生速率HCNMHG應(yīng)用,但是周看人們研究物L(fēng)(C,)第i種化合物的去除速率。質(zhì)量模式中所取的網(wǎng)絡(luò)數(shù)成干上萬,一般又要模擬較長收稿日期:1999-09-14段時間的反應(yīng)因而對計算速度的要求也在提高,即需要作者簡介劉峻峰(1974-)男北京市人碩士研究生種快速且較為精確的計算方法。GEAR法逐漸成為人們對新研究的快速方法精確度的檢驗手段。1977年Y&B法161第1期劉峻峰等大氣化學(xué)模式的數(shù)值計算方法45和1978年QSsA"擬穩(wěn)態(tài)假設(shè))法的提出標(biāo)志著快速計算較18-201,并提出了一批新的或是改進(jìn)的方法,如法的誕生。此后直到現(xiàn)在又有許多科學(xué)家對大氣化學(xué)模式 LSODEL211,RADM21, Gong andTWOSTEP 241中剛性O(shè)DEs數(shù)值解法的效率和精度進(jìn)行了研究和比 SMVGEAR25] MQSSA和MY&B201表1幾種常見大氣化學(xué)模式及其計算方法Table 1 Atmospheric chemistry models and their numerical solvers機(jī)理縮寫名參考物種數(shù)反應(yīng)方程數(shù)計算法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資料主要來源 Kuhn et al.41998筆者在前面簡要介紹了GEAR, LSODE, SMVGEARhF(Cm)Cm-1為前一步的C值,CQSSA, MQSSA,Y&B,MY&B,小參數(shù)法, RADM Gong and為本步的C值Cho及 TWOSTEP等計算方法。在下面結(jié)合各方面的研究對這些方法進(jìn)行了分類并對重要方法的精確性和計算效率n——計算步數(shù)做了比較和分析每步校正次數(shù)1方法介紹積分開始時,Z0=(CC′o),C0為初值,C'o=dC/d1.1GEAR法和改進(jìn)GEAR法(C)第一步步長取值較小P取1以后根據(jù)計算的誤GEAR法及與GEAR法類似的 LSODE, LSODES,差大小可按一定的方法自動調(diào)整h和P。SMVGEAR等代表了一類較精確的數(shù)值計算方1.1.2 LSODE與 LSODES1.1.1GEAR法LSODE( Livermore solver for ordinary differential equaGEAR數(shù)值計算法151具有適用范圍廣、可以自動調(diào)節(jié)ons521是建立在GEAR法基礎(chǔ)上,被廣為應(yīng)用計算ODEs時間步長與誤差,以及計算結(jié)果準(zhǔn)確度高等特點,是一種的較精確方法。此方法允許使用者根據(jù)ODEs的剛性程度為理想的計算方法。但其最大的缺點是計算步驟繁瑣要頻選擇Adam多步計算公式(對于非剛性體系)或BDF公式繁計算 Jacobian行列式,當(dāng)計算的物種比較多時會耗費(fèi)大 backward difference multistep formulas對于剛性體系)大量的機(jī)時和內(nèi)存。GEAR法的數(shù)值計算公式如下氣化學(xué)模式的ODEs體系剛性較大故多用BDF法。BD)F的a.預(yù)測公式計算公式為(2)C)+△tbf(Cb.校正公式式中:C"+1—當(dāng)前t"1的濃度Zm+1=Zm+1-[(l1-lnh)(xZm+1)](3)中國煤化工C.最終公式(4)CNMHG邇以刊酊炬陣結(jié)構(gòu)稀少的不同來源(是式中zn—(k)×(P+1矩陣來自于用戶直接輸入,還是來自于調(diào)用 Jacobian子程序或是A-(P+1)×(P+1)矩陣l≤P≤6多次調(diào)用方程評估子程序)分別采用不同方法簡化計算步Ⅰ——單位矩陣驟提高計算效率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%。似性但也有其特點1.3Y&B和MY&B其最大的特點是利用了矢量式計算機(jī)在計算內(nèi)循環(huán)時1.3.1Y&B可以顯著提高浮點運(yùn)算效率這一現(xiàn)象設(shè)計計算代碼的。矢Y&B數(shù)值計算法又稱 Hybrid法l6,由 Young和Bori量化是指計算機(jī)在計算內(nèi)部DO…LOOP循環(huán)時,可以提速于1977年提出其目的是建立一種先預(yù)測再迭代校正的積到其最大速度的90%以上。當(dāng)然若循環(huán)中包含一些特殊分方法。此法后被 McRae odman等采用并改進(jìn)。 Yamani作如遞歸、矩陣運(yùn)算時)矢量化也會減慢因而此法也減1o等將其應(yīng)用于 CALGRID大尺度模式計算中271少了不必要的矩陣操作以提高速度Y&B法開始時按反應(yīng)物的特征時間與積分時間步長的對于化學(xué)模式計算中的ODEs可以認(rèn)為是二維的即關(guān)系將反應(yīng)物分成2類:①剛性物質(zhì)(τ<△t),其在大氣中維是化學(xué)物種數(shù),一維是網(wǎng)格數(shù)。 Jacobson等研究認(rèn)為油于的壽命較短非剛性物質(zhì)(r;>△t)在大氣中的壽命較長般情況下網(wǎng)格數(shù)會遠(yuǎn)遠(yuǎn)大于物種數(shù)故以網(wǎng)格數(shù)矢量化的公式如下速度比以化學(xué)物種數(shù)矢量化的速度快。計算每步的階和時a.剛性(x;≥△tn)間步長的方法同于 LSODE。減少矩陣運(yùn)算的措施是引入預(yù)測值JSPARSE法。此法先對物種進(jìn)行排序,將出現(xiàn)頻率較少、和其他物種關(guān)系簡單的物種排在矩陣的頂部而出現(xiàn)頻率C:= Ci+(P(10)多、關(guān)系復(fù)雜的物種排在矩陣的底部這樣處理可以最大限修正值度地提高矩陣和零相乘的幾率洏而后 JSPARSE將所有零和矩陣相乘的項刪除。這是其相對于 LSODE法的主要差別。(11)當(dāng)然有一點SMVGEAR運(yùn)行在非矢量化的計算b.非剛性(r<△tn)機(jī)上時和 LSODE法在性能上是相差不大的預(yù)測值1.2QSsA和 MQSSAC(2r"-4tn)+2P?r"△1.2.1 QSSAQSSA數(shù)值計算法17又稱擬穩(wěn)態(tài)近似法是 Hesstwedu修正值等于1978年提出的。其主要思想是按照積分時間步長△t△t(P"+P;r"+r;y2與壽命τ的比將計算公式分成3種r+r;-△tna.Δt/τ>10,此時化合物的壽命很短按穩(wěn)態(tài)假設(shè)處(13)理即將P;和L;看成是常數(shù)1.3.2MY&BPMathur20等人指出倘若對n的物種應(yīng)用穩(wěn)態(tài)假(6)設(shè)則會大大提高計算效率。他們提出了將反應(yīng)物分成3類b.0.01<△t/r<10,此時化合物的壽命與所取步長的方法即快速物質(zhì)τ;<0.2△n中速物質(zhì)0.2△tn≤r;接近5△tn和慢速物質(zhì)r2>5△tn。對于中速和慢速物質(zhì)應(yīng)用公式(7)(10~13)快速物質(zhì)則應(yīng)用公式(6) Mathur等同樣將自動c,△t/x<0.01此時化合物的壽命較長時間步長控制機(jī)制引入到了Y&B法見式(9+1=CLn+△(P-LC)(8)14小參數(shù)法(SPM)公式(7是建立在P和L;是常數(shù)的基礎(chǔ)上的公式6廂(8)小參數(shù)法由我國韓天敏提出281適用于中等剛性的是為了增加計算速度。當(dāng)△t/r>10時用公式(6算所引ODEs。計算公式為起的誤差小于公式(7)所引起誤差的4.5×10-3%而公式M= p[ef(tn+I Ch+1)+Ch+I-C] (14)(8X簡單 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;。每次計算出的結(jié)果用上式校正后將P和L;代入1.5CNMHG公式(6~8)中繼續(xù)迭代RADM法22是針對RADM( Regional acid depositon mod此外 Mathur等還引入了時間步長控制機(jī)制:el)模式提出。其對18種自由基作了穩(wěn)態(tài)假設(shè),與QSSA法類似對時間步長也規(guī)定了上下界。計算公式如下a.OH自由基濃度由解下式方程得到第1期劉峻峰等大氣化學(xué)模式的數(shù)值計算方法47al oh]+ b[ OH]+c=16)各不相同但最終目的都是為了提高計算速度和精確程度對于其他自由基應(yīng)用穩(wěn)態(tài)假設(shè)見公式13然而就象測不準(zhǔn)原理中描述的時間和空間不可能同時測準(zhǔn)C,H,O,與HO計算公式一樣剛性O(shè)DEs的數(shù)值解法也存在著計算速度和計算精度(L:)(17)間難以同時滿足的矛盾。計算速度的提高是建立在犧牲計算精度的基礎(chǔ)上的而高精度的計算方法必然伴隨著占有大d.其他物種量的計算機(jī)資源。在模式的實際應(yīng)用過程中有時對計算速C(1-Ltn2)+P△t度要求很高而計算精度只要滿足一定的誤差范圍即可故L"△tn/2(18)采用一定的方法簡化計算過程可以大大提高計算速度1.6 Gong and Cho筆者介紹的方法基本上可以分成3類:第1類是GEARng and Cho法23將時間步長定為0.5h根據(jù)物質(zhì)的類包括 GEAR LSODE LSODES) SMVGEAR等方法其基壽命將其分為2類:①r≥104s稱為慢物質(zhì)熜τ<104s稱本思路來源于GEAR法;第2類是穩(wěn)態(tài)假設(shè)類,此類方法的為快物質(zhì)。計算公式為基本思路是將物種按其壽命大間步長之比分類,對不C=(C1C2)同類的物種分別采用不同的計算公式,包括QSSA(MQSSA),Y&BMY&B) RADM Gong and Cho等方法第3類是(19)d迭代法包括小參數(shù)法和 TWOSTEP法慢物質(zhì)在策(n+1)s濃度用式(20)算在各種數(shù)值計算方法中GEAR法是一種較為精確的計Cr= Ci +f(C )At(20)算方法也是研究其他快速方法的基礎(chǔ)。在研究模式機(jī)理階對于快物質(zhì)則用牛頓迭代公式計算出n-(n+1)s的濃度段對計箅結(jié)果的精確程度要求較高,般使用該方法這點增量δC可從表1看出。但是該法的計算速度較慢不適于大尺度的C2+-1+C2+f(Cm1+1,C2+-1)氣質(zhì)模式計算。 LSODE和 SMVGEAR兩者為了提高計算速(21)度都減少了計算 Jacobian矩陣的量,即只有在收斂失敗的情況下才計算。 SMVGEAR則又利用了矢量計算機(jī)的特點來其中C2=C21+8C2而送代計算出的C2又可以利用加快計算速度?？傮w而言第1類的3種方法基本上可以滿式(22增加C1的精度足較高的計算精度,其時間步長可根據(jù)計算過程動態(tài)改變。C+=C+(C1)(22)穩(wěn)態(tài)假設(shè)法代表了一類高速計算方法。由于大氣體系中物1.7 TWOSTEP種的壽命相差很大是造成ODEs剛性的來源故其核心內(nèi)容WOSTEP法2428是專為解大氣化學(xué)中的ODEs設(shè)計是對壽命很小和很大的物種采用較為簡單的公式而對壽命的。其用 Gauss-Seidel迭代代替了牛頓迭代公式如下適中的物種進(jìn)行復(fù)雜計算這樣可以提高計算速度。但是此C+l= Y+ yf( tn+I Cn+I)(23)類方法的精度依賴于時間步長取值大小,當(dāng)Δt的取值較大式中+5=時計算精度較差而Δt取值較小雖然提高了精度但計算c+2速度較低。QSSA法和Y&B法從形式上可以看到是比較相似的,尤其經(jīng)過改進(jìn)后的 MQSSA和MY&B法可近似認(rèn)為是一種方法。( ong and Chou法介于 LSODE法和Y&B法之間其對壽命較長的物種采用一般的積分方法,而對壽命較利用 Gauss-Seidel迭代式(23)可寫成短物質(zhì)的計算方法類似于 LSODE法。RADM法建立在C=C)Y+mxtt11、(24)RADM模式的基礎(chǔ)上,它雖然針對特定模式提出但也有其普遍的適用意義特別是它提出了18種自由基的計算順序TWOSTEP的迭代計算公式為眾所周知在數(shù)值計算過程中不同的計算順序會產(chǎn)生不同的結(jié)果,該算法給出的合理計算順序?qū)笕说难芯繒硪欢ǖ膸椭?。?類方法主要是用迭代法來求解,其中(k=12(25) TWOSTEP法假定了本步所用濃度是通過前2次的濃度計式中:m—物種數(shù)。算得汁汗田工內(nèi)竿剛性的體系當(dāng)然式(25)中物種的計算順序?qū)k的計算結(jié)果有影中國煤化工因而有必要對其進(jìn)行誤響但試驗證明采用式25廂(26)對結(jié)果的影響是不大的。差和CNMHG交方法是計算各物種的均C=F(C1C2…Ck=12,Ck…,,Cm)(26)方根(RMS)對于第k種物質(zhì)其RMS定義為2各算法間的比較RISk(27)以上所列方法雖然在基本原理、公式結(jié)構(gòu)及條件參數(shù)上式中C",C"—分別代表精確值和計算值環(huán)境科學(xué)研究第13卷總步數(shù),m=12,,M[5] Stern R M. Entwichlung und Anwendung eines ells mit ver在滿足了一定的誤差范圍前提下看其對相同計算機(jī)所schiedenen chemischen Mechanismen[ A ] In: Meteorologische用時間來比較速度的高低。Say19等對 LSODE Gong&Abhandlungen Vol. 8. Institut Fur Meteorologie der Freien UrCho,Y& B QSSA4種方法作了比較。其模式選用了CBMersitat Berlin[ C ] Verlag von Dietrich Reimer, Berlin, ISSNⅣ機(jī)理(包含33種化合物和81個反應(yīng)),對10幾種物質(zhì)采0342-4324,1994用不同方法時的RMS與CPU耗時情況進(jìn)行了比較?？蓺w[6] Milford J B Gao D Russell A G et al. Use of sensitivity analysis納如下:①隨著計算誤差的降低,各種方法的CPU耗時也在to compare chemical mechanisms for air-quality modeling J ]. Environ Sci Technol 1992 26: 1179-1 189成倍增加。② LSODE法的RMS隨著CPU耗時的增加降低[ 7] Simpson D Andersson-Skold Y Jenkin M E. Updating the chemi-最快而QSSA法降低最慢。③當(dāng)對CPU耗時要求較少時4cal scheme for the EMEP MSC-w oxidant model: current status種算法的RMS近似QSA法的性能較高。 ODMAN等研[R]. EMEP MSC-W Note 2/93, Norwegian Meteorological Insti究認(rèn)為使用Y&B和改進(jìn)后的QSA法可以使計算精度得到tute Oslo ,Norway 1993保證的前提下大幅提高計算速度。他們同樣采用CBM-Ⅳ[8] Simpson D. Biogenic emissions in Europe2: implication for ozone機(jī)理(一些條件不同于上者),以Y&B和改進(jìn)的QSSA法分control strategie J ] Journal of Geophysical Rescarch 1995 .100別與 LSODE法作了比較結(jié)果發(fā)現(xiàn)在與 LSODE的誤差控22891-22906在1%之內(nèi)時,Y&B和QSSA法分別是 LSODE法的40與20[9] Pilling M Saunders s,Jenkin M,etal. Tropospheric Chemistry低計算過程中對QSSA法使用了質(zhì)量守恒的控制技術(shù),此Modelling[Eb/oLI.Downloadedfromhttp:/www.chem技術(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了比較結(jié)論較為相似。for use in long-range transport/ acid deposition computer modeling結(jié)論[J]. J Geophys Res 1986 91: 10905--1093611] Strand A Hov O. A two-dimensional global study of tropospheric通過對以上各種方法的介紹和討論可以看出大氣化學(xué)ozone production[ J]. j Geophys Res 1994 99 22877-2289模式中剛性O(shè)DFs的計算問題的解決方法是多種多樣的但[12] Andersson-Skold Y, Updationg the chemical scheme for the IVL大多都圍繞在計算精度和計算時間的取舍上建立的可分為tochemical trajectory model[ R]. IVL B-1151, IVL,Gote-2類,一類在可以接受的誤差范圍內(nèi)盡可能地提高計算速borg sweden 1995另一類則追求在可以忍受的計算速度下提高計算精度。顯[13] Ruhnke r. Ein Verfabren zur Analyse von chemischen Reak然隨著三維模式和大尺度模式的廣泛應(yīng)用,研究前者的人占tionssystemen und zur Erstellung von Chemienodulen fur atmo-多數(shù)spharische Modelld Ph-D]. Thesis University Essen germany在各種方法中 GEAR LSODE QSSA和Y&B法一直是人們研究的熱點。其中GEAR和 LSODE等方法代表了一類[14] Kuhn M Builtjes P J H Poppe D et al. Intercomparison for the精度較高的算法,是人們研究其他方法時精確度的參照法gas-phase chemistry in several chemistry and transport models而 LSODE因其在計算速度上的優(yōu)勢已被廣泛應(yīng)用。QSSA[J]. Atmos Environ 1998 32 593--709和Y&B法作為快速計算方法其計箅精度與系統(tǒng)穩(wěn)定性能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參考文獻(xiàn)[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é)模式的數(shù)值計算方法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)治陸金甫.數(shù)值分析基砒[M].北京高等教育出版社kinetic J ]. SIAM J Sci Comput 1994, 15: 1243-12501998.445-45125] Jacobson M Z, Turco R P. SMVGEAR: A sparse-matrix ve