第一節(jié)海洋數(shù)值模型的發(fā)展概況觀測(cè)，試驗(yàn)理論分析數(shù)值模型海洋的研究手段海洋數(shù)值模型發(fā)展的客觀要求海水運(yùn)動(dòng)方程推導(dǎo)過(guò)程中的動(dòng)力學(xué)假設(shè)，比如海水靜力學(xué)假設(shè)，都是建立在大部分空間尺度大于10km，時(shí)間尺度大于慣性周期（2π/f）的情況下的海水運(yùn)動(dòng)方程不存在精確解：（1）偏微分方程本身是非線性的，（2）將一個(gè)實(shí)際問(wèn)題具體化時(shí)所需要的環(huán)境場(chǎng)如海底地形、海岸線幾何形狀、表面驅(qū)動(dòng)力等都無(wú)法給出解析函數(shù)表達(dá)式，只能在離散化的空間和時(shí)間區(qū)域取得因此，求解海水運(yùn)動(dòng)方程只能通過(guò)近似求解的辦法來(lái)確定，而這些數(shù)值解又可能引入明顯的近似誤差海洋數(shù)值模型發(fā)展的歷史海洋模型到按其水平網(wǎng)格的離散方式以及所使用的垂向坐標(biāo)系的不同大致經(jīng)歷了如下幾個(gè)發(fā)展階段最早出現(xiàn)并且還在使用的海洋模型是Bryan等人開發(fā)的基于原始方程的低階精度的有限差分模型，它在水深方向采用z坐標(biāo)系目前常見的HOPS（HarvardOceanPredictionSystem），MOM（

GFDL

ModularOceanModel），POP

(ParallelOceanProgram)，NCOM（NCARCommunityOceanModel）模型在某種程度上都可以認(rèn)為是該模型改進(jìn)版海洋數(shù)值模型發(fā)展的歷史上個(gè)世紀(jì)70年代，sigma坐標(biāo)系開始應(yīng)用于海洋模型在水深方向，比如目前被廣泛使用的POM（PrincetonOceanModel）、ECOM（EstuarineCoastalandOceanModel）、ROMS（RegionalOceanModelingSystem）模型都屬于這種類型的模型與傳統(tǒng)的

z

坐標(biāo)系相比，sigma坐標(biāo)可以更好地貼合海底地形變化海洋數(shù)值模型發(fā)展的歷史自上世紀(jì)90年代以來(lái)，非結(jié)構(gòu)化網(wǎng)格技術(shù)開始在海洋模型中得到應(yīng)用，相應(yīng)的也出現(xiàn)了基于非結(jié)構(gòu)化網(wǎng)格離散方式的有限元模型，如SEOM模型（SpectralFiniteElement

Ocean

Model），和有限體積模型，如FVCOM雖然與傳統(tǒng)的結(jié)構(gòu)化網(wǎng)格相比,非結(jié)構(gòu)化網(wǎng)格可以更好地?cái)M合陸地邊界，但是代碼實(shí)現(xiàn)上的困難以及計(jì)算穩(wěn)定性的問(wèn)題使其迄今還沒有得到非常廣泛的應(yīng)用新一代的海洋數(shù)值模型新一代的海洋模型廣泛采用隨地坐標(biāo)系（terrain-followingcoordinates），進(jìn)而促進(jìn)了有關(guān)時(shí)間步長(zhǎng)，對(duì)流項(xiàng)和壓力梯度項(xiàng)等數(shù)值算法的改進(jìn)進(jìn)入新世紀(jì)以來(lái)，下一代的海洋數(shù)值動(dòng)力模型正在緊鑼密鼓的研制中，代表性的是TOMS(Terrain-followingOceanModelingSystem)，它融合了目前最先進(jìn)的物理知識(shí)、數(shù)值方法和數(shù)據(jù)同化技術(shù)

/WWWPUBLIC/htdocs.pom/TOMS.htm

海洋模型的商業(yè)軟件就海流的仿真建模研究而言，現(xiàn)階段國(guó)外已形成了不少具有代表性的、具有強(qiáng)大前后處理功能及核心仿真建模技術(shù)的軟件或系統(tǒng)（商業(yè)軟件）如加拿大DalhousieUniversity研究的DALCOAST河口海岸預(yù)報(bào)系統(tǒng)；丹麥的DHIWater&Environment機(jī)構(gòu)開發(fā)的的MIKE系列軟件系統(tǒng)；荷蘭的WLDelftHydraulics開發(fā)的的Delft－3D軟件第二節(jié)

基于POM模型的一般介紹水動(dòng)力模型基本特點(diǎn)垂直方向的

??坐標(biāo)變換水平網(wǎng)格采用正交曲線坐標(biāo)和ArakawaC差分格式水平時(shí)間差分采用顯格式，而垂直差分為隱格式自由表面可以模擬水位變化垂直和水平方向的混合擴(kuò)散分別采用2.5階的Mellor-Yamada湍流閉合模式和Smagorinski模式內(nèi)外模態(tài)分別處理速度較慢的內(nèi)重力波和速度較快的外重力波以提高整個(gè)模式計(jì)算效率包含了海水的熱動(dòng)力過(guò)程Sigma坐標(biāo)變換：hs=-1z=0z=H(x,y)

s=0

V(i,j+1)

V(i,j)

U(i+1,j)

U(i,j)

??

(i,j)模型采用（a）垂向sigma坐標(biāo)系和（b）水平ArakawaC網(wǎng)格(1)控制方程連續(xù)方程海水運(yùn)動(dòng)方程(2)(3)(4)溫度方程湍流動(dòng)能方程鹽度方程湍流混合長(zhǎng)度方程(5)(6)(7)(8)其中，ω是基于sigma坐標(biāo)系的垂向流速矢量，它與笛卡兒坐標(biāo)系下的垂向流速W之間的關(guān)系可表示為另外，方程中分別表示海水的實(shí)際密度，Boussinesq近似密度以及密度擾動(dòng)其中，是實(shí)際水深(9)M—Y湍流閉合模型（profq）在湍流動(dòng)能及混合長(zhǎng)度方程中分別表示湍流混合系數(shù)，熱擴(kuò)散系數(shù)和湍流動(dòng)能的垂直擴(kuò)散系數(shù)；墻近似函數(shù)，其中，κ=0.4是von

Karman常數(shù)湍流動(dòng)能和混合長(zhǎng)度方程組由下列等式關(guān)系閉合其中，是穩(wěn)定函數(shù)(10)在2.5階M—Y湍流閉合模型中，可表示為其中，。在不穩(wěn)定層次條件下，的上限值為0.023，對(duì)應(yīng)的值分別為2.0145和2.4401；在穩(wěn)定層次條件下，它下限值為-0.28，對(duì)應(yīng)的值分別為0.0470和0.0461。

是2.5階M-Y湍流閉合模式參數(shù)，由實(shí)驗(yàn)測(cè)得(11)M—Y湍流模型的適用條件M-Y湍流模型已被廣泛應(yīng)用于淺海潮汐，風(fēng)生表面混合層的模擬以及底邊界層的研究；該模型在混合較弱的層化流體中對(duì)湍流混合系數(shù)計(jì)算效果不佳，應(yīng)用于河口的可行性也有待于探討；雖然該模型無(wú)論是在物理上還是在數(shù)學(xué)上都存在著明顯的缺陷，但是目前依然得到廣泛應(yīng)用。只有深入了解海水的微細(xì)結(jié)構(gòu)以及海洋湍流結(jié)構(gòu)和性質(zhì)后，才有可能構(gòu)建出更為完善的湍流閉合模型海水運(yùn)動(dòng)方程的水平擴(kuò)散項(xiàng)其中水平切應(yīng)力是水平渦粘性系數(shù)，由Smagorinsky模式計(jì)算(12)(13)Smagorinsky水平擴(kuò)散模式其中參數(shù)C

為無(wú)量綱值，其取值為0.1或0.2，如果計(jì)算網(wǎng)格劃分得足夠細(xì)的話，C可以取值為0

由Smagorinsky模式計(jì)算的水平渦粘性系數(shù)隨著網(wǎng)格精度的改善和速度梯度的減小而減小(14)溫度鹽度的水平擴(kuò)散項(xiàng)其中是水平熱擴(kuò)散系數(shù)，一般地，是一個(gè)小值，取0.1或0.2，甚至在某些情況下可以取為0(15)(16)垂向邊界條件1.海面(17)2.海底其中，分別是湍流閉合模型參數(shù)和摩擦速度；分別是海面和海底的摩擦系數(shù)；(18)模型的外模態(tài)控制近岸環(huán)流的海水運(yùn)動(dòng)方程包含了運(yùn)動(dòng)速度較快的外重力波和速度較慢的內(nèi)重力波，因此進(jìn)行模式分裂可以大大提高模型的計(jì)算效率模型的外模態(tài)就是計(jì)算自由水位和垂向平均的流速變化，需要較小的時(shí)間步長(zhǎng)內(nèi)模態(tài)主要模擬流速矢量、溫度、鹽度等三維結(jié)構(gòu)變化，對(duì)時(shí)間步長(zhǎng)的要求相對(duì)寬松采用模式分裂技巧可以有效地減少因外模態(tài)計(jì)算所消耗的計(jì)算時(shí)間，從而提高計(jì)算效率垂向積分的海水運(yùn)動(dòng)方程其中：連續(xù)方程海水運(yùn)動(dòng)方程(19)(20)(21)水平擴(kuò)散項(xiàng)水平數(shù)值耗散項(xiàng)(22)(23)模型代碼中的一些符號(hào)im,jm,kb,imm1,jmm1,kbm1,mode,isplitdte,dti,days,umolz,zz,dz,dzzaam2d,art,aru,arv,dum,dvm,fsmdx,dy,h,el,d,ua,va,ut,vt,et,corswradwusurf,wvsurf,wubot,wvbot,wtsurf,wssurf,radu,v,ω,t,s,rho,l,q2km,kh,kq,aam,aahrmean,tclim,sclim數(shù)值積分

流

程

圖外模態(tài)和內(nèi)模態(tài)的相互作用內(nèi)模態(tài)的計(jì)算方法三維變量的計(jì)算（以T為例）可以分解為垂直擴(kuò)散項(xiàng)的計(jì)算和水平對(duì)流擴(kuò)散項(xiàng)的計(jì)算，以溫度方程為例方程寫為求解過(guò)程分為兩步，第一步計(jì)算水平對(duì)流擴(kuò)散項(xiàng)，采用中心差分格式，由advt子程序?qū)崿F(xiàn)(24)(25)水平對(duì)流擴(kuò)散項(xiàng)垂直擴(kuò)散項(xiàng)第二步計(jì)算垂向擴(kuò)散項(xiàng)，數(shù)值方法采用“蛙跳”格式，由proft子程序?qū)崿F(xiàn)由于“蛙跳”格式在奇數(shù)時(shí)間步時(shí)的解與在偶數(shù)時(shí)間步時(shí)的解會(huì)發(fā)生偏離，因此每一時(shí)間步的計(jì)算結(jié)束后還需對(duì)上述解進(jìn)行平滑處理，即α為一個(gè)小值，一般取0.05。處理后，(26)(27)計(jì)算網(wǎng)格的安排二維外模態(tài)網(wǎng)格分布三維內(nèi)模態(tài)網(wǎng)格分布水平對(duì)流項(xiàng)的計(jì)算■雖然模型采用有限差分計(jì)算格式，但是對(duì)于每個(gè)網(wǎng)格對(duì)流項(xiàng)的計(jì)算都是按照有限容積的方法進(jìn)行處理，即溫度對(duì)流過(guò)程可表示為■速度的對(duì)流過(guò)程與溫度相似，可寫為其中，是由z坐標(biāo)轉(zhuǎn)化為σ坐標(biāo)后產(chǎn)生的彎曲項(xiàng)(28)(29)垂直擴(kuò)散項(xiàng)的計(jì)算■垂直擴(kuò)散項(xiàng)的計(jì)算公式（第k層，1<k<kb-1）■可以改寫為(30)(31)■展開后合并同類項(xiàng)得到其中(32)(33)■現(xiàn)在假定溫度解的形式為■由(34)得到的代入到(32)就可以得到其中的值由(33)求得，的值由上一層的值確定(34)(35)■當(dāng)k=1時(shí)，即表層海水，海水溫度主要取決于海表面溫度通量。公式（34）的解可以近似寫為■這樣，表層海水的溫度只能根據(jù)式（31）求解，即■

上式可以進(jìn)一步表示為(36)(37)■式（37）還可以進(jìn)一步寫為■再與溫度的通解（34）比較■

可以得到(34)(38)(39)■當(dāng)短波輻射通量的計(jì)算公式■

r，ad1，ad2均為光輻射常數(shù)，根據(jù)不同水質(zhì)取值如下(40)ntp

12345JerlovtypeIIaIbIIIIIr

0.580.620.670.700.78ad1(m)0.350.60

1.01.51.4ad2(m)

23.020.017.014.07.9來(lái)源：Jerlov，1976；PaulsonandSimpson，1977■當(dāng)k=kb-1時(shí)，即底層海水，假定海底的熱通量為0，根據(jù)前面的推導(dǎo)方式同樣可以得到底層海水溫度的計(jì)算公式■

對(duì)于鹽度方程垂向擴(kuò)散項(xiàng)的計(jì)算與溫度方程相同，只是不考慮太陽(yáng)短波輻射這一項(xiàng)(41)時(shí)間步長(zhǎng)的CFL限制條件●由于模型水平方向采用顯格式，因此時(shí)間步長(zhǎng)的選取必須要滿足CFL穩(wěn)定條件?！駥?duì)于外模態(tài)，時(shí)間步長(zhǎng)的限制條件為其中，可能預(yù)見到的最大流速●對(duì)于內(nèi)模態(tài)，時(shí)間步長(zhǎng)的限制條件較外模態(tài)的情形寬松很多，主要是速度較快的外重力波已經(jīng)在外模態(tài)中考慮了。一般取值30~50即可

(42)時(shí)間步長(zhǎng)的其它限制條件●對(duì)于動(dòng)量或標(biāo)量還有其它一些時(shí)間限制條件其中或●以及

其中分別為地球角速率和地理緯度(43)(44)側(cè)開邊界條件（bcond）陸地及岸線是由dum、dvm、fsm控制的，在陸地上這些變量的值設(shè)為0，有水的地方設(shè)為1子程序bcond（idx），idx=1對(duì)應(yīng)的水位邊界條件；idx=2對(duì)應(yīng)垂向平均流速邊界條件；idx=3對(duì)應(yīng)三維水平流速邊界條件；idx=4對(duì)應(yīng)溫鹽邊界條件；idx=5對(duì)應(yīng)垂向流速邊界條件；idx=6對(duì)應(yīng)湍流動(dòng)能和湍流混合長(zhǎng)度邊界條件模型對(duì)開邊界條件的要求很高，而開邊界條件本身又具有很大的不確定性，因此有必要對(duì)模型的外模態(tài)和內(nèi)模態(tài)的開邊界分別進(jìn)行處理數(shù)值積分

流

程

圖FormulaBoundaryCodeInflowcondition

EASTuaf(im,j)=2*bc(j)/(h(im,j)+elf(im,j)+h(imm1,j+elf(imm1,j))elf(im,j)=elf(imm1,j)vaf(im,j)=setWESTuaf(2,j)=2*bc(j)/(h(1,j)+elf(1,j)+h(2,j)+elf(2,j))elf(1,j)=elf(2,j)vaf(1,j)=setNORTHvaf(i,jm)=2*bc(i)/(h(i,jm)+elf(i,jm)+h(i,jmm1)+elf(i,jmm1))elf(i,jm)=elf(i,jmm1)uaf(i,jm)=setSOUTHvaf(i,2)=2*bc(i)/(h(i,1)+elf(i,1)+h(i,2)+elf(i,2))elf(i,1)=elf(i,2)uaf(i,1)=setElevationconditionh

=BC

EASTelf(imm1,j)=bc(j)elf(im,j)=elf(imm1,j)cosmeticuaf(im,j)=uaf(imm1,j)vaf(im,j)=setWESTelf(2,j)=bc(j)uaf(2,j)=uaf(3,j)vaf(1,j)=setNORTHelf(i,jmm1)=bc(i)elf(i,jm)=elf(i,jmm1)cosmeticvaf(i,jm)=vaf(i,jmm1)uaf(i,jm)=setSOUTHelf(i,2)=bc(i)vaf(i,2)=vaf(i,3)uaf(i,1)=set

外模態(tài)

開

邊

界

條

件

一FormulaBoundaryCodeRadiationEASTuaf(im,j)=sqrt(grav/h(imm1,j))*el(imm1,j)+bc(j)elf(im,j)=elf(imm1,j)vaf(im,j)=setWESTuaf(2,j)=-sqrt(grav/h(2,j))*el(2,j)+bc(j)elf(1,j)=elf(2,j)vaf(1,j)=setNORTHvaf(i,jm)=sqrt(grav/h(i,jmm1))*el(i,jmm1)+bc(i)elf(i,jm)=elf(i,jmm1)uaf(i,jm)=setSOUTHvaf(i,2)=-sqrt(grav/h(i,2))*el(i,2)+bc(i)elf(i,1)=elf(i,2)uaf(i,1)=setRadiation

EASTgae=dte*sqrt(grav*h(im,j))/dx(im,j)uaf(im,j)=gae*ua(imm1,j)+(1.-gae)*ua(im,j)elf(im,j)=elf(imm1,j)vaf(im,j)=setWESTgae=dte*sqrt(grav*h(2,j))/dx(2,j)uaf(2,j)=gae*ua(3,j)+(1.-gae)*ua(2,j)elf(1,j)=elf(2,j)vaf(1,j)=setNORTHgae=dte*sqrt(grav*h(i,jm))/dy(i,jm)vaf(i,jm)=gae*va(i,jmm1)+(1.-gae)*va(i,jm)elf(i,jm)=elf(i,jmm1)uaf(i,jm)=setSOUTHgae=dte*sqrt(grav*h(i,2))/dy(i,2)vaf(i,2)=gae*va(i,3)+(1.-gae)*va(i,2)elf(i,1)=elf(i,2)uaf(i,1)=set

外模態(tài)

開

邊

界

條

件

二流速水位流速FormulaBoundaryCodeRadiationEASTgae=dte*sqrt(grav*h(imm1,j))/dx(imm1,j)elf(imm1,j)=gae*el(imm2,j)+(1.-gae)*el(imm1,j)elf(im,j)=elf(imm1,j)uaf(im,j)=uaf(imm1,j)vaf(im,j)=setWESTgae=dte*sqrt(grav*h(2,j))/dx(1,j)elf(2,j)=gae*el(3,j)+(1.-gae)*el(2,j)uaf(2,j)=uaf(3,j)vaf(2,jm)=setNORTHgae=dte*sqrt(grav*h(i,jmm1))/dy(i,jmm1)elf(i,jmm1)=gae*el(i,jmm2)+(1.-gae)*el(i,jmm1)elf(i,jm)=elf(i,jmm1)vaf(i,jm)=vaf(i,jmm1)uaf(i,jm)=setSOUTHgae=dte*sqrt(grav*h(i,2))/dy(i,2)elf(i,2)=gae*el(i,3)+(1.-gae)*el(i,2)vaf(i,2)=vaf(i,3)，uaf(i,1)=setCyclicEAST(I=IM)elf(im,j)=elf(3,j)uaf(im,j)=uaf(3,j)，vaf(im,j)=vaf(3,j)WEST(I=1)elf(1,j)=elf(imm2,j)，elf(2,j)=elf(imm1,j)uaf(2,j)=uaf(imm1,j)，vaf(2,j)=vaf(imm1,j)NORTH(J=JM)elf(i,jm)=elf(i,3)uaf(i,jm)=uaf(i,3)vaf(i,jm)=vaf(i,3)SOUTH(J=1)elf(i,1)=elf(i,jmm2)，elf(i,2)=elf(i,jmm1)uaf(i,2)=uaf(i,jmm1)vaf(i,2)=vaf(i,jmm1)

外模態(tài)

開

邊

界

條

件

三水位FormulaBoundaryCodeInflowconditionEASTuf(im,j,k)=bc(j,k)vf(im,j,k)=setU=BC

WESTuf(2,j,k)=bc(j,k)vf(1,j,k)=set

NORTHvf(i,jm,k)=bc(i,k)uf(i,jm,k)=setSOUTHvf(i,2,k)=bc(i,k)uf(i,1,k)=setRadiation:EASTgai=sqrt(h(im,j)/hmax)uf(im,j,k)=gai*u(imm1,j,k)+(1.-gai)*u(im,j,k)vf(im,j,k)=setWESTgai=sqrt(h(2,j)/hmax)uf(2,j,k)=gai*u(3,j,k)+(1.-gai)*u(2,j,k)vf(1,j,k)=setNORTHgai=sqrt(h(i,jm)/hmax)vf(i,jm,k)=gai*v(i,jmm1,k)+(1.-gai)*v(i,jm,k)uf(i,jm,k)=setSOUTHgai=sqrt(h(i,2)/hmax)vf(i,2,k)=gai*v(i,3,k)+(1.-gai)*v(i,2,k)uf(i,1,k)=set

內(nèi)模態(tài)開邊界條件（一）FormulaBoundaryCodeUpstreamadvectiononTorSEASTuf(im,j,k)=t(im,j,k)-dti/(dx(im,j)+dx(imm1,j))*((u(im,j,k)+abs(u(im,j,k)))*(t(im,j,k)-t(imm1,j,k))+(u(im,j,k)-abs(u(im,j,k)))*(tbe(j,k)-t(im,j,k)))WESTuf(1,j,k)=t(1,j,k)-dti/(dx(1,j)+dx(2,j))*((u(1,j,k)+abs(u(1,j,k)))*(t(1,j,k)-tbw(j,k))+(u(1,j,k)-abs(u(1,j,k)))*(t(2,j,k)-t(1,j,k)))