二维潮波数值模拟

合集下载

长江感潮河段平面二维潮流数值模拟

长江感潮河段平面二维潮流数值模拟曾小辉;李国杰;姜昱【摘要】基于正交曲线网格下的有限体积法,建立河口潮汐的水深平均二维模型,运用Simple算法求解河口潮汐的二维浅水方程.以长江下游江阴水道为计算实例,用“露滩冻结”动边界技术处理方法模拟江阴水道二维潮流场,较好地解决了边滩、浅滩边界随水位的变动问题.实例和验证结果表明:计算的潮位和潮流过程与实测过程吻合较好,计算的流场合理.对计算结果的前后可视化处理,采用AutoCAD的DXF 数据接口生成流场和网格图,以便更好地进行直观分析.【期刊名称】《水运工程》【年(卷),期】2012(000)004【总页数】5页(P12-16)【关键词】二维潮流场;数值模拟;SIMPLER算法;动边界技术【作者】曾小辉;李国杰;姜昱【作者单位】中交第二航务工程勘察设计院有限公司,湖北武汉430071;中交第二航务工程勘察设计院有限公司,湖北武汉430071;中交第二航务工程勘察设计院有限公司,湖北武汉430071【正文语种】中文【中图分类】TP391有限差分法的原理比较简单，若边界不复杂，数学推理和程序编制都比较简单，并且计算内存占用比较少，计算速度快，在非恒定性比较强的问题中应用比较多。

但有限差分方法在处理复杂边界时灵活性比较差，计算精度不高，在复杂边界上常常出现虚假流动的现象，有时会出现水量不守恒现象。

相对而言，有限元法数学推理和程序编制复杂些，且因这种方法的网格可以为三角形或四边形等任意网格，虽然这种网格结构易于处理边界和地形比较复杂的问题，但需要数据结构定义计算单元之间的位置关系，占用内存比较大，计算速度相对较慢。

有限体积法是将计算区域划分成若干不规则形状的单元或控制体，对每个控制体分别进行水量和动量平衡计算，在计算出通过每个控制体边界沿法向输入（出）的流量和动量后，便得到计算时段末各控制体平均水深和流速。

有限体积法与有限单元法和有限差分法的数值逼近相比，其物理意义更直接、更明晰。

浙江近海潮汐潮流的数值模拟-海洋学报

第25卷第5期海洋学报Vol.25,N o.5 2003年9月ACTA OCEANOLOGICA SINICA September,2003浙江近海潮汐潮流的数值模拟陈倩1,2,3,黄大吉2,章本照1(11浙江大学力学系,浙江杭州310027;21国家海洋局第二海洋研究所浙江杭州310012;31国家海洋局海洋动力过程与卫星海洋学重点实验室,浙江杭州310012)摘要:用三维陆架海模式(HAM SOM)对浙江近海的潮汐、潮流进行了数值模拟,并采用网格嵌套和动边界技术对原模式作了改进,以提高计算的精度,改进后的模式在浙江近海的应用中被证明是成功的.沿岸50个潮位站计算与实测值的比较表明,加入动边界以后的小区域细网格计算较之粗网格以及未加动边界以前精度普遍提高,比较的均方差结果为:M2分潮振幅差416cm,相角差7114b;S2分潮振幅差510cm,相角差514b;K1分潮振幅差2125cm,相角差5176b;O1分潮振幅差1156cm,相角差515b,可见计算与实测符合良好.另外,选取了105个实测潮流点,比较了表层M2和K1分潮流调和常数分量U cos N,U sin N,V cos G,V sin G的实测值与计算值的偏差,结果表明计算与实测的符合程度较好.在此基础上,给出了各主要分潮的潮位同潮图、潮流同潮图、潮汐性质、潮流性质、潮流椭圆和潮流的运动形式等,发现4个主要分潮M2,S2,K1,O1在本区内均未出现无潮点;M2分潮流在29b18c N,122b46c E处有一个圆流点.此外还得到了一些有意义的结论,都与实测情况符合良好,从而对整个浙江沿海区域的潮汐潮流特性有了一个全面认识.关键词:数值模拟;网格嵌套;动边界;潮汐潮流;浙江近海中图分类号:P731123文献标识码:A文章编号:0253-4193(2003)05-0009-121引言浙江近海岸线曲折,地形复杂,港湾众多、岛屿星罗棋布.由北往南的主要港湾有杭州湾、宁波-舟山深水港、象山港、三门湾、乐清湾、温州湾等,均属强潮海区.沿岸和诸岛屿上设有一些验潮站和潮流测点.关于浙江沿岸各海区内潮汐潮流实测资料的研究已有不少成果,本文作者[1]也曾以多年来沿岸各潮位站观测资料以及海岸带和海岛调查的实测海流资料为依收稿日期:2002-09-25;修订日期:2002-12-10.基金项目:浙江省自然科学基金资助项目(499018);国家自然科学基金资助项目(40076010).作者简介:陈倩(1975)),女,浙江省舟山市人,国家海洋局第二研究所和浙江大学力学系联合培养博士,从事海洋动力过程研究.E-mail:dajih2001@10海洋学报25卷据,分析研究了整个浙江沿海区域的潮汐潮流特征,得到了一些规律性结论,但是实际观测既昂贵又费时,而且实测点有限,该海域的流场分布又较为复杂,因此要由点到面准确地掌握大面积潮位和潮流的分布变化规律,数值模拟计算仍不失为一种廉价而有效的研究手段.有关浙江近海各海区的数值模拟,前人已做过不少工作,曹德明等[2,3]用有限差分法对杭州湾的潮汐、潮流进行了二维数值模拟;李身铎等[4]采用R坐标下的三维数值模式模拟了杭州湾三维潮波运动;曹欣中等[5]做了宁波、舟山内海域的潮流场数值模拟计算;董礼先等[6]数值模拟了象山港水域的潮波运动;周大成[7]采用平面二维三角形单元显性有限元潮流模型对椒江河口大潮潮流特性进行了数值模拟;李孟国等[8]建立了时间二次插值的三角形网格显式差分数学模型,对瓯江口海区的潮流场进行了成功的模拟.但是,由于该区地形和岸线条件复杂,为提高数值计算的精度,绝大多数的计算工作是针对某局部海区进行的,极少把浙江近海作为一个整体进行数值计算,而且大多数的研究成果局限于某几点的潮位验证或某时刻的流场分布,未给出整个研究海域的潮分布规律(如各分潮的等振幅线和同潮时线等).也有许多模拟东海或渤海、黄海、东海的数值计算文章[9~11],但由于浙江近海只是其中很小一块区域,故其所提供的该区的潮汐潮流分布信息就不够细致.因此,把浙江近海作为一个整体,在较高精度下进行数值计算,并由此得到一系列规律性成果,是很有意义的,而且为了弄清潮汐潮流对物质输送和扩散以及浙江近海温度、盐度分布的影响,也需要把浙江近海作为一个整体来进行计算.本文采用三维陆架海模式对整个浙江海区的潮汐、潮流进行了数值模拟计算,并引入网格嵌套技术使得小区域内的细网格计算精度提高,引入动边界技术来处理露滩问题,使该区潮模拟更加真实可靠.本文用这两项技术对原有的HAMSOM模式进行了改进,改进后的模式被应用于浙江近海潮汐、潮流的三维数值计算中,经验证效果良好.在验证计算值与实测值符合良好的基础上,给出了各主要分潮的潮位等振幅线和同潮时线、潮汐性质、潮流等振幅线和同潮时线、潮流椭圆、潮流性质及潮流的运动形式.2数值模式及其应用211三维陆架海模式(HAMSOM)简介HAMSOM(H amburg shelf ocean model)是由德国汉堡大学海洋研究所Backhaus和他的同事们发展的三维斜压陆架海模式.自20世纪80年代初发展至今,模式已有了不少改进[12~15].下面对HAM SOM的框架及其特点作简要说明.HAMSOM是一垂向分层模式,控制方程建立在任一垂向层上.这是为简化数值计算,通过对原始三维运动方程和连续方程组进行层积分处理(即把相应的方程对模式的垂向某一层积分),得到层积分的连续方程、运动方程和温度、盐度方程[13],从而将三维问题转化为二维问题,层与层之间通过垂向动量交换联系起来.此外,针对限制时间步长的线性不稳定因子,采取相应的措施,使模式的时间步长不受稳定性限制:(1)科氏力项,在运动方程中通过对科氏力项引入一个稳定的二阶旋转矩阵来克服它在时间迭代过程中产生的线性不稳定[13];(2)外重力波项,对运动方程中的正压梯度力项和连续方程中的水平散度项,采用半隐式的差分格式,以克服由外重力波所引起的对稳定性的限制;(3)垂向黏性(扩散)项,对运动方程中的垂向黏性项和温度、盐度方程中的垂向扩散项,采用半隐式的差分格式,以克服该两项对稳定性的限制.通过以上处理并作进一步推导后,可以得到关于水位的椭圆型方程及输运量的三对角方程.用超松弛法可求解关于水位增量的椭圆型方程,从而得到下一时间步长的水位值,用追赶法可求解三对角方程得到水平输运量和水平流速分量,然后求解层积分方程,确定垂向流速场.由于本文仅探讨浙江近海的主要正压动力过程:潮汐和潮流,故采用HAM SOM 模式的正压模拟部分,对温度、盐度方程暂不考虑.212 模式的应用21211 计算海区和模式安排本计算海域为27b ~31b N,12011b ~12311b E 的整个浙江沿岸水域,区域地形如图1所示.经向和纬向水平网格均取1c ,纬线(x 方向)上网格的大小由纬度来校正,垂向取8层,这是为反映岛屿间水道或海峡通道地形的剧烈变化而设置.从上向下各层的厚度分别为10,20,30,40,50,75,100,150m.时间步长取5m in.陆地边界取法向流速为0;开边界潮位给定,包含5个主要分潮(M 2,S 2,N 2,K 1和O 1),由大区域模式(渤海、黄海、东海,如图2)提供.开边界处流速的边界条件根据Orlanski 辐射条件来确定.图1 区域地形(等值线单位:m)及实测比较点的位置图2 大、小计算海区及网格示意图在模式应用过程中,由于研究海域的特殊性,将引入两项数值技术,即网格嵌套和动边界技术.21212 引入网格嵌套技术为小区提供开边界条件浙江近海是个相对于渤海、黄海、东海等较小的海域,且沿岸地形复杂,岸线曲折.为了得到反映该区特征的数值计算结果,更好地刻划岸线以及尽可能反映影响因素(如江河入海径流),都希望把网格划分得足够细,但是小区域的开边界条件往往难以给定,因为研究区域小,115期陈倩等:浙江近海潮汐潮流的数值模拟开边界可能引入误差的影响就比较大,为减小此误差,又希望把计算海域尽量扩大.如果在扩大海域下采用单一的细网格,将会增加计算量,因此我们引入网格嵌套技术来解决这一矛盾,即对外部大区域用粗网格,对内部小区域用细网格分别进行计算,其中大区的计算结果为小区提供所需的开边界条件,这样既保证了小区域计算精度,又可避免开边界上过大的误差影响.本文取24b~41b N,117b~131b E的整个渤海、黄海、东海为大区域,经向和纬向水平网格都取10c,小区域即浙江沿岸海域的大小计算海区如图2所示.开边界上潮位边界条件为F=E5i=1H i cos(R i t-g i),(1)式中,i从1到5分别表示5个主要分潮M2,S2,N2,K1和O1.小区域计算中,开边界上的调和常数由大区的计算结果经线性插值得到[16~18].图3浙江近海的潮滩示意图21213引入动边界模拟潮滩涨落浙江沿岸河口、港湾众多,滩地广泛分布.图3为浙江近海的潮滩示意图,其中黑色区域为潮滩.潮滩沙脊等随着潮位的涨落时而淹没时露出,相应水域的面积也随之增大或减小.为正确模拟这类变化的水域,我们引入动边界技术.本文将采用水位判别法[17~23].在计算过程中,随时对滩地网格点作状态判别,根据该时刻某点的瞬时水深判断其是否淹没,以确定该点是否参与计算.对于所有的潮滩点,由于最大可能的瞬时水深都不超过10m(实测最大潮差为8193m),故垂向仅有一层,我们可以把它当作二维问题来讨论.对于平面上任一潮滩点(i,k),瞬时水深h i,k= D i,k+F i,k,其中D i,k,F i,k分别为(i,k)点的滩地高程及潮位.由于算法的限制,对动边界的判别指标往往不能以单元水深为0来确定,而需要引入一个富裕水深D0(本文中取D0=011m)来保证求解的稳定性.落潮过程中,F i,k<0,当h i,k[D0时,认为该点干出,不参与计算,且令其流速为0;反之,则认为该点淹没,参与计算.涨潮过程中,水边界线随潮位的上升向高滩推进,则计算网格点增多.因新增网格点原无潮位值,故取其周围4点(i+1,k),(i-1,k),(i,k+1),(i,k-1)中已淹水的诸点潮位的平均值:F i,k=(F i+1,k E i+1,k+F i-1,k E i-1,k+F i,k+1E i,k+1+F i,k-1E i,k-1)(E i+1,k+E i-1,k+E i,k+1+E i,k-1),(2)式中,E i,k=1水点,0陆点.涨潮时,F i,k>0,当h i,k\D0时,认为该点淹没,参与计算,反之不参与计算.12海洋学报25卷3 计算结果及与实测的比较根据潮汐调和分析理论,分离主要分潮调和常数所需的最短潮位时间序列为15d,在实际应用中,常用1个月的潮位时间序列以获得较好的结果.计算的初始条件为:当t =0时,u =v =w =F =0,其中u,v 分别为水平流速的东分量和北分量,w 为垂向流速,F 为未扰动海面上的潮高.潮汐和潮流模拟时,模式共运行了32d,在2d 内模式已由零初始状态完全建立起来,模式产生的水位和流速数据每隔1h 进行储存,后30d 的逐时数据用于分析计算潮汐和潮流的调和常数,这里我们用最小二乘法来计算7个主要分潮(M 2,S 2,N 2,K 1,O 1,M 4和MS 3)的调和常数.为了验证计算结果的可靠性,我们从历史观测资料中挑选出50个有代表性的沿岸水位站,其位置如图1中空心五角星所示.此外,选出105个潮流比较点,其位置如图1中实心圆点所示.通过比较模拟计算与实测点的调和常数来验证计算的可靠性.表1为50个潮位站上4个主要分潮调和常数之差值(该值为实测减模拟所得)的均方差比较.为验证网格嵌套和动边界技术应用的效果,我们分别列出了10c 固定边界、1c 固定边界、1c 动边界3种情况下的均方差.由表1可见,计算与实测都符合较好.M 2分潮振幅差值中有82%的站点1c 细网格的计算精度优于10c 粗网格,1c 动边界的计算精度又优于1c 固定边界;10%的站点1c 固定边界效果最好,1c 动边界其次,10c 粗网格最差;另外8%的站点1c 动边界计表1 50个潮位站上4个主要分潮调和常数的模拟与观测之差值的均方差比较M 2S 2K 1O 1$H /cm$g /(b )$H /cm $g /(b )$H /cm $g /(b )$H /cm $g /(b )10c812110160817010167310010100210781701c615081786120719221688100118071151c 动边界41607114510051402125517611565150算结果最好,10c 粗网格次之,1c 固定边界最差.其他分潮的振幅、相位差统计结果与此皆相差不大.因此,综合来看,应用了网格嵌套和动边界技术之后,各点的计算精度普遍得到提高,效果明显.在下面的讨论中,我们将仅以1c 动边界情况下所得的模拟计算结果作为研究对象.图4a,b 分别为50个潮位站上4个主要分潮(M 2,S 2,K 1,O 1)调和常数的振幅和相角的实测值(x 轴)与计算值(y 轴)的对比结果,图中的点均分布在第一象限的角平分线附近,表明计算值与观测值比较一致.浙江近海潮流的数值模拟虽然是在三维的HAMSOM 模式下进行,但除近海底外,潮流的垂向变化都较小,因此仅以表层的潮流特征来讨论.图5a,b 分别给出了105个测流站表层M 2,K 1分潮流调和常数分量U cos N ,U sin N ,V cos G ,V sin G 的实测值(x 轴)和计算值(y 轴)的对比结果.测点均匀分布于计算海区.由图可见,计算值与实测值符合较好.误差统计分析表明,对M 2分潮流来说,两者偏差绝对值小于10cm/s 者占93%,最大偏差为1617cm/s;对K 1分潮流来说,两者偏差绝对值小于3cm/s 者占97%,最大偏差为4159cm/s.135期陈倩等:浙江近海潮汐潮流的数值模拟图4 模拟与实测各主要分潮M 2,S 2,K 1,O 1调和常数之比较图5 模拟与实测分潮流分量U cos N ,U sin N ,V cos G ,V sin G 之比较a.M 2分潮流,b.K 1分潮流4 计算结果的讨论411 同潮时线和等振幅图图6a,b,c,d 分别为主要半日分潮M 2,S 2和主要全日分潮K 1,O 1的同潮时线和等振幅线图,由图可见,M 2,S 2分潮同潮时线的走向基本一致.在M 2同潮图中,250b 同潮时线以三门湾为顶点,呈八字形向两旁伸展.图中S 2的300b 同潮时线类似分布.同潮时线的这一分布特征表明:半日潮波进入陆架后,由东南向西北挺进,首先在三门湾口附近达到高潮,然后分南北两支传播.K 1,O 1分潮的等振幅线图大致相似,两者的振幅都是由东向西略有增大,但增幅很小.同潮时线分布表明,K 1分潮由东北向西南传播;O 1分潮以西北-东南向传入本海域,在三门湾附近传向渐变为东北-西南.在本海域内,各分潮都没有出现无潮点.14海洋学报 25卷图6潮位的同潮时线和等振幅线a1M2分潮, b.S2分潮, c.K1分潮, d.O1分潮等振幅线/cm;同潮时线/(b)155期陈倩等:浙江近海潮汐潮流的数值模拟412 潮汐类型潮汐类型是根据潮型数F =(H K 1+H Q 1)/H M 2来划分的,它反映了日分潮与半日分潮的相对重要性.根据值的大小,一般可把潮汐分为4种类型,即规则半日潮(010<F [015)、不规则半日潮(015<F [210)、不规则日潮(210<F [410)和规则日潮(F >410).由图7a 可见,浙江近海以半日分潮为主.大部分海域内F [015,为正规半日潮,如浙北的杭州湾以及浙中、浙南的沿岸海域;小部分地区015<F [210,为非正规半日潮,主要分布在甬江两侧并连及舟山群岛部分海区.其产生原因是由于潮波变形,H M 2减小,而H Q 1和H K 1基本不变,故比值增大.图7 潮型数(F )a.潮位,b.表层潮流413 潮流的同潮时线和等振幅图图8a,b 分别给出了M 2,K 1分潮流表层合成流等振幅线和同潮时线的分布.由图可见,M 2,K 1分潮流等振幅线的分布显示,在本海域内存在几个明显的强流区.浙北的杭州湾和舟山群岛诸水道内潮流最强,M 2分潮流的振幅可达120cm/s 以上,这是由于潮流受地形影响显著.由于杭州湾的喇叭口地形,定海、岱山、嵊泗等海区的一些狭长水道,潮流通道的截面较小,故潮流速很强.另外,在三门湾、温州湾等港湾区,潮流振幅随近岸距离的减小而有明显增强趋势,这种变化趋势主要是由海区的水深条件、岸线和地形等因素造成的.M 2分潮流在29b 18c N,122b 46c E 处有一个圆流点,这与文献[11]中所得结果一致.同潮时线的旋转方向为反时针旋转.16海洋学报 25卷图8 潮流的同潮时线和等振幅线(表层)a 1M 2分潮, b.K 1分潮等振幅线/cm,等潮时线/(b )图9 M 2分潮流旋转率K 的分布(表层)414 潮流类型和潮流的运动形式潮流类型的划分标准与潮汐类型相类似.根据我国5港口工程技术规范6[24]的规定,采用(W O 1+W K 1)/W M 2作为指标,其中W O 1,W K 1,W M 2分别为O 1,K 1,M 2分潮流的椭圆长轴.由图7b 可见,绝大部分海域表层潮流F 值在012左右,均小于015,故浙江近海基本都属于正规半日潮流区.潮流的运动形式由潮流的椭圆旋转率K 来描述.K 值为潮流椭圆的短轴与长轴之比.当K 值大于0125时,潮流表现出较强的旋转性,而当其小于0125时,潮流表现为往复流.K 值前面的正负号表示潮流的旋转方向,正号为左旋,负号为右旋.由于本区域内半日潮流具有支配地位,因此我们给出M 2分潮流旋转率K 的分布来表征潮流的旋转特征,如图9所示.由图可见,在沿岸的港湾、河口水域及潮汐通道等处,潮流运动由于受地形、边界条件的制约,往复流的性质非常明显,K 的绝对值多小于012.外海或离岸较远且较宽敞的海区K 175期陈倩等:浙江近海潮汐潮流的数值模拟的绝对值大于014,呈旋转流态,如大目洋、三门湾、温州湾外海等.从旋转方向上看,大多数站点的K 为负号,呈右旋,这是由于地球自转效应产生的结果.另外,杭州湾口、浙闽交界水域各有一明显左旋区,这是潮波干涉区的影响所致.415 潮流椭圆图10a,b 分别给出了5c 为间隔的计算点上主要半日分潮M 2和主要全日分潮K 1的潮流椭圆长短轴分布.潮流的椭圆长轴指示了最大流速和最大流速方向.由图可见M 2,K 1分潮流椭圆长短轴的分布较为相似.从南往北,流速渐增.远岸区域旋转性较强,流速较弱;近岸区域则多为往复流,流速较强.M 2分潮流椭圆长轴的分布比较规则,其变化与地形密切相关.本区东南海域椭圆长轴的走向为西北-东南,这显示了潮波的传入方向.潮波传入后分为两支,往北和往南的传播方向各不相同,这在椭圆长轴的走向上有所反映.往南的一支表现为椭圆长轴在浙南海区多为东-西走向,而往北一支潮波则表现为椭圆长轴由浙中的东北-西南走向渐变为浙北的西北-东南走向,进入杭州湾后,基本为东-西走向.在港湾区或河汊水道处,椭圆长轴的方向一般与岸线或港湾水道走向相一致.K 1分潮流椭圆长轴之走向,在三门湾以北区域多为西北-东南,从三门湾往南渐变为东北-西南走向.在近岸区域,椭圆长轴的方向与岸线平行.图10 潮流椭圆长短轴分布图(表层)a.M 2分潮流,b.K 1分潮流5 小结本文用三维陆架海模式对整个浙江近海的潮汐和潮流进行了三维数值模拟,并针对浙江近海岸线曲折、潮滩广泛分布等地形特点,引入网格嵌套和动边界技术以提高计算的精度.改18海洋学报 25卷进后的HAM SOM 模式被成功地应用于该区的潮流数值模拟中,通过与实测点调和常数的比较,验证了这两项技术引入原模式后的良好效果,计算精度普遍得到提高.这是HAMSOM 模式首次应用于浙江近海的潮汐、潮流计算,该模式在研究陆架海动力学上有着独特的优越性.模式在应用过程中根据实际需要所作的改进也被证明是成功的,这使得该模式得到进一步的完善.改进后的模式可以嵌套用于局部更小区域的精细计算,也可以处理边界变动的问题.除模式的改善和成功应用外,本文将浙江近海作为一个整体进行数值计算并得到了一批反映规律的研究成果.本文给出了全区范围内4个主要分潮M 2,S 2,K 1,Q 1的潮汐同潮图,发现在本海域内,这4个分潮都没有出现无潮点;给出了M 2,K 1分潮潮流同潮图,发现M 2分潮流在29b 18c N,122b 46c E 处有一个圆流点,K 1分潮流在该区内无圆流点;此外还得到了潮型数F ,M 2分潮流椭圆率K 和潮流椭圆的分布规律.以上模拟结果都与实测情况[1]符合良好,而模拟结果由于不受实测点的限制,因此更全面、完整.参考文献:[1] 陈倩,黄大吉,章本照.浙江近海潮汐特征的研究[J].东海海洋,2003,21(2):1)12.[2] 曹德明,方国洪.杭州湾潮汐潮流的数值计算[J].海洋和湖沼,1986,17(2):93)101.[3] 曹德明,方国洪.杭州湾和钱塘江潮波的联合数值模型[J].海洋学报,1988,10(5):521)530.[4] 李身铎,顾思美.杭州湾潮波三维数值模拟[J].海洋与湖沼,1993,24(1):7)15.[5] 曹欣中,唐龙妹,张月秀.宁波、舟山内海域实测海流分析及潮流场的数值模拟[J].东海海洋,1996,14(2):1)9.[6] 董礼先,苏纪兰.象山港潮波响应和变形研究.Ò.象山港潮波数值研究[J].海洋学报,1999,21(2):1)8.[7] 周大成.椒江河口大潮潮流特性的数值分析[J].浙江水利水电专科学校学报,2001,13(3):13)14,26.[8] 李孟国,王正林.瓯江口潮流数值模拟[J].长江科学院院报,2002,19(2):19)22.[9] 赵保仁,方国洪,曹德明.渤、黄、东海潮汐潮流的数值模拟[J].海洋学报,1994,16(5):1)10.[10] 万振文,乔方利,袁业立.渤、黄、东海三维潮波运动数值模拟[J].海洋与湖沼,1998,29(6):611)616.[11] 王凯,方国洪,冯士.渤海、黄海、东海M 2潮汐潮流的三维数值模拟[J].海洋学报,1999,21(4):1)13.[12] BACKHAUS J O.A sem-i implicit scheme for the shallow water equations for applications to shelf sea modelling[J ].ContS helf Res,1983,(2):243)254.[13] BACKHAUS J O.A three -di m ensional model for the simulation of shelf sea dynamics [J ].Deutsche HydrographischeZeitschrift,1985,38:165)178.[14] BACKHAUS J O,HAINBUCHER D.A finite difference general ci rculation model for shelf seas and i ts application to lowfrequency variability on th e North European shelf[A].NIHOUL J C ,JAM ART B M .Three -Dimensional M odels of M arine and Estuarine Dynamics[M ].Amsterdam:Elsevier Science Publishing B.V.,1987.221)244.[15] 黄大吉,陈宗镛,苏纪兰.三维陆架海模式在渤海中的应用[J].海洋学报,1996,18(5):1)13.[16] 辛文杰.差分模型网格嵌套边界技术在工程潮流计算中的应用[J].水利水运科学研究,1999,12:355)360.[17] 辛文杰.河口、海湾平面潮流数值计算中的几个问题[J].水动力学研究与进展,1993,8(3):348)354.[18] FULT ON S R.An adaptive multigrid barotropic tropical cyclone track model[J].M onthly Weather Review ,2001,129(1):138)151.[19] 曹德明.胶州湾潮汐潮流的数值计算[J].海洋科学集刊(21)[C].北京:科学出版社,1984.157)164.[20] 孙英兰,张越美.胶州湾三维变动边界的潮流数值模拟[J].海洋与湖沼,2001,32(4):355)361.[21] 李燕初,蔡文理.ADI 潮汐模型的活动边界方法及其效应[J].海洋学报,1993,15(2):115)120.[22] 张存智,杨连武,窦振兴.具有潮滩移动边界的浅海环流有限元模型[J].海洋学报,1990,12(1):1)13.[23] 韩康,吴冠,张存智.普兰店湾潮流场数值模拟[J].海洋环境科学,2001,20(1):42)46.[24] 中华人民共和国交通部.港口工程技术规范(上册)[M ].北京:人民交通出版社,1987.195期陈倩等:浙江近海潮汐潮流的数值模拟20海洋学报25卷Nu merical simulation of tide and tidal currentsin the seas adjacent to ZhejiangCH EN Qian1,2,3,HUANG Da-ji2,ZHANG Ben-zhao1(1.Department of M echa nics,Zhe j iang Univ iersity,Hangz hou310012,China;2.Second I nstitute of Oceanography,State Oceanic A d ministration,Hangz hou310012,China;3.K ey L aboratory of Ocean Dynamic Proc esses and Satellite Oceanogra-p hy of S tate Oce a nic A d ministration,Hangz hou310012,China)Abstract:By means of three-dimensional baroclinic primitive equation model)))Hamburg shelf o cean model (HA M SOM),the tide and tidal cur rents in the seas adjacent to Zhejiang are simulated.Fur thermore,the original model is impro ved by two numerical technologies)))nested g rid and mov ing boundary method,which are intro-duced to increase the computat ional precision.T he impro ved mode is proved to be successful while it is applied to the seas adjacent to paring the computed values w ith t hose of50tidal observ ator ies,it is found that the computat ional precision w ith fine gr ids and moving boundar y is gener ally higher t han that with coarse gr ids or fix ed boundary.T he root-mean-square values of comparativ e results show that t he difference betw een the simulated and t he observed amplitudes of M2constituent is only4.6cm,the differ ence of phase-lags is7.14b;the difference of amplitudes and phase-lags of S2constituent are5.0cm and5.4b;the difference of amplitudes and phase-lags of K1 constituent ar e2.25cm and5.76b;the difference of amplitudes and phase-lags of O1constituent are1.56cm and 5.5b.T hese indicate that the computational r esults agree with the observed ones very well.I n addition,105current stat ions are chosen,and the difference between the calculated and the observed harmonic co nstant,U cos N,U sin N, V cos G,V sin G of M2and K1component curr ents at surface layer is compared,and the results also show a goo d a-gr eement.Based on these results,the co-amplitude and co-phase lag lines of the main co mponent tides,the type of t ide,tidal current ellipse,the type and the mov ing mode of tidal current ar e given.It is found that the four main constituents M2,S2,K1,O1have no tide-fr ee points in this ar ea;and the M2co mponent current has o ne current-amphidromic point at29b18c N,122b46c E.In addition,some meaning ful results are concluded,and ag ree w ith the observed ones w ell.T hereby,a thorough kno wledg e of the character istics of tides and t idal currents is got in the w ho le coastal zone of Z hejiang Prov ince.Key words:numerical simulation;nested gr id;mo ving boundary;tides and tidal cur rents;seas adjacent to Zhejiang。

基于潮波数值模拟技术的潮位推算及其精度验证

基于潮波数值模拟技术的潮位推算及其精度验证缪锦根;刘雷;董玉磊【摘要】基于潮波数值模拟技术的潮位推算算法可以推算测区内任一区域的潮位,减少临时验潮站,尤其是海上定点验潮站的布设,避免潮位改正模型的误差.文章利用实测潮位资料对基于数值模拟技术的推算潮位进行了验证,取得了较好的结果.【期刊名称】《水道港口》【年(卷),期】2010(031)003【总页数】4页(P195-198)【关键词】数值模拟;潮位推算;精度验证【作者】缪锦根;刘雷;董玉磊【作者单位】天津海事局海测大队,天津,300222;天津海事局海测大队,天津,300222;天津海事局海测大队,天津,300222【正文语种】中文【中图分类】O242.1中国沿海潮波主要是太平洋潮波传入后形成的协振潮，由天体引潮力产生的强迫潮很小，但受地形岸线、海底摩擦、地转偏向力以及浅水效应等因素影响，动力变化巨大，潮汐性质复杂，潮差空间分布悬殊，因此在水深测量的同时必须进行潮位观测，以提供必要的水位改正资料。

传统的水位改正方法是在岸上或海上设立验潮站，实行双站分带或三（多）站分区水位改正，有效降低了潮差对测量结果的影响，提高了测深精度，但也存在如下问题：（1）海上及孤岛验潮站设立成本较高，且无精度及完整性保障；（2）要求各验潮站的潮位数据为水深测量期间的同步数据，不能充分利用测区内已有的潮位资料；（3）原则上要求测量区域内潮时、潮差在验潮站之间或控制范围内线性变化。

事实上，实际的潮波运动要复杂得多，没有一定密度的验潮站无法真实地反映测区内的潮波特点。

基于余水位空间相关性的潮位推算方法有效地解决了具有历史验潮资料海区存在的前2个问题，而日趋成熟和完善的潮波数值模拟技术使第三个问题的解决成为了可能［1-9］。

1 潮位推算基本原理潮位推算的基本思想是假设水位由潮汐和非潮汐水位（余水位）2部分构成，并且非潮汐水位具有良好的空间一致性，采用长期验潮站的实测水位剥离潮汐和非潮汐水位部分，并用非潮汐部分改正预报站的潮位。

中国近海潮波运动数值模拟

中国近海潮波运动数值模拟本文基于球面坐标系下的二维垂线平均潮波运动方程建立中国近海潮波数学模型，模型区域包括渤海、黄海、东海、南海、泰国湾和环台湾岛海域，网格尺寸2′×2′，网格数1175×955。

在考虑引潮力情况下，计算模拟了中国近海的复合潮波运动；并对分布于各个海域的281个潮位站的4个主要分潮（M2、S2、K1、O1）潮位调和常数以及13个海洋预报站的潮流资料进行了验证，验证结果基本合理。

针对计算结果绘制了主要分潮包括浅水分潮（以M4为例）的潮汐同潮图和潮流同潮图，对中国近海潮汐和潮流分布即潮波运动进行了分析，并和前人的结果进行比较，结果基本吻合。

整个东中国海的潮波主要是太平洋潮波经台湾和九州之间的水道传入的协振波，南海的潮波主要是太平洋潮波经吕宋海峡传入的协振波。

东海和南海主要通过台湾海峡进行水量和潮能交换。

在东中国海基本以半日潮为主，尤其是M2占优，而在南海基本以全日潮为主。

由于受到地形影响、边界的反射、地转偏向力和陆架浅海的摩阻作用，潮波在各海区或以前进波或以无潮点和圆流点为主要特征的旋转潮波系统组成了复杂的潮波系统。

在模型计算的基础上对台湾海峡的M2分潮的潮汐分布特征和传播机制进行了探讨。

认为台湾海峡的M2分潮主要是有北部的前进波和南部的前进驻波系统组成，由吕宋海峡进入的太平洋潮波和广东、福建沿海岸线的相互作用形成了南端前进驻波现象和北部地形边界的放大效应产生的潮能幅聚现象是台湾海峡M2分潮分布特征的主要原因。

二维潮流泥沙数值模拟及可视化研究

及深入地系统研究二维水沙运动基本规律．是建
象．并利用数值计算进行近似求解以复演某种自然或人工过程这些在严世强．熊德琪Ｉｌｌ对潮流数值模拟研究进展中均有所提及。与物理模型相比。
一
ｌ０一
吉林水利
二维潮流泥沙数值模拟及可视化研究
朱嘉玲等２０１３年１Ｏ月
一
Ｏｈ＋生＋生：０
ＯｔＯｘａｖ
计算区域内的任意位置及相应模拟时间段以内的任意时刻进行潮流及泥沙运动复演．同时还可对港口、航道、围垦、码头等水利及海洋工程实施之后的潮流及泥沙运动进行模拟研究。主界面分为
第１０期（总第３７７期）
［文章编号］１００９ — ２８４６（２０１３）１０ — ００１０ — ０６
吉林水利
２０１３年ｌ０月
二维潮流泥沙数值模拟及可视化研究
朱嘉玲．石然２１０００９）（河海大学港口海岸与近海工程学院，江苏南京
是二维泥沙数学模型问题至今尚未能很好地解决有两点重要原因．其一是由于泥沙运动规律极为复杂．对冲积河道的水流和泥沙运动规律认识不
足，很多问题，诸如输沙率、阻力等问题，均没有较好地解决现有的二维泥沙数值模拟大都是直接引用一维计算公式．对于这样做是否适合很少涉

波—流共同作用下长江口二维悬沙数值模拟

波—流共同作用下长江口二维悬沙数值模拟波—流共同作用下长江口二维悬沙数值模拟摘要：本文通过数值模拟方法，研究了波—流共同作用下长江口内二维悬沙运移的情况。

通过建立数学模型和计算方法，模拟了长江口内悬沙的输运过程，并对波浪和流场的相互影响进行了探讨。

结果表明，波浪和流场的共同作用对悬沙输运过程具有显著影响，进一步深化了对长江口悬沙运动规律的认识。

一、引言长江口是我国最重要的河口之一，其水动力和悬沙运动规律的研究对于河口水文、水沙动力学的理解具有重要的意义。

波浪和流场作为长江口内两个主要的力学驱动因素，其相互作用对悬沙的输运过程有着重要影响。

因此，基于数值模拟方法，研究波流共同作用下长江口内二维悬沙运移的规律具有重要的理论和实际意义。

二、模型建立在本文研究中，建立了二维悬沙运移数值模型，考虑了波浪、流场以及悬沙之间的相互作用。

模型采用了数值解的方法，通过离散化长江口内区域，并运用流动方程、输运方程和释放方程来模拟波流共同作用下的悬沙运动过程。

模型的输入数据包括入口流速、入口浓度和入口波浪参数等，模拟的输出结果为悬沙在长江口内的输运情况和浓度分布。

三、数值模拟结果通过数值模拟得出的结果显示，长江口内的悬沙输运过程受到波浪和流场的共同作用。

在波浪作用下，悬沙会随着波动而上升和下降，这会影响悬沙的输运速度和方向。

而流场则会决定悬沙的平均输运速度和方向。

模拟结果还显示，在大波浪和强流场共同作用的情况下，悬沙会有较大的浓度差异和输运速度。

四、悬沙输运规律分析根据数值模拟结果，我们进一步分析了悬沙输运的规律。

首先，当波浪和流场方向相一致时，悬沙的输运速度较快；当方向相反时，输运速度较慢。

其次，波浪的幅度和周期对悬沙的上升和下降速度有较大影响，较大的波浪幅度和周期会引起更剧烈的悬沙运动。

此外，在大流量的情况下，悬沙的输运速度较快，浓度较高。

五、结论通过本文的研究，我们深入了解了长江口内波流共同作用下的悬沙运移规律。

瓯江口波流耦合作用下二维流场数值模拟研究

近岸海区的水动力特性。瓯江口海区岛屿林立、沙滩密布、水道纵横、
（１）
地形复杂，波浪潮流相互作用显著。近年来有很多关于该海区水动力环境的数值研究＿３４。，其中以李孟国等＿５刮的研究最具代表性。但以上研究一般都是将波浪和潮流分开来处理，考虑波流耦合作用的少；即使是耦合计算，也是单向耦合的多，双向耦合的少。因此建立双向波流耦合作用数值
Ｅ－ｍａｉｌ：ｗａｃｋｙＯ１＠１２６．ｃｏｉｎ
代表使方程能量守恒的源项，．ｓｉ指风输入的能
量；５指波与波之问的非线性作用引起的能量损耗；Ｓ。指由白帽引起的能量损耗；Ｓｈｏｔ指底摩阻引起的能量损耗；Ｓ指由于水深变化引起的波浪破
一
Ｃｄ
＝
式中：ｆ —— 时间；
、
．
１．２５５×１０一。
ｌ｛＜７ｍ／ｓ
ｙ、ｚ —— 右手Ｃａｒｔｅｓｉａｎ坐标系；
ｄ — —静止水深；
— —
｛（Ｑ８＋０．１３６５ｆ）ｘ１０７ｍ／ｓ ≤ Ｉ ≤ ２５ｍ／ｓ
１３０
瓯江口波流耦合作用下二维流场数值模拟研究—— 王留洋，陈晓亮
碎产生的能量损耗。
１．２潮流模型

河网二维水流数值模拟

第二章
曲线正交网格变换与水流控制方程 ............................... 7
2.1 概述................................................................................................................................... 7 2.2 Possion 方程的构造 .......................................................................................................... 7 2.3 方程的求解..................................................................................................................... 12 2.4 网格计算实例................................................................................................................. 12 2.5 基本水流控制方程 ......................................................................................................... 14 2.6 正交曲线坐标下的方程 ................................................................................................. 14 2.7 小结................................................................................................................................. 16

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

! 二维潮波数值模拟与预报! 计算海域包括渤海和黄海北部（117.50 -126.67 E,34.00 - 41.00 N）! 空间分辨率是10'x10'! 开边界处的h和g由TMD TOOLBOX求得!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!program hw4implicit noneinterfacesubroutine h_get(h,dx,dy,dt,u,v,d,cd)integer::i,jreal,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cdreal,dimension(43,56),intent(in)::d,u,vreal,dimension(43,56),intent(inout)::hreal,dimension(43),intent(in)::dxend subroutine h_getsubroutine v_get(v,u,cu,cv,dx,dy,dt,d,h,f)integer::i,jreal::kk=1.5e-3,gg=9.8real,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cu,cvreal,dimension(43,56),intent(in)::d,h,ureal,dimension(43,56),intent(inout)::vreal,dimension(43),intent(in)::dx,freal,dimension(43,56)::uu,ssend subroutine v_getsubroutine u_get(u,v,cu,cv,dx,dy,dt,d,h,f)integer::i,jreal::kk=1.5e-3,gg=9.8real,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cu,cvreal,dimension(43,56),intent(in)::d,h,vreal,dimension(43,56),intent(inout)::ureal,dimension(43),intent(in)::dx,freal,dimension(43,56)::vv,rrend subroutine u_getend interface! < PART 1 >! dx, dy : 网格宽度，单位（m）! dt : 时间步长，取为M2分潮周期的1/240，单位（s）! d : 各个网格点的水深! h : 海水相对静止海面的高度! u, v : 海水水平流速! uu, vv : u、v的平均值! f : 科氏参量! cd : 水深控制场! cu, cv : 速度控制场! h0, g0 : 开边界处的振幅和迟角! h_check : 用于检验迭代是否达到要求的精度! eps : 迭代精度!< PART 2 >! h_3d : 迭代完成后得到的水位! a,b,g,hh : 调和常数! l1,l2,l3,l4,y1,y2 :最小二乘法正规方程组的系数!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!real,parameter::pi=3.14159625,dy=6371004*2*pi/360/6real,parameter::dt=24.8412/2*3600/240,w=360/(240*dt),eps=0.01 !M2 !real,parameter::dt=24.8412*3600/240,w=360/(240*dt),eps=0.01 !K1 integer::i,j,k,m,n,tinteger,dimension(43,56)::cd,cu,cvreal,dimension(43,56)::d,h,u,v,uu,vv,h_checkreal,dimension(43)::dx,freal,dimension(19:54)::h0,g0real,dimension(43,56,240)::h_3dreal::l1,l2,l3,l4,cos_g,sin_greal,dimension(43,56)::a,b,g,hh,y1,y2real,dimension(240)::temp1,temp2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! STEP1:读取水深数据、控制场和开边界点的振幅与迟角!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!open(1,file='BHhai.dat')read(1,'(56f5.0)') (d(i,:),i=1,43)close(1)open(2,file='Bctrlh.dat')read(2,'(56i1)') (cd(i,:),i=1,43)close(2)open(3,file='data.out_M2')!open(3,file='data.out_K1')do i=19,54read(3,'(42x,f6.4,f11.2)') h0(i),g0(i)end doclose(3)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! STEP2: 计算各个纬度的dx和f,计算各个网格点的流速控制场cu和cv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!do i=1,43dx(i)=6371004*cosd(41.0-(i-1)/6)*2*pi/360/6f(i)=2*7.292115e-5*sind(41.0-(i-1)/6)end dodo i=2,42do j=2,55cu(i,j)=cd(i,j)*cd(i,j+1)end doend dodo i=2,42do j=2,55cv(i,j)=cd(i,j)*cd(i-1,j)end doend do!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! STEP3:迭代循环求解!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!h_check=0do j=19,54h(42,j)=h0(j)*cosd(-g0(j))end doprint*,'迭代控制误差',epsprint*,'迭代次数'do while(maxval(abs(h-h_check))>eps)t=1h_check=hdo while(t<=240)do j=19,54h(42,j)=h0(j)*cosd(w*t*dt-g0(j))end docall h_get(h,dx,dy,dt,u,v,d,cd)call v_get(v,u,cu,cv,dx,dy,dt,d,h,f)call u_get(u,v,cu,cv,dx,dy,dt,d,h,f)t=t+1do j=19,54h(42,j)=h0(j)*cosd(w*t*dt-g0(j))end docall h_get(h,dx,dy,dt,u,v,d,cd)call u_get(u,v,cu,cv,dx,dy,dt,d,h,f)call v_get(v,u,cu,cv,dx,dy,dt,d,h,f)t=t+1end dok=k+1print*,kend do!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !STEP 4:迭代收敛后输出h，u，v!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!t=1do while(t<=240)do j=19,54h(42,j)=h0(j)*cosd(w*t*dt-g0(j))end docall h_get(h,dx,dy,dt,u,v,d,cd)h_3d(:,:,t)=hcall v_get(v,u,cu,cv,dx,dy,dt,d,h,f)call u_get(u,v,cu,cv,dx,dy,dt,d,h,f)t=t+1do j=19,54h(42,j)=h0(j)*cosd(w*t*dt-g0(j))end docall h_get(h,dx,dy,dt,u,v,d,cd)h_3d(:,:,t)=hcall u_get(u,v,cu,cv,dx,dy,dt,d,h,f)call v_get(v,u,cu,cv,dx,dy,dt,d,h,f)t=t+1end do!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! STEP 5 : 用最小二乘法求解各个网格点的调和常数hh和g !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 最小二乘计算原理：《计算方法引论》P52! 调和分析原理：《海洋水文环境要素的分析方法和预报》P97! 求正规方程组的系数l1=0do t=1,240l1=l1+cosd(w*t*dt)**2end dol2=0do t=1,240l2=l2+sind(w*t*dt)*cosd(w*t*dt)end dol3=l2l4=0do t=1,240l4=l4+sind(w*t*dt)**2end dodo i=1,43do j=1,56do t=1,240temp1(t)=cosd(w*t*dt)*h_3d(i,j,t)temp2(t)=sind(w*t*dt)*h_3d(i,j,t)end doy1(i,j)=sum(temp1)y2(i,j)=sum(temp2)end doend do!求解调和常数do i=1,43do j=1,56if(cd(i,j)==0)cyclea(i,j)=(l4*y1(i,j)-l2*y2(i,j))/(l1*l4-l2*l3)b(i,j)=(l3*y1(i,j)-l1*y2(i,j))/(l2*l3-l1*l4)hh(i,j)=sqrt(a(i,j)**2+b(i,j)**2)if(hh(i,j)==0)cyclecos_g=a(i,j)/hh(i,j)sin_g=b(i,j)/hh(i,j)if(cos_g==1.and.sin_g==0)theng(i,j)=0else if(cos_g>0.and.sin_g>0)theng(i,j)=atand(b(i,j)/a(i,j))else if(cos_g==0.and.sin_g==1)theng(i,j)=90else if(cos_g<0.and.sin_g>0)theng(i,j)=atand(b(i,j)/a(i,j))+180else if(cos_g==-1.and.sin_g==0)theng(i,j)=180else if(cos_g<0.and.sin_g<0)theng(i,j)=atand(b(i,j)/a(i,j))+180else if(cos_g==0.and.sin_g==-1)theng(i,j)=270elseg(i,j)=atand(b(i,j)/a(i,j))+360end ifend doend doopen(4,file='H_M2.dat')!open(4,file='H_K1.dat')write(4,'(56f8.3)') ((hh(i,j),j=1,56),i=1,43) close(4)open(5,file='g_M2.dat')!open(5,file='g_K1.dat')write(5,'(56f8.3)') ((g(i,j),j=1,56),i=1,43) close(5)end program hw4!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !子程序!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!计算海水相对于静止海面的高度hsubroutine h_get(h,dx,dy,dt,u,v,d,cd)implicit noneinteger::i,jreal,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cdreal,dimension(43,56),intent(in)::d,u,vreal,dimension(43),intent(in)::dxreal,dimension(43,56),intent(inout)::hdo i=2,41do j=2,55if(cd(i,j)/=0)thenh(i,j)=h(i,j)&&-dt/dx(i)*((d(i,j)+h(i,j)+d(i,j+1)+h(i,j+1))/2*u(i,j)-&&(d(i,j)+h(i,j)+d(i,j-1)+h(i,j-1))/2*u(i,j-1))&&-dt/dy*((d(i,j)+h(i,j)+d(i-1,j)+h(i-1,j))/2*v(i,j)-&&(d(i,j)+h(i,j)+d(i+1,j)+h(i+1,j))/2*v(i+1,j))elsecycleend ifend doend doend subroutine h_get! 计算vsubroutine v_get(v,u,cu,cv,dx,dy,dt,d,h,f)implicit noneinteger::i,jreal::kk=1.5e-3,gg=9.8real,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cu,cvreal,dimension(43,56),intent(in)::d,h,ureal,dimension(43,56),intent(inout)::vreal,dimension(43),intent(in)::dx,freal,dimension(43,56)::uu,ssdo i=2,42do j=2,55if(cu(i,j)/=0)thenif(i==42)thenuu(i,j)=(uu(i,j-1)+uu(i,j+1))/2elseuu(i,j)=(u(i,j)+u(i,j-1)+u(i-1,j)+u(i-1,j-1))/4end ifelsecycleend ifend doend dodo i=2,42do j=2,55ss(i,j)=sqrt(uu(i,j)**2+v(i,j)**2)if(cv(i,j)/=0)thenv(i,j)=&&1.0/(1.0/dt+0.5*kk*ss(i,j)/((d(i,j)+d(i-1,j))/2+(h(i,j)+h(i-1,j))/2))& &*(v(i,j)/dt&&-v(i,j)*0.5*kk*ss(i,j)/((d(i,j)+d(i-1,j))/2+(h(i,j)+h(i-1,j))/2)&&-0.5*uu(i,j)*(v(i,j+1)-v(i,j-1))/dx(i)&&-0.5*v(i,j)*(v(i-1,j)-v(i+1,j))/dy&&-f(i)*uu(i,j)&&-gg*(h(i-1,j)-h(i,j))/dy)elsecycleend ifend doend doend subroutine v_get!计算usubroutine u_get(u,v,cu,cv,dx,dy,dt,d,h,f)implicit noneinteger::i,jreal::kk=1.5e-3,gg=9.8real,intent(in)::dt,dyinteger,dimension(43,56),intent(in)::cu,cvreal,dimension(43,56),intent(in)::d,h,vreal,dimension(43,56),intent(inout)::ureal,dimension(43),intent(in)::dx,freal,dimension(43,56)::vv,rrdo i=2,42do j=2,55if(cv(i,j)/=0)thenif(i==42)thenvv(i,j)=v(i,j)elsevv(i,j)=(v(i,j)+v(i,j+1)+v(i+1,j)+v(i+1,j+1))/4end ifelsecycleend ifend doend dodo i=2,42do j=2,55rr(i,j)=sqrt(u(i,j)**2+vv(i,j)**2)if(cu(i,j)/=0)thenu(i,j)=&&1.0/(1.0/dt+0.5*kk*rr(i,j)/((d(i,j)+d(i,j+1))/2+(h(i,j)+h(i,j+1))/2))& &*(u(i,j)/dt&&-u(i,j)*0.5*kk*rr(i,j)/((d(i,j)+d(i,j+1))/2+(h(i,j)+h(i,j+1))/2)&&-0.5*u(i,j)*(u(i,j+1)-u(i,j-1))/dx(i)&&-0.5*vv(i,j)*(u(i-1,j)-u(i+1,j))/dy&&+f(i)*vv(i,j)&&-gg*(h(i,j+1)-h(i,j))/dx(i))elsecycleend ifend doend doend subroutine u_get。