1.POM模式概况 前文用了两篇文章的篇幅讲述了如何从0到1实现一个浅水方程,涉及到了交错网格、差分离散化和初边值条件的处理等等。本文就来探讨一下海洋模式中最经典的模式之一,POM模式。, m$ c. n1 n* I+ l! E" V
POM模式的全名为Princeton Ocean Model,在1970年代由G.Mello和Alan Blumberg所开发。经过发展和维护,逐渐成为了可以胜任数值实验和业务化应用的经典模式。尽管从2021年的今天来看,这个模式可能略微跟不上时代,但其经典型和代表性是模式学习者所绕不开的。后续很多海洋模式都是从POM中修改而得到的。POM是一个串行模式,所有代码都写在一个Fortran文件之中。不涉及多文件编译,而且代码结构清晰,是模式学习者初学的首选。除此之外,对于模式的高性能计算的学习者来说,优化POM模式也是很好的实战案例。倘若能用MPI把POM模式改写成并行代码,对代码能力的锻炼是很显著的。
" F. {) |% k& A6 l5 J6 B" } POM模式的原始控制方程如下。' D# x; ]2 |# j5 t' G+ |
) r9 Z. _) i* {' d" q0 B! j( Q$ A g
9 \1 h6 o( ]0 q4 h2.Sigma坐标系 前文讲述的浅水模式,介绍了蛙跳格式和交错网格。由于浅水方程组对NS方程做了垂向平均,因此前面提到的网格都是水平网格。在真实的海洋模拟中,水平尺度大于垂直尺度。海底地形起伏较大,所模拟的海区水深可能从几十米到几千米深。如果使用传统的笛卡尔正交坐标系,会出现垂直步长dz不论怎么取都不能满足所有需求。假如近岸水深50m,远洋水深10000m,如果dz取5m,近岸则有10层的网格,而远洋则会出现2000层的网格,造成了极大的计算资源浪费。而如果dz取的比较大,在浅海地区的层数就少的可怜。除此之外,笛卡尔正交坐标系划出的锯齿状很难贴合边界,由下图可以看出,Z坐标系中被底地形横切的网格,不论当做海洋还是当做底地形都会影响精度。
& _/ _3 L$ B: E POM模式给出的解决方案是采用sigma坐标系,该坐标系也被称为地形追随坐标系。有图中可以看出,该坐标系能把海洋各个位置均等的划分同样的层数,在边界上也能很好的贴合地形。因此,在推导POM的方程时,要做的第一步就是将上述控制方程一一进行sigma变换,得到在sigma坐标下的控制方程。* V% ?2 ^/ o" H5 b. l+ |* G3 S
9 @+ U. x* |3 t/ E9 S# m& b! t7 S: ~
根据链式法则,就可以得到每个导数项的关系。& `' \2 c" w2 p- G# M' M) n. V1 i3 E8 `
" O4 I# k) ]0 a, g
1 u- w( \3 C; v( A
用s代表x,y,t的任一项,D海底到自由表面的高度,即
5 O H( t4 q/ d0 L ,可得到如下表达式。
5 t# i( E: N: L( }4 m1 A& @3 N5 G) X. f; n
+ \, N- {" O$ K4 s5 I( W- y8 ^2 Z- y5 w" j 由此推导下去即可得到sigma坐标下的控制方程,推导过程极其繁琐,再此省略了推导的中间过程,直接给出结果。为表示方便,后文sigma坐标系的变量中省略其右上角的星号。若对推导的详细过程感兴趣,见文末参考资料。
! ~* ^ `5 Y3 a O# Y* I
/ J v4 Y) }* C3 T6 V' l3 F7 ?- r2 Q+ v- s& A! I, R
3.内外模态分离 首先,再回顾一下第一篇文章所讨论过的CFL条件,上次是从数学的角度理解CFL条件为什么能确保线性偏微分方程稳定,这次从波动的角度理解一下CFL条件的意义。由于海洋和大气的动力框架系统为高度非线性系统,因此其稳定性变得更加难以控制。CFL条件是一种很好的参考,而无法绝对确保稳定。
- {# b# g$ P( K( I9 q- ~
( H; {. }% I* F4 E
1 P' L, O. D: p5 {) y CFL条件中c的物理意义是波速。假设 ,那么此时 。可以看出网格的步长比和波的传播速度相同,意味着这样的网格分辨率是无法分辨这个波的。而当. L6 ? O9 U8 z! U# @
时,波速比步长比要大,同样是无法解得这个波的运动状态。这样描述或许不够严谨,但是有助于理解CFL条件的物理意义。结合海洋的实际情况来看,在表层的表面重力波的波速约为200-300m/s,而在海洋内部的重力内波波速远小于表面重力波,大概是在5m/s左右。可以看出,海洋内部的运动过程和海洋表层的运动过程时间尺度相差较大,表层明显快于内部。再回看CFL条件,可以看出如果要想同时满足海洋表层和海洋内部的稳定,表层就需要迁就内部。而POM模式的做法是将表层和内部分离。把表层的正压过程和内部的斜压过程分别称为外模态和内模态,分别设置时间步长。 & H" J, K }9 P: y& G! f5 g
先来看外模态(即正压模态),该部分也被称为快过程,时间步长较小。处理方式类似于浅水方程的推导,对其所在区域做垂向积分,忽略了水平扩散项。在sigma坐标系下的方程组如下所示。# I5 ?5 h9 }* c9 r0 _
" e4 a3 P2 U, R; j0 L7 f
: g5 m! ~! o- Y
对于内模态,则方程形式和第二部分列的形式一样。由于外模态时间步长短,内模态时间步长较长。在POM模式中,内模态的时间步长通常是外模态的数十倍。如果将POM模式的整体结构写成伪代码的话,可以写成如下形式。内模态的时间步长是外模态的isplit倍,这样外模态就可以嵌套在内模态的循环里写。+ A9 J7 S* f/ l9 R
program POM, I5 N; U3 Q( W" S# u, [* z$ [4 m
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. E5 n' J$ T* H+ I, T9 c
call advt !计算T和S call proft
. i3 I1 F. F5 u" E9 h4 m% \8 k call dens !计算密度 call advu !计算u call profu
) A1 l/ Q0 f; s: ?8 b. @ call advv !计算v call profv
, A5 }8 p3 Y: ~) ` print !将结果输出 end do !内模态结束end4.湍流闭合方案& J; P% B8 ^' K1 t
! G0 p) Q. `+ K9 ?5 o
通过观察可以发现,本文最开端给出的POM原始方程的运动方程和温盐方程都有 或 。而在这些方程的末尾,也有 , , 或 这些项。这些项的存在使得方程的未知量多于待求解的变量,而如果忽略这些项则会对模拟结果大打折扣。因此,需要解决这些参数的设置问题,而POM模式选择了使用Mellor-Yamada方案,具体形式如下。
# P: I0 u3 g/ D0 ~8 g+ [
6 s! Q# K( m! v+ A- S- c; c
1 n7 Q4 L% g" D7 l% a& u2 x' J6 S4 P' k' _5 Q
湍流是一个十分复杂的现象,如果想理解湍流参数化方案,就需要理解什么是湍流。下文将从湍流的本质讲起,逐渐引出湍流参数化方案的全貌。当对模式的动力框架有了比较明确的理解之后,再去看模式代码甚至修改模式代码,就会容易很多。
. G# {# C" `4 [; x/ D( ]版权声明 本文创作的初衷是用于帮助数值模式的学习者。欢迎转载,转载请私信并注明作者和出处,请勿用于任何商业用途。
4 t; x- W- S* f8 ?! Q7 x- _参考资料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).
( `/ `) M% x* O3 J y/ A) J* i: M L$ o: v z
|