数值海洋与大气模式（四）：POM模式框架

liangbingquan · 发表于海洋技术专业 2021-11-26 12:58:46

1.POM模式概况　　前文用了两篇文章的篇幅讲述了如何从0到1实现一个浅水方程，涉及到了交错网格、差分离散化和初边值条件的处理等等。本文就来探讨一下海洋模式中最经典的模式之一，POM模式。
　　POM模式的全名为Princeton Ocean Model，在1970年代由G.Mello和Alan Blumberg所开发。经过发展和维护，逐渐成为了可以胜任数值实验和业务化应用的经典模式。尽管从2021年的今天来看，这个模式可能略微跟不上时代，但其经典型和代表性是模式学习者所绕不开的。后续很多海洋模式都是从POM中修改而得到的。POM是一个串行模式，所有代码都写在一个Fortran文件之中。不涉及多文件编译，而且代码结构清晰，是模式学习者初学的首选。除此之外，对于模式的高性能计算的学习者来说，优化POM模式也是很好的实战案例。倘若能用MPI把POM模式改写成并行代码，对代码能力的锻炼是很显著的。
　　POM模式的原始控制方程如下。

$\begin{aligned} \frac{\partial U}{\partial t}+U \frac{\partial U}{\partial x}+V \frac{\partial U}{\partial y}+W \frac{\partial U}{\partial z}-f V &=-\frac{1}{\rho_{0}} \frac{\partial p}{\partial x}+\frac{\partial}{\partial z}\left(K_{M} \frac{\partial U}{\partial z}\right)+F_{x} \\ \frac{\partial V}{\partial t}+U \frac{\partial V}{\partial x}+V \frac{\partial V}{\partial y}+W \frac{\partial V}{\partial z}+f U &=-\frac{1}{\rho_{0}} \frac{\partial p}{\partial y}+\frac{\partial}{\partial z}\left(K_{M} \frac{\partial V}{\partial z}\right)+F_{y} \\ \frac{\partial U}{\partial z} &=-\rho g \\ \frac{\partial T}{\partial t}+U \frac{\partial T}{\partial x}+V \frac{\partial T}{\partial y}+W \frac{\partial T}{\partial z} &=\frac{\partial}{\partial z}\left(K_{H} \frac{\partial T}{\partial z}\right)+F_{T} \\ \frac{\partial S}{\partial t}+U \frac{\partial S}{\partial x}+V \frac{\partial S}{\partial y}+W \frac{\partial S}{\partial z} &=\frac{\partial}{\partial z}\left(K_{H} \frac{\partial S}{\partial z}\right)+F_{S} \\ \rho &=\rho(T, S, p) \end{aligned}\\$

登录/注册后可看大图

2.Sigma坐标系　　前文讲述的浅水模式，介绍了蛙跳格式和交错网格。由于浅水方程组对NS方程做了垂向平均，因此前面提到的网格都是水平网格。在真实的海洋模拟中，水平尺度大于垂直尺度。海底地形起伏较大，所模拟的海区水深可能从几十米到几千米深。如果使用传统的笛卡尔正交坐标系，会出现垂直步长dz不论怎么取都不能满足所有需求。假如近岸水深50m，远洋水深10000m，如果dz取5m，近岸则有10层的网格，而远洋则会出现2000层的网格，造成了极大的计算资源浪费。而如果dz取的比较大，在浅海地区的层数就少的可怜。除此之外，笛卡尔正交坐标系划出的锯齿状很难贴合边界，由下图可以看出，Z坐标系中被底地形横切的网格，不论当做海洋还是当做底地形都会影响精度。

登录/注册后可看大图

　　POM模式给出的解决方案是采用sigma坐标系，该坐标系也被称为地形追随坐标系。有图中可以看出，该坐标系能把海洋各个位置均等的划分同样的层数，在边界上也能很好的贴合地形。因此，在推导POM的方程时，要做的第一步就是将上述控制方程一一进行sigma变换，得到在sigma坐标下的控制方程。

$x^{*}=x, \quad y^{*}=y, \quad \sigma=\frac{z-\eta}{H+\eta}, \quad t^{*}=t\\$

登录/注册后可看大图

　　根据链式法则，就可以得到每个导数项的关系。

$\begin{aligned} \frac{\partial}{\partial x} &=\frac{\partial x^{*}}{\partial x} \frac{\partial}{\partial x^{*}}+\frac{\partial y^{*}}{\partial x} \frac{\partial}{\partial y^{*}}+\frac{\partial \sigma}{\partial x} \frac{\partial}{\partial \sigma}+\frac{\partial t^{*}}{\partial x} \frac{\partial}{\partial t^{*}} \\ \frac{\partial}{\partial y} &=\frac{\partial x^{*}}{\partial y} \frac{\partial}{\partial x^{*}}+\frac{\partial y^{*}}{\partial y} \frac{\partial}{\partial y^{*}}+\frac{\partial \sigma}{\partial y} \frac{\partial}{\partial \sigma}+\frac{\partial t^{*}}{\partial y} \frac{\partial}{\partial t^{*}} \\ \frac{\partial}{\partial z} &=\frac{\partial x^{*}}{\partial z} \frac{\partial}{\partial x^{*}}+\frac{\partial y^{*}}{\partial z} \frac{\partial}{\partial y^{*}}+\frac{\partial \sigma}{\partial z} \frac{\partial}{\partial \sigma}+\frac{\partial t^{*}}{\partial z} \frac{\partial}{\partial t^{*}} \\ \frac{\partial}{\partial t} &=\frac{\partial x^{*}}{\partial t} \frac{\partial}{\partial x^{*}}+\frac{\partial y^{*}}{\partial t} \frac{\partial}{\partial y^{*}}+\frac{\partial \sigma}{\partial t} \frac{\partial}{\partial \sigma}+\frac{\partial t^{*}}{\partial t} \frac{\partial}{\partial t^{*}} \end{aligned}\\$

登录/注册后可看大图

　　用s代表x,y,t的任一项，D海底到自由表面的高度，即

$D=H+\eta$

登录/注册后可看大图

，可得到如下表达式。

$\begin{aligned} \frac{\partial \sigma}{\partial s} &=\frac{\partial}{\partial s}\left(\frac{z-\eta}{H+\eta}\right) \\ &=-\frac{(z-\eta)}{(H+\eta)^{2}} \frac{\partial}{\partial s}(H+\eta)-\frac{1}{H+\eta} \frac{\partial \eta}{\partial s} \\ &=-\frac{\sigma}{(H+\eta)} \frac{\partial}{\partial s}(H+\eta)-\frac{1}{H+\eta} \frac{\partial \eta}{\partial s} \\ &=-\frac{\sigma}{D} \frac{\partial D}{\partial s}-\frac{1}{D} \frac{\partial \eta}{\partial s} \end{aligned}\\$

登录/注册后可看大图

　　由此推导下去即可得到sigma坐标下的控制方程，推导过程极其繁琐，再此省略了推导的中间过程，直接给出结果。为表示方便，后文sigma坐标系的变量中省略其右上角的星号。若对推导的详细过程感兴趣，见文末参考资料。

登录/注册后可看大图

3.内外模态分离　　首先，再回顾一下第一篇文章所讨论过的CFL条件，上次是从数学的角度理解CFL条件为什么能确保线性偏微分方程稳定，这次从波动的角度理解一下CFL条件的意义。由于海洋和大气的动力框架系统为高度非线性系统，因此其稳定性变得更加难以控制。CFL条件是一种很好的参考，而无法绝对确保稳定。

$\mu=c \frac{\Delta t}{\Delta x} \leq 1\\$

登录/注册后可看大图

　　CFL条件中c的物理意义是波速。假设

$\mu=1$

登录/注册后可看大图

，那么此时

$c=\frac{\Delta x}{\Delta t}$

登录/注册后可看大图

。可以看出网格的步长比和波的传播速度相同，意味着这样的网格分辨率是无法分辨这个波的。而当

$\mu>1$

登录/注册后可看大图

时，波速比步长比要大，同样是无法解得这个波的运动状态。这样描述或许不够严谨，但是有助于理解CFL条件的物理意义。结合海洋的实际情况来看，在表层的表面重力波的波速约为200-300m/s，而在海洋内部的重力内波波速远小于表面重力波，大概是在5m/s左右。可以看出，海洋内部的运动过程和海洋表层的运动过程时间尺度相差较大，表层明显快于内部。再回看CFL条件，可以看出如果要想同时满足海洋表层和海洋内部的稳定，表层就需要迁就内部。而POM模式的做法是将表层和内部分离。把表层的正压过程和内部的斜压过程分别称为外模态和内模态，分别设置时间步长。

　　先来看外模态（即正压模态），该部分也被称为快过程，时间步长较小。处理方式类似于浅水方程的推导，对其所在区域做垂向积分，忽略了水平扩散项。在sigma坐标系下的方程组如下所示。

$\begin{array}{l} \frac{\partial \eta}{\partial t}+\frac{\partial \bar{U} D}{\partial x}+\frac{\partial \bar{V} D}{\partial y}=0 \\ \frac{\partial \bar{U} D}{\partial t}+\frac{\partial \bar{U}^{2} D}{\partial x}+\frac{\partial \overline{U V} D}{\partial y}-\tilde{F}_{x}-f \bar{V} D+g D \frac{\partial \eta}{\partial x}=-<w u(0)>+<w u(-1)>+G_{x}- \\ \frac{g D}{\rho_{0}} \int_{-1}^{0} \int_{\sigma}^{0}\left[D \frac{\partial \rho^{\prime}}{\partial x}-\frac{\partial D}{\partial x} \sigma^{\prime} \frac{\partial \rho^{\prime}}{\partial \sigma}\right] d \sigma^{\prime} d \sigma \\ \frac{\partial \bar{V} D}{\partial t}+\frac{\partial \overline{U V} D}{\partial y}+\frac{\partial \bar{V}^{2} D}{\partial y}-\tilde{F}_{y}+f \bar{U} D+g D \frac{\partial \eta}{\partial y}=-<w v(0)>+<w v(-1)>+G_{y}- \\ \frac{g D}{\rho_{0}} \int_{-1}^{0} \int_{\sigma}^{0}\left[D \frac{\partial \rho^{\prime}}{\partial y}-\frac{\partial D}{\partial y} \sigma^{\prime} \frac{\partial \rho^{\prime}}{\partial \sigma}\right] d \sigma^{\prime} d \sigma \end{array}\\$

登录/注册后可看大图

　　对于内模态，则方程形式和第二部分列的形式一样。由于外模态时间步长短，内模态时间步长较长。在POM模式中，内模态的时间步长通常是外模态的数十倍。如果将POM模式的整体结构写成伪代码的话，可以写成如下形式。内模态的时间步长是外模态的isplit倍，这样外模态就可以嵌套在内模态的循环里写。
program POM
Init Paramter       !初始化各种参数，如im,jm等 Init Variable       !初始化T,S,U,V,W等 do iint=1,iend    !内模态循环       call advct()    !计算平流项       call baropg() !计算压强梯度力       do iext=1,isplit !外模态循环          compute el !计算eta          compute ua !计算正压ua          compute va !计算正压va          compute ut,vt!计算正压平均速度       end do          !外模态循环结束       adjust u,v    !       call vertvl    !计算垂直速度       call advq       !计算km,kh       call profq
      call advt       !计算T和S       call proft
      call dens       !计算密度       call advu       !计算u       call profu
      call advv       !计算v       call profv
      print          !将结果输出 end do             !内模态结束end4.湍流闭合方案

$\begin{array}{l} F_{x}=\frac{\partial}{\partial x}\left(2 A_{M} \frac{\partial U}{\partial x}\right)+\frac{\partial}{\partial y}\left[A_{M}\left(\frac{\partial U}{\partial y}+\frac{\partial V}{\partial x}\right)\right] \\ F_{y}=\frac{\partial}{\partial y}\left(2 A_{M} \frac{\partial V}{\partial y}\right)+\frac{\partial}{\partial x}\left[A_{M}\left(\frac{\partial U}{\partial y}+\frac{\partial V}{\partial x}\right)\right] \\ F_{T}=\frac{\partial}{\partial x}\left(A_{H} \frac{\partial T}{\partial x}\right)+\frac{\partial}{\partial y}\left(A_{H} \frac{\partial T}{\partial y}\right) \\ F_{S}=\frac{\partial}{\partial x}\left(A_{H} \frac{\partial S}{\partial x}\right)+\frac{\partial}{\partial y}\left(A_{H} \frac{\partial S}{\partial y}\right) \end{array}\\$

登录/注册后可看大图

　　通过观察可以发现，本文最开端给出的POM原始方程的运动方程和温盐方程都有

$K_{M}$

登录/注册后可看大图

或

$K_{H}$

登录/注册后可看大图

。而在这些方程的末尾，也有

$F_{x}$

登录/注册后可看大图

,

$F_{y}$

登录/注册后可看大图

,

$F_{T}$

登录/注册后可看大图

或

$F_{S}$

登录/注册后可看大图

这些项。这些项的存在使得方程的未知量多于待求解的变量，而如果忽略这些项则会对模拟结果大打折扣。因此，需要解决这些参数的设置问题，而POM模式选择了使用Mellor-Yamada方案，具体形式如下。

$\begin{array}{c} \frac{\partial q^{2}}{\partial t}+U \frac{\partial q^{2}}{\partial x}+V \frac{\partial q^{2}}{\partial y}+W \frac{\partial q^{2}}{\partial z}=\frac{\partial}{\partial z}\left(K_{q} \frac{\partial q^{2}}{\partial z}\right)+2 K_{M}\left[\left(\frac{\partial U}{\partial Z}\right)^{2}+\left(\frac{\partial V}{\partial Z}\right)^{2}\right]+\frac{2 g}{\rho_{0}} K_{H} \frac{\partial \rho}{\partial z}-\frac{2 q^{3}}{B_{1} l}+F_{q} \\ \frac{\partial q^{2} l}{\partial t}+U \frac{\partial q^{2} l}{\partial x}+V \frac{\partial q^{2} l}{\partial y}+W \frac{\partial q^{2} l}{\partial z}=\frac{\partial}{\partial z}\left(K_{q} \frac{\partial q^{2} l}{\partial z}\right)+l E_{1} K_{M}\left[\left(\frac{\partial U}{\partial Z}\right)^{2}+\left(\frac{\partial V}{\partial Z}\right)^{2}\right]+\frac{l E_{1} g}{\rho_{0}} K_{H} \frac{\partial \rho}{\partial z}-\frac{q^{3}}{B_{1}} \widetilde{W}+F_{l} \end{array}\\$

登录/注册后可看大图

　　湍流是一个十分复杂的现象，如果想理解湍流参数化方案，就需要理解什么是湍流。下文将从湍流的本质讲起，逐渐引出湍流参数化方案的全貌。当对模式的动力框架有了比较明确的理解之后，再去看模式代码甚至修改模式代码，就会容易很多。
版权声明　　本文创作的初衷是用于帮助数值模式的学习者。欢迎转载，转载请私信并注明作者和出处，请勿用于任何商业用途。
参考资料A Manual for POM and GOMO. Xiaomeng Huang, Xing Huang. Users Guide For A Three-Dimensional Numerical Ocean Model. George L. Mellor. CEE262c Lecture 8: Sigma coordinates and mode splitting: The Princeton Ocean Model (POM).

数值海洋与大气模式（四）：POM模式框架

相关帖子

您看了很久哦，登陆下吧！