0%

MITgcm入门指南

1. MITgcm模式介绍

MITgcm是一个用于研究大气、海洋与气候的数值模式,它是一个基于有限体积法的非静力平衡模式。MITgcm巧妙地通过流体动力学方程的同构,使得其动力内核既可以进行大气运动的模拟,也可以进行海洋流动的模拟。

图1 海洋与大气运动的 z-p 同构

如果对MITgcm模式的方程和数值算法感兴趣,可以参考MITgcm的官方文档,里面对于模式的kernel、参数化方案等有较为详细的介绍。

2. 安装与编译

首先我们需要获取MITgcm的代码,可以通过git进行:

1
git clone https://github.com/MITgcm/MITgcm.git

如果git不可用,也可以直接访问MITgcm的github页面手动下载解压。

下载之后,cd进入MITgcm的根目录,应该包含以下内容:

1
2
3
4
5
6
7
8
9
10
model  #包含主程序源码
eesupp #包含运行环境的源码
pkg #包含各package的源码
doc #包含rst格式的MITgcm文档
tools #包含各种有用的工具,比如genmake2,用于生成makefile的bash脚本
utils #包含各种功能,比如与matlab和python的间交互,utils/scripts包含将不同处理器或pile的输出文件合并的csh脚本
verification #包含样例实验
jobs #包含样例任务脚本
lsopt #Line search优化代码
optim #MITgcm和line search代码的接口

接下来我们需要选定一个想要进行的实验,编译将在此基础上进行。所有的样例实验都存放在verification目录下,本文将以tutorial_baroclinic_gyre为例进行演示,进行其他实验的步骤类似。我们首先进入目标实验的编译目录:

1
cd verification/tutorial_baroclinic_gyre/build

接着,需要生成Makefile:

1
../../../tools/genmake2 -mods ../code -optfile ../../../tools/build_options/linux_amd64_gfortran

其中-mods命令告诉genmake使用../code中的内容覆盖模型默认的源码,包括SIZE.h(配置模式的区域大小,后面介绍);-optfile命令告诉genmake在执行期间运行指定的optfile,这是一个bash脚本,包含环境变量、路径、编译器选项等,MITgcm提供了非常多预置的optfile,位于tools/build_options目录下,可以根据自己的机器进行选择。

生成Makefile之后,我们通过以下命令创建依赖关系:

1
make depend

最后,可以使用以下命令编译代码:

1
make

如果成功编译,会生成一个名为mitgcmuv的可执行文件,因为编译后build目录内文件较多,可以使用stat mitgcmuv来方便地检查编译是否成功进行。

3. 实例演示:以tutorial_baroclinic_gyre为例

3.1 试运行

我们从MITgcm的根目录进入运行目录:

1
cd run

目前这是一个空目录,我们需要将可执行文件和输入文件放置进来:

1
2
cp ../build/mitgcmuv .
cp ../input/* .

这样,我们就做好了运行前的全部准备,可以直接运行

1
./mitgcmuv

样例实验通常默认了非常少的时间步数,运行后立刻可以得到结果,可以将输出与../results/output.txt进行对比,如果结果相同,则说明程序执行正常,我们接着可以开始对实验参数进行更具体的配置。

图2 baroclinic ocean gyre

本文运行的实验是Baroclinic Ocean Gyre,模拟风和浮力强迫下的斜压双涡漩海洋环流,位于verification/tutorial_baroclinic_gyre目录下。本实验使用球坐标网格并设有15个垂直层,横跨热带到中纬度,并受到恒定的东西风应力强迫,并且风应力的强度在南北方向上成正弦变化。

3.2 编译配置(Compile-time Configuration)

接下来我们对实验中的参数进行配置,配置分为两部分:编译配置和运行配置,编译配置需要在编译前进行,每次修改之后都需要重新编译,而运行配置则不需要。

3.2.1 code/packages.conf
1
2
3
4
# packages.conf文件应该包含以下内容:
gfd
diagnostics
mnc

在这个文件中我们指定了我们希望包含的库。gfd是一个预定义的库组,包含了大多数设置所需要的标准库。mnc库用于将输出以netCDF格式存储;diagnostics库允许用户从用于模型诊断的众多变量中选择输出,并指定输出频率,提供多种时间平均和快照选项,如果没有启用这个库,输出将局限于少量的快照。

3.2.2 code/SIZE.h

这个文件指定了模型的几何。以下代码计算了整个模型域的水平大小:

1
2
&           Nx  = sNx*nSx*nPx,
& Ny = sNy*nSy*nPy,

整个模型域的大小是$62\times 62$,我们在这里将模型在水平方向上划分成四个区块,每个区块的大小为$31\times 31$,也就是

1
2
3
4
5
&           sNx =  31,
& sNy = 31,

& nSx = 2,
& nSy = 2,

为了正确地处理边界条件,模型的区块需要具有一定的重叠区域(overlap),MITgcm许可的最小的overlap是2,对于更加复杂的对流方案或参数化方案可能需要更大的overlap,我们在这里保持2即可:

1
2
&           OLx =   2,
& OLy = 2,

以下代码设置了nPxnPy,分别表示在$x$和$y$方向上使用的进程数量,我们暂时使用单进程,因此这两个参数设置为1:

1
2
&           nPx =   1,
& nPy = 1,

最后,我们指定模型具有15个垂直层:

1
&           Nr  =   15)

修改完上述参数后,别忘记重新编译。

3.3 运行配置(Run-time Configuration)

input/data文件里包含了模型的大部分参数,可以根据自己的需要进行修改。

第一部分主要是物理参数设置:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# Model parameters
# Continuous equation parameters
&PARM01
viscAh=5000., #水平粘度系数
viscAr=1.E-2, #竖直粘度系数
no_slip_sides=.TRUE., #侧边界无滑移边界条件
no_slip_bottom=.FALSE., #底部无滑移边界条件
diffKhT=1000., #水平温度扩散系数
diffKrT=1.E-5, #竖直温度扩散系数
#静力不稳定混合参数化
ivdc_kappa=1.,
implicitDiffusion=.TRUE.,
eosType='LINEAR', #使用线性状态方程
tRef=30.,27.,24.,21.,18.,15.,13.,11.,9.,7.,6.,5.,4.,3.,2., #参考温度
tAlpha=2.E-4, #热膨胀系数
sBeta=0., #盐度收缩系数
rhoNil=999.8, #表面参考密度
gravity=9.81, #重力加速度
# 自由表面公式
rigidLid=.FALSE.,
implicitFreeSurface=.TRUE.,
exactConserv=.TRUE.,
saltStepping=.FALSE., #盐度时间积分
# globalFiles=.TRUE.,
&

另一部分比较重要的是对时间参数的设置:

1
2
3
4
5
6
7
8
9
startTime=0.,  #开始时间
endTime=12000., #结束时间
deltaT=1200., #时间步长
pChkptFreq=622080000., #检查点文件保存频率
chkptFreq=155520000., #检查点文件保存频率
dumpFreq=31104000., #快照输出频率
#监视器输出选项
monitorFreq=1200.,
monitorSelect=2,

input/bathy.bininput/windx_cosy.bin文件称分别指定了海底的地形以及风应力,这两个文件可以通过input/gendata.minput/gendata.py生成。

完成对所有参数的的修改之后,就可以正式进行模式的运行了,可以将程序提交后台,并将日志输出到文件中:

1
./mitgcmuv > output.log &

3.4 输出格式

程序运行结束之后,可以查看日志文件的末尾,如果显示

1
PROGRAM MAIN: Execution ended Normally

则说明程序顺利完成。此时run目录下会有很多pickup.*文件,这些文件可以用来restart,也就是从某个时间点开始运行,这样即使我们的程序中途退出了,我们也不需要再次从$t=0$开始运行。

由于我们开启了mnc功能,我们可以获得netCDF格式的输出文件,这些文件被存放在run/mnc_test_[iter]子目录下,iter从0001开始,每次重新运行程序都会递增1。

由于我们在运行的时候将模型分为了4个区块,因此数据也是分4个区块存储的,MITgcm提供了便捷的工具将结果合并为完整的netCDF文件,我们可以在运行目录run下执行:

1
2
3
4
5
6
7
8
% ln -s ../../../utils/python/MITgcmutils/scripts/gluemncbig .
% ./gluemncbig -o grid.nc mnc_test_*/grid.t*.nc
% ./gluemncbig -o state.nc mnc_test_*/state*.t*.nc
% ./gluemncbig -o dynDiag.nc mnc_test_*/dynDiag*.t*.nc
% ./gluemncbig -o surfDiag.nc mnc_test_*/surfDiag*.t*.nc
% ./gluemncbig -o phiHyd.nc mnc_test_*/phiHyd*.t*.nc
% ./gluemncbig -o phiHydLow.nc mnc_test_*/phiHydLow*.t*.nc
% ln -s mnc_test_0001/dynStDiag.0000000000.t001.nc dynStDiag.nc

主要的netCDF输出文件包括:

  • grid.nc:包含模型网格变量
  • state.nc:包含状态变量($U, V, W, Temp, S, Eta$)的快照
  • surfDig.nc:包含data.diagnostics的列表1指定的二维诊断量(ETAN, TRELAX, MXDEPTH等)
  • dynDiag.c:包含列表2指定的三维诊断量(THETA, PHIHYD, UVEL, VVEL, WVEL等)
  • dynStDiag.nc:包含统计-动力学诊断量
  • phiHyd.nc, phiHydLow.nc:分别为静水压势异常的三维快照场和底部静水压势异常的二维快照场

4. 结果分析

我们将模型运行100年,查看系统达到稳态之后的状态。由于系统运行所需的时间较长,且占用存储空间较大,我们在本文中就不改变参数进行实验了,我们将对单次运行之后得到的结果进行分析。

图3 不同深度处的平均位温
图4 不同深度处的位温标准差

模型表层的温度最开始为$30^\circ C$,在外界的forcing下表层温度迅速降低,从图3中可以看到在最开始的5年内,表层温度差不多就达到了平衡,而深度的海水温度变化比较缓慢。从图4中可以看到,表层海水的梯度较大,水平方向上分布不均匀,而深层海水的水平方向比较均匀。

图5 海平面的高度分布
图6 温度弛豫的分布

模型在运行到100年的时候已经很接近稳态了,我们取这一年的年平均数据进行分析。图5显示了海平面的高度,可以看到在风应力和科里奥利力的作用下,海平面的高度呈现南高北低。图6显示了温度弛豫的分布,我们可以看到海水与外界的热量交换在水平方向上是怎么分布的,表层海水在内部输送和表面交换的作用下达到热量平衡。

图7 位温在400m深度的水平分布
图8 位温在北纬$33.5^\circ$的分布

我们还可以观察海水内部的温度是如何分布的,图7和图8分别显示了在400m深度和北纬$33.5^\circ$处的温度分布。

5. 总结

本文介绍了MITgcm模式的安装、编译与运行过程。MITgcm功能强大、使用方便,如果想要更深入地学习,可以访问官方文档,里面不仅包含对本文所有内容的介绍,还有对所有tutorial实验的讲解,以及更加灵活的物理过程和参数化方案的配置。