GRAPES (Global and Regional Assimilation Prediction System) 是第1个完全由我国自主开发，拥有持续自主创新能力的多尺度数值天气预报系统^[1-3]。其四维变分同化系统 (4DVAR) 是在三维变分同化系统 (3DVAR)^[4-7]的基础上通过发展数值预报模式的切线性/伴随模式而来。它可将不同时刻、不同地区、不同性质的气象观测资料，包括最优插值客观分析难以应用的卫星、雷达等非常规观测资料作为一个整体同时进行考虑，得到与预报模式协调一致的初始场^[8-10]。

由于四维变分资料同化具有计算过程复杂、计算量大、内存使用量大等特点，因此对其主要组成部分切线性/伴随模式进行优化是必要的。由于伴随模式的计算复杂性，选择一种更优的实现方法对伴随模式的内存性能和计算性能的提高有重要影响。本文主要描述了在采用极限断点存储技术优化伴随模式的过程中所必需的一种中间数据管理支持技术。该技术在原有的栈结构基础上进行优化设计，有效促进了伴随模式优化，极大提高了其运行效率。

GRAPES伴随模式是在GRAPES并行计算软件框架^[11-13]的支持下，先依据非线性模式开发其切线性模式，得到切线算子，而后转置得到伴随算子，从而得到非线性模式的伴随模式。

设非线性变换F可表示成一系列中间变换f_k的集合，则预报模式可记为

(1)

切线模式记为

(2)

伴随模式记为

(3)

从式 (3) 可以看出，与非线性模式相比较，伴随模式在形式上类似非线性模式的反向计算。计算中的任一子段，需要对应该段非线性模式的基态。考虑到空间和时间效率，这些基态采取不同的实现方式会产生出不同形式的伴随模式，同时也有着不同的计算性能。

目前，得到基态较常用的方式有3种：全计算、全存储和断点存储策略^[14]。全计算与全存储策略在存储空间和计算时间上各具优势，而在另一方面却具有最差的性能。断点存储策略平衡了二者的劣势，是一种折衷手段，其存储空间和计算时间决定于程序中设置的断点数量。

极限断点存储策略的基本思想是在全存储策略实现下，由于基态的保存和恢复同样耗费CPU时间，如果该时间比从前一基态计算所需要的时间多，则采用从前一基态计算而不是保存来实现该基态。相比全存储策略有效减小内存开销，该方法比全存储技术内存占用要小，也比全存储技术执行时间要短。

在极限断点存储实现中，基态的实现方法有计算和存储两种。具体于每一个基态的实现方法决定于前一基态到自身的计算时间。这样，如果将基态保存在栈结构中，它们之间存在先进后出和先进先出的关系，并不符合栈中数据先进后出的这一特点。因而，必须有一种新的数据栈技术来支持极限断点存储策略的实现。

如图 1所示，使用极限断点存储技术实现时，子程序中内部断点X_i+1在X_i后保存，在X_i恢复后恢复。两者的关系与栈的结构形成了错位，而不能用普通的栈来管理这些数据。以下介绍一种复合的数据栈管理结构，可以有效满足使用极限断点存储技术实现GRAPES伴随模式所存在的混合数据管理需求。

图 1

图 1. 断点间先进先出关系示意图 Fig 1. The FIFO (First In First Out) relationship of checking-point

压、弹栈现场：指在压栈或弹栈过程中，某一时刻的栈顶地址和栈顶所在结点的结点地址。

接口：作为承上启下的纽带，应该具备如下几个功能：① 拥有多条子链的能力；② 可识别出每一条链的父、子链性质；③ 可获取父链的压栈或弹栈地址；④ 可获取每一条子链的压栈或弹栈地址。

节点：作为链的基本构件，除了具备可链接关系外，还需具备指示出父、子接口的功能。

基本链结构：是指通过接口构成父、子链关系的两条链结构。

多子链：是指在某个节点的同一接口上同时链接多个基本链结构的复合结构。

嵌套链：是指在基本链结构的子链节点再通过接口级联基本链结构所形成的复合结构。

嵌套多链栈：是指以多子链、嵌套链为基础链接形成的混合结构。

在实际运行中，一个接口可能含有多条子链或利用多个接口形成多级、多层次的链结构，总体上是一个嵌套的多链栈。

运行分为压栈过程和弹栈过程。

压栈时，将需先进先出的数据块压栈到当前链中，在当前节点创建接口构造子链，后进后出的数据块压栈进入到子链中；压栈完毕返回到接口，在接口保存子链的压栈地址；最后返回到父链；嵌套子链的压栈方式同上。

弹栈时，由于先进先出的数据块在父链中，因此可以首先获得该块数据，而后进入到子链中，并通过接口找到正确的子链及其弹栈位置，获取到后进后出的数据块。

嵌套多链栈是基于原数据栈基本结构，采取增加子链的方式，将不同关系的数据块分开，从而实现先进先出或先进后出综合功能。上一级链和下级链之间通过设置专用的接口来构成连接关系。图 2为采用接口方式实现的可嵌套的多链栈 (不同时刻进栈的数据块用不同灰度区分)。单独从某一条链来看，依然是一个基本的栈结构，具有后进先出的特点。

图 2

图 2. 嵌套多链栈结构示意图 Fig 2. The schematic diagram of nested and double chained stack

嵌套多链栈通过一个简单的接口实现了多子链和嵌套子链的统一管理，对上层调用屏蔽了复杂的实现过程。实现中保证了原程序中PUSH，LOOK，POP函数接口的一致性，整个栈接口简单，使用方便，容易理解，给应用程序开发者提供了统一的上层调用接口，减轻了编程和调试的负担，且有效支持了GRAPES伴随模式从断点存储技术到极限断点存储技术的转换。由于父、子链可以较为随意的压栈、弹栈，打乱了基本栈固有的先进后出的观念，使得功能更加灵活和强大。

优化的正确性可采用伴随模式的正确性检验方法，也可采用结果比较法，如图 3所示，采用结果比较可以看出优化前后的结果没有差异。

图 3

图 3. GRAPES模式6 h预报的850 hPa位势高度 (单位: gpm) (a) 优化前结果，(b) 优化后结果 Fig 3. An example of 6-hour 850 hPa geo-potential height forecast of GRAPES (unit: gpm) (a) the result before optimization, (b) the result after optimization

试验程序使用GRAPES有限区域切线伴随模式软件包，试验数据采用的初始场是2005年8月28日00：00(世界时，下同) 国家气象中心全球业务预报系统T213的分析场，网格为37×31×15，水平格距为1.0°，积分步长为600 s，在IBM-cluster 1600上运行。从图 3a与图 3b可以看出，优化前后的结果完全相同，没有差异，故说明优化方法是正确的。

为了检验在嵌套多链栈支持下极限断点存储技术的优化效果，对优化前、后的GRAPES伴随模式进行了测试对比。试验数据采用的初始场是2005年8月28日00：00国家气象中心全球业务预报系统T213的分析场，网格为37×31×15，水平格距为1.0°，积分步长为600 s。试验平台为Intel至强64位处理器架构的集群平台，Intel Fortran编译器，编译优化级别为3。

通过图 4与表 1的测试结果来看，采用极限断点存储技术实现伴随模式相比较断点存储技术虽造成堆栈最大值增加了1倍，但增加的内存在程序总使用内存中所占比例不高于30%，而程序的执行时间却下降了约50%，即达到200%的加速。

图 4

图 4. GRAPES伴随模式运行内存占用动态图 (a) 优化前后临时数据保存总量动态图，(b) 优化前后伴随模式占用内存总量动态图 Fig 4. The dynamic memory consumption of GRAPES adjiont mode (a) total amount of temporary data storage before and after optimization, (b) total memory consumption of adjoint mode before and after optimization

表 1

表 1 优化前后对比 Table 1 The result before and after optimization 时间 堆栈最大值/MB 程序占用总内存/MB ad_integrat执行时间/s 优化前 约55 约200 49.79 优化后 约110 约255 24.10

表 1 优化前后对比 Table 1 The result before and after optimization

本文主要介绍了一种基于基本链表栈基础上建立的一种灵活的数据栈存取管理结构——嵌套多链栈结构，并基于该结构使用极限断点存储技术对GARPES伴随代码进行优化改造。从测试的结果看，在整个栈的使用中，最长的两条子链在GCR求解^{[8, 11]}和插值运算^[15]中形成。每条子链保存的断点总量相对于整个程序占用内存来讲并不多，这种累积不会造成内存使用的大幅度增加，在付出较小代价的前提下获得较大的运行性能。

GRAPES伴随模式的性能提高，使得整个四维变分同化模式的运行时间也相应缩短，业务运行的实时性更强。同时由于区域模式和全球模式的一致性，这种优化可以直接移植到全球模式上，接口简单，优化方便，可维护性较强，更好地促进GRAPES产品全面投入业务运行。