前言
随着研究的更加深入，如水汽收支方程这样简单的方程已难以满足我们的研究需求，我们将向着更复杂的数学问题迈进。
气象学和海洋学研究通常会遇到需要数值求解的反演问题。经典反演问题之一是求解流函数给定涡度的垂直分量和适当的边界条件：
如今， xarray成为大数据地球科学中常用的数据结构。由于整个 4D 数据以及坐标信息都被组合到xarray中，因此解决反演问题变得非常简单，唯一的输入只是xarray.DataArray的一个。在球形地球上的反演，就像一些气象问题一样，可以利用像Windspharm这样的球谐函数(或球面谐波，spherical harmonics) ,使用 快速傅里叶变化(FFT) 比这里使用的逐次超松弛法(SOR) 更有效。然而，在海洋的情况下，在存在陆地/海洋掩模的情况下，SOR方法绝对是更好的选择。
Tips：
逐次超松弛法 (Successive Overrelaxation Method，SOR) 是一种求解线性方程组的方法，通过外推 Gauss-Seidel 方法得出。该外推法采用前一个迭代与每个分量连续计算的Gauss-Seidel迭代之间的加权平均值：
表示 Gauss-Seidel 迭代， 是松弛因子。选择一个值将加快迭代收敛到解的速度。
>1，称为逐次超松弛迭代。 =1，退化为Gauss-Seidel 迭代。 <1，低松弛迭代。
这个方法是众所周知的。
Xinvert
有了前面的介绍，求解椭圆型偏微分方程就十分简单了，xinvert是一个强大到不可思议的Python库，旨在求解以下抽象方程：
是二阶偏微分算子, 是要反演的未知函数, 是强迫函数。
对于二维情况，方程的一般形式：(1) 为：
其中系数 都是已知变量。当条件在定义域内处处满足，则上式为椭圆型方程。在这种情况下，可以使用SOR迭代方法进行反演。当或 ，它是一个抛物线或双曲线方程(高中内容，忘记了的回去狠狠地复习一百万遍)。无论哪种情况，SOR 都无法收敛到最终解。
使用SOR求解的一个问题是，在 Python 中使用显式循环的迭代速度将非常……慢！这里一个非常合适的解决方案是使用numba。或使用GPU来加快求解速度。
有时方程的一般形式。 (2) 可以转化为标准形式（即标准化）：
在这种情况下， 应满足以确保其椭圆度。椭圆条件在不同的问题中有其自己的物理意义。也就是说，系统处于稳定状态，对任何小的扰动都是稳定的。
气象学和海洋学中的许多问题都可以转化为方程（1）的形式。 (2) 或等式。 （3）。然而，其中一些是在3D情况下的（如准地转Omega方程）：
或在四阶情况下（Munk模型）：
因此，我们实现了四个基本求解器来考虑上述四个方程。 (2-5)或更多。如果问题不属于这四种类型中的一种，将为此类问题添加一个求解器。但是我们并不希望是这样，因为我们希望使求解器尽可能最小和通用。
安装
通过 pip 安装
pip install xinvert
从github安装
git clone https://github.com/miniufo/xinvert.git
cd xinvert
python setup.py install
经典的物理问题
这里我们将气象学和海洋学中的所有反演问题总结为下表。如果发现更多符合式（1）的问题，则可以进一步扩展该表。
迭代收敛过程
SOR迭代的整个收敛过程可以看成：
from xinvert import animate_iteration
# output has 1 more dimension (iter) than input, which could be animated over.
# Here 40 frames and loop 1 per frame (final state is after 40 iterations) is used.
psi = animate_iteration(invert_Poisson, vor, iParams=iParams,
loop_per_frame=1, max_frames=40)
绘图为了方便，我们使用proplot库来绘制，其作者是来自英国东英吉利大学的博士后研究员Luke Davis。
下面我们将拿出一个大气领域一个海洋领域的经典方程来反演，大气领域我们选择Gill-Matsuno模型，海洋领域选择Stommel-Munk模型，这两个模型也是教科书上的经典例子。
Gill-Matsuno模型
1.简介
Gill-Matsuno 模型（Matsuno 1966 ；Gill 1980）是研究热带大气对指定热源响应的经典模型。模型的反演意味着在给定任何类型的加热模式的情况下，我们可以获得大气流量和质量方面的稳态响应。球坐标上的模型可以写为：
, , 和 是相对于参考状态的异常风和质量响应, , 和 .和假设为零（静止参考状态）。给定非绝热加热场，就会得到风场(, ) 和质量场 。
2.理论
Gill-Matsuno 模型关于三个未知数是线性的, (水平风场), 和 (位势)。使用简单的替换，可以获得它们的显式表达式：
其中，, , 和 ， 方程（4）是一个椭圆方程，只要 , , 和 都为正且 不远小于 . 以分量形式，方程(4)可以写成:
可以使用 SOR 方法在给定 的情况下反演，从而利用方程 (5) 和 (6) 计算风场 和 。该方法可以稍作修改，以适应一般的二维求解器，具体如下：
虽然方程 (7) 和 (8) 是等价的，但目前尚不清楚哪种形式更适合离散化和迭代。在这里，我们选择方程 (8) 的形式用于 SOR 方法。
3.Python Example
import sys
sys.path.append('../../../')
import numpy as np
import xarray as xr
lon = xr.DataArray(np.linspace(0, 360, 144), dims='lon', coords={'lon':np.linspace(0, 360, 144)})
lat = xr.DataArray(np.linspace(-90, 90, 73), dims='lat', coords={'lat':np.linspace(-90, 90, 73)})
lat, lon = xr.broadcast(lat, lon)
# three patterns of heating
Q1 = 0.05*np.exp(-((lat-0)**2+(lon-120)**2)/100.0)
Q2 = 0.05*np.exp(-((lat-10)**2+(lon-120)**2)/100.0) \
- 0.05*np.exp(-((lat+10)**2+(lon-120)**2)/100.0)
Q3 = 0.05*np.exp(-((lat-10)**2+(lon-120)**2)/100.0)
反演在纬度/经度平面上以质量 和风场 表示的大气响应非常简单：
from xinvert import invert_GillMatsuno, cal_flow
iParams = {
'BCs' : ['fixed', 'periodic'],
'mxLoop' : 600,
'tolerance': 1e-5,
'optArg' : 1.4,
}
mParams = {
'epsilon': 1e-5,
'Phi' : 5000,
}
h1 = invert_GillMatsuno(Q1, dims=['lat','lon'], mParams=mParams, iParams=iParams)
h2 = invert_GillMatsuno(Q2, dims=['lat','lon'], mParams=mParams, iParams=iParams)
h3 = invert_GillMatsuno(Q3, dims=['lat','lon'], mParams=mParams, iParams=iParams)
u1, v1 = cal_flow(h1, dims=['lat','lon'], BCs=['fixed','periodic'], mParams=mParams, vtype='GillMatsuno')
u2, v2 = cal_flow(h2, dims=['lat','lon'], BCs=['fixed','periodic'], mParams=mParams, vtype='GillMatsuno')
u3, v3 = cal_flow(h3, dims=['lat','lon'], BCs=['fixed','periodic'], mParams=mParams, vtype='GillMatsuno')
请注意，这里演示的是一个全球的情景，与原始的 平面情况不同。结果可以可视化为：
#%% plot wind and streamfunction
import proplot as pplt
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
def maskout_cmap(cmap, clevs, maskout):
for msk in maskout:
if msk not in clevs:
raise Exception('maskout {0} not in clevs {1}'.format(maskout, clevs))
N = len(clevs) + 1
colors = plt.cm.get_cmap(cmap, N) # define the colormap (even number)
colors = [colors(i) for i in range(N)]
if type(clevs) is list:
clevsLst = clevs
else:
clevsLst = clevs.tolist()
strIdx = clevsLst.index(maskout[0]) + 1
endIdx = clevsLst.index(maskout[1]) + 1
for i in range(strIdx, endIdx):
# print(colors[i])
colors[i] = np.array((1, 1, 1, 1))
if len(colors) != N:
print(len(colors))
print(N)
raise Exception('lengths not equal')
# create the new map
return mcolors.LinearSegmentedColormap.from_list('msk', colors, N)
lat, lon = xr.broadcast(u1.lat, u1.lon)
fig, axes = pplt.subplots(nrows=3, ncols=1, figsize=(7, 6.9), sharex=3, sharey=3,
proj=pplt.Proj('cyl', lon_0=180))
skip = 1
fontsize = 11
cmap = maskout_cmap('bwr', [-0.05, -0.04, -0.03, -0.02, -0.01, 0.,
0.01, 0.02, 0.03, 0.04, 0.05], [-0.01, 0.01])
ax = axes[0]
ax.contourf(Q1, cmap=cmap, levels=np.linspace(-0.05, 0.05, 12))
ax.contour(h1, cmap='jet')
ax.set_title('Gill-Matsuno response - type 1 ($\epsilon=10^{-5}, \Phi=5\\times 10^3$)', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip+1], lat.values[::skip,::skip+1],
u1.values[::skip,::skip+1], v1.values[::skip,::skip+1],
width=0.0012, headwidth=10., headlength=12., scale=250)
# headwidth=1, headlength=3, width=0.002)
ax.set_ylim([-40, 40])
ax.set_xlim([-150, 110])
ax = axes[1]
ax.contourf(Q2, cmap=cmap, levels=np.linspace(-0.05, 0.05, 12))
ax.contour(h2, cmap='jet')
ax.set_title('Gill-Matsuno response - type 2 ($\epsilon=10^{-5}, \Phi=5\\times 10^3$)', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip+1], lat.values[::skip,::skip+1],
u2.values[::skip,::skip+1], v2.values[::skip,::skip+1],
width=0.0012, headwidth=10., headlength=12., scale=250)
ax.set_ylim([-40, 40])
ax.set_xlim([-150, 110])
ax = axes[2]
p=ax.contourf(Q3, cmap=cmap, levels=np.linspace(-0.05, 0.05, 12))
ax.contour(h3, cmap='jet')
ax.set_title('Gill-Matsuno response - type 3 ($\epsilon=10^{-5}, \Phi=5\\times 10^3$)', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip+1], lat.values[::skip,::skip+1],
u3.values[::skip,::skip+1], v3.values[::skip,::skip+1],
width=0.0012, headwidth=10., headlength=12., scale=250)
ax.set_ylim([-40, 40])
ax.set_xlim([-150, 110])
fig.colorbar(p, loc='r', label='', rows=(1,3), ticks=0.01)
axes.format(abc='(a)', coast=True, lonlines=40, latlines=15, lonlabels='b', latlabels='l',
grid=False, labels=False)
在这种情况下，加热模式是从30-60天带通滤波的OLR（向外长波辐射）异常中得到的。我们只是想看看基于Gill-Matsuno模型的大气响应是否接近观测结果。从数据文件中读取数据：
import numpy as np
import xarray as xr
ds = xr.open_dataset('../../../Data/MJO.nc')
# observed h and u,v anomalies (30-60day filtered)
h_ob = ds.hl
u_ob = ds.ul
v_ob = ds.vl
# observed heating anomaly (30-60day filtered)
Q = (ds.ol*-0.0015).where(np.abs(ds.lat)<60, 0)
iParams = {
'BCs' : ['fixed', 'periodic'],
'mxLoop' : 800,
'tolerance': 1e-5,
'optArg' : 1.4,
}
mParams1 = {'epsilon': 1e-5, 'Phi': 5000}
mParams2 = {'epsilon': 7e-6, 'Phi': 8000}
mParams3 = {'epsilon': 7e-6, 'Phi': 10000}
h1 = invert_GillMatsuno(Q, dims=['lat','lon'], mParams=mParams1, iParams=iParams)
h2 = invert_GillMatsuno(Q, dims=['lat','lon'], mParams=mParams2, iParams=iParams)
h3 = invert_GillMatsuno(Q, dims=['lat','lon'], mParams=mParams3, iParams=iParams)
u1, v1 = cal_flow(h1, dims=['lat','lon'], BCs=['fixed','periodic'], mParams=mParams1, vtype='GillMatsuno')
u2, v2 = cal_flow(h2, dims=['lat','lon'], BCs=['fixed','periodic'], mParams=mParams2, vtype='GillMatsuno')
u3, v3 = cal_flow(h3, dims=['lat','lon'], BCs=['fixed','periodic'], mParams=mParams3, vtype='GillMatsuno')
lat, lon = xr.broadcast(ds.lat, ds.lon)
fig, axes = pplt.subplots(nrows=2, ncols=2, figsize=(7, 4.4), sharex=3, sharey=3,
proj=pplt.Proj('cyl', lon_0=180))
skip = 1
fontsize = 11
cmap = maskout_cmap('bwr', [-0.1, -0.09, -0.08, -0.07, -0.06, -0.05,
-0.04, -0.03, -0.02, -0.01, 0., 0.01, 0.02,
0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1],
[-0.01, 0.01])
ax = axes[0,0]
ax.contourf(Q, cmap=cmap, levels=np.linspace(-0.1, 0.1, 21))
ax.contour(h_ob, color='black', linewidth=1.2,
levels=[-25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25])
ax.set_title('observed mass and wind anomalies', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip+1], lat.values[::skip,::skip+1],
u_ob.values[::skip,::skip+1]*5, v_ob.values[::skip,::skip+1]*5,
width=0.0016, headwidth=10., headlength=12., scale=300)
# headwidth=1, headlength=3, width=0.002)
ax.set_ylim([-30, 30])
ax.set_xlim([-160, -20])
ax = axes[0,1]
ax.contourf(Q, cmap=cmap, levels=np.linspace(-0.1, 0.1, 21))
ax.contour(h1, color='black', linewidth=1.2, levels=11)
ax.set_title('response ($\epsilon=10^{-5}, \Phi=5\\times 10^3$)', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip+1], lat.values[::skip,::skip+1],
u1.values[::skip,::skip+1], v1.values[::skip,::skip+1],
width=0.0016, headwidth=10., headlength=12., scale=300)
# headwidth=1, headlength=3, width=0.002)
ax.set_ylim([-30, 30])
ax.set_xlim([-160, -20])
ax = axes[1,0]
ax.contourf(Q, cmap=cmap, levels=np.linspace(-0.1, 0.1, 21))
ax.contour(h2, color='black', linewidth=1.2, levels=11)
ax.set_title('response ($\epsilon=7\\times 10^{-6}, \Phi=5\\times 10^3$)', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip+1], lat.values[::skip,::skip+1],
u2.values[::skip,::skip+1], v2.values[::skip,::skip+1],
width=0.0016, headwidth=10., headlength=12., scale=300)
ax.set_ylim([-30, 30])
ax.set_xlim([-160, -20])
ax = axes[1,1]
p=ax.contourf(Q, cmap=cmap, levels=np.linspace(-0.1, 0.1, 21))
ax.contour(h3, color='black', linewidth=1.2, levels=11)
ax.set_title('response ($\epsilon=7\\times 10^{-6}, \Phi=10^4$)', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip+1], lat.values[::skip,::skip+1],
u3.values[::skip,::skip+1], v3.values[::skip,::skip+1],
width=0.0016, headwidth=10., headlength=12., scale=300)
ax.set_ylim([-30, 30])
ax.set_xlim([-160, -20])
fig.colorbar(p, loc='b', label='heating Q', ticks=0.02, length=1)
axes.format(abc='(a)', coast=True, grid=False, labels=False,
lonlines=20, latlines=10, lonlabels='b', latlabels='l')
很明显，这样的简单模型已经再现了大部分观测到的信号，尽管响应中心存在一些偏移。
Stommel-Munk 模型
1.简介
经典的风驱动海洋环流模型是 Stommel-Munk 模型（Stommel 1948; Munk 1950）。该模型研究的是表层海洋环流的水平流函数 对风应力旋度 的响应（其中），可以写为：
如果 且 ，那么方程 (1) 退化为 Stommel (1948) 模型，这是一个关于 的二阶偏微分方程；如果 且 ，方程 (1) 退化为 Munk (1950) 模型，这是一个四阶偏微分方程。因此，方程 (1) 也被称为 Stommel-Munk 模型。给定风应力旋度 或风应力 本身，可以使用 SOR 迭代反演水平流函数。
2.理论
这里简要回顾了 Stommel-Munk 模型的推导过程。起点是以矢量形式表示的水平动量方程：
其中 是应力。引入 运算符作为二维向量的逆时针旋转（例如，如果 ，那么 ）。对方程 (2) 的两边取旋度 可以得到涡度方程 ：
假设处于稳态（去掉第一项）、线性（去掉第二项）且无辐散（去掉第三项），则可以得到：
假设在混合层内流动是垂直均匀的，并且在表面 ，在混合层底部 （线性阻力或牛顿冷却），则方程 (4) 变为：
经过一些整理，并引入流函数 ，它就成为了 Stommel-Munk 模型：
或者等价于：
这里， 是流函数的双调和扩散系数（即动量的调和扩散系数）， 是线性阻尼系数， 是混合层深度， 是海水的常数参考密度。虽然方程 (6) 和 (7) 是等价的，但目前尚不清楚哪种形式更适合离散化和迭代。在这里，我们选择方程 (6) 的形式用于 SOR 方法。
3.Python Example
3.1 经典的Stommel例子
Stommel (1948) 已经证明，给定一个随经度变化的风应力分布，在海盆中部的应力为零的情况下，可以在 平面上得到对称的海洋环流，而在 平面上则会得到西边界强化的环流。以下代码可以再现这一现象。
### classical cases ###
import sys
sys.path.append('../../../')
import numpy as np
import xarray as xr
xnum = 201
ynum = 151
Lx = 1e7 # 10,000 km
Ly = 2 * np.pi * 1e6 # 6249 km
R = 0.0008 # Rayleigh friction
depth = 200 # fluid depth 200
beta = 1.8e-11
F = 0.3
xdef = xr.DataArray(np.linspace(0, Lx, xnum), dims='xdef', coords={'xdef':np.linspace(0, Lx, xnum)})
ydef = xr.DataArray(np.linspace(0, Ly, ynum), dims='ydef', coords={'ydef':np.linspace(0, Ly, ynum)})
ygrid, xgrid = xr.broadcast(ydef, xdef)
tau_ideal = xr.DataArray(-F * np.cos(np.pi * ygrid / Ly), dims=['ydef','xdef'], coords={'ydef':ydef, 'xdef':xdef})
# finite difference for curl
curl_tau = - tau_ideal.differentiate('ydef')
在纬度/经度平面内，以流函数 和风场 表示的海洋环流反演过程如下：
from xinvert import invert_Stommel, cal_flow
import proplot as pplt
iParams = {
'BCs' : ['fixed', 'fixed'],
'mxLoop' : 3000,
'tolerance': 1e-8,
'optArg' : 1.9,
'undef' : np.nan,
}
mParams1 = {'beta': 0 , 'R': R, 'D': depth}
mParams2 = {'beta': beta, 'R': R, 'D': depth}
S1 = invert_Stommel(curl_tau, dims=['ydef','xdef'], coords='cartesian', mParams=mParams1, iParams=iParams)
S2 = invert_Stommel(curl_tau, dims=['ydef','xdef'], coords='cartesian', mParams=mParams2, iParams=iParams)
u1, v1 = cal_flow(S1, dims=['ydef','xdef'], coords='cartesian')
u2, v2 = cal_flow(S2, dims=['ydef','xdef'], coords='cartesian')
array = [
[1, 2, 2, 3, 3,],
]
fig, axes = pplt.subplots(array, figsize=(7,3.5), sharex=0, sharey=3)
skip = 3
fontsize = 11
ax = axes[0]
ax.plot(tau_ideal[:,0], tau_ideal.ydef)
ax.plot(tau_ideal[:,0]-tau_ideal[:,0], tau_ideal.ydef, color='k', linewidth=0.3)
ax.set_ylabel('y-coordinate (m)', fontsize=fontsize-2)
ax.set_title('wind stress', fontsize=fontsize)
ax = axes[1]
m=ax.contourf(S1/1e6*depth, cmap='rainbow', levels=np.linspace(0, 140, 15))
ax.set_title('$f$-plane Stommel solution', fontsize=fontsize)
p=ax.quiver(xgrid.values[::skip+1,::skip], ygrid.values[::skip+1,::skip],
u1.values[::skip+1,::skip], v1.values[::skip+1,::skip],
width=0.0014, headwidth=8., headlength=12., scale=10)
ax.colorbar(m, loc='b')
ax.set_ylim([0, Ly])
ax.set_xlim([0, Lx])
ax.set_xlabel('x-coordinate (m)', fontsize=fontsize-2)
ax.set_ylabel('y-coordinate (m)', fontsize=fontsize-2)
ax = axes[2]
m=ax.contourf(S2/1e6*depth, cmap='rainbow', levels=np.linspace(0, 90, 19))
ax.set_title('$\\beta$-plane Stommel solution', fontsize=fontsize)
ax.quiver(xgrid.values[::skip+1,::skip], ygrid.values[::skip+1,::skip],
u2.values[::skip+1,::skip], v2.values[::skip+1,::skip],
width=0.0014, headwidth=8., headlength=12., scale=8)
# headwidth=1, headlength=3, width=0.002)
ax.colorbar(m, loc='b')
ax.set_ylim([0, Ly])
ax.set_xlim([0, Lx])
ax.set_xlabel('x-coordinate (m)', fontsize=fontsize-2)
ax.set_ylabel('y-coordinate (m)', fontsize=fontsize-2)
axes.format(abc='(a)', grid=False, ylabel='y-coordinate (m)')
很明显，在 平面上，环流得到了增强。
3.2 经典的Munk例子
Munk (1950) 通过添加侧向扩散项并去除底部线性阻力摩擦，修改了 Stommel 的模型。我们将通过添加一个双调和系数 并将底部线性阻力系数 设为零，使用 Munk 模型来再现上述示例。这很简单：
from xinvert import invert_StommelMunk, cal_flow
iParams = {
'BCs' : ['fixed', 'fixed'],
'mxLoop' : 4000,
'tolerance': 1e-8,
'optArg' : 1.0,
}
mParams1 = {'R': 0, 'D': depth, 'A4':1e4, 'beta':0}
mParams2 = {'R': 0, 'D': depth, 'A4':1e4, 'beta':beta}
h1 = invert_StommelMunk(curl_tau, dims=['ydef','xdef'], coords='cartesian', mParams=mParams1, iParams=iParams)
h2 = invert_StommelMunk(curl_tau, dims=['ydef','xdef'], coords='cartesian', mParams=mParams2, iParams=iParams)
u1, v1 = cal_flow(h1, dims=['ydef','xdef'], coords='cartesian')
u2, v2 = cal_flow(h2, dims=['ydef','xdef'], coords='cartesian')
array = [[1, 2, 2, 3, 3,],]
fig, axes = pplt.subplots(array, figsize=(7,3.5), sharex=0, sharey=3)
skip = 3
fontsize = 11
ax = axes[0]
ax.plot(tau_ideal[:,0], tau_ideal.ydef)
ax.plot(tau_ideal[:,0]-tau_ideal[:,0], tau_ideal.ydef, color='k', linewidth=0.3)
ax.set_ylabel('y-coordinate (m)', fontsize=fontsize-2)
ax.set_title('wind stress', fontsize=fontsize)
ax = axes[1]
m=ax.contourf(h1/1e6*depth, cmap='rainbow', levels=np.linspace(0, 50, 11))
ax.set_title('$f$-plane Munk solution', fontsize=fontsize)
p=ax.quiver(xgrid.values[::skip+1,::skip], ygrid.values[::skip+1,::skip],
u1.values[::skip+1,::skip], v1.values[::skip+1,::skip],
width=0.0014, headwidth=8., headlength=12., scale=4)
ax.colorbar(m, loc='b')
ax.set_ylim([0, Ly])
ax.set_xlim([0, Lx])
ax.set_xlabel('x-coordinate (m)', fontsize=fontsize-2)
ax.set_ylabel('y-coordinate (m)', fontsize=fontsize-2)
ax = axes[2]
m=ax.contourf(h2/1e6*depth, cmap='rainbow', levels=np.linspace(0, 50, 11))
ax.set_title('$\\beta$-plane Munk solution', fontsize=fontsize)
ax.quiver(xgrid.values[::skip+1,::skip], ygrid.values[::skip+1,::skip],
u2.values[::skip+1,::skip], v2.values[::skip+1,::skip],
width=0.0014, headwidth=8., headlength=12., scale=4)
# headwidth=1, headlength=3, width=0.002)
ax.colorbar(m, loc='b')
ax.set_ylim([0, Ly])
ax.set_xlim([0, Lx])
ax.set_xlabel('x-coordinate (m)', fontsize=fontsize-2)
ax.set_ylabel('y-coordinate (m)', fontsize=fontsize-2)
axes.format(abc='(a)', grid=False, ylabel='y-coordinate (m)')
3.3 使用完整 Stommel-Munk 模型
现在我们要进入一个更接近真实的案例：使用气候学风应力和地形数据来计算全球海洋中的风驱动环流。首先加载数据，这里我们使用SODA的风应力旋度数据：
#%% real cases
import xarray as xr
# load in wind stress from SODA product
ds = xr.open_dataset('../../../Data/SODA_curl.nc')
curl = ds.curl
使用 Stommel-Munk 模型进行反演：
from xinvert import invert_StommelMunk, cal_flow
iParams = {
'BCs' : ['fixed', 'periodic'],
'mxLoop' : 4000,
'tolerance': 1e-14,
'optArg' : 1.0,
'undef' : np.nan,
}
mParams = {'R': 1e-4, 'D': 100, 'A4':3e3, 'beta':beta, 'rho0':1027}
# inversion for streamfunction and flow vector
h1 = invert_StommelMunk(curl[0], dims=['lat','lon'], mParams=mParams, iParams=iParams)
u1, v1 = cal_flow(h1, dims=['lat','lon'], BCs=['extend','periodic'])
lat, lon = xr.broadcast(u1.lat, u1.lon)
fig, axes = pplt.subplots(nrows=1, ncols=1, figsize=(7, 3.5), sharex=3, sharey=3, proj=pplt.Proj('cyl', lon_0=180))
skip = 3
fontsize = 12
ax = axes[0]
p=ax.contourf(h1/1e6*100, cmap='bwr', levels=np.linspace(-80,80,33))
ax.set_title('Stommel-Munk solution to wind stress curl', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip], lat.values[::skip,::skip],
u1.values[::skip,::skip], v1.values[::skip,::skip],
width=0.0012, headwidth=10., headlength=12., scale=70)
ax.set_ylim([-50, 75])
ax.set_xlim([-180, 180])
ax.colorbar(p, loc='b', label='', ticks=15, length=1)
axes.format(abc='(a)', coast=True, lonlines=40, latlines=15, lonlabels='b', latlabels='l', grid=False, labels=False)
很明显，许多观测到的特征都被这个简单模型再现了，包括西部边界流、双旋环流、南极绕极流（ACC）等。我们放大西北太平洋部分，以获得更清晰的视图：
fig, axes = pplt.subplots(nrows=1, ncols=1, figsize=(7, 4.3), sharex=3, sharey=3, proj=pplt.Proj('cyl', lon_0=180))
skip = 2
fontsize = 12
ax = axes[0]
p=ax.contourf(h1/1e6*100, cmap='bwr', levels=np.linspace(-80,80,33))
ax.set_title('Stommel-Munk solution to wind stress curl', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip], lat.values[::skip,::skip],
u1.values[::skip,::skip], v1.values[::skip,::skip],
width=0.0012, headwidth=10., headlength=12., scale=80)
ax.set_ylim([-0, 65])
ax.set_xlim([-80, 60])
ax.colorbar(p, loc='b', label='', ticks=15, length=1)
axes.format(abc='(a)', coast=True, lonlines=20, latlines=10, lonlabels='b', latlabels='l', grid=False, labels=False)
我们还可以使用不同季节的风应力，以及在高纬度地区较大的空间变化底部摩擦 进行计算（以使得 ACC 的输送有所减少）。
# loading climatological January
curl_Jan = curl[0].load()
# loading climatological July
curl_Jul = curl[6].load()
iParams = {
'BCs' : ['fixed', 'periodic'],
'mxLoop' : 4000,
'tolerance': 1e-14,
'optArg' : 1.0,
'undef' : np.nan,
}
# use spatial-varying bottom friction
mParams = {'R':2e-4 * (2 - 1*np.cos(np.deg2rad(lat))), 'rho0':1027, 'A4':3e3, 'beta':beta, 'D':depth}
# inversion
h1 = invert_StommelMunk(curl_Jan, dims=['lat','lon'], mParams=mParams, iParams=iParams)
h2 = invert_StommelMunk(curl_Jul, dims=['lat','lon'], mParams=mParams, iParams=iParams)
u1, v1 = cal_flow(h1, dims=['lat','lon'], BCs=['extend','periodic'])
u2, v2 = cal_flow(h2, dims=['lat','lon'], BCs=['extend','periodic'])
lat, lon = xr.broadcast(u1.lat, u1.lon)
fig, axes = pplt.subplots(nrows=2, ncols=1, figsize=(7, 6), sharex=3, sharey=3, proj=pplt.Proj('cyl', lon_0=180))
skip = 3
fontsize = 11
ax = axes[0,0]
p=ax.contourf(h1/1e6*depth, cmap='bwr', levels=np.linspace(-60,60,13), extend='both')
ax.set_title('Stommel-Munk response (January)', fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip+2], lat.values[::skip,::skip+2],
u1.values[::skip,::skip+2], v1.values[::skip,::skip+2],
width=0.001, headwidth=10., headlength=12., scale=30)
ax.set_ylim([-76, 75])
ax.set_xlim([-180, 180])
ax.colorbar(p, loc='r', label='', ticks=10)
ax = axes[1,0]
p=ax.contourf(h2/1e6*depth, cmap='bwr', levels=np.linspace(-80,80,17), extend='both')
ax.set_title('Stommel-Munk response (July)',
fontsize=fontsize)
ax.quiver(lon.values[::skip,::skip+2], lat.values[::skip,::skip+2],
u2.values[::skip,::skip+2], v2.values[::skip,::skip+2],
width=0.001, headwidth=10., headlength=12., scale=30)
ax.set_ylim([-76, 75])
ax.set_xlim([-180, 180])
ax.colorbar(p, loc='r', label='', ticks=10)
axes.format(abc='(a)', coast=True, lonlines=40, latlines=15, lonlabels='b', latlabels='l', grid=False, labels=False)
显著的特征是北半球在7月的环流减弱，以及索马里急流在印度洋上空的增强。
API接口
References
Munk, W., 1950: On the wind-driven ocean circulation. Journal of Meteorology, 7, 79-93.
Stommel, H., 1948: The westward intensification of wind‐driven ocean currents. Eos, Transactions American Geophysical Union, 29, 202-206.
Matsuno, T., 1966: Quasi-Geostrophic Motions in the Equatorial Area. Journal of the Meteorological Society of Japan, 44, 25-43.
Gill, A. E., 1980: Some Simple Solutions for Heat-Induced Tropical Circulation. Quarterly Journal of the Royal Meteorological Society, 106, 447-462.
https://github.com/miniufo/xinvert
往期回顾
Python | SST | 高通滤波 | EOF | 论文复现 Python | 解决cmaps库报错的问题 Python | 北大西洋涛动 | NAO指数 | EOF Python | HadSLP2 ASCII转NetCDF Python | 半球间经度转换 Python | 大气科学 | 偏相关 37个CMIP6模型对流参数化方案
