此次作业首先用fortran实现了最简单的欧拉显式差分格式，调用csv_file模块输出为csv文件，并用python进行画图；整个过程用make全自动进行。之后的每次模拟都只在这个最简单的模拟的基础上做一处改动。我们首先检验了不同的时间分辨率造成的影响。之后我们改变$\kappa,f_0,\beta$的值，对这三个参数进行了比较初步的敏感性分析。最后选用其他的时间积分格式，包括Leap-frog格式，松野格式，对时间积分格式的性能也做了一些分析。
向前差显格式的核心步骤是： $$ \frac{\zeta^{n+1}-\zeta^n}{\tau} =-u^n\frac{d_h \zeta^n}{d_hx}-v^n\frac{d_h \zeta^n}{d_hy}+\kappa\nabla_h^2\zeta^n \tag{1} $$ 其中$\tau$为时间分辨率，$\zeta^n,u^n,v^n$为$t = t_n$时刻的涡度、纬向速度和经向速度，$\kappa$为涡度扩散系数，$\frac{d_h}{d_hx},\frac{d_h}{d_hx},\nabla_h$为离散的空间梯度（在边界上取单边差分，在内部取中央差分）。这里时间步长选取一小时，也就是3600秒。后面的实验都只在此实验基础上做一处改动。
Fortran执行的动力学部分的流程图如下：
st=>start: 开始
iox=>inputoutput: 读入位势高度的初始值
op=>operation: 计算PSI,U,V.ZETA
op2=>operation: 计算HT,U,V
dump=>inputoutput: 用csv格式输出HT,PSI,U,V,ZETA
cond=>condition: 是否积分满24小时
sub1=>subroutine: 用公式(1)更新ZETA
sub2=>subroutine: 求解泊松方程得到PSI
#io=>inputoutput: 输入输出框
io=>inputoutput: 输出末状态
e=>end: 结束框
st->iox->op->cond
cond(yes)->io(bottom)->e
cond(no)->dump->sub1(bottom)->sub2(right)->op2(top)->cond
因为Fortran没有画图模块，在实际操作中我使用了make，达到了类似python, fortran, shell 三种语言混编的效果。Makefile的内容如下：
flow.gif ht.gif: mass_plot.py u%.csv v%.csv ht%.csv
/home/mace/c/Anaconda3/python.exe mass_plot.py
rm *.csv
rm flow*.png
rm dvor*.png
rm vor*.png
rm ht*.png
u%.csv v%.csv ht%.csv: mine.f90
gfortran -o mine mine.f90 csv_file.o
./mine它所实现的整个运行和画图流程如下：
st=>start: 开始
fortran=>subroutine: 用Fortran进行模拟，并输出csv
os=>inputoutput: 得到目录下的文件列表
rd=>inputoutput: 读入下一个文件
cond=>condition: 读完了？
u=>condition: 是U吗？
ht=>condition: 是HT吗？
zeta=>condition: 是ZETA吗？
getv=>inputoutput: 读入对应时刻的V
quiver=>operation: 画箭头图
contour=>operation: 画等值线图
gif=>operation: 将图片排序并制作gif
delete=>end: 删除过程文件，结束
e=>end: 结束
st->fortran->os->rd(right)->u(no)->zeta(no)->ht(no)->rd
cond(yes)->gif(right)->delete(right)-e
cond(no)->rd
u(yes)->getv->quiver->cond
zeta(yes)->contour->cond
ht(yes)->contour->cond
这两个常数的表达式分别为 $$ f_o = 2\Omega\sin\varphi\ \beta = \frac{d f_0}{dy} = \frac{d\varphi}{dy}\frac{d f_0}{d\varphi}=\frac{2\Omega\cos\varphi}R $$ 期中$\varphi,\Omega,R$分别代表纬度，地球自转频率和地球半径。模拟的中心纬度$45^\circ N$作为$\varphi$，取$\Omega = 7.2921\times10^{-5}s^{-1},R = 6.371\times10^6m$。求得 $$ f_0 = 1.031\times10^{-4}\ \beta = 1.619\times10^{-11} $$
涡度的扩散系数在实验上是一个非常难以计算的数字，为了选取合适的$\kappa$值，我参考了1 。文中提到的多个模型$\kappa$从$7\times10^4$至$2\times10^6$不等。作为基础款的模拟，我们选用了$\kappa = 0$。这一方面是为了减少主观因素，另一方面也是为了使得对照组可以更好的凸显出平滑和其他积分格式的作用。$\kappa$的其他取值将在敏感性分析中给出。
解法的核心思想是用快速傅里叶变换后的乘法代替求导/积分运算。这一部分要解决的问题是：
$$
\Delta \psi =\zeta\
$$
最后我们还需要考虑函数的边界条件。由于左边的拉普拉斯算子是二次的，对于任意的a,b有： $$ \Delta(ay+b) = 0 $$ 由此可以满足y方向上的固定边界条件。在x方向上则直接将西边界上的值设定成东边界上的值即可。
程序中调用的是Numerical Recipe的多维FFT算法(正向逆向各一次)，程序读入n维数组变形而成的一维数组，相邻的两位代表复数的实部和虚部。普通的FFT算法可以通过下面的式子简单的拓展到高维（以二维为例）： $$ \begin{align} H(n_1,n_2) &= FFT_1(FFT_2[h(k,1k_2)])\ &= FFT_2(FFT_1[h(k,1k_2)]) \end{align} $$ 其中$FFT_1,FFT_2$分别代表对第1，2个坐标轴做一维的快速傅里叶变换。接下来我将介绍一下FFT的实现方法，以及程序中其他的细节。FFT算法的核心是Danielson-Lanczos引理： $$ \begin{align} F_k &= \sum{j=0}^{N-1}e^{2\pi ijk/N}f_j\ &= \sum{j=0}^{N/2-1}e^{2\pi i(2j)k/N}f_j+\sum_{j=0}^{N/2-1}e^{2\pi i(2j+1)k/N}f_j \ &= \sum_{j=0}^{N/2-1}e^{2\pi ijk/(N/2)}f_j+W^k\sum_{j=0}^{N/2-1}e^{2\pi ijk/(N/2)}f_j \ &=F_k^{偶}+W^kF_k^{奇} \end{align} $$ 其中$W = e^{2\pi i/N}$。每一次操作都将要完成的任务分成了两份，只考虑奇数或只考虑偶数。这个过程可以重复的进行1分2，2分4将矩阵逐步细分。程序要求N是2的整次幂，因此最终分到的每个部分就只有一个数，此时FFT的计算就没有难度了。 $$ F_k^{奇偶偶\cdots奇偶} = f_j $$ 如果将上标倒序，将每个“偶”变成0，将每个“奇”变成1我们就得到了j的二进制表示。这也就是为什么在程序的一开始要将数组颠倒一下。
除此之外另一个有趣的地方是三角函数值的计算，这个算法在写的时候非常细腻，为了避免重复计算三角函数值（比较昂贵），使用了递推的方法： $$ \cos(\theta+\delta) = \cos\theta-[a\cos\theta+b\sin\theta]\ \sin(\theta+\delta) = \sin\theta-[a\sin\theta-b\cos\theta] $$ 其中$a = 2\sin^2\left(\frac\delta2\right),b = \sin\delta$，$\delta$为每次增加的角度。虽然这个公式很简单，可以用和角公式快速得到，但确实很大程度上加快了计算。
24小时预报结果如下：
风场：
位势高度：
涡度场：
虽然流场和高度场还比较光滑，但从涡度场已经可以看到一些高频的信号被明显的放大了。这说明向前差格式基本可以胜任24小时预报，但是还是有进步的空间。
起初在实现向前差格式的时候我犯了这样的一个错误，正常的流程本来应该是
psi=>start: PSI
uv=>operation: get U,V
nuv=>operation: get new U,V from PSI
zeta=>operation: get ZETA
ahead=>operation: find ZETA of the next moment
psolve=>operation: get PSI by solving Poisson equation
psi(right)->uv(right)->zeta(right)->ahead(right)->psolve(right)->nuv(top)->ahead
我在设计的时候，为了使每个时刻的ZETA都对应当时的U,V，每次更新PSI的时候都重新更新一次ZETA，具体流程如下图。
psi=>start: PSI
uv=>operation: get U,V
zeta=>operation: get ZETA
ahead=>operation: find ZETA of the next moment
psolve=>operation: get PSI by solving Poisson equation
psi(right)->uv(right)->zeta(right)->ahead(right)->psolve(top)->uv
而这导致了极大的误差，甚至淹没了整个物理过程。这背后的原因是由ZETA求解离散的泊松方程得到PSI和用PSI的中央差得到ZETA并不是彼此的逆运算。
在查阅资料的过程中我看到了另一种求解泊松方程的方法2，这个方法是在随机的三角形网格上的有限元方法。这个方法给了我一定的启发。如果把PSI变形成一个NX*NY维的列向量$\vec{\psi}$，那么求解ZETA的过程可以写成：
$$
\vec\zeta = A\vec\psi\tag{2}
$$
其中$A$代表内部中央差分，边界上单侧差分求得ZETA的算符，它可以表示成：
$$
A = D_x\otimes I_y+I_x\otimes D_y\tag{3}
$$
值得注意的是虽然求解$A$的逆需要的时间是$O(NX^3NY^3)$，但矩阵$A$并不随时间变化，因此求一次逆就可以一劳永逸，每次求解PSI只需要用$O(NX\cdot NY)$的时间计算一个内积即可。相较之下二维FFT算法每一次的时间复杂度是$O(NX\log(NX)NY\log(NY))$，因此在积分的步数比较多的时候，直接求解线性方程所用的时间是比FFT算法更快的。
比较遗憾的是，我安装的gfortran的lapack包在链接的时候遇到了一些问题，因此没能自己实现这一想法。
同为向前差显格式，选时间分辨率分别为3600s,1800s,900s，我们做了如下三次模拟：
容易看出时间分辨率的增加抑制了高频信号的增长，这和预期非常一致。
五点平滑的过程可以写为 $$ f_{i,j} = f_{i,j}-\frac s4(f_{i+1,j}+f_{i-1,j}+f_{i,j+1}+f_{i,j-1}-4f_{i,j}) $$
我们在实际操作时在每次更新过ZETA后，都对ZETA做一次5点平滑，这样理论上可以滤掉一部分高频波动。其中$s$是一个有一定主观色彩的参数，这个值直接决定了模拟过程中各种波长的信号随时间的衰减。参考了1 中s的选取，我采用了$s = 0.03,s=0.1$这两种模拟。
在$s=0.03时$，滤波有一定效果但并不完全令人满意：
而在取$s=0.1$时，从位势高度上来看，有一部分有用的信息也被滤掉了：
这说明用好平滑方法需要一些技术和经验，并不是想象中那么容易。
1: 张耀科. 正压涡度方程中考虑粘性项与平滑过程的一些意见[J]. 气象学报, 1965 (1): 3.
在所有的常数之中$f_0$是唯一一个没有在核心的时间积分步骤中出现的。它只出现在位势高度和流函数PSI的相互转换之中，理想情况下对PSI,ZETA,U,V的影响都是线性的，而对位势高度的预报没有影响。例如如果我们选取$65^\circ N$对应的地转频率，得到的位势高度预报几乎和没有改变$f_0$之前完全相同：
然而由于$f_0$较小是U,V较大，涡度方程的平流部分可能会发生不稳定，在选取$25^\circ N$对应的$f_0$带入模式时，确实有明显的不稳定出现：
简而言之，模型在保证计算稳定性的情况下对$f_0$不敏感。
从物理原理上来分析$\beta$可以影响Rossby波的速度，所以如果模型准确的话，预测结果就应该和$\beta$有关。我们对$\beta$为原本值的一半、为原本值的两倍分别进行了数值实验，结果依次为：
以上两组结果和对照组比较可以发现$\beta$越大高压向西移动越慢。这是因为$\beta$越大，高压东侧的北风因为$\beta$效应产生的正涡度越强，阻碍了高压向东移动。这个现象也可以用Rossby波的相速度理论来解释。
整体来看，模型比较好的捕捉到了$\beta$效应。模拟表明，$\beta$的大小对系统有一定的影响。但是在24小时内，即使$\beta$成倍的变化，对预报准确性的影响也不是很大。
这一系列实验中我们分别选取$\kappa = 7\times10^4,2\times10^5,2\times10^6$。最后一个量明显是偏大的，引入这一组实验主要是因为前两组的滤波结果并不十分理想。
虽然最后一组的滤波效果比较理想，但类似于平滑算法中$s = 0.1$那一组，这一组也丢失了一部分有用的信息。这说明$\kappa$的选取也是很需要知识和经验的。
总而言之，因为前两组得到的结果和对照组都十分接近，我们认为模型对$\kappa$并不十分敏感。
松野格式可以表示为 $$ \frac{\zeta^{n+1}*-\zeta^n}{\tau} =-u^n\frac{d_h \zeta^n}{d_hx}-v^n\frac{d_h \zeta^n}{d_hy}+\kappa\nabla_h^2\zeta^n \ \frac{\zeta^{n+1}-\zeta^n}{\tau} =-u^n\frac{d_h \zeta^n_}{d_hx}-v_^n\frac{d_h \zeta^n_}{d_hy}+\kappa\nabla_h^2\zeta_*^n \tag{4} $$ 不难看出松野格式直接使得计算量翻倍了，但也让格式更接近隐式格式了。流程图只需要更改时间积分这一处。最终预报的结果如下：
在24小时内没有高频信号的急剧增加，可以说模拟的挺不错。
Leap-frog格式就是把空间差分也写成中心差分的形式，具体表示为 $$ \frac{\zeta^{n+1}-\zeta^{n-1}}{2\tau} =-u^n\frac{d_h \zeta^n}{d_hx}-v^n\frac{d_h \zeta^n}{d_hy}+\kappa\nabla_h^2\zeta^n\tag{5} $$ 虽然计算量没有增加，但是保存涡度的内存需求大约增加了一倍。流程图也需要做一些调整。
st=>start: 开始
iox=>inputoutput: 读入位势高度的初始值
op=>operation: 计算PSI,U,V.ZETA
op2=>operation: 计算HT,U,V
odzt=>operation: 保存OLDZETA
sub=>subroutine: 用公式(1)更新ZETA
zc=>operation: 在TEMP里暂存ZETA
qc=>operation: 用TEMP给OLDZETA赋值
dump=>inputoutput: 用csv格式输出HT,PSI,U,V,ZETA
cond=>condition: 是否积分满24小时
sub1=>subroutine: 用公式(5)更新ZETA
sub2=>subroutine: 求解泊松方程得到PSI
#io=>inputoutput: 输入输出框
io=>inputoutput: 输出末状态
e=>end: 结束框
st->iox->op->odzt->sub->cond
cond(yes)->io(bottom)->e
cond(no)->dump->zc->sub1(bottom)->qc->sub2(right)->op2(top)->cond
最终预报的涡度场如下图所示
和松野格式不分伯仲。
见动画，使用了最为可靠的几个模型。依次为向前差，向前差+平滑（s=0.03)，向前差（dt=900)，向前差（k=7e4)，松野格式，Leap-frog格式。
此次作业实现了正压涡度方程的向前差求解，预报出了高压向东的移动。在此基础上证明了模式对$\kappa,\beta,f_0$都不是很敏感。为了使计算更加稳定，我们选用了增加时间分辨率、平滑、加大粘性项、改变时间积分格式这几种方法。其中平滑和加大粘性项的尝试并不十分成功