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