热带太平洋的1.5层C网格海洋模式数值模拟

一般而言我们认为，海洋中的各个要素除了在水平层面上有较大的分布差异外，在垂直方向上也存在一定的分布差异，特别是在温跃层上下存在明显的差异。但是另一方面，海洋要素在垂直方向上的分布随时间的变化主要存在于温跃层以浅。因此海洋可以简化为由温、盐、流、表面高度剧烈变化的活跃的上层，与温、盐基本不变，流速极小的被动的下层两个部分组成。这两个部的分界面就是温跃层深度面。一般我们将这样的抽象海洋模式称为“1.5层海洋模式”。通常1.5层海洋模式被应用于温跃层深度、风生环流的模拟等领域。

1.5层海洋模式也是由最基础的不可压缩流体的动力学方程组得到的，这个方程组又是由纳维-斯托克斯方程（N-S公式）以及连续性方程组成。再引入静力平衡的假设，并对方程组进行垂直积分，得到最终的方程。接下来就要进行差分运算了。

本文采用的C网格是将xy方向上的速度场和海面高度场交错放置的一种形式（图 1）。与传统的重叠放置相比，C网格最大的优势就在于空间差分天然地为中央差分格式，有效避免了因为空间差分而导致发散的情况。

图1  C网格xy方向速度场与海面高度场交错配置示意

在空间上采用C网格形式配置的同时，为保证时间上的稳定，这里也采用中央差分格式，即蛙跃格式。由于蛙跃格式需要起步项，本文采用了前差起步（欧拉格式）。

此外还需要设置其他一些常量。首先将区域限定在30N-30S的热带太平洋，经向剖分为171个格点，纬向为61个格点，空间步长为10^5米；时间上假设步长为5400秒。初始速度场、混合层厚度变化量均为0。外部强迫场仅考虑12个月的气候态风应力强迫，由于其时间步长为1个月，与模式中的时间步长不符，需插分到每一个的时间格点上。科氏力参数本文取f平面近似，其中beta取2*10^-11；边界为无滑侧边界。约化重力g'取9*10^-2，水平湍流粘性系数Ah取2*10^4，平均海水密度rho取1025。陆地全部设为干单元，海洋全部为湿单元。

---

以下是运算的主程序

close all; clear; clc;
load( 'E:\海洋数值模拟\实习\pacific_winds.mat' );
load( 'E:\海洋数值模拟\实习\pacific_mask.mat' );
delta_x = 10 ^ 5;
delta_y = 10 ^ 5;
g = - 9 * 10 ^ ( - 2 );
A_h = 2 * 10 ^ ( 4 );
rho = 1025;
U =