海洋数值模型的理论及应用-课件

第一第一节海洋数海洋数值模模型的型的发展概况展概况观测，试验观测，试验理论分析理论分析数值模型数值模型海洋的研究手段海洋的研究手段海洋数海洋数值模型模型发展的客展的客观要求要求海水运动方程推导过程中的动力学假设，比如海海水运动方程推导过程中的动力学假设，比如海水静力学假设，都是建立在大部分空间尺度大于水静力学假设，都是建立在大部分空间尺度大于10km10km，时间尺度大于惯性周期（，时间尺度大于惯性周期（2 2/f/f）的情况）的情况下的下的海水运动方程不存在精确解：（海水运动方程不存在精确解：（1 1）偏微分方程）偏微分方程本身是非线性的，（本身是非线性的，（2 2）将一个实际问题具体化）将一个实际问题具体化时所需要的环境场如海底地形、海岸线几何形状、时所需要的环境场如海底地形、海岸线几何形状、表面驱动力等都无法给出解析函数表达式，只能表面驱动力等都无法给出解析函数表达式，只能在离散化的空间和时间区域取得在离散化的空间和时间区域取得因此，求解海水运动方程只能通过近似求解的办因此，求解海水运动方程只能通过近似求解的办法来确定，而这些数值解又可能引入明显的近似法来确定，而这些数值解又可能引入明显的近似误差误差海洋数海洋数值模型模型发展的展的历史史海洋模型到按其水平网格的离散方式以及所使用海洋模型到按其水平网格的离散方式以及所使用的垂向坐标系的不同大致经历了如下几个发展阶的垂向坐标系的不同大致经历了如下几个发展阶段段最早最早出现并且还在使用的海洋模型是出现并且还在使用的海洋模型是BryanBryan等人等人开发的基于原始方程的低阶精度的有限差分模型，开发的基于原始方程的低阶精度的有限差分模型，它在水深方向采用它在水深方向采用z z坐标系坐标系目前常见的目前常见的HOPS（Harvard Ocean Prediction System），MOM（GFDL Modular Ocean Model），POP(Parallel Ocean Program)，NCOM（NCAR Community Ocean Model）模模型在某种程度上都可以认为是该模型改进版型在某种程度上都可以认为是该模型改进版海洋数海洋数值模型模型发展的展的历史史上个世纪上个世纪7070年代年代，sigma坐坐标系开始系开始应用于用于海洋模型在水深方向，比如目前被广泛使用海洋模型在水深方向，比如目前被广泛使用的的POM（Princeton Ocean Model）、ECOM（Estuarine Coastal and Ocean Model）、ROMS（Regional Ocean Modeling System）模型都属于模型都属于这种种类型的型的模型模型与传统的与传统的 z z 坐标系相比，坐标系相比，sigmasigma坐标可以更坐标可以更好地贴合海底地形变化好地贴合海底地形变化海洋数海洋数值模型模型发展的展的历史史自上世纪自上世纪9090年代以来年代以来，非，非结构化网格技构化网格技术开始开始在海洋模型中得到在海洋模型中得到应用，相用，相应的也出的也出现了基于了基于非非结构化网格离散方式的有限元模型构化网格离散方式的有限元模型，如如SEOM模型模型（Spectral Finite Element Ocean Model），和有限体），和有限体积模型模型，如如FVCOM虽然与传统的结构化网格相比虽然与传统的结构化网格相比,非结构化网格非结构化网格可以更好地拟合陆地边界，但是代码实现上的可以更好地拟合陆地边界，但是代码实现上的困难以及计算稳定性的问题使其迄今还没有得困难以及计算稳定性的问题使其迄今还没有得到非常广泛到非常广泛的的应用应用新一代的海洋数新一代的海洋数值模型模型新一代的海洋模型广泛采用随地坐新一代的海洋模型广泛采用随地坐标系（系（terrain-following coordinates），），进而促而促进了有了有关关时间步步长，对流流项和和压力梯度力梯度项等数等数值算法算法的改的改进进入新世纪以来进入新世纪以来，下一代的海洋数值动力模型，下一代的海洋数值动力模型正在紧锣密鼓的研制中，代表性的是正在紧锣密鼓的研制中，代表性的是TOMS(Terrain-following Ocean Modeling System)，它，它融合了目前最先进的物理知识、数值方法和数融合了目前最先进的物理知识、数值方法和数据同化技术据同化技术 海洋模型的商海洋模型的商业软件件就海流的仿真建模研究而言，现阶段国外已形就海流的仿真建模研究而言，现阶段国外已形成了不少具有代表性的、具有强大前后处理功成了不少具有代表性的、具有强大前后处理功能及核心仿真建模技术的软件或系统能及核心仿真建模技术的软件或系统（商业软（商业软件）件）如加拿大如加拿大Dalhousie University研究研究的的DALCOAST河口海岸河口海岸预报系系统；丹麦的丹麦的DHI Water&Environment 机构开机构开发的的的的MIKE系系列列软件系件系统；荷荷兰的的WL Delft Hydraulics开开发的的的的Delft3D软件件第二第二节基于基于POM模型的一般介模型的一般介绍水水动力模型基本特点力模型基本特点n垂直方向的垂直方向的 坐标坐标变换变换n水平网格采用正交曲线坐标和水平网格采用正交曲线坐标和ArakawaCArakawaC差分格式差分格式n水平时间差分采用显格式，而垂直差分为隐格式水平时间差分采用显格式，而垂直差分为隐格式n自由表面可以模拟水位变化自由表面可以模拟水位变化n垂直和水平方向的混合扩散分别采用垂直和水平方向的混合扩散分别采用2.52.5阶的阶的Mellor-YamadaMellor-Yamada湍流闭合模式和湍流闭合模式和SmagorinskiSmagorinski模式模式n内外内外模态分别处理速度较慢的内重力波和速度较模态分别处理速度较慢的内重力波和速度较快的外重力波以提高整个模式快的外重力波以提高整个模式计算效率计算效率n包含了海水的热动力过程包含了海水的热动力过程Sigma坐坐标变换：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 a）垂向）垂向sigmasigma坐标系和（坐标系和（b b）水平）水平Arakawa CArakawa C网格网格(1)(1)控制方程控制方程连续方程连续方程海水运动方程海水运动方程(2)(2)(3)(3)(4)(4)温度方程温度方程湍流动能方程湍流动能方程盐度方程盐度方程湍流混合长度方程湍流混合长度方程(5)(5)(6)(6)(7)(7)(8)(8)其中，其中，是基于是基于sigma坐坐标系的垂向流速矢量，系的垂向流速矢量，它与笛卡儿坐它与笛卡儿坐标系下的垂向流速系下的垂向流速W之之间的关系的关系可表示可表示为另外，方程中另外，方程中另外，方程中另外，方程中 分别表示海水的实际密度，分别表示海水的实际密度，分别表示海水的实际密度，分别表示海水的实际密度，BoussinesqBoussinesqBoussinesqBoussinesq 近似密度以及密度扰动近似密度以及密度扰动近似密度以及密度扰动近似密度以及密度扰动 其中，其中，是实际水深是实际水深(9)(9)MY湍流湍流闭合模型（合模型（profq）在湍流动能及混合长度方程中在湍流动能及混合长度方程中在湍流动能及混合长度方程中在湍流动能及混合长度方程中 分别分别表表示湍流混合系数示湍流混合系数，热热扩散系数和湍流动能扩散系数和湍流动能的垂直的垂直扩扩散系数散系数；墙墙近似近似函数函数 ，其中，其中 ，=0.4=0.4 是是vonvon KarmanKarman 常数常数湍流动能和混合长度方程组由下列等式关系闭合湍流动能和混合长度方程组由下列等式关系闭合 其中，其中，是稳定函数是稳定函数(10)(10)在在2.52.5阶阶MYMY湍流闭合模型中，湍流闭合模型中，可表示为可表示为其中，其中，。在不稳定层次。在不稳定层次 条件下，条件下，的上限值为的上限值为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湍流闭合模式参数湍流闭合模式参数，由实验测得，由实验测得(11)(11)MY湍流模型的适用条件湍流模型的适用条件M-YM-Y湍流模型已被广泛应用于浅海潮汐，风生表湍流模型已被广泛应用于浅海潮汐，风生表面混合层的模拟以及底边界层的研究；面混合层的模拟以及底边界层的研究；该模型在混合较弱的层化流体中对湍流混合系该模型在混合较弱的层化流体中对湍流混合系数计算效果不佳，应用于河口的可行性也有待数计算效果不佳，应用于河口的可行性也有待于探讨；于探讨；虽然该模型无论是在物理上还是在数学上都存虽然该模型无论是在物理上还是在数学上都存在着明显的缺陷，但是目前依然得到广泛应用。在着明显的缺陷，但是目前依然得到广泛应用。只有深入了解海水的微细结构以及海洋湍流结只有深入了解海水的微细结构以及海洋湍流结构和性质后，才有可能构建出更为完善的湍流构和性质后，才有可能构建出更为完善的湍流闭合模型闭合模型 海水运海水运动方程的水平方程的水平扩散散项其中水平切应力其中水平切应力 是水平涡粘性系数，由是水平涡粘性系数，由Smagorinsky模式模式计算算(12)(12)(13)(13)Smagorinsky水平水平扩散散模式模式其中参数其中参数C C 为无量纲值，其取值为为无量纲值，其取值为0.10.1或或0.20.2，如果计算网格划分得足够细的话，如果计算网格划分得足够细的话，C C 可以可以取值为取值为0 0 由由SmagorinskySmagorinsky模式计算的水平涡粘性系数随模式计算的水平涡粘性系数随着网格精度的改善和速度梯度的减小而减小着网格精度的改善和速度梯度的减小而减小(14)(14)温度温度盐度的水平度的水平扩散散项其中其中 是水平热扩散系数，一般地，是水平热扩散系数，一般地，是一个小是一个小值，取值，取0.10.1或或0.20.2，甚至在某些情况下可以取为，甚至在某些情况下可以取为0 0(15)(15)(16)(16)垂向垂向边界条件界条件1.1.海面海面(17)(17)2.2.海底海底其中，其中，分别是湍流闭合模型参数和摩擦速度；分别是湍流闭合模型参数和摩擦速度；分别是海面和海底的摩擦系数；分别是海面和海底的摩擦系数；(18)(18)模型的外模模型的外模态控制近岸环流的海水运动方程包含了运动速度较快的控制近岸环流的海水运动方程包含了运动速度较快的外重力波和速度较慢的内重力波，因此进行模式分裂外重力波和速度较慢的内重力波，因此进行模式分裂可以大大提高模型的计算效率可以大大提高模型的计算效率模型的外模态就是计算自由水位和垂向平均的流速变模型的外模态就是计算自由水位和垂向平均的流速变化，需要较小的时间步长化，需要较小的时间步长内模态主要模拟流速矢量、温度、盐度等三维结构变内模态主要模拟流速矢量、温度、盐度等三维结构变化，对时间步长的要求相对宽松化，对时间步长的要求相对宽松采用模式分裂技巧可以有效地减少因外模态计算所消采用模式分裂技巧可以有效地减少因外模态计算所消耗的计算时间，从而提高计算效率耗的计算时间，从而提高计算效率垂向垂向积分的海水运分的海水运动方程方程其中：其中：连续方程连续方程海水运动方程海水运动方程(19)(19)(20)(20)(21)(21)水平水平扩散散项水平数水平数值耗散耗散项(22)(22)(23)(23)模型代模型代码中的一些符号中的一些符号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数数值积分分流流程程图外模外模态和内模和内模态的相互作用的相互作用内模内模态的的计算方法算方法三维变量的计算（三维变量的计算（以以T为例为例）可以分解为垂直扩散项的计）可以分解为垂直扩散项的计算和水平对流扩散项的计算，以温度方程为例方程写为算和水平对流扩散项的计算，以温度方程为例方程写为求解过程分为两步，第一步计算求解过程分为两步，第一步计算水平对流扩散项水平对流扩散项，采用，采用中心差分中心差分格式，由格式，由advt子程序实现子程序实现(24)(24)(25)(25)水平对流扩散项水平对流扩散项垂直扩散项垂直扩散项第二步计算第二步计算垂向扩散项垂向扩散项，数值方法采用，数值方法采用“蛙跳蛙跳”格式格式，由由proft子程序实现子程序实现由于由于“蛙跳蛙跳”格式在奇数时间步时的解与在偶数时间步时格式在奇数时间步时的解与在偶数时间步时的解会发生偏离，因此每一时间步的计算结束后还需对上的解会发生偏离，因此每一时间步的计算结束后还需对上述解进行平滑处理，即述解进行平滑处理，即为一个小值，一般取为一个小值，一般取0.05。处理后。处理后，(26)(26)(27)(27)计算网格的安排算网格的安排二维外模态网格分布二维外模态网格分布三维内模态三维内模态网格分布网格分布水平水平对流流项的的计算算 虽然模型采用虽然模型采用有限差分有限差分计算格式，但是对于每个网格计算格式，但是对于每个网格对流项的计算都是按照对流项的计算都是按照有限容积有限容积的方法进行处理，即的方法进行处理，即温度对流过程可表示为温度对流过程可表示为 速度的对流过程与温度相似，可写为速度的对流过程与温度相似，可写为 其中，其中，是由是由z z坐标转化为坐标转化为坐标后坐标后产生的弯曲项产生的弯曲项 (28)(28)(29)(29)垂直垂直扩散散项的的计算算 垂直扩散项的计算公式（垂直扩散项的计算公式（第第k k层，层，1kkb-11kkb-1）可以改写为可以改写为(30)(30)(31)(31)展开后合并同类项得到展开后合并同类项得到 其中其中(32)(32)(33)(33)现在假定温度解的形式为现在假定温度解的形式为 由由(34)(34)得到的得到的 代入到代入到(32)(32)就可以得到就可以得到其中其中 的值由的值由(33)(33)求得，求得，的值由上的值由上一层的一层的 值确定值确定(34)(34)(35)(35)当当k=1k=1时时，即表层海水，海水温度主要取决于海表面，即表层海水，海水温度主要取决于海表面温度通量。公式（温度通量。公式（3434）的解可以近似写为）的解可以近似写为 这样，表层海水的温度只能根据式（这样，表层海水的温度只能根据式（3131）求解，即）求解，即 上式可以进一步表示为上式可以进一步表示为(36)(36)(37)(37)式（式（3737）还可以进一步写为）还可以进一步写为 再与温度的通解（再与温度的通解（3434）比较）比较 可以得到可以得到(34)(34)(38)(38)(39)(39)当短波辐射通量当短波辐射通量 的计算公式的计算公式 r r，ad1ad1，ad2ad2 均为光辐射常数，根据不同水质取值如下均为光辐射常数，根据不同水质取值如下(40)(40)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，1977 当当k=kb-1k=kb-1时时，即底层海水，假定海底的热通量为，即底层海水，假定海底的热通量为0 0，根据前面的推导方式同样可以得到底层海水温度的计根据前面的推导方式同样可以得到底层海水温度的计算公式算公式 对于盐度方程垂向扩散项的计算与温度方程相同，对于盐度方程垂向扩散项的计算与温度方程相同，只是不考虑太阳短波辐射这一项只是不考虑太阳短波辐射这一项(41)(41)时间步步长的的CFL限制条件限制条件 由于模型水平方向采用显格式，因此时间步长的选取由于模型水平方向采用显格式，因此时间步长的选取必须要满足必须要满足CFLCFL稳定条件。稳定条件。对于外模态，时间步长的限制条件为对于外模态，时间步长的限制条件为 其中其中 ，可能预见到的最大流速可能预见到的最大流速 对于内模态，时间步长的限制条件较外模态的情形对于内模态，时间步长的限制条件较外模态的情形宽松很多，主要是速度较快的外重力波已经在外模态宽松很多，主要是速度较快的外重力波已经在外模态中考虑了。一般中考虑了。一般 取值取值30305050即可即可 (42)(42)时间步步长的其它限制条件的其它限制条件 对于动量或标量还有其它一些时间限制条件对于动量或标量还有其它一些时间限制条件 其中其中 或或 以及以及 其中其中 分别为地球角速率和地理纬度分别为地球角速率和地理纬度 (43)(43)(44)(44)侧开开边界条件（界条件（bcond）陆地及岸线是由陆地及岸线是由dumdum、dvmdvm、fsmfsm控制的，在陆地控制的，在陆地上这些变量的值设为上这些变量的值设为0 0，有水的地方设为，有水的地方设为1 1子程序子程序bcondbcond（idxidx），），idx=1idx=1对应的水位边界条对应的水位边界条件；件；idx=2idx=2对应垂向平均流速边界条件；对应垂向平均流速边界条件；idx=3idx=3对应三维水平流速边界条件；对应三维水平流速边界条件；idx=4idx=4对应温盐边对应温盐边界条件；界条件；idx=5idx=5对应垂向流速边界条件；对应垂向流速边界条件；idx=6idx=6对应湍流动能和湍流混合长度边界条件对应湍流动能和湍流混合长度边界条件 模型对开边界条件的要求很高，而开边界条件模型对开边界条件的要求很高，而开边界条件本身又具有很大的不确定性，因此有必要对模本身又具有很大的不确定性，因此有必要对模型的外模态和内模态的开边界分别进行处理型的外模态和内模态的开边界分别进行处理 数数值积分分流流程程图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)+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(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外外模模态开开边界界条条件件一一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外外模模态开开边界界条条件件二二流速水位流速水位流速流速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)=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(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)外外模模态开开边界界条条件件三三水位水位FormulaBoundaryCodeInflow conditionEASTuf(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内模内模态开开边界条件（一）界条件（一）FormulaBoundaryCodeUpstream advection on T or SEASTuf(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)NORTHuf(i,jm,k)=t(i,jm,k)-dti/(dy(i,jm)+dy(i,jmm1)*(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 replace uaf with uf,etc.and t,s,q2 and q2l are handled similar to elf内模内模态开开边界条件（二）界条件（二）第三第三节POM模型的模型的实际应用用渤海河流羽状流的数值模拟渤海河流羽状流的数值模拟渤海及入海河流渤海及入海河流河流径流河流径流输入模型入模型将河流视为点源，而点源的模拟采用的是一个水将河流视为点源，而点源的模拟采用的是一个水平网格的面积平网格的面积 径流动量对动量方程的贡献可以忽略不计径流动量对动量方程的贡献可以忽略不计受河流径流影响的二维连续方程受河流径流影响的二维连续方程 三维连续方程可表示为三维连续方程可表示为 内外模式相容内外模式相容30天时纯浮力羽流条件下渤海的盐度分布，等值线间隔为天时纯浮力羽流条件下渤海的盐度分布，等值线间隔为0.5 psu，最外层等值线为最外层等值线为30 psu。（。（a）海水表层；（）海水表层；（b）海水底层）海水底层 30天时潮致羽流条件下渤海的盐度分布，等值线间隔为天时潮致羽流条件下渤海的盐度分布，等值线间隔为0.5 psu，最外层等值线为最外层等值线为30 psu。（。（a）海水表层；（）海水表层；（b）海水底层）海水底层 东风东风南风南风风场加入风场加入1010天后渤海表层海水流场天后渤海表层海水流场东风东风南风南风风场加入风场加入10天后渤海表层盐度的水平分布状况天后渤海表层盐度的水平分布状况。等值。等值线间隔为线间隔为0.5 psu。海河口附近站点（海河口附近站点（1171174545E E，38 38 5555N N）盐度变）盐度变化的时间序列化的时间序列海水表层海水表层海水底层海水底层7 7月份渤海表层盐度的水平分布状况，等值线间隔为月份渤海表层盐度的水平分布状况，等值线间隔为1.0 psu1.0 psu。计算结果。计算结果(左图）；实测结果（右图）左图）；实测结果（右图）