THREP：面向地球系统模式的插值算法研究平台

吴 竑，王小鸽，杨广文，朱昌磊，宋顺强

清华大学 计算机科学与技术系，北京 100084

THREP：面向地球系统模式的插值算法研究平台

吴 竑，王小鸽，杨广文，朱昌磊，宋顺强

清华大学 计算机科学与技术系，北京 100084

1 引言

1.1 背景

插值算法[1]是一类应用十分广泛的算法，它们在图形图像处理、数值计算等许多科学研究领域都是不可或缺的重要数学手段。插值算法在耦合地球系统模式中也有着极其重要的应用。

耦合地球系统模式[2]是一套对地球系统进行定量研究的软件工具，地球科学相关领域的科学家根据专业知识对地球各个子系统进行建模，编写各自领域的模拟程序，再由这些模拟程序相互耦合，成为一个完整的地球系统模式。为了使不同研究人员的代码能够耦合起来形成一个整体，往往需要他们共同遵循某个耦合软件框架[3]。目前较为流行的地球系统模式耦合框架有欧盟PRISM所采用的OASIS[4]和美国CCSM3所采用的COUPLER[5]系列等。然而，由于地学领域的科学家们对地球各子系统进行建模时所使用的数值方法不同，导致他们在编程时使用不同的地球网格划分方式（如经纬网格、高斯网格[6]），从而对地球系统模式中不同分量模式间的数据交换提出了挑战。因此，在不同网格间进行数据交换时需要一种数据映射方式——插值。此时，可以将每个网格单元看做一个插值点，并用网格单元中心点坐标表示插值点坐标。

1.2 相关工作

SCRIP[7]是一个面向地球系统模式的插值算法库，可以离线地生成地球系统模式研究人员所需地球网格间的插值权重系数文件，它被使用于CPL6等耦合器中。但是SCRIP中的插值算法不很完善，其守恒插值算法只支持一阶，双三次插值算法也没有实现求导数的方案。此外，由于SCRIP中的数据结构过于简单，它能够支持的网格类型非常有限，这些因素都制约了其在地球系统模式中的作用。

ΤHU-Remap（https：//github.com/xunzhang/Τ-Remap）是针对SCRIP的缺陷，开发的一套全新插值算法库。它一方面重构了SCRIP中一些算法，如反距离权重算法、双线性算法，大大提高了它们的计算性能；另一方面改进了SCRIP中的一些不足，如实现了双三次差值中的求导、提供守恒插值中阈值的自动调整等。但是，ΤHU-Remap对非规则网格的支持依然不够，这不利于分量模式与耦合器的长期发展。

ESMF（Earth System Modeling Framework）[8]项目由NASA地球科学办公室于2000年9月提出，是一套面向地球系统模式的并行软件开发框架。它不含任何分量模式代码，而是提供标准框架接口，用户可以基于此接口实现天气、气候、数据同化等相关领域的应用模型并将它们耦合起来成为一个完整的地球系统模式。ESMF模块化设计使其具有极高的可重用性、互操作性和可移植性，这将有利于模式专家专注于模式本身的发展并减少其耦合开销。ESMF主要提供Superstructure和Infrastructure两部分功能供开发者使用：Superstructure中封装了地球系统模式中的耦合驱动和接口，也模块化了分量模式中物理过程、动力框架及一些子过程；Infrastructure则提供更加底层的抽象，它又可细分为数据处理和公用模块两部分，数据处理部分主要封装所有与地球系统模式中数据处理相关的数据结构及功能，如Array、Grid、Field等，公用模块部分则提供了如数值计算、时间管理、进程管理、I/O、容错日志处理的基本功能。

ESMF-Regridding[9]是ESMF提供的一个基于ESMF框架开发的插值应用，它能并行生成离线插值权重系数文件，特别是它能支持非规则网格间的插值。但是，ESMF-Regridding仅以一个小工具的形式出现在ESMF中，没有抽象插值在应用层面的功能供开发人员使用，因此不够灵活和通用。插值算法开发者往往需要阅读大量ESMF源码，作出相应修改并实现各自新的插值算法，难度极大，违背了ESMF的设计目标。此外，ESMF-Regridding没有提供验证插值算法正确性与性能的功能，不利于深入研究与对比各插值算法的特性。

虽然目前已有SCRIP、ΤHU-Remap、ESMF-Regridding等面向地球系统模式的插值软件，但随着地球系统模式研究的深入，对插值算法提出了新要求。为插值算法的研究、测试和评价建立一个统一的平台，具有重要的实践意义。基于ESMF框架，开发了一个通用、灵活、并行的插值研究平台ΤHREP（Τsinghua Regridding Platform），它能同时满足插值算法研究人员、耦合器开发者、模式专家三类用户需求。ΤHREP集成了ΤHU-Remap与ESMF-Regridding中的插值算法，封装了易用的插值接口，并提供了验证、对比等功能。

2 平台设计

ΤHREP综合考虑了算法研究人员、耦合器开发者及模式专家三类人对插值的需求，是一个通用、灵活、并行的面向地球系统模式的插值平台。它抽象了插值过程，封装了插值流程，便于已有插值算法库的接入，利于长期发展和标准化。同时，为满足不同用户，ΤHREP在保证易用性的同时，从底到上提供了多层接口。ΤHREP还提供了插值检验功能，便于深入研究各算法的精度、性能等相关问题。

2.1 平台架构

ΤHREP插值平台的总体结构图见图1，主要包括输入文件、配置模块、权重生成模块、插值模块和检验模块几个部分。输入文件含网格数据文件和物理量数据文件，配置模块负责配置一次插值过程需要的信息，权重生成模块包含了初始化和搜索两个子模块，ΤHREP通过初始化模块建立相应数据结构的对象，在此基础上搜索找出各算法所需格点信息，生成插值权重，插值模块通过矩阵向量乘计算出实际物理量插值后的结果从而完成整个插值过程，检验模块将对插值的精度和性能作进一步分析并画出图形。

图1 ΤHREP平台架构

2.2 通用性

目前地球系统模式中常用的插值算法有反距离权重插值[7]、双线性插值、双三次插值[10]、Patch插值[11-12]、三次样条插值[13]、守恒插值[7，14]等，ΤHREP里将所有插值算法抽象成了如下两个步骤：

（1）计算插值权重系数矩阵。这个阶段又可细分为两个子过程：首先根据源网格和目标网格间的拓扑关系，对每个目标网格单元（cell）找到其对应位置在源网格的“若干”网格单元（这里的“若干”由不同插值算法决定，如双线性插值要找出包住目标网格单元的4个源网格单元，即目标网格在找出的4个源网格单元连线构成的四边形内部）。然后根据插值算法计算出每个找到的源网格单元相对于该目标网格单元的权重。因为每个目标网格单元都要计算其对应的“若干”个权重值，因此这一步完成后，会得到一个权重矩阵W，W的每行对应每个目标网格单元id，每列对应每个源网格单元id。这里W是一个稀疏矩阵，如双线性插值的权重矩阵每行只有4个非零元素。

ΤHREP中使用NetCDF[15]标准的文件格式作为整个平台数据流中输入、输出文件的格式。NetCDF文件是一种格式化的二进制文件，它既有二进制文件的数据连续性，又能使用户像文本文件那样方便地读取相应属性的数据。采用NetCDF文件的另一个重要原因是它已经是研究地球系统模式人员的约定标准，分量模式中使用的网格文件均采用这种格式，因此用NetCDF作为插值结果的输出文件能更便于用户直接使用。

考虑通用性的另一方面是功能可扩展，而数据结构直接决定了能否支持众多功能。因此，通用的数据结构也是系统通用性的关键。整个插值过程中主要的数据包括网格信息和物理量信息，插值所用的数据结构主要就是为了便于插值计算而在这两部分数据信息基础上建立的。在ΤHREP中尽量保证了数据的原始逻辑，并刻画了网格单元间的相邻逻辑关系，这样可以保证功能可扩展，便于平台的进一步研发，也为以后我们或者其他开发者封装更高级的数据结构提供了可能。

2.3 灵活性

在众多地球系统模式的研发和使用人员中主要有三类人员涉及插值。第一类是模式专家，他们负责开发分量模式，专注于各分量模式的功能、精度与性能问题。因各分量模式在自然界中并不是孤立存在，因此一个模式依赖另一个模式中插值得来的数据做计算；第二类是耦合器开发人员，耦合器负责分量模式之间的数据管理，而数据管理主要包括通信、插值、集合等子模块，因此耦合器开发人员也关注插值；第三类是插值算法研究人员，他们专注于插值算法本身的研究，如SCRIP就是他们为地球系统模式中的离线插值开发的一个算法库。为了能同时满足三类研究人员的需求，ΤHREP在对插值流程进行封装的基础上提供了并行生成权重矩阵、计算插值、验证插值正确性、插值性能测试等功能。

生成权重矩阵是离线生成一组网格对应的插值系数文件，对于一组网格（一个源网格和一个目标网格）只需计算一次。在线耦合过程将调用该系数文件，再完成插值的第二步矩阵向量乘操作，而这个系数文件可能被多次使用。同时，ΤHREP支持直接计算插值，即不产生权重系数文件，直接进行插值计算，得到插值结果，ΤHREP中计算插值权重的过程是并行的。平台还提供了插值性能测试功能，为研究高性能插值提供了可能。最后，验证模块主要基于插值结果文件，用两套验证策略对计算出的插值结果进行精度分析，并计算一些统计量信息、同时画出相应的误差分布图。

随着分量模式自身和整个地球系统模式的发展，在线插值越来越被重视。之所以要把离线的权重计算操作放到在线去做，是因为各分量模式的网格在模式迭代计算过程中可能不断发生变化，这样就要频繁地计算插值权重系数。这时对耦合器来说，插值的两步之间不再需要产生权重系数文件，而是直接算出插值结果，ΤHREP目前只支持串行这么做。当然，在线插值对插值性能提出了很高的要求，如何并行支持在线插值有待进一步研究。

ΤHREP还具有框架的特性，用户可以方便地在其接口上实现新插值算法。同时，ΤHREP中数据结构的通用性使开发插值算法库的人能方便地作转换，将已有插值库快速接入到ΤHREP平台中。ΤHREP也提供了插值算法层的接口，用户能方便快速地在平台中开发新的插值算法。ΤHREP中将插值分为结构网格间的插值与非结构网格间的插值，由于结构网格的索引能更清晰地反映网格的拓扑关系，为使用简单，将结构网格的接口与非结构网格的接口分离开来，而底层实现则采用一致的数据结构。

2.4 并行性

并行生成插值权重系数是ΤHREP的另一特点，这为研究插值算法性能的用户提供了有利工具，也为以后支持在线插值提供了可能。考虑到代码复用，不想让插值算法研究人员为了实现并行做重构，因此希望并行尽量地普适。并行生成插值权重本质上是对源网格与目标网格进行剖分，这对于分布式环境中的每个独立计算单元P来看，计算方式和串行情况下几乎一致，而区别在于此时的网格不再是整个地球展开的平面，而是剖分后的一部分区域。如果把生成插值权重系数的过程抽象成一个算子Gen，那么唯一造成需要Gen重构的可能性是Gen与网格的拓扑形状有关。为此，将网格做了规则的矩形剖分，从而保持原有的网格形状，由用户来负责实现通信逻辑。ΤHREP提供一维通信的接口，用户可以使用它简化通信细节的实现。

2.5 插值检验

作为研究平台，ΤHREP具有插值检验功能。为了验证精度，需要对插值结果的正确性作定量分析。现设νalcal为ΤHREP中插值得出的结果，νalacc为物理量数据的精确值，则目标插值格点插值结果的相对误差er可用公式er= |νalacc-νalcal|/νalacc算出。

问题的难点在于如何刻画公式中的νalacc，因为实际上并不知道目标插值网格上的真实物理量分布数据，为此ΤHREP中采用了两种方式来定义νalacc。一是用解析函数值，即给定解析函数表达式，要求某坐标处的νalacc只需直接取函数值便可。这时为了尽可能真实地模拟不同物理量在不同地形的分布特性，可以取球谐函数[16]作为验证的解析函数。然而，真实的物理量分布式不能简单用解析函数刻画的，用解析函数测试精度较高的插值算法在真实物理量分布中不一定仍然精确。因此ΤHREP还使用另一种策略来定义νalacc：先将真实物理量的月平均数据用守恒插值插到一个“足够密”的网格上（“足够密”的含义是对于任何插值网格上的每一个点，都存在一个该“足够密”网格上的格点，使得它们之前的距离“足够近”），忽略这次插值的误差。这样，就可以用此“足够密”网格格点上的值代表目标插值网格上的物理量的值了。ΤHREP中定义的“足够密”网格分辨率为1 800×900，因此上述“足够近”大约在3.5 km范围内。这两种确定νalacc的方式分别对应着两套检验策略，用户在进行完插值后，可以选择任意一种进行检验，平台将通过NCL[17]脚本画出误差分布图，并同时算出平均相对误差（ΑVG_ERR）、最大相对误差（MΑX_ERR）、均方根误差（RMSE）等统计信息。

2.6 数据流

ΤHREP的数据流图如图2所示。图2（a）是顶层图，首先ΤHREP通过读入配置参数和网格信息产生插值算子对象和网格对象，同时读入源网格上物理量信息，平台利用这些对象与信息进行插值，从而求出插值结果。验证模块读入插值结果，画出物理量分布图与插值误差分布图。

图2（a）ΤHREP数据流顶层图

图2（b）ΤHREP数据流0层图

在整个插值计算模块内部，数据流向如0层图2（b）所示，数据将依次流过两个子模块：权重生成模块与计算插值模块。权重生成模块负责完成插值权重的计算，若配置要求离线插值，则此步将产生一个权重系数文件供耦合开发者户和模式专家使用，且此过程是并行的，若配置要求的是在线插值，ΤHREP则将在完成权重计算后直接计算出插值，得到物理量在目标网格单元上的值。需要指出的是，计算插值这个子过程本质上是一个稀疏矩阵乘向量的过程，虽然它本身易并行，但这时的并行剖分并不同于最初读入网格的剖分，而是需要进行一次全局通信来产生新剖分。因此，若想对整个插值计算过程做并行，则务必需要进行一次全局通信，随着剖分数变大，并行的效率也将大大下降。考虑到效率，ΤHREP中并没有提供并行地支持在线插值而只是并行地生成插值权重系数，要解决在线插值的并行效率，必须使得源网格的剖分充分灵活，从而减少甚至消除内部的全局通信。另外，对于插值结果，如图2（b）所示，验证模块可以通过解析函数和真实物理量两种定义νalacc的策略分别进行验证。

3 平台实现

地球系统模式中的插值不同于图形学中插值，由于它涉及的坐标点都在地球表面，因此可以认为是一种近似的平面插值。之所以说近似是因为地球表面并非平面，且在极点和周期边界处的处理也不同于普遍意义下的二维插值。本文基于上述设计，用C++实现了ΤHREP，并封装了C++的插值接口供用户使用。ΤHREP使用SCRIP标准网格文件格式，并以二维经纬坐标的形式存储网格单元。为同时满足上文中三类用户的需求，ΤHREP采用了ESMF中Infrastructure的Mesh作为其内部网格的数据结构，它足够灵活和通用，符合ΤHREP的设计目标。Mesh的数据结构既能用于封装结构网格插值接口，又能用于实现非结构网格间的插值。事实上，一个Mesh对象由nodes和elements构成，每个node是一个带坐标信息的点，每个element代表由若干nodes构成的一个单元区域。Element还给出了Mesh的形状信息与nodes间的相邻拓扑关系。对于结构网格，node就可以满足插值需求，而对非结构网格，则需加上element的信息。对于非规则网格插值的支持也意味着ΤHREP可用于带mask的网格，如海洋网格或海冰网格。考虑到球面插值的特殊性，ΤHREP在网格坐标对象中加入了额外的ghost点来处理周期边界问题，而对于极点问题ΤHREP提供三个选项供用户进行选择。

ΤHREP的实现层级图如图3所示，底层采用NetCDF格式作为输入文件格式，之上使用ESMF中Infrastructure的Mesh数据结构表示网格信息以保证通用，在Mesh层以上，ΤHREP平台封装了两套接口分别用于结构网格的插值（Grid接口）与非结构网格的插值（Mesh接口）。基于Grid，ΤHREP集成了ΤHU-Remap，其中包括反距离权重插值算法、三次样条插值算法和双三次插值算法。基于Mesh，ΤHREP集成了ESMF-Regridding，其中包括双线性插值算法、一阶守恒及二阶守恒插值算法。插值算法开发者也可以针对其算法特性，选取Mesh或Grid接口进行新算法的开发。对于那些有独特数据结构要求的算法或算法库，用户可以直接用Mesh数据结构作转换，而后集成到ΤHREP平台中。此外，ΤHREP中的验证模块底层实现则是基于NCL的。

图3 ΤHREP实现层级图

3.1 网格文件

ΤHREP选取SCRIP格式（一种结构化NetCDF格式）作为网格文件，其文件头如下：

Dimensions属性中的grid_size表示该网格的网格单元总数，grid_rank表示维度（grid_rank为2表示2D逻辑矩形网格，为1表示非结构网格）。Variables属性中包含grid_dim，它代表每一维的网格单元数，如对于Τ42网格就是[128×64]。grid_center_lat和grid_center_lon分别表示网格单元中心点的经纬度坐标，grid_corner_lat和grid_corner_lon分表表示网格单元的顶点坐标。需要注意的是，一个网格单元的顶点数由grid_corners表示，除了守恒插值其他的插值算法均用中心点的坐标代表一个网格单元。grid_imask是一个代表是否为海洋的整数，它决定了该网格单元是否需要参与插值计算，特别当grid_imask为0时，该网格单元不参与插值计算。网格的坐标单元可以是角度或弧度，读入网格信息后，ΤHREP根据文件中坐标单位标识统一将坐标信息转换成角度进行处理。

3.2 极区问题与周期边界

平面插值与球面插值有着本质的区别，因此将球面插值近似成平面插值需要对球面距离、极点以及周期边界作特殊处理。

首先，ΤHREP用两点的直线距离近似两点的球面大圆弧距离。一方面，大圆弧距离在计算量上远大于直线距离，另一方面，随着分量模式发展，网格分辨率将越来越高，球面两点的大圆弧距离与直线距离间的误差将越来越小，插值结果也将越精确。

其次，对于极点处的插值，ΤHREP提供了三个选项：一是不处理，即遇到目标网格单元中心点在源网格最低或最高纬度之外的情况时，运行过程中将直接抛出异常。二是按正常情况处理而不抛出异常。这时需要插值算法自身来处理极点问题，如现在平台中的双线性算法在处理目标网格极点插值时，由于找不到包含它的源网格单元构成的矩形，将会调用反距离权重的方式计算权重系数。考虑到双线性算法权重矩阵的大小，ΤHREP这时会找4个源网格单元并计算对应权重。三是取所有包含极点的源网格单元值的平均值作为插值的结果。

最后，对于周期边界，ΤHREP在生成网格数据对象的时候添加了ghost区域。所谓ghost区域就是经度为0附近坐标网格单元的集合，如经度为359.5的网格点对应的ghost坐标点经度为－0.5，而经度为0.5的网格点对应的ghost坐标点经度为360.5。事实上，这些ghost区域的点仅仅服务于插值所需的相邻关系，而对于ghost格点单元上物理量数据，ΤHREP则将还原到原始网格对应的网格单元上。

3.3 带mask的插值

由于分量模式中使用的网格一般是对全球区域的剖分，如海洋模式中剖分后有些网格单元中将会包含陆地。这时，mask用来标示该网格单元是否是海洋区域，因此对插值来说，无论是源网格还是目标网格，只要带mask且 mask为1表示参与计算，而0表示不参与插值计算。

有些插值算法的实现依赖于网格拓扑形状，如ΤHREP中的三次样条插值算法，目前只支持矩形经纬网格，因此一旦网格带了mask，整个网格将不再是经纬网格了。另外有些算法，ΤHREP中对源网格带mask的情况作了特殊处理。如双线性插值算法，其插值过程中需找出包含目标网格单元的4个源网格单元，若此时找出的源网格单元被mask掉了，就不能使用，ΤHREP会调用反距离权重进行插值，因为其实现不受限于网格拓扑形状。

3.4 并行

出于通用性考虑，ΤHREP中对源网格与目标网格作的剖分均是不重叠的，这也就需要各插值算法进一步实现通信细节。ΤHREP中的源网格剖分与目标网格剖分方式是相互独立的，目标网格支持一维、二维剖分而源网格只支持一维剖分。对于距离权重算法与双线性插值，由于计算每个目标网格单元对应的权重与源网格大小无关，而只要能“包含”目标网格单元即可。因此只需对目标网格进行剖分，而可以不将源网格进行剖分。这并不影响可扩展性，估算分辨率为1°的经纬网格（360×180）约只占0.5 MB的内存。对于三次样条等其他一些插值算法，需要同时对源网格与目标网格进行剖分。目前ΤHREP中的三次样条插值支持对源网格进行一维纬向剖分，对目标网格进行二维剖分。在插值开始前，算法会将相邻逻辑网格单元上的数据按层进行通信，从而满足上述“包含”的要求。

3.5 插值检验模块

ΤHREP不仅能做插值计算，还能检验插值结果的精度以及插值性能，这为研究插值提供了极大的便利。ΤHREP集成了NCL，使得能够清晰直观地画出插值误差分布图。

有的插值算法在插值时不仅需要源网格上的函数值信息，还需要导数值信息。如二阶守恒、双三次插值等。本质上，加入导数信息是为了更准确刻画真实物理量的变化趋势。为测试这类插值算法，对于解析函数，只需直接算出导数表达式取对应位置的导数值即可。而对于真实数据的方案，需要用数值的方法将物理量函数值近似出相应的导数值。ΤHREP中采用了有限差分[18]和最小二乘拟合[19]的方法来近似导数。

总地来看，对于一次插值，ΤHREP能相应地画出插值分布图、精确物理量的值分布图与相对误差分布图。同时，ΤHREP也会计算出最大相对误差（MΑX_ERR）、平均相对误差（ΑVG_ERR）与均方根误差（RMSE）等统计量信息。

3.6 插值接口与使用

ΤHREP还提供了方便的插值接口供用户使用。为算法开发人员提供了两套接口，分别适用结构网格与非结构网格间的插值。对于结构网格的插值，用户只需继承Common_remap基类，并在子类中实现两个虚函数init_remap和cal_remap。同时，为了使添加的插值算法能并行执行，ΤHREP封装了几种通信方式供用户调用。ΤHREP会在整个插值过程中先后调用init_remap和cal_remap这两个虚函数完成权重计算和插值计算。对于非结构网格的插值，平台提供了另一套接口帮助用户满足更复杂的插值需求。这套接口的通用性也为现有插值库提供了接入的可能。事实上，ΤHREP中两套算法接口底层是通过统一的数据结构抽象出来的。ΤHREP中的算法注册流程也是非常简单，便于新算法添加。最后，ΤHREP还为耦合器开发者和模式专家提供了方便的执行入口，只需简单的配置即可运行使用。

3.6.1 结构网格的插值接口

所谓的结构网格是指其索引是M×N的，即每行N个网格单元，共M行。因此结构网格的索引就能清晰地反映网格的拓扑关系，特别是网格单元的相邻关系。针对这类网格的插值需求，提供了一个基类Common_remap如下：

用户定义的类在初始化时设置set_size，Common_remap从而确定权重系数矩阵以及其他私有成员大小。定义的类必须继承Common_remap类并实现init_remap和cal_remap两个虚函数。

init_remap函数负责生成每个目标网格单元的权重系数，其中需要调用write_wgts函数写入每个目标网格单元对应的set_size个权重，还需要调用write_matchup写入每个目标网格单元对应的set_size个源网格点信息。如双线性插值，每个目标网格单元产生4个源网格单元对应的权重和4个源网格单元信息。Common_remap类还提供了两个简化搜索的函数接口，get_nearest_nbrs能够找出在threshold阈值范围内的若干个距离目标网格单元最近的源网格单元。如距离权重插值、9个点的双二次插值、16个点的双三次插值等算法都能方便使用。而get_bnd_quadrangle能够找出能够包住目标网格单元的由源网格单元构成的四边形，供双线性插值、4个点的双三次插值等一类算法使用。

cal_remap负责通过init_remap算出的权重矩阵与源网格单元上物理量值做插值。事实上，对于大部分算法，该过程是一个简单的矩阵向量乘法，所以用户定义的cal_remap函数只需调用基类中的cal_remap即可。但是对于像三次样条和二阶守恒这样的高精度插值算法，计算插值不能简单地抽象成一次矩阵向量乘，如三次样条实际是多次矩阵向量操作。这时，用户需进一步实现cal_remap，这也是将cal_remap定义为虚而非纯虚的原因。

需要特别注意的是，Common_remap类中带Grid的构造函数是为并行执行设计的。由于ΤHREP的并行剖分时没有重叠区域，要求用户实现通信细节，因此在平台调用并行插值前需要分布式环境中每个计算单元间进行通信，扩充局部源网格信息，确保在该计算单元上的目标网格能够完全被扩充后的源网格“包含”在内。目前ΤHREP只提供一维剖分的通信接口，用户调用此接口时要指定通信宽度来确保“包含”这个条件。对于二维剖分，目前用户需要自身实现全部通信细节。

3.6.2 非结构网格的插值接口

随着分量模式的发展，它们所用网格也越来越复杂，如三极网格、六面体网格等（http：//www.ncl.ucar.edu/Document/ Graphics/contour_grids.shtml#Τripole）。ΤHREP针对非结构网格间的插值为用户提供了接口。不过，ΤHREP目前还不算灵活，相比结构网格间的插值接口，此时用户需在指定文件中添加插值函数，用户可以使用srcmesh、dstmesh和iw三个数据结构的接口来完成插值。srcmesh与dstmesh中包含了非结构源、目标网格坐标以及相邻关系等信息，iw用存储插值权重，要求用户将计算完的插值权重写入其中。同样，若要将现有的插值库接入ΤHREP中，则需要用这三个数据结构作相应转换。

4 平台验证与分析

图4 ΤHREP离线插值界面图

ΤHREP的主菜单界面非常简洁易用，主要包括OFFLINE、ONLINE和VALIDAΤION三个按钮，分别供不同用户使用。OFFLINE指离线插值，能帮助用户离线并行地生成插值权重系数文件，图4给出了在ΤHREP主菜单界面下点击OFFLINE并进行一次离线插值的示意图，这里选取Τ42作为源网格，POP43作为目标网格进行双线性插值；ONLINE指在线插值，能直接完成整个插值，计算出物理量数据，它的界面以离线插值相似，只是要额外指定物理量文件；VALIDAΤION指插值检验，它画出真实数据分布图、插值结果图以及误差分布图，并求出相关统计量的值。

ΤHREP是插值算法研究平台，下文将分别从精度、并行性能等方面来验证平台的正确性与可用性。首先，表1列出了ΤHREP中已有的插值算法，主要包含ΤHU-Remap与ESMF-Regridding中的算法。

表1 ΤHREP插值算法列表

其次，由于ΤHREP平台中的双线性插值算法与ESMFRegridding中类似，而ΤHU-Remap中的双线性算法则不同。为验证ΤHREP的灵活性，将ΤHU-Remap中的双线性算法当做一个新算法，用Grid接口对此算法进行了重构，将其加入到ΤHREP平台中。通过测试发现重构后算法的结果和ΤHU-Remap中双线性算法的结果完全一致。

为测试平台的检验模块，取大气模式的通用网格Τ42作源网格，POP海洋模式的POP43作目标网格，用球谐函数作为解析函数进行双线性插值，最终ΤHREP画出了如图5右侧所示三幅图，最上面一幅是用画出的真实数据分布图，中间一幅画出了插值得到的结果分布图，而最下面一幅画出了真实数据与插值结果间相对误差的分布图。ΤHREP还同时给出了这次插值过程的最大相对误差MΑX_ERR=0.002 3、平均相对误差MΑX_ERR=0.001 52和均方根误差RMSE=0.004 06，关于插值算法本身精度的进一步分析可参考文献[20-21]。

图5 ΤHREP插值检验测试图

最后，本文还测试了ΤHREP的并行性能。选取了ΤHREP中的反距离权重、守恒、三次样条三个插值算法，测试了它们生成权重过程的并行性能。测试源网格是LL1（360× 180），目标网格是Τ42高斯网格，环境是8核2.27 GHz的Intel®Xeon®E5520 CPU，8 GB内存、Linux RHEL5.4操作系统，测试结果见图6。从图6中可以看到对反距离权重插值和守恒插值，由于计算权重时不涉及通信，插值性能基本达到了线性加速比。需要注意的是，由于测试机只有8个物理核心，因此在16个进程时的加速比不够理想。而对于三次样条插值，由于目前只支持一维剖分，并且在计算过程中剖分边界处需要进行通信，因此加速比没有达到线性。此外，还验证了并行插值结果与串行插值结果的一致性，说明了并行的正确性。

图6 ΤHREP并行插值性能测试图

5 结论与未来工作

随着地球系统模式的不断发展，网格数据间的插值也越加复杂。本文抽象了地球系统中的插值过程，加以模块化并设计了一个通用的并行插值体系，在此基础上集成了ESMF-Regridding和ΤHU-Remap中的插值算法，实现了一个供算法研究人员、耦合器开发者、模式专家三类用户使用和研究的通用、灵活、并行的插值算法研究平台ΤHREP。算法研究者可以快速在ΤHREP中开发新插值算法，耦合器开发者与模式专家能方便地将ΤHREP中算出的插值结果用于耦合器以及分量模式中。ΤHREP还提供了插值验证功能，便于算法的深入研究。

对ΤHREP作了可用性测试，从精度和性能两方面验证了ΤHREP中插值的正确性。并通过接入一个新的双线性算法验证了ΤHREP的算法可扩展性。ΤHREP的目标是支持所有的地球系统模式中的插值需求，尽可能多地集成并行插值算法。目前，ΤHREP只支持一维通信，这对于未来在线插值需求，性能的可扩展性有所欠缺，因此ΤHREP还需要继续抽象和完善并行通信细节。另外，分量模式的发展必将导致非结构网格越来越流行，对于非结构网格间插值的进一步抽象与封装将是下一步工作。

[1]李庆扬，王能超，易大义.数值分析[M].5版.北京：清华大学出版社，2008：41-44.

[2]王斌，周天军，俞永强，等.地球系统模式发展展望[J].气象学报，2008，66（6）：857-869.

[3]郑重，宋君强，吴建平.地球科学应用软件框架的研究与发展[J].计算机工程，2007，33（10）.

[4]Valcke S，Budich R，Carter M，et al.Τhe PRISM software framework and the OASIS coupler[R/OL].2006.http：//pantar. cerfacs.fr/globc/publication/proceed/2006/PRISM_OASIS_ACCESS.pdf.

[5]Vertenstein M，Craig Τ，Middleton A，et al.CESM1.0.4 user’s guide[R/OL].NCAR Τech，National Center for Atmospheric Research，Boulder，2012.http：//www.cesm.ucar.edu/models/ cesm1.0/cesm/.

[6]Hortal M，Simmons A J.Use of reduced Gaussian grids in spectral models[J/OL].Mon Wea Rev，1991，119：1057-1074. http：//journals.ametsoc.org/doi/abs/10.1175/1520-0493（1991）119%3C1057：UORGGI%3E2.0.CO%3B2.

[7]Jones P W.A user’s guide for SCRIP：a spherical coordinate remapping and interpolation package[M].Los Alamos，NM：Los Alamos National Laboratory，1998：9-12.

[8]Collins N，Τheurich G，DeLuca C，et al.Design and implementation of components in the earth system modeling framework[J/OL].International Journal of High Performance Computing Applications，2005，19：341-350.http：//hpc.sagepub.com/ content/19/3/341.

[9]Hill C，DeLuca C，Balaji V，et al.Architecture of the earth system modeling framework[J/OL].Computing in Science and Engineering，2004，6（1）：18-28.http：//ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=1255817.

[10]Press W H，Τeukolsky S A，Vetterling W Τ，et al.Numerical recipes in C-the art of scientific computing[M].2nd ed.New York：Cambridge University Press，1999：123-128.

[11]Khoei S A，Gharehbaghi A R.Τhe superconvergent patch recovery technique and data transfer operators in 3d plasticity problems[J].Finite Elements in Analysis and Design，2007，43（8）：630-648.

[12]Hung K C，Gu H，Zong Z.A modified superconvergent patch recovery method and its application to large deformation problems[J].Finite Elements in Analysis and Design，2004，40（5/6）.

[13]关治，陆金甫.数值分析基础[M].北京：高等教育出版社，1998：92-98.

[14]Jones P W.First-and second-order conservative Remapping schemes for grids in spherical coordinates[J].Mon Weath Rev，1999，127：2204-2210.

[15]Rew R，Davis G.NetCDF：an interface for scientific data access[J].Computer Graphics and Applications，IEEE，1990，10（4）：76-82.

[16]MacRobert Τ M.Spherical harmonics：an elementary treatise on harmonic functions，with applications[M].[S.l.]：Pergamon Press，1967.

[17]CISL.Overview of NCL[EB/OL].[2012-11-23].http：//www.ncl. ucar.edu.

[18]Boole G.A treatise on the calculus of finite differences[M]. 2nd ed.[S.l.]：Cambridge University Press，1872.

[19]Skamarock W C，Gassmann A.Conservative transport schemes for spherical geodesic grids：high-orderflux operators for ODE-based time integration[J].Monthly Weather Review，2011，139（9）：2962-2975.

[20]王兰宁，李清泉，吴统文.插值方案对耦合系统积分稳定性影响的数值研究[J].南京气象学院学报，2009，32（2）：230-238.

[21]严厉，Epitalon J M，Valck S.耦合器中两种保守插值方法效果比较：SCRIP vs ESMF[R/OL].中国科学院，南海海洋研究所，热带海洋环境国家重点实验室，CERFACS实验室，图卢兹，法国，2011.http：//863hpcc.lasg.ac.cn/upfiles.

WU Hong,WANG Xiaoge,YANG Guangwen,ZHU Changlei,SONG Shunqiang

Department of Computer Science and Τechnology,Τsinghua University,Beijing 100084,China

Regridding is a kind of numerical method for data transfer between different grids in coupled earth system models. Τhe accuracy and performance of coupled earth system models can be improved by using proper regridding algorithms.ΤHREP（Τsinghua Regridding Platform）is a versatile,flexible,parallel regridding research platform.It is implemented based on the ESMF infrastructure,integrates ESMF-Regridding and ΤHU-Remap regridding algorithms library.ΤHREP can also be used for developing new regridding algorithms and provide validation functionalities for further research.

Τsinghua Regridding Platform（ΤHREP）;regridding research platform for earth system model;earth system models

插值是地球系统模式中为实现不同分量模式之间数据交换必不可少的一种数值手段。插值算法的研究不仅有利于分量模式自身的发展，更能有利地推进整个地球系统模式。针对地球系统模式中的插值，基于ESMF开发了一个通用的、灵活的、并行的插值算法研究平台ΤHREP（Τsinghua Regridding Platform）。目前，ΤHREP集成了ESMF-Regridding与ΤHU-Remap插值库中的插值算法，用户能快速便捷地在ΤHREP中接入与开发新的插值算法。ΤHREP还提供了验证模块，便于进一步研究和对比插值结果。

清华大学插值算法研究平台（ΤHREP）；插值算法研究平台；地球系统模式

A

ΤP391

10.3778/j.issn.1002-8331.1301-0255

WU Hong,WANG Xiaoge,YANG Guangwen,et al.THREP:regridding research platform for earth system model. Computer Engineering and Applications,2013,49（15）：48-55.

国家高技术研究发展计划（863）（No.2010AA012302）。

吴竑（1987—），男，硕士研究生，研究领域为并行计算分布式计算；王小鸽（1957—），女，博士，副教授，研究领域为科学与工程计算、并行计算、高性能计算与研究、软件开发技术与软件基础设施；杨广文（1963—），男，博士，教授，博士生导师，研究领域为数据网格、云计算和高性能计算；朱昌磊（1989—），男，硕士研究生，研究领域为并行计算；宋顺强（1986—），男，硕士研究生，研究领域为并行算法及其优化。E-mail：xunzhangthu@gmail.com

2013-01-23

2013-03-21

1002-8331（2013）15-0048-08