1、第一节海洋数值模型的第一节海洋数值模型的发展概况发展概况观测，试验观测，试验理论分析理论分析数值模型数值模型海洋的研究手段海洋的研究手段1ppt课件海洋数值模型发展的客观要求海洋数值模型发展的客观要求海水运动方程推导过程中的动力学假设，比如海水海水运动方程推导过程中的动力学假设，比如海水静力学假设，都是建立在大部分空间尺度大于静力学假设，都是建立在大部分空间尺度大于10km10km，时间尺度大于惯性周期（时间尺度大于惯性周期（2 2/f/f）的情况下的）的情况下的海水运动方程不存在精确解：（海水运动方程不存在精确解：（1 1）偏微分方程本身）偏微分方程本身是非线性的，（是非线性的，（2 2）将
2、一个实际问题具体化时所需要）将一个实际问题具体化时所需要的环境场如海底地形、海岸线几何形状、表面驱动的环境场如海底地形、海岸线几何形状、表面驱动力等都无法给出解析函数表达式，只能在离散化的力等都无法给出解析函数表达式，只能在离散化的空间和时间区域取得空间和时间区域取得因此，求解海水运动方程只能通过近似求解的办法因此，求解海水运动方程只能通过近似求解的办法来确定，而这些数值解又可能引入明显的近似误差来确定，而这些数值解又可能引入明显的近似误差2ppt课件海洋数值模型发展的历史海洋数值模型发展的历史海洋模型到按其水平网格的离散方式以及所使用的海洋模型到按其水平网格的离散方式以及所使用的垂向坐标系的
3、不同大致经历了如下几个发展阶段垂向坐标系的不同大致经历了如下几个发展阶段最早最早出现并且还在使用的海洋模型是出现并且还在使用的海洋模型是BryanBryan等人开发等人开发的基于原始方程的低阶精度的有限差分模型，它在的基于原始方程的低阶精度的有限差分模型，它在水深方向采用水深方向采用z z坐标系坐标系目前常见的目前常见的HOPS（ Harvard Ocean Prediction System ），MOM（ GFDL Modular Ocean Model），POP (Parallel Ocean Program)，NCOM（NCAR Community Ocean Model）模型在某种模型
4、在某种程度上都可以认为是该模型改进版程度上都可以认为是该模型改进版3ppt课件海洋数值模型发展的历史海洋数值模型发展的历史上个世纪上个世纪7070年代年代，sigma坐标系开始应用于海洋坐标系开始应用于海洋模型在水深方向，比如目前被广泛使用的模型在水深方向，比如目前被广泛使用的POM（Princeton Ocean Model）、ECOM （Estuarine Coastal and Ocean Model） 、ROMS（Regional Ocean Modeling System）模型都属于这种类型的模型模型都属于这种类型的模型与传统的与传统的 z z 坐标系相比，坐标系相比，sigmasi
5、gma坐标可以更好地坐标可以更好地贴合海底地形变化贴合海底地形变化4ppt课件海洋数值模型发展的历史海洋数值模型发展的历史自上世纪自上世纪9090年代以来年代以来，非结构化网格技术开始在，非结构化网格技术开始在海洋模型中得到应用，相应的也出现了基于非结海洋模型中得到应用，相应的也出现了基于非结构化网格离散方式的有限元模型构化网格离散方式的有限元模型，如如SEOM模型模型（ Spectral Finite Element Ocean Model），和有），和有限体积模型限体积模型，如如FVCOM虽然与传统的结构化网格相比虽然与传统的结构化网格相比, ,非结构化网格可以非结构化网格可以更好地拟合陆
6、地边界，但是代码实现上的困难以更好地拟合陆地边界，但是代码实现上的困难以及计算稳定性的问题使其迄今还没有得到非常广及计算稳定性的问题使其迄今还没有得到非常广泛泛的的应用应用5ppt课件新一代的海洋数值模型新一代的海洋数值模型新一代的海洋模型广泛采用随地坐标系（新一代的海洋模型广泛采用随地坐标系（ terrain-following coordinates），进而促进了有关时间步），进而促进了有关时间步长，对流项和压力梯度项等数值算法的改进长，对流项和压力梯度项等数值算法的改进进入新世纪以来进入新世纪以来，下一代的海洋数值动力模型正在，下一代的海洋数值动力模型正在紧锣密鼓的研制中，代表性的是紧锣
7、密鼓的研制中，代表性的是TOMS (Terrain-following Ocean Modeling System)，它融合了目前，它融合了目前最先进的物理知识、数值方法和数据同化技术最先进的物理知识、数值方法和数据同化技术 http:/www.aos.princeton.edu/WWWPUBLIC/htdocs.pom/TOMS.htm 6ppt课件海洋模型的商业软件海洋模型的商业软件就海流的仿真建模研究而言，现阶段国外已形成就海流的仿真建模研究而言，现阶段国外已形成了不少具有代表性的、具有强大前后处理功能及了不少具有代表性的、具有强大前后处理功能及核心仿真建模技术的软件或系统核心仿真建模技
8、术的软件或系统（商业软件）（商业软件）如加拿大如加拿大Dalhousie University研究研究的的DALCOAST河口海岸预报系统河口海岸预报系统；丹麦的丹麦的DHI Water & Environment 机构开发的的机构开发的的MIKE系列系列软件系统软件系统；荷兰的荷兰的WL Delft Hydraulics开发的的开发的的Delft3D软件软件7ppt课件8ppt课件水动力模型基本特点水动力模型基本特点n垂直方向的垂直方向的 坐标坐标变换变换n水平网格采用正交曲线坐标和水平网格采用正交曲线坐标和ArakawaCArakawaC差分格式差分格式n水平时间差分采用显格式，而垂直差分
9、为隐格式水平时间差分采用显格式，而垂直差分为隐格式n自由表面可以模拟水位变化自由表面可以模拟水位变化n垂直和水平方向的混合扩散分别采用垂直和水平方向的混合扩散分别采用2.52.5阶的阶的Mellor-YamadaMellor-Yamada湍流闭合模式和湍流闭合模式和SmagorinskiSmagorinski模式模式n内外模态分别处理速度较慢的内重力波和速度较内外模态分别处理速度较慢的内重力波和速度较快的外重力波以提高整个模式计算效率快的外重力波以提高整个模式计算效率n包含了海水的热动力过程包含了海水的热动力过程9ppt课件Sigma坐标变换：坐标变换：hs = -1z = 0z = H(x,
10、 y) (a)s = 0 V (i, j+1) V (i, j) U (i+1, j) U (i, j) (i, j)(b)模型采用（模型采用（a a）垂向）垂向sigmasigma坐标系和（坐标系和（b b）水平）水平Arakawa CArakawa C网格网格zHhsh-=(1)(1)10ppt课件控制方程控制方程0DUDVxyths=2UDU DUVDUfVDgDtxyxhs-20MxoKgDDUdFxDxDssssss-=2VDUVDV DVfUDgDtxyyhs20MyoKgDDVdFyDyDssssss-=连续方程连续方程海水运动方程海水运动方程(2)(2)(3)(3)(4)(4)
11、11ppt课件温度方程温度方程湍流动能方程湍流动能方程HTKTDTUDTVDTTRFtxyDssss=-HSKSDSUDSVDSSFtxyDsss=盐度方程盐度方程22222qKq DUq DVq DqqtxyDsss=3221222()()MHqoKUVgDqKFDBlsss-22222qKq lDUq lDVq lDq lq ltxyDsss=322131()()MHloKUVgDqEEKWFDBsss-湍流混合长度方程湍流混合长度方程(5)(5)(6)(6)(7)(7)(8)(8)12ppt课件 dzzzzWUVdttxyDDDUVttxxyyhhhsss=o=其中，其中， 是实际水深是
12、实际水深DHh=(9)(9)13ppt课件,MHqKKK21( /)WE l kL= 分别分别表表示湍流混合系数示湍流混合系数，热热扩散系数和湍流动能扩散系数和湍流动能的垂直的垂直扩扩散系数散系数；墙墙近似近似函数函数 ，其中，其中 ，=0.4=0.4 是是vonvon KarmanKarman 常数常数湍流动能和混合长度方程组由下列等式关系闭合湍流动能和混合长度方程组由下列等式关系闭合 111()()LzHzh-=-21( /)WE lL= 0.2MMHHqKlqSKlqSKlq=其中，其中， 是稳定函数是稳定函数,MHSS(10)(10)14ppt课件在在2.52.5阶阶MYMY湍流闭合模
13、型中，湍流闭合模型中， 可表示为可表示为其中，其中， 。在不稳定层次。在不稳定层次 条件下，条件下， 的上限值为的上限值为0.0230.023，对应的，对应的 值分别为值分别为2.01452.0145和和2.44012.4401；在稳定层次；在稳定层次 条件下，它下限值为条件下，它下限值为-0.28-0.28，对应的，对应的 值分别为值分别为 0.04700.0470和和 0.04610.0461。 是是2.52.5阶阶M-YM-Y湍流闭合模式参数湍流闭合模式参数，由实验测得，由实验测得,MHSS0.42753.3541 34.6761 6.1270.4941 34.676hMhhHhGSGG
14、SG-=-=-2201hglGq Ds=10DshG10Ds1231,E E E B(11)(11),MHSS,MHSS15ppt课件M-YM-Y湍流模型已被广泛应用于浅海潮汐，风生湍流模型已被广泛应用于浅海潮汐，风生表面混合层的模拟以及底边界层的研究；表面混合层的模拟以及底边界层的研究；该模型在混合较弱的层化流体中对湍流混合系数计该模型在混合较弱的层化流体中对湍流混合系数计算效果不佳，应用于河口的可行性也有待于探讨；算效果不佳，应用于河口的可行性也有待于探讨；虽然该模型无论是在物理上还是在数学上都存在着虽然该模型无论是在物理上还是在数学上都存在着明显的缺陷，但是目前依然得到广泛应用。只有深明
15、显的缺陷，但是目前依然得到广泛应用。只有深入了解海水的微细结构以及海洋湍流结构和性质后，入了解海水的微细结构以及海洋湍流结构和性质后，才有可能构建出更为完善的湍流闭合模型才有可能构建出更为完善的湍流闭合模型 16ppt课件其中水平切应力其中水平切应力 是水平涡粘性系数，由是水平涡粘性系数，由Smagorinsky模式计算模式计算()()()()xxxxyyxyyyFHHxyFHHxy=2()2xxMxyyxMyyMUAxUVAyxVAy=MA(12)(12)(13)(13)17ppt课件Smagorinsky模式模式其中参数其中参数C C 为无量纲值，其取值为为无量纲值，其取值为0.10.1或
16、或0.20.2，如，如果计算网格划分得足够细的话，果计算网格划分得足够细的话，C C 可以取值为可以取值为0 0 由由SmagorinskySmagorinsky模式计算的水平涡粘性系数随着模式计算的水平涡粘性系数随着网格精度的改善和速度梯度的减小而减小网格精度的改善和速度梯度的减小而减小1/222211 = 22MuuvvAC x yxxyy (14)(14)18ppt课件其中其中 是水平热扩散系数，一般地，是水平热扩散系数，一般地， 是一个小值，是一个小值，取取0.10.1或或0.20.2，甚至在某些情况下可以取为，甚至在某些情况下可以取为0 0HA22()(), , ,xyFHqHqT
17、S q q lxy=22xHyHqAxqAy=HMAA(15)(15)(16)(16)19ppt课件垂向边界条件垂向边界条件1/2221010101022212/3,(,)(0)0, (0) ,(0)(0),(0) = (0), 0MasHKuvCUVUVDKTSwwsDqqBussss= -= - 1. 1. 海面海面0s(17)(17)20ppt课件2. 2. 海底海底1s -1/222222/321,(,)( 1)0,0( 1), (1) = ( 1), 0MbbbbbHKuvCuvu vDKTSDqqBussss=-=-其中，其中， 分别是湍流闭合模型参数和摩擦速度；分别是湍流闭合模型
18、参数和摩擦速度；1,B u,sbC C分别是海面和海底的摩擦系数；分别是海面和海底的摩擦系数；(18)(18)221 = max, 0.0025ln1+/bkboCHzs-21ppt课件控制近岸环流的海水运动方程包含了运动速度较快的控制近岸环流的海水运动方程包含了运动速度较快的外重力波和速度较慢的内重力波，因此进行模式分裂外重力波和速度较慢的内重力波，因此进行模式分裂可以大大提高模型的计算效率可以大大提高模型的计算效率模型的外模态就是计算自由水位和垂向平均的流速变模型的外模态就是计算自由水位和垂向平均的流速变化，需要较小的时间步长化，需要较小的时间步长内模态主要模拟流速矢量、温度、盐度等三维结
19、构变内模态主要模拟流速矢量、温度、盐度等三维结构变化，对时间步长的要求相对宽松化，对时间步长的要求相对宽松采用模式分裂技巧可以有效地减少因外模态计算所消采用模式分裂技巧可以有效地减少因外模态计算所消耗的计算时间，从而提高计算效率耗的计算时间，从而提高计算效率22ppt课件垂向积分的海水运动方程垂向积分的海水运动方程0DUDVxyth=2+xUDU DUVDFfVDgDwuwutxyxh-= -1+ ooxogDDGDddxxsssss- 2+yVDUVDV DFfUDgDwvwvtxyyh-= -1+ ooyogDDGDddyysssss- -1-1, ooUUdVVdss=其中：其中：连续方
20、程连续方程海水运动方程海水运动方程(19)(19)(20)(20)(21)(21)23ppt课件= 2+ + = + +2xMMyMMUUVFH AHAxxyyxUVVFHAH Axyxyy2222xxxyyyU DUVDU DUVDGFFxyxyUVDV DUVDV DGFFxyxy=-=-(22)(22)(23)(23)24ppt课件im, jm, kb, imm1, jmm1, kbm1, mode, isplitdte, dti, days, umolz, zz, dz, dzzaam2d, art, aru, arv, dum, dvm, fsmdx, dy, h, el, d, u
21、a, va, ut, vt, et, corswradwusurf, wvsurf, wubot, wvbot, wtsurf, wssurf, radu, v, , t, s, rho, l, q2km, kh, kq, aam, aahrmean, tclim, sclim25ppt课件 START9000 IINT =1,IENDPrint STOPSet Parmeters, Initial ValuesSet Parameters Initial Values BAROPG8000 Adjust integral of U,V to match UT, VTVERTVL BCOND(
22、5)ADVQ(Q2) ADVQ(Q2L) PROFQ BCOND(6) ADVT(T) ADVT(S) PROFT(T) PROFT(S) BCOND(4) ADVU ADVV PROFU PROFV BCOND(3) 9000 Compute EL BCOND(1) ADVAVECompute UA,VACompute UT,VT for use in Internal ModeBCOND(2)8000Figure 2. Flow diagram of the code. The boxes with sidebars contain subroutines. IEXT = 1,ISPLIT
23、ADVCT26ppt课件Time n-1tntn+1tDTEExternal ModeETBETETFUTB VTBUTF VTFooooFeedbackDTIInternal Mode27ppt课件三维变量的计算（三维变量的计算（以以T为例为例）可以分解为垂直扩散项的计）可以分解为垂直扩散项的计算和水平对流扩散项的计算，以温度方程为例方程写为算和水平对流扩散项的计算，以温度方程为例方程写为1+( )( ) =HDTTRAdv TDif TKtDsss-求解过程分为两步，第一步计算求解过程分为两步，第一步计算水平对流扩散项水平对流扩散项，采用，采用中心差分中心差分格式，由格式，由advt子程序
24、实现子程序实现1111 =() + ()2nnnnnDT DTAdv TDif Tt-(24)(24)(25)(25)水平对流扩散项水平对流扩散项垂直扩散项垂直扩散项28ppt课件第二步计算第二步计算垂向扩散项垂向扩散项，数值方法采用，数值方法采用“蛙跳蛙跳”格式格式，由由proft子程序实现子程序实现1111+11= 2 tnnn+nHnDT DTTRKDsss-由于由于“蛙跳蛙跳”格式在奇数时间步时的解与在偶数时间步时格式在奇数时间步时的解与在偶数时间步时的解会发生偏离，因此每一时间步的计算结束后还需对上的解会发生偏离，因此每一时间步的计算结束后还需对上述解进行平滑处理，即述解进行平滑处理
25、，即11 = + 2 + 2nnnnsTTTTT-为一个小值，一般取为一个小值，一般取0.05。处理后。处理后，11 ,nnnsTTTT-(26)(26)(27)(27)29ppt课件VA(I,J+1)UA(I,J)UA(I+1,J)VA(I,J)yxh(I,J) V(I,J+1)U(I,J,K)U(I+1,J,K)V(I,J)yxT(I,J,K) Q(I,J,K)plan view W(I,J,K) Q(I,J,K)U(I,J,K)U(I+1,J,K)sx T(I,J,K)elevation view W(I,J,K+1) Q(I,J,K+1)Z(K)ZZ(K)Z(K+1)二维外模态网格分布
26、二维外模态网格分布三维内模态三维内模态网格分布网格分布30ppt课件 虽然模型采用虽然模型采用有限差分有限差分计算格式，但是对于每个网格计算格式，但是对于每个网格对流项的计算都是按照对流项的计算都是按照有限容积有限容积的方法进行处理，即的方法进行处理，即温度对流过程可表示为温度对流过程可表示为 速度的对流过程与温度相似，可写为速度的对流过程与温度相似，可写为 ()( ) = + + xyxyyxxyTAdv T h hDh UTDh VTh hss-()( ) = + () + xyxyyxxyxyUAdv U h hDh UUDh UVh hfVDh hss-其中，其中， 是由是由z z坐标
27、转化为坐标转化为坐标后坐标后产生的弯曲项产生的弯曲项 ()() = xyyxxyxyVhUhfh hh h-(28)(28)(29)(29)31ppt课件 垂直扩散项的计算公式（垂直扩散项的计算公式（第第k k层，层，1kkb-11kkb-1） 可以改写为可以改写为1= 2 tkkkHDT DTTRKDsss-111212 t 1= HHkkkkkkkkkkkKKT TTTTTDdzdzzdzz-121radradkkktD dz-(30)(30)(31)(31)32ppt课件 展开后合并同类项得到展开后合并同类项得到 其中其中1kkk1kka + a + c1c = dkkkTTT-1k2k
28、21k12a2c2d =radradHkkkHkkkkkkktKDdzdzztKDdzdzztTD dz- = - = -(32)(32)(33)(33)33ppt课件 现在假定温度解的形式为现在假定温度解的形式为 由由(34)(34)得到的得到的 代入到代入到(32)(32)就可以得到就可以得到其中其中 的值由的值由(33)(33)求得，求得， 的值由上的值由上一层的一层的 值确定值确定k1k ee +ggkkTT=(34)(34)1kT-kkkkk 1kk-1kkkkk 1aee = a + c1ee1cggdgg = a + c1ee1-(35)(35)kkka , c ,dkkee ,
29、 ggk 1k 1ee, gg-34ppt课件 当当k=1k=1时时，即表层海水，海水温度主要取决于海表面，即表层海水，海水温度主要取决于海表面温度通量。公式（温度通量。公式（3434）的解可以近似写为）的解可以近似写为 这样，表层海水的温度只能根据式（这样，表层海水的温度只能根据式（3131）求解，即）求解，即 上式可以进一步表示为上式可以进一步表示为(36)(36)111ee =0, ggT=12110112121012= radradHHKKtT TTTTTDdzDdzzDdzz-wtsurf-2111212211122=wtsurfradradHt KtT TTTDdzDdzdzz-1
30、a-(37)(37)35ppt课件 式（式（3737）还可以进一步写为）还可以进一步写为 再与温度的通解（再与温度的通解（3434）比较）比较 可以得到可以得到111 212112a1=awtsurfradradtTT TDdz-(34)(34)1121 ee +ggTT=(38)(38)11111211aeea112ggwtsurfradrada1t TDdz=-=-(39)(39)36ppt课件 当短波辐射通量当短波辐射通量 的计算公式的计算公式 r r，ad1ad1，ad2ad2 均为光辐射常数，根据不同水质取值如下均为光辐射常数，根据不同水质取值如下(40)(40)kradkkD zD
31、 zad1ad2krad = swradre(1r)e-ntp 1 2 3 4 5Jerlov type I Ia Ib II IIIr 0.58 0.62 0.67 0.70 0.78ad1 (m) 0.35 0.60 1.0 1.5 1.4ad2 (m) 23.0 20.0 17.0 14.0 7.9来源：来源：Jerlov，1976；Paulson and Simpson，197737ppt课件 当当k=kb-1k=kb-1时时，即底层海水，假定海底的热通量为，即底层海水，假定海底的热通量为0 0，根据前面的推导方式同样可以得到底层海水温度的计根据前面的推导方式同样可以得到底层海水温度的
32、计算公式算公式 对于盐度方程垂向扩散项的计算与温度方程相同，对于盐度方程垂向扩散项的计算与温度方程相同，只是不考虑太阳短波辐射这一项只是不考虑太阳短波辐射这一项(41)(41)kb 1kb 1kb 2121kb-1kb 1kb 22cgg radrad = c(1ee)1tTDdzT-38ppt课件 由于模型水平方向采用显格式，因此时间步长的选取由于模型水平方向采用显格式，因此时间步长的选取必须要满足必须要满足CFLCFL稳定条件。稳定条件。 对于外模态，时间步长的限制条件为对于外模态，时间步长的限制条件为 其中其中 ， 可能预见到的最大流速可能预见到的最大流速1/222111 +tEtCxy
33、-max2tCgHU=maxU 对于内模态，时间步长的限制条件较外模态的情形对于内模态，时间步长的限制条件较外模态的情形宽松很多，主要是速度较快的外重力波已经在外模态宽松很多，主要是速度较快的外重力波已经在外模态中考虑了。一般中考虑了。一般 取值取值30305050即可即可 /EItt(42)(42)39ppt课件 对于动量或标量还有其它一些时间限制条件对于动量或标量还有其它一些时间限制条件 其中其中 或或122111 +4ItAxy- 以及以及 其中其中 分别为地球角速率和地理纬度分别为地球角速率和地理纬度 11 =2sinItfMAA=HAA=, (43)(43)(44)(44)40ppt
34、课件陆地及岸线是由陆地及岸线是由dumdum、dvmdvm、fsmfsm控制的，在陆地上控制的，在陆地上这些变量的值设为这些变量的值设为0 0，有水的地方设为，有水的地方设为1 1子程序子程序bcondbcond（idxidx），），idxidx=1=1对应的水位边界条件；对应的水位边界条件；idxidx=2=2对应垂向平均流速边界条件；对应垂向平均流速边界条件；idxidx=3=3对应三对应三维水平流速边界条件；维水平流速边界条件；idxidx=4=4对应温盐边界条件；对应温盐边界条件；idxidx=5=5对应垂向流速边界条件；对应垂向流速边界条件；idxidx=6=6对应湍流动对应湍流动能
35、和湍流混合长度边界条件能和湍流混合长度边界条件 模型对开边界条件的要求很高，而开边界条件本身模型对开边界条件的要求很高，而开边界条件本身又具有很大的不确定性，因此有必要对模型的外模又具有很大的不确定性，因此有必要对模型的外模态和内模态的开边界分别进行处理态和内模态的开边界分别进行处理 41ppt课件 START9000 IINT =1,IENDPrint STOPSet Parmeters, Initial ValuesSet Parameters Initial Values BAROPG8000 Adjust integral of U,V to match UT, VTVERTVL BC
36、OND(5)ADVQ(Q2) ADVQ(Q2L) PROFQ BCOND(6) ADVT(T) ADVT(S) PROFT(T) PROFT(S) BCOND(4) ADVU ADVV PROFU PROFV BCOND(3) 9000 Compute EL BCOND(1) ADVAVECompute UA,VACompute UT,VT for use in Internal ModeBCOND(2)8000Figure 2. Flow diagram of the code. The boxes with sidebars contain subroutines. IEXT = 1,IS
37、PLITADVCT42ppt课件FormulaBoundaryCode Inflow condition 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
38、) + 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) = set Elevation condition h = 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(
39、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) = setDU = BC43ppt课件FormulaBoundaryCodeRadiationEASTuaf(im,j) = sqrt(grav/h(imm1,j)* el(imm1,j) + bc(j)elf(im,j) = elf(imm1,j)vaf(im,j) = se
40、tWESTuaf(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)/
41、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)
42、 = 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) = setHU ceh= BC0eUUctx=ce=gH流速水位流速水位流速流速44ppt课件FormulaBoundaryCodeRadiationEASTgae = dte*sqrt(grav*h(imm1,j)/dx(imm1,j)elf(imm1,j) = gae*el(imm2,j) + (1.-gae) *el(
43、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)
44、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) = set Cyclic EAST(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(im
45、m1,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)0ectxhh=水位水位45ppt课件FormulaBoundaryCodeInflow conditionEASTuf(im,j,k) = bc(j,k)vf(im,j,k) = setU = BC WESTuf(
46、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) = s
47、etNORTHgai = 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) = set0iUUctx=46ppt课件FormulaBoundaryCodeUpstream advection on T or SEASTuf(im,j,k) = t(im,j,k) -dti/(dx(im,j)+dx(imm1,j) *
48、(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) NORTHuf(i,jm,k) = t(i,jm,k) -dti/(dy(i,jm)+dy(i,jmm1) *
49、 (v(i,jm,k) + abs(v(i,jm,k) * (t(i,jm,k)-t(i,jmm1,k) + (v(i,jm,k) - abs(v(i,jm,k) * (tbn(i,k)-t(i,jm,k)SOUTHuf(i,1,k) = t(i,1,k) -dti/(dy(i,1)+dy(i,2) * (v(i,1,k) + abs(v(i,1,k) * (t(i,1,k)-tbs(i,k) + (v(i,1,k) - abs(v(i,1,k) * (t(i,2,k)-t(i,1,k)CyclicMuch the same as the external case except replac
50、e uaf with uf, etc. and t, s, q2 and q2l are handled similar to elf0TTUtx=47ppt课件48ppt课件渤海及入海河流渤海及入海河流0100200300050100Time (d)Runoff (108 m3) YRHRDR49ppt课件河流径流输入模型河流径流输入模型将河流视为点源，而点源的模拟采用的是一个水将河流视为点源，而点源的模拟采用的是一个水平网格的面积平网格的面积 径流动量对动量方程的贡献可以忽略不计径流动量对动量方程的贡献可以忽略不计受河流径流影响的二维连续方程受河流径流影响的二维连续方程 ( , , )Du
