已有 11165 次阅读 2014-3-28 21:53 |个人分类:海洋模式|系统分类:科研笔记| POM模式, 入门指南

开学一个多月以来，自学了POM的使用。因为在学的时候遇到好多困难走了许多弯路，想求助互联网时候发现资源很少，现在就想把自己的学习心得写一写，希望能够帮助像我一样的初学者解决一些基本的问题。由于时间有限，篇幅有限，我就写一部份贴出一部分。本人才疏学浅，不足之处、疏漏之处望各位前辈批评指正，有什么问题可以直接给我留言。

1 下载POM的源代码以及grid.f程序

      在POM官网上下载Pom2k.f以及Pom2k.c，第一次模拟海洋的运动确实很激动，但是激动过后发现离自己想要的结果还需要做一些工作。相信有数值计算基础你的已经想到了，要求解描述海洋运动的原始方程组，原始方程组的数值差分求解是最关键的，这部分正是POM提供给我们的，能够简化我们的许多工作。但是仅仅有差分求解方程组的程序是不够的，对于微分方程组，我们需要初始条件、边界条件来确定方程组的数值解，要正确使用POM模拟海洋运动，我们需要提供的正是这些边界条件。我们要做的工作主要包括提供驱动模式上下边界条件，侧边界条件，以及初始的海洋的3维温盐场。

另外，grid.f正是一个提供这些所需的数据的程序，但是我不推荐使用grid.f这个程序。第一，这个程序的代码很老，用新的编译器总是出现问题；第二，为了了解POM模型（我认为这是一个新手最应该做的），应该自己尝试多编程，多实践，只是使用别人写好的代码，对自己的学习进步没有好处；第三，grid.f中的插值方式，对各种条件的处理方法对于使用者来说是确定的，没法修改，而驱动POM模式通常是需要做很多数值实验的。综上所述，我建议新手应该自己编程给POM模式建立input文件，当然grid.f这个程序可以用来给自己作为参考。

其次，在这里我推荐使用Matlab制作POM输入文件，我认为Matlab提供的插值函数很成熟，自己没必要再为了运行POM而去再编写插值函数，而且Matlab提供了图形输出功能，能让你随时检查自己制作的input数据是否有误。在接下来的部分，我都是基于Matlab来制作驱动POM的input文件。

2 建立模式网格

网格的概念这里有两个，一个是POM模式的网格，也就是模拟的区域内微分方程离散出的网格点，这些网格点上的经纬度、温度、盐度等数据是我们需要输入到POM模式里的。这个网格叫做计算网格，是我们人为设定的。第二个网格，就是我们使用的原始数据的网格，比如我们从网上下载的WOA09数据，它的网格是全球1度乘1度的，WOA09的数据都是分布在这些网格点上的。如果我们要模拟0.5度乘0.5度的精度，那么我们要做的就是将原始数据网格点上的数据通过特定的插值方法插值到我们的计算网格点上去。

所以第一步要建立计算网格，确定计算网格的经纬度。建立计算网格等间距的计算网格在Matlab中很容易，只需要三个命令即可，比如要模拟南海，即可采用：

Lon=100：0.5：125;

Lat =0：0.5：25;

[LON LAT]=meshgrid(Lon,Lat);

这样计算网格的经纬度便存储在了LON，LAT中了。接着，要把这个数据从Matlab中写出到文件中，便于POM读取。如何写出以及POM如何读我就不啰嗦了，值得注意的这是水位的计算点的经纬度，POM中要把LON,LAT数据读给变量east_e和north_e。

如果要建立更精确的网格，或者曲线网格，可以使用一些其他的工具，比如seagrid。

在建立了计算网格之后，下一步工作是将需要的数据都插值到计算网格上去。主要包括地形数据、风场数据、辐射数据、初始的温盐数据

3 地形数据

地形数据对应着求解方程组的下表面边界条件，通过实验以及阅读前人的文献，发现地形会对模拟结果产生很重要的结果，所以这是第一步，也是很重要的一部。

首先，确定自己的模拟区域后，我们需要的自然是描述研究区域的地形数据了，通常用的是etopo5地形数据，这在网上可以自行下载。另外在高精度的模拟中，也可以使用etopo2和etopo1数据。下载地形数据后，读如到Matlab中（如何处理原始数据我就不啰嗦了，每种数据都有标准读取程序的）。接着，为了进行插值，我们需要在Matlab中按照第一步的方法，再把原始地形数据的经纬度网格也建立起来。比如地形数据是1/12度的，便可以使用如下命令

Lon_e=100：1/12：125;

Lat_e =0：1/12：25;

[LON_e LAT_e]=meshgrid(Lon_e,Lat_e);

注意：上述Lon_e, Lat_e的范围是假设的，具体要根据你的原始数据的范围确定。现在有了计算网格的经纬度LON、LAT，原始数据网格的经纬度LON_e、LAT_e，原始数据的值data_n，便可以使用Matlab命令interp2（…）来得到计算网格点上数据值了。

有时候处理地形复杂的区域为了保持计算结果的稳定需要对地形进行平滑处理，对地形的处理就在此处进行。平滑有很多方法，也有很多现成的程序，可以自己写，也可将POM中的slpmax函数移植到matlab中使用。

接着要对地形进行修改，由于使用etopo数据陆地为正值，海洋为负值，这里需要将正值（陆地）设置为0，将海洋里的水深设置为正（POM中水深为正值）。接着为了保持计算稳定，要设置地形数据的最大最小值，一般最大值根据实际情况确定，最小值确定为1或者5或者10米。

对地形数据做的其他修改也在此处进行，比如人为地修改边界，特定点的地形等等。  

将上述处理后得到的计算网格点上的地形数据存储，读如到POM中的h(im.jm)变量中即可。地形数据输入设定完毕。