楊 飛，楊文俊，楊 森

（1.長(zhǎng)江科學(xué)院水利部江湖治理與防洪重點(diǎn)實(shí)驗(yàn)室，武漢 430010；2.重慶市武隆縣水務(wù)局，重慶 408500）

ELCIRC源程序代碼分析

楊 飛1，楊文俊1，楊 森2

（1.長(zhǎng)江科學(xué)院水利部江湖治理與防洪重點(diǎn)實(shí)驗(yàn)室，武漢 430010；2.重慶市武隆縣水務(wù)局，重慶 408500）

開放源程序代碼ELCIRC是采用基于水平向無結(jié)構(gòu)網(wǎng)格、垂向z坐標(biāo)體系和半隱格式的歐拉－拉格朗日有限體積／有限差分方法來解淺水方程。研究了ELCIRC中除控制方程組求解以外的部分，分析了源程序在插值計(jì)算、拓?fù)浣Y(jié)構(gòu)、定解條件、分層信息和干濕法等技術(shù)上的具體操作，尤其是對(duì)歐拉拉格朗日模型特有的逆向追蹤算法做了詳細(xì)闡述。同時(shí)討論了一些細(xì)節(jié)問題，在不恰當(dāng)?shù)牡胤浇o出一些參考和建議，避免小數(shù)做分母、相近數(shù)相減引起的較大誤差。然而，這些經(jīng)過實(shí)際應(yīng)用的算法原理，可以為我們開發(fā)水動(dòng)力學(xué)模型提供有價(jià)值的參考。

ELCIRC；逆向追蹤；源代碼；三維水動(dòng)力學(xué)模型；算法

1 研究背景

ELCIRC［1］模型（Eulerian-Lagrangian CIRCulation）是美國(guó)俄勒岡健康與科學(xué)大學(xué)海岸陸地邊緣陸緣研究中心（Oregon Health＆Science University OHSU’s Center for Coastal and Land-Margin Research）基于對(duì)哥倫比亞河的研究，開發(fā)的哥倫比亞環(huán)流模型CORIE的一部分，可以單獨(dú)用于模擬計(jì)算。模型已進(jìn)行基準(zhǔn)測(cè)試，并應(yīng)用到哥倫比亞河的物理過程的模擬。國(guó)內(nèi)外應(yīng)用逐漸增加，Myers和Aikman2003年研究St.Johns河流域洪泛區(qū)淹沒過程［2］。楊金艷等2006年利用ELCIRC計(jì)算了長(zhǎng)江口的潮流場(chǎng)［3］，該模型的中文介紹可以參考楊金艷［3］碩士畢業(yè)論文。Gong等2008年將ELCIRC應(yīng)用到海南新村瀉湖研究中［4］，Gong等2009年使用ELCIRC檢驗(yàn)了Chesapeake海灣水位對(duì)選定的3組東北風(fēng)響應(yīng)及正壓模式下潮下帶海灣河口水交換［5］。Dias等2009年利用水動(dòng)力ELCIRC和輸移VELA，VELApart模型在葡萄牙福爾摩挲河口瀉湖的應(yīng)用，研究入口位置改變的影響［6］。Picado等2010年以葡萄牙阿威羅瀉湖為例研究了局部地形調(diào)整對(duì)潮汐可能產(chǎn)生的影響，采用ELCIRC模型研究在過水面積增加、瀉湖的水動(dòng)力條件減弱下可能的潮汐變化［7］。

ELCIRC采用了基于水平向無結(jié)構(gòu)網(wǎng)格、垂向z坐標(biāo)體系下，半隱格式的歐拉－拉格朗日有限體積／有限差分方法來解淺水方程。算法上考慮了多種紊流閉合模式，也包括了潮汐勢(shì)，大氣壓梯度項(xiàng)及水氣交換，可以模擬多種物理過程［1］。圖1為以3D斜壓哥倫比亞河環(huán)流模擬為例的ELCIRC求解過程結(jié)構(gòu)圖，網(wǎng)格水平節(jié)點(diǎn)33 634個(gè)，水平單元50 389個(gè)，垂向分層62個(gè)，對(duì)應(yīng)總的活動(dòng)柱面約為2.2×1012個(gè)［1］。本文主要涉及的是動(dòng)量方程組求解以外的參數(shù)處理等，文中大部分程序語言來自文獻(xiàn)［8］中的ELCIRC5＿4c。

2 拓?fù)潢P(guān)系

程序首先使用平面的單元、邊、頂點(diǎn)來表示實(shí)際區(qū)域，然后進(jìn)行垂向劃分。實(shí)際的計(jì)算三維網(wǎng)格為柱體，實(shí)際中用到的變量定義位置如圖2。字母t和s（或tem和sal）表示溫度和鹽度，q2和xl表示紊流模型中的紊動(dòng)能和混合長(zhǎng)，字母p（或nd）表示節(jié)點(diǎn)處，字母e（或el）和c表示單元中心，字母s（或sd）表示側(cè)面中心，字母z表示垂向分層，0表示初始，diff表示擴(kuò)散，tmp表示過程量，half表示相鄰垂向分層平均，tot表示垂向累積，rho表示密度，bt表示逆向追蹤。對(duì)于流速，下標(biāo)n和t表示法向和切向，下標(biāo)u，v和w表示節(jié)點(diǎn)處流速，后面的數(shù)字1和2對(duì)應(yīng)迭代計(jì)算過程不同時(shí)段。圖中d表示厚度，但程序中多數(shù)變量會(huì)將d作為首字母，表示雙精度數(shù)值，將n作為首字母表示整數(shù)值，如單元、節(jié)點(diǎn)和邊的序號(hào)。

圖1 ELCIRC計(jì)算流程圖［1］Fig.1 Flow chart of ELCIRC［1］

圖2 程序常用變量定義位置Fig.2 Position of program variables in grid

2.1 平面位置關(guān)系數(shù)組

nx（4，4，3），考慮三角形和四邊形2種單元，見式

（1）。數(shù)組實(shí)現(xiàn)循環(huán)三角形1－2－3循環(huán)和四邊形的

1－2－3－4循環(huán)。注意，在平面網(wǎng)格生成文件hgrid.in中，單元對(duì)應(yīng)的節(jié)點(diǎn)就是按照空間位置逆時(shí)針的順序表述，數(shù)組可以循環(huán)尋址，起點(diǎn)任意。實(shí)際中調(diào)用非零元素來得出邊與節(jié)點(diǎn)在單元中的相對(duì)位置關(guān)系。

2.2 邊正向判定

離散網(wǎng)格的每條邊其實(shí)是矢量，有（正）方向，端點(diǎn)定義為起止節(jié)點(diǎn)，邊的正向和單元編號(hào)關(guān)系緊密。一般情況下單元編號(hào)是按照空間位置由下到上、由左到右順序編排的，這給邊的定義帶來極大的方便。對(duì)于內(nèi)部邊，對(duì)應(yīng)2個(gè)相鄰單元，邊的法線方向?yàn)橛删幪?hào)低的垂直邊指向編號(hào)高的；對(duì)于組成區(qū)域邊界的邊，其只對(duì)應(yīng)一個(gè)單元，即邊界單元，邊的法線方向?yàn)橛蓡卧獌?nèi)部垂直邊指向外。邊的切向均為法線方向逆時(shí)針旋轉(zhuǎn)90度。假定邊為j，存儲(chǔ)時(shí)采用x軸正向到邊夾角的正弦值snx（j），余弦值sny（j），邊的法向單位向量表示為（snx（j），sny（j）），切向單位向量為（－sny（j），snx（j））。

在計(jì)算單元i的水位時(shí)，其中j邊的法向正負(fù)符號(hào)Ssign（i，j）表達(dá)式具體為

3 逆向追蹤

逆向追蹤從邊的中點(diǎn)（側(cè)面中心x0，y0，z0）出發(fā)，根據(jù)該點(diǎn)t時(shí)刻的速度，顯式倒推判斷Δt以前該質(zhì)點(diǎn)位置（xt，yt，zt）。網(wǎng)格無水和流速小于10－6m／s，不需要逆向追蹤計(jì)算，垂向流速為側(cè)面四頂點(diǎn)垂向流速平均值。由于方程建立在邊上，每條邊都要考慮逆向追蹤方程，上面是依據(jù)邊來搜尋逆向追蹤的需要網(wǎng)格，邊界以外的邊都對(duì)應(yīng)2個(gè)單元，一般優(yōu)先考慮編號(hào)小的單元。

3.1 定位（quick-search）

（x0，y0）在nel單元內(nèi)，倒推得到的（xt，yt）也在nel內(nèi)時(shí)，確定zt所在層jlev即可，但不能溢出邊界；（xt，yt）不在nel內(nèi)時(shí)，找出nel與（x0，y0）、（xt，yt）連線（這時(shí)（x0，y0）可以先微調(diào)內(nèi)移，使其不在邊上）相交的一邊及交點(diǎn)（用到下面的intersect2函數(shù)），求出交點(diǎn)的高程，并用交點(diǎn)代替（x0，y0）作為起點(diǎn)，共用該邊的單元作為新的nel單元，繼續(xù)迭代直到（x0，y0）、（xt，yt）在同一單元內(nèi)。（x0，y0）的速度值是取所屬邊速度（兩端點(diǎn)速度均值，求解動(dòng)量方程不考慮時(shí)間因素，流場(chǎng)過程知道時(shí)求解物質(zhì)輸移方程時(shí)要考慮速度隨時(shí)間變化），這樣簡(jiǎn)化默認(rèn)相鄰節(jié)點(diǎn)間的速度差異在誤差允許范圍內(nèi)。筆者認(rèn)為可以根據(jù)端點(diǎn)線性插值交點(diǎn)速度，保證一階精度。每次確定新起點(diǎn)后，（xt，yt）根據(jù)剩余時(shí)間和起點(diǎn)速度重新確立。整個(gè)過程如圖3所示，該圖時(shí)間步長(zhǎng)dtb不同于文獻(xiàn)1中圖4所示整個(gè)時(shí)間步長(zhǎng)dt。從t到t－dtb的迭代次數(shù)不能預(yù)先知道，因此要限定以避免進(jìn)入死循環(huán)。

圖3 一次逆向追蹤迭代示意圖Fig.3 Back-track iteration for one sub step

整個(gè)過程中并沒有考慮z向速度的逆向追蹤，這與動(dòng)量積分方程建立在水平法向和切向是一致的。逆向追蹤遇到邊界時(shí)，只取切向速度。時(shí)間步長(zhǎng)Δt劃分越小，誤差就越小，因此把原來的時(shí)間步長(zhǎng)在計(jì)算有效性（整個(gè)程序的計(jì)算時(shí)間長(zhǎng)短和計(jì)算量的大?。┣疤嵯逻m當(dāng)細(xì)分ndelt個(gè)子步長(zhǎng)dtb，進(jìn)行迭代ndelt次逆向追蹤迭代（dtb＝Δt／ndelt），每次迭代得到的（xt，yt）作為下次的起點(diǎn)。

3.2 交叉點(diǎn)（intersect2）

判斷端點(diǎn)為（x1，y1）、（x2，y2）和（x3，y3）、（x4，y4）兩線段有無交點(diǎn)，若有交點(diǎn)，則存在0－1之間的2個(gè)數(shù)tt1，tt2，滿足：

判斷tt1，tt2是否在0和1之間，即可確定有無交點(diǎn)，若有，則可求出交點(diǎn)。

3.3 逆向追蹤（back-track）

利用quick-search函數(shù)找到（x0，y0）對(duì)應(yīng)的（xt，yt）以及nel后，確定（xt，yt）所在的平面，垂向插值nel單元在該平面頂點(diǎn)值，然后根據(jù)在單元的相對(duì)位置，平面插值（xt，yt）。

單元為三角形時(shí)，（xt，yt）與頂點(diǎn)連線將其分為3個(gè)小三角形，（xt，yt）的值為3頂點(diǎn)值的加權(quán)平均，權(quán)重為頂點(diǎn)所對(duì)應(yīng)的小三角形與單元面積之比。單元為四邊形時(shí)，需要雙線性投影變換，中心對(duì)應(yīng)新的坐標(biāo)原點(diǎn)，原來的點(diǎn)1，2，3，4或點(diǎn)A，B，C，D映射坐標(biāo)為（－1，－1），（1，－1），（1，1），（－1，1），如圖4。分別討論平行四邊形、梯形、不規(guī)則四邊形，任意一點(diǎn)P的投影變換后的坐標(biāo)（ξ，η）。

圖4 四邊形單元雙線性投影Fig.4 Bilinear projection of quadrilateral unit

平行四邊形時(shí)，如圖5（a）所示。

式中：X，E分別是2條邊的中點(diǎn)，O為平行四邊形的中心，（xt，yt）即為P點(diǎn)。

梯形時(shí)，以其中一組對(duì)邊AB和DC平行為例，另一種情況一樣。對(duì)邊中點(diǎn)連線將梯形分為4個(gè)區(qū)域。如圖5（b），M，N為對(duì)角線BD和AC的中點(diǎn)，E，F(xiàn)，G，H為AB，BC，CD，DA的中點(diǎn)，O為HF和EG的交點(diǎn)。SP平行底邊，延長(zhǎng)線交BC于I，ξ大小為SP／SI，η大小為OS／OG。

圖5 圖形變換Fig.5 Graphic transformation

一般對(duì)不規(guī)則四邊形，如圖5（c）所示，其算法為

式中：α＝2 H→F×2M→N；β＝2 H→F×2→EG－4 O→P×2 M→N；γ＝－4 O→P×2→EG。求解一元二次方程式（9）得到ξ的2個(gè)根，選擇絕對(duì)值小于1的根，代入式（10）中求解η。

3.4 算法誤差

程序原有的計(jì)算不規(guī)則四邊形（圖5（c））使用到向量M→N和→AD×C→B（記為d）。數(shù)值計(jì)算要求數(shù)

eta值穩(wěn)定，在設(shè)計(jì)算法時(shí)還應(yīng)盡量避免誤差危害，防止有效數(shù)字損失，要避免兩相近數(shù)相減和用絕對(duì)值很小的數(shù)做除數(shù)，還要注意運(yùn)算次序和減少運(yùn)算次數(shù)。在程序判斷語句中使用了M→N的分量數(shù)值和d是否

eta大于小量small3＝10－5作為判斷依據(jù)，對(duì)于動(dòng)量方程求解器計(jì)算精度為7位有效數(shù)字的程序來說，繼續(xù)將其作為計(jì)算乘數(shù)和除數(shù)，筆者認(rèn)為會(huì)帶來不可忽略的誤差?？紤]到在網(wǎng)格剖分基本優(yōu)化為正交網(wǎng)格，將會(huì)有更多的單元格處在臨界條件附近，盡管采用雙倍字長(zhǎng)計(jì)算，相近數(shù)相減仍會(huì)帶來有效位數(shù)的減少。因此需要對(duì)算法進(jìn)行修改。

先利用式（7）求出近似值ξ0，令ξ＝ξ0＋Δξ，Δξ是一個(gè)未知校正量，然后將ξ0＋Δξ代入式（9）得到式（11），參照式（9）表示成式（12）。

這樣由（13）或迭代計(jì)算式（14）求解會(huì)使得計(jì)算的精度有所改善。對(duì)求得的小量Δξ加以判斷，Δξ太大（如大于0.1）則不采用。

2

2

2

程序中包含判斷ξ＞1.1，η＞2.0是否成立的語句，顯示出了求解結(jié)果的誤差，當(dāng)然，由于速度梯度不大，這種求解誤差給結(jié)果帶來的影響不是很明顯，但是較大的計(jì)算量和不可忽略的誤差確實(shí)暴露出這種插值方法的缺陷，現(xiàn)實(shí)中存在一些較為成熟的插值算法。為了追求更高的計(jì)算精度和求解可靠性，有必要詳細(xì)根據(jù)流體特征分析速度的插值算法。

最后，由所得的ξ，η來確定加權(quán)系數(shù)，進(jìn)行插值。

對(duì)于物質(zhì)輸移的逆向追蹤，要將單元細(xì)分為4個(gè)子區(qū)，四邊形是采用對(duì)邊中點(diǎn)連線劃分成4個(gè)子四邊形，三角形則是3條中位線劃分成4個(gè)子三角形，判斷逆向追蹤后質(zhì)點(diǎn)所屬子區(qū)，在子區(qū)內(nèi)部進(jìn)行插值計(jì)算（圖6、圖7）。圖6，圖7中areas（1），areas（2），areas（3），areas（4）表示所在子三角形面積，作為插值計(jì)算的權(quán)重。

圖6 物質(zhì)輸移的四邊形插值Fig.6 Interpolation on quadrilateral unit for transport equation

圖7 物質(zhì)輸移的三角形插值Fig.7 Interpolation on triangle unit for transport equation

4 定解條件

初始場(chǎng)分布條件是在節(jié)點(diǎn)處給定的，溫度和鹽度都要進(jìn)行邊、單元初始化賦值。陸地邊界和開邊界在計(jì)算中有很大的差別，陸地邊界實(shí)際限定了流場(chǎng)的范圍，干濕判別標(biāo)準(zhǔn)是根據(jù)水深是否大于某一臨界值確定。在高出水面的陸地區(qū)域，也要參與網(wǎng)格剖分，所有剖分的網(wǎng)格都會(huì)參與計(jì)算；垂向網(wǎng)格剖分時(shí)，水體以下部分的柱體也參與計(jì)算，相同水平位置的無水柱體，均賦予同一水平位置上最底部有水單元的值，這樣處理是為了網(wǎng)格整齊，方便調(diào)用求解。

對(duì)于水位，一般會(huì)有全區(qū)域的單一初始值，在計(jì)算過程中會(huì)逐漸調(diào)整為合理值。恒定水流求解問題水流只需給出流量的上邊界（進(jìn)口）條件和水位的下邊界（出口）條件，鹽度和溫度場(chǎng)則根據(jù)需要給出進(jìn)口處的變化過程，非恒定水流求解問題要給出流量的上邊界過程和水位的下邊界過程。要將給出的進(jìn)口處流量轉(zhuǎn)化為流速，只是簡(jiǎn)單的利用開邊界寬度來除流量，得到單寬流量，再除以邊界水深，因此研究區(qū)域和進(jìn)口邊界有相當(dāng)距離，避免失真流速的影響。

5 網(wǎng)格前處理

在預(yù)先分層的網(wǎng)格上確定實(shí)際水體的頂?shù)讓樱ㄔ闯绦驑?biāo)注中M表示頂層號(hào)，m表示底層號(hào)），初始化單元、節(jié)點(diǎn)、邊的分層信息。每次迭代求解出水位值后，都要重新計(jì)算單元、節(jié)點(diǎn)、邊的分層信息。以單元為例進(jìn)行說明，實(shí)際水深小于某一特定值時(shí)，單元為無水單元，頂層號(hào)kfe＝0；否則計(jì)算單元新的底層號(hào)kbe，頂層號(hào)kfe，各層厚度dc和上下鄰層平均厚度dchalf。如圖8，ztot指的是層面的垂直坐標(biāo)，nvrt個(gè)層對(duì)應(yīng)nvrt＋1個(gè)層面，jlev表示層號(hào)。逆向追蹤中用到的zt表示某質(zhì)點(diǎn)的實(shí)際高度，zup指該點(diǎn)到所在層層頂?shù)木嚯x。頂?shù)淄瑢訒r(shí)，dc＝dchalf＝水深；不同層時(shí)頂?shù)讓雍馾c和dchalf，如圖9所示。

圖9 頂?shù)讓邮疽鈭DFig.9 Schem atic diagram of top and bottom elements

圖8 垂向分層示意圖Fig.8 Schematic diagram of verticalmeshes

dzhalf，dz要在動(dòng)量方程做分母（dzhalf（i，kbs－1）、dzhalf（i，kfs）不做分母，kbs，kfs表示邊的底層號(hào)和頂層號(hào)），所以要檢查兩者是否小于某限定值denom（略小于判定網(wǎng)格干濕時(shí)的臨界水深），初始水位線應(yīng)盡量靠近分層的中上部，作為簡(jiǎn)化，程序中的初始水位設(shè)定要求水位線不能靠近層面線。迭代過程中，這種檢查只針對(duì)頂?shù)淄瑢印g，vg表示j邊上的第k層面上的流速可通過式（15）、（16）線性插值得到

圖10 四邊形單元節(jié)點(diǎn)插值Fig.10 Interpolation w ith nodes on quadrilateral unit

其中utmp（k），vtmp（k）在柱體第k層側(cè)面（平面網(wǎng)格中仍對(duì)應(yīng)j邊）面心上的橫向和縱向流速。

如圖10所示根據(jù)距離倒數(shù)加權(quán)平均求得4條邊nd1，nd2，nd3，nd4的公共節(jié)點(diǎn)i的x方向流速uu2（i，k）為

式中u1，u2，u3，u4和d1，d2，d3，d4分別對(duì)應(yīng)nd1，nd2，nd3，nd4四邊中點(diǎn)x方向流速和四邊長(zhǎng)度。Y方向流速vv2（i，k）為

v

g（nd1，k），vg（nd2，k），vg（nd3，k），vg（nd4，k）分別對(duì)應(yīng)nd1，nd2，nd3，nd4四邊中點(diǎn)y方向流速，distj（nd1），distj（nd2），distj（nd3），distj（nd4），即式（17）中的d1，d2，d3，d4，表示所在邊的長(zhǎng)度。節(jié)點(diǎn)垂向流速ww2（i，k）則是按周圍4個(gè)單元垂向流速面積加權(quán)平均求得。

其中areas（1），areas（2），areas（3），areas（4）為4單元面積，we1，we2，we3，we4為4單元垂向流速。節(jié)點(diǎn)垂向流速則是按周圍4個(gè)單元垂向流速面積加權(quán)平均求得。

6 結(jié) 語

ELCIRC作為水動(dòng)力學(xué)開放源代碼，在河口海岸數(shù)值模擬中取得一定的成功。文中給出了部分源代碼的具體分析，可以為水動(dòng)力學(xué)模型開發(fā)人員提供方法實(shí)現(xiàn)的參考。源代碼凝結(jié)了美國(guó)俄勒岡健康與科學(xué)大學(xué)海岸陸地邊緣陸緣研究中心科研人員的辛勤勞動(dòng)，在此表示感謝，筆者只是對(duì)其優(yōu)秀的程序進(jìn)行部分粗略的學(xué)習(xí)，重點(diǎn)著眼于其核心部分的歐拉拉格朗日模型逆向追蹤過程，發(fā)現(xiàn)了其中的不完善的地方，希望能給源程序在改進(jìn)過程中盡一點(diǎn)努力。ELCIRC的開發(fā)人員已經(jīng)在此代碼的基礎(chǔ)上進(jìn)行改進(jìn)，開發(fā)出SELFE程序，本文的大部分都可以為讀者學(xué)習(xí)SELFE提供幫助。

致謝：非常感謝清華大學(xué)研究生李健在ELCIRC源代碼學(xué)習(xí)上給予悉心指導(dǎo)和莫大的幫助！

［1］ ZHANG Y L，BAPTISTA A M，MYERSE P.A Crossscale Model for 3D Baroclinic Circulation in Estuary-Plume-Shelf Systems：I.Formulation and Skill Assessment［J］.Continental Shelf Research，2004，24（18）：2174－2214.

［2］ MYERSE P，AIKMAN F，ZHANG A J.A Forecast Circulation Model of the St.Johns River，F(xiàn)lorida［C］∥Proceedings of the Eighth International Conference on Estuarine and Coastal Modeling，November 3－5，2003， Monterey，California，USA：144－156.

［3］ 楊金艷.ELCIRC模型在長(zhǎng)江口的應(yīng)用［D］.南京：河海大學(xué)，2006.（YANG Jin-yan.Application of ELCRIC Model in Changjiang Estuary［D］.Nanjing：Hohai University，2006.（in Chinese））

［4］ GONG W P，SHEN J，WANG D R.Mean Water Level Setup／Setdown in the Inlet-Lagoon System Induced by Tidal Action：A Case Study of Xincun Inlet，Hainan Island in China［J］.Acta Oceanologica Sinica，2008，27（5）：63－80

［5］ GONG W P，SHEN J，CHO K H，et al.A Numerical Model Study of Barotropic Subtidal Water Exchange Between Estuary and Subestuaries（Tributaries）in the Chesapeake Bay During Northeaster Events［J］.Ocean Modelling，2009，26：170－189

［6］ DIAS JM，SOUSA M C，BERTIN X，et al.Numerical Modeling of the Impactof the Anc~ao Inlet Relocation（Ria Formosa，Portugal）［J］.Environmental Modelling＆Software，2009，24（6）：711－725.

［7］ PICADO A，DIAS JM，F(xiàn)ORTUNATO A B.Tidal Changes in Estuarine Systems Induced by LocalGeomorphologic Modifications［J］.Continental Shelf Research，2010，30（17）：1854－1864.

［8］ OHSU’s Center for Coastal and Land-Margin Research.ELCIRC Description：3D Baroclinic Circulation Model［EB／OL］.（2006－05－01）［2011－05－14］，http：／／www.stccmop.org／CORIE／software／elcirc／r ecent.html

（編輯：劉運(yùn)飛）

Algorithm of ELCIRC Source Code

YANG Fei1，YANGWen-jun1，YANG Sen2
（1.Key Laboratory of River Regulation and Flood Control of MWR，Yangtze River Scientific Research Institute，Wuhan 430010，China；2.Water Affairs Bureau ofWulong County，Chongqing City，Chongqing 408500，China）

Open-source code ELCIRC（Eulerian-Lagrangian CIRCulation）solves the shallow water equations using a semi-implicit Eulerian-Lagrangian finite volume／finite difference method with horizontally unstructured grids and vertically unstretched z-coordinates.ELCIRC source code aside from the governing equation solution is analyzed in this paper.The operation of interpolation，topological structure，definite condition，hierarchy information，and wetting and dryingmethod are described in detail.Back-track，as the key technology of ELCIRC，is expounded comprehensively.Moreover，some detailed problems are discussed，and a few references and suggestions are given to avoid the errors caused by similar number subtraction and by employing small number as the denominator.ELCIRC has been applied in practice and could serve as a valuable reference for developing hydrodynamic models.

ELCIRC；back-track；source code；3-D hydrodynamic model；algorithm

TV148

A

1001－5485（2013）05－0097－06

2013，30（05）：97－102

10.3969／j.issn.1001－5485.2013.05.021

2012－04－25；

2012－06－13

國(guó)家自然科學(xué)基金面上項(xiàng)目（51079008）

楊 飛（1985－），男，山東泰安人，碩士研究生，從事水力學(xué)及河流動(dòng)力學(xué)研究，（電話）13810609717（電子信箱）yangfeihaoyun＠163.com。