在过去的30年里,切线性和伴随技术在气象和海洋等领域中被广泛应用,它是解决数据同化、集合预报和敏感性分析等问题的非常强大的工具[1-2]。目前切线性和伴随技术在欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts, ECMWF)的研发和业务中得到了非常成功的应用,如ECMWF的四维变分同化(four-dimensional variational data assimilation, 4DVar)系统[3],基于奇异向量技术的集合预报系统[4]和观测诊断工具的预报观测敏感性(forecast sensitivity to observation, FSO)分析系统[5],使ECMWF在全球数值预报领域长期占据领先地位。
GRAPES (Global and Regional Assimilation PrEdiction System)是第1个完全由我国自主开发、拥有持续自主创新能力的多尺度数值天气预报系统,其资料同化系统主要发展目标是通过建立对卫星观测资料的同化能力,提高数值预报的初值质量[6-7]。而4DVar方法是同化大量观测资料,尤其是卫星资料的最佳办法之一,它的核心是计算观测和模式背景场距离代价函数的极小化方法,而伴随模式是计算代价函数梯度的一个有效工具。所以选择基于伴随技术的4DVar系统是发展GRAPES资料同化系统的一种有效方法。
2010年国防科学技术大学基于GRAPES全球非线性模式(non-linear model,NLM)1.0版本,开发了GRAPES全球切线性模式(tangent linear model,TLM)和伴随模式(adjoint model,ADM)1.0版本[8]。但目前GRAPES全球NLM 2.0版本发展较快,采用新的并行框架和上游点插值方案提高计算效率[9],并且GRAPES资料同化系统也有了较大的改变[10],GRAPES全球TLM和ADM 1.0版本已经不能与GRAPES全球NLM 2.0版本兼容,且运行效率低,制约了GRAPES全球4DVar系统的发展。本文以面向GRAPES全球4DVar系统和GRAPES奇异向量集合预报系统[11]业务化发展为需求,基于GRAPES全球TLM和ADM 1.0版本,结合目前的GRAPES全球NLM 2.0版本重新设计和优化了TLM和ADM,开发高效、稳定、简洁的GRAPES全球TLM和ADM 2.0版本。
1 切线性和伴随模式简介4DVar中的梯度计算以及敏感性分析都需要伴随技术,对于随着时间向前演变的大气模式而言,ADM则是反向随时间演变。ADM的开发均针对TLM进行,要构建一个TLM,就需要NLM的线性化和轨迹处理[12]。
1.1 TLM介绍假设NLM用算子M表示,则大气状态x随时间ti的演变为
|
(1) |
小扰动dx随时间的切线性演变可以通过M的一阶近似为
|
(2) |
TLM算子是L=dM/dx。依据经典的泰勒展开公式,切线性程序的近似性检验公式为
|
(3) |
其中, α为一个渐变系数,|| ||表示矩阵的Frobonius范数,F(α)的值接近于1的有效位数就是切线性程序的近似精度。这是一个当扰动dx渐进逼近0时的一个理论测试方法,用来验证切线性代码的数值正确性。如果切线性程序正确,那么当α从1.0(此时dx的量级大概为x的0.1倍左右)开始,以0.1倍的速度变小时,F(α)的精度逐渐提高,但当α变得更小时,受到计算机的舍入误差影响,F(α)的精度逐渐变差[12]。
然而对于一个真实的数值预报模式系统,如4DVar系统中,一般扰动的大小为分析增量dxa=(xa-xb)(xa为分析场,xb为背景场)的幅度,分析扰动dxa在TLM中随时间的演变的近似性更加具有实际意义。定量上可以通过分析比较dxa在TLM下的演变和在NLM扰动演变之间的绝对平均偏差ε检验TLM的切线性近似程度[12]:
|
(4) |
切线性算子L的伴随算子L*,可以通过内积〈, 〉的方式定义:
|
(5) |
ADM是从时间终点反向积分到初始时刻,得到扰动变量在初始时刻的梯度和敏感性,它可以用来计算4DVar目标函数的梯度。一般伴随代码是相对于切线性代码逐句开发的,依据式(5),伴随代码正确性可通过扰动dx严格按照下面公式测试:
|
(6) |
即计算演化时刻切线性变量的内积及初始时刻变量和伴随变量的内积,如果两次内积的差异在计算误差范围内,则伴随代码正确,同样式(6)可用来测试整个ADM的正确性[12]。
切线性和伴随代码的开发是一项艰巨的工程,虽然有很多自动伴随工具可以提供切线性和伴随代码的自动生成[13],但事实上它们产生的代码可读性差、难以维护、运行效率较低,并且有很多不足。GRAPES全球TLM和ADM在借助TAPENADE工具[14]自动生成那些逻辑结构简单的程序,并进行手动调优,而那些逻辑结构复杂的代码为手动逐句开发。
2 GRAPES切线性和伴随模式设计GRAPES全球模式是一个非静力的、半隐式半拉格朗日的均匀经纬格点的数值预报模式[7]。为了开发一个适合4DVar和奇异向量计算的高效稳定运行的GRAPES全球TLM和ADM,本文对GRAPES全球NLM的程序进行了梳理,对那些简单的线性化程序进行归类,在其切线性和伴随开发时,不需要考虑基态问题,其切线性和伴随易于实现,在GRPAES全球NLM中有很多辅助性程序,如极区滤波、水平扩散和变量变换等均为该类程序。
对另外一些复杂的程序,如半隐式半拉格朗日上游点的计算和插值,以及Helmholts方程求解等,它们是NLM的核心代码,程序中包含着大量的循环迭代和控制结构,在开发其切线性和伴随程序时,需要考虑基态如何处理,为了开发高效运行的伴随程序,即需要从整体上来设计和考虑,也需要从细节上逐句调优。如GRPAES模式中Helmholts方程模块,虽然Helmholts求解器从理论上是线性的,但由于其在GRAPES模式中采用的是预条件的广义共轭余差(GCR)算法[15],Helmholts计算模块具有多次迭代、非线性、结构复杂且耗时的特点,其切线性和伴随的开发需要考虑基态的处理方法,使其计算效率达到最佳,下面会以其为例说明GRAPES全球TLM和ADM开发如何调优。
2.1 基态设计这里基态是指切线性和伴随程序中所需部分非线性程序变量,它是计算切线性和伴随程序的关键信息,必须保持基态在切线性和伴随相对应的程序中一致。在设计非线性程序的切线性和伴随代码时,基态处理是其中非常关键的部分,一般有以下两种方式:①计算法。切线性和伴随程序所需的基态通过计算非线性程序获取,简单直观,但需要耗费较多的计算时间。②存储法。切线性和伴随程序所需的基态由非线性程序预先计算存储起来,需要耗费更多的内存,但可以减少计算时间。
一般对单个程序,切线性程序运行的时候会同时运行非线性程序,所以其基态一般通过计算法得到,但伴随程序的计算顺序与切线性反向,需要采用保存基态的办法。而在真实的模式中,包含了大量的子程序,基态的处理较为复杂。
在考虑到4DVar的极小化和奇异向量求解计算中,需要多次反复迭代调用TLM和ADM,计算效率是业务化发展必须考虑的问题。GRAPES全球TLM和ADM的设计和开发中以存储法为主,即在运行TLM时,通过数据压栈方法将ADM中程序所需的基态保存到内存中,然后在运行ADM时通过数据出栈方法中取出相应的基态,并辅以计算法达到最佳计算效率。如针对Helmhots模块的伴随设计,在GRAPES全球ADM 1.0版本中,在每次ADM积分时,每计算1次Helmhots伴随模块,需要先计算1次Helmhots非线性模块来保存基态供伴随模块使用;而在GRAPES全球ADM 2.0版本里,在TLM计算Helmhots切线性模块时,将所有积分步的基态预存到内存中,提供给ADM的Helmhots伴随模块使用。这种通过基态存储法虽然需要更多的内存,但可以在ADM中减少1次Helmhots模块计算,最终可以为TLM和ADM计算共节省10%左右的时间(TLM略微增加,ADM减少20%左右)。
2.2 轨迹设计这里轨迹指NLM在积分过程中产生的每个积分步的模式状态变量。一般4DVar系统采用高低分辨率结合的方式,即外循环采用高分辨率NLM计算新息向量,内循环采用低分辨率的TLM和ADM进行极小化计算,如ECMWF内外循环分为两步:首先计算外循环后,同时将轨迹保存到硬盘,然后在进行低分辨率的内循环时读取轨迹并插值到低分辨率上给TLM和ADM使用[3]。而在GRAPES全球4DVar系统的设计,为了节约轨迹从硬盘读取时间,实现了内外循环的一体化,即在外循环计算新息向量时,将高分辨率NLM的轨迹插值并保存到内存,供给内循环极小化的TLM和ADM使用,这样在同一个执行程序中需要考虑模式从高分辨率到低分辨率的切换和保存轨迹,需要更多的内存,但可以大大提高4DVar的计算效率。
如图 1所示,GRAPES全球4DVar系统结合本节的轨迹设计和2.1节中基态设计。首先运行外循环的高分辨率GRAPES全球NLM计算新息向量,将轨迹插值到内循环的低分辨率并保存到内存中,供TLM和ADM使用;其次在内循环的极小化每次迭代中,读取轨迹并运行TLM,同时将TLM中的切线性程序计算的基态压栈到内存中(如前面提到的Helmhots求解),运行ADM时将相应的基态取出供相应的伴随模块使用。
|
|
| 图1 GRAPES全球4DVar中的轨迹和基态设计 Fig.1 The design of trajectory and base state in GRAPES Global 4DVar | |
2.3 切线性物理过程设计
切线性物理过程是GRAPES全球4DVar和奇异向量求解计算中一个非常重要的模块,为了减少NLM中物理过程计算效率低和非线性强的影响,一般采用简化的切线性物理过程。在GRAPES全球TLM和ADM引入切线性和伴随物理过程时,采用将物理过程轨迹和基态计算和扰动计算分离的方法。如图 2所示,在GRAPES全球NLM的积分过程中采用完备的物理过程方案[16],并保存相对应的TLM和ADM中切线性物理过程所需的轨迹和基态,如NLM采用完备的MRF方案[7],而TLM和ADM中采用简化的切线性MRF方案[17],NLM采用复杂的SASS积云对流方案[7],而TLM和ADM中采用简化的切线性CUDU积云对流方案[7],NLM采用双参数云物理过程方案[18],而TLM和ADM中采用简化的切线性LCEC87大尺度凝结方案[7]。
|
|
| 图2 切线性物理过程的设计方案 Fig.2 The design of tangent linear physics | |
如在4DVar计算中,在外循环NLM计算MRF方案时,将内循环TLM和ADM中边界层方案所需要的轨迹基态保存到内存并插值后供切线性MRF方案使用。显然这样的设计有以下3个优点:①TLM和ADM中切线性物理过程的轨迹基态为NLM的全物理过程提供,保证了切线性物理过程的高质量基态;②更高的计算效率,TLM和ADM只需要计算简化的物理过程,其计算效率不会受制于NLM中物理过程的低计算效率影响(如辐射和路面过程计算耗时);③非线性物理过程和简化的切线性物理过程松散耦合,易于独立发展和维护。
3 切线性和伴随检验GRAPES全球TLM和ADM的检验的试验设置:起报时间为2013年5月26日06:00(世界时,下同),模式水平分辨率为1°×1°,垂直为60层,积分时间为6 h,积分步长为1200 s,扰动向量采用GRAPES全球4DVar的分析增量,双精度编译。
3.1 切线性近似检验下面以GRPAES全球模式中的Helmhots求解模块为例,利用式(3)检验其切线性近似性。
表 1给出了GRAPES全球模式中Helmhots求解模块在模式积分第1步(h(0))和最后1步(h(6))的切线性近似测试结果,当α=10-6时,切线性近似精度可以达到7位精度,这符合式(3)中当α趋近于零时,F(α)趋近于1的结论,但随着α越来越小,计算精舍入误差开始影响近似结果,测试结果表明,GRAPES全球TLM的Helmhots模块的切线性版本是正确的。
|
|
表 1 Helmhots求解模块的切线性近似测试 Table 1 The tangent test of Helmhots subroutine |
GRAPES全球模式是一个非常复杂的系统,由多个非线性系统组成,且程序包含大量非线性和判断(if)开关结构[19],这些不连续性导致了切线性程序的近似性问题。利用切线性检验式(3)对一个真实数值预报模式进行整体的切线性近似测试,很难得到理论上的完美测试结果,如GRAPES模式中的半隐式半拉格朗日模块受到判断(if)开关结构的影响就比较大:①上游点计算中关于东西边界周期性变换(0和2π时)的判断(if)开关结构导致的不连续性,非线性测试时所引入的不同量级小扰动可能会产生不同上游点在周期性边界上的跳跃变化(如未加扰动时某一格点变量的上游点可能小于2π,如果加小扰动时会导致这个格点变量的上游点可能大于2π,需要减去2π),导致了不同量级情况下小扰动的非线性近似的差异,会导致切线性近似测试的失真;②上游点插值时,垂直方向上的边界层和中间层的判断(if)开关结构(边界层和中间层的插值方案不一样),测试时不同量级的小扰动导致其边界层的选择的不连续性,会导致切线性近似测试失真。
在GRAPES全球TLM的开发中,对子程序利用式(3)进行单元测试时,除半隐式半拉格朗日模块外,其他子程序可得到完美的切线性近似检验结果。但半隐式半拉格朗日模块的不连续问题会直接影响和反映到GRAPES全球TLM的切线性整体测试结果中,从而影响到整个TLM模式切线性近似性测试,不过在小扰动为分析增量量级时,TLM的测试结果仍近似为1,这说明TLM能够适用于真实4DVar的极小化计算。
3.2 切线性物理过程检验在4DVar和GRAPES奇异向量计算的实际应用中,应使用切线性物理过程,2.3节已介绍切线性物理过程与非线性物理过程的差异,所以对增加了切线性物理过程的GRAPES全球TLM无法按照式(3)测试。这里通过比较扰动随时间的切线性演变与非线性演变的近似程度分析切线性物理过程下TLM的合理性。
如图 3所示,模式第5层位温扰动在无切线性物理过程TLM在的6 h演化结果(图 3c)与无物理过程的NLM的演化结果(图 3a)非常相似,这说明在动力框架下TLM的切线性近似是合理的;但6 h演化结果与全物理过程下NLM的演化结果(图 3b)差异较大:其一是无物理过程下扰动增长幅度较大,其二是全物理过程下NLM的演化结果具有更多的小尺度信息;而增加了简化的切线性物理过程后的TLM的扰动幅度得到一定抑制(图 3d),大小与全物理过程下NLM的扰动演化相当,但依然未能模拟出更多的小尺度信息。
|
|
| 图3 位温扰动在模式第5层的6 h演化(a)无物理过程的非线性演变,(b)全物理过程的非线性演变,(c)无物理过程的切线性演变,(d)简化切线性物理过程的切线性演变 Fig.3 6 h evolution of the potential temperature perturbation at the 5th level of model (a) non-linear evolution of no physical process, (b) non-linear evolution of all physical processes, (c) tangent evolution of no physical process, (d) tangent evolution of simple tangent physical process | |
基于式(4), 图 4给出了位温扰动在有无切线性物理过程的TLM的6 h演化结果相对于全物理过程NLM下的平均绝对误差的纬向平均分布,由图 4可以看出, 没有切线性物理过程时,位温扰动的切线性演化在中高纬度地区底层相对于非线性扰动演化增长迅速(图 4a),增加了简化切线性物理过程后(图 4b),对位温扰动演化的增长具有很好的抑制作用,使其与非线性扰动演化更加接近。所以增加简化的切线性物理过程可能为GRAPES全球TLM模式带来好的效果,尤其对边界层的影响较为显著。
|
|
| 图4 位温扰动的平均绝对偏差(a)无物理过程的TLM,(b)简化切线性物理过程的TLM Fig.4 The mean absolute error of the potential temperature perturbation (a) TLM with no tangent physical process, (b) TLM with simple tangent physical process | |
综上可知,GRAPES全球TLM增加的几个简化切线性物理过程具有一定正效果,但模拟效果与全物理过程NLM还有较大差距,这也说明GRAPES全球TLM中需要增加更多与NLM一致的切线性物理过程方案。
3.3 伴随正确性测试GRAPES全球ADM以及其物理过程的伴随程序是严格按照GRAPES全球TLM对代码的逐行开发的,其正确性可以采用式(6)测试。下面是切线性物理过程下的GRAPES全球ADM的测试结果:
|
(7) |
可以看出,其有效的相等位数为15位,在计算精度范围内,表明简化切线性物理过程下的GRAPES全球ADM相对于TLM是正确的。
4 TLM和ADM效率分析一般在数值预报模式中,物理过程计算比较耗时,为了节约4DVar和奇异向量的计算时间,在第2章GRAPES全球TLM和ADM的结构设计基础上,采用了2.3节的物理过程轨迹基态计算与扰动计算分离的方法。为了更好地分析GRAPES全球TLM和ADM在不同情形下的计算效率,下面分两种情况进行性能测试:一是无物理过程的纯动力框架下的并行效率测试,为了分析具有相同程序结构的TLM和ADM相对于NLM的效率;二是在物理过程情况下的并行效率测试,为了测试真实情况下TLM和ADM的计算效率。
测试设置:模式垂直为60层,水平分辨率为1°×1°,预报时间为6 h,积分步长为1200 s,双精度编译,在中国气象局的IBM-Cluster1600机器上进行测试,每个计算节点为32个核。由于在GRAPES全球TLM和ADM的设计中需要保存大量轨迹和基态,需要更多的内存和节点,这里从4个计算节点开始测试,并且针对确定的计算核分成两种情况(情况1:小内存,32核/节点;情况2:大内存,16核/节点)测试,目的是测试内存对计算效率的影响。
4.1 无物理过程的性能测试表 2给出了无物理过程下的GRAPES全球NLM,TLM和ADM计算效率测试结果,由表 2可知:①GRAPES全球TLM和ADM的计算效率高。由表 2不同核在两种情况下(小内存和大内存)的计算时间上统计可以看出,小内存情况下,TLM的计算时间是NLM的1.1~1.5倍,ADM的计算时间是NLM的1.5~1.9倍;大内存情况下, TLM的计算时间是NLM的1.1~1.3倍,ADM的计算时间是NLM的1.6~2倍。一般认为TLM的计算时间是NLM的1.5倍,而ADM的计算时间是NLM的2~3倍,这说明GRAPES全球TLM和ADM的研发在计算效率上是成功的。②增加内存可以提高计算效率。将表 2中的小内存的计算时间(如NLM中的128计算核下的4节点×32核/节点, 26.51 s)减去大内存(8节点×16核/节点, 20.75 s)的计算时间,然后除以小内存的计算时间,得到图 5的结果。可以看出,同样的核数下,更多的内存可以显著提高模式计算效率,NLM平均可以提高24%左右的计算效率,TLM平均可以提高29%左右的计算效率,ADM平均可以提高18%左右的计算效率。由于TLM积分中会将ADM需要的基态压栈到内存中,所以它比NLM需要更多的内存,增加内存可以提高更多的计算效率;而在ADM计算中通过数据出栈的方式得到每积分步基态,会慢慢释放TLM中占用的动态内存,所以增加内存得到的效率提高略低一些。③GRAPES全球TLM和ADM的加速效率高。图 6给出了NLM,TLM和ADM在大内存情况的并行计算效率,可以看出ADM的并行加速效率最高,其次是TLM,而NLM较差,并且当计算核增加到1024(64节点×16核/节点)时,TLM和ADM的计算时间与NLM的计算时间几乎相等,这是由于ADM模式所需的内存最大,当节点增加时其加速效果高。
|
|
表 2 无物理过程下TLM和ADM的计算效率(单位:s) Table 2 The parallel efficiency of TLM and ADM without physical processes (unit:s) |
|
|
| 图5 增加内存的计算效率提高率 Fig.5 Parallel efficiency rates by increasing memory | |
|
|
| 图6 GRAPES全球NLM, TLM和ADM的加速效率 Fig.6 The speedup efficiency of GRAPES Global NLM, TLM and ADM | |
4.2 有物理过程的性能测试
采用与4.1节同样的方法对有物理过程的GRAPES全球NLM,TLM和ADM进行并行效率测试,其中NLM的物理过程设置与TLM,ADM的切线性和伴随物理过程设置参考2.3节。
由表 3可知,NLM中非线性物理过程耗时大,比较表 2和表 3,可以看出,GRAPES全球NLM中增加了物理过程后,计算效率大大降低,并行加速效率更低,其计算时间超过无物理过程时的2倍,最大可达到5倍以上(1024核(64节点×16核/节点)),从512核(32节点×16核/节点)到1024核(64节点×16核/节点)的加速效果为0,主要原因是GRAPES全球NLM中的所采用的辐射RRTM和边界层CoLM方案都很耗时,且它们采用的是模块自带的并行方案,与模式交互时存在效率瓶颈。切线性和伴随物理过程耗时少,比较表 2和表 3,可以看出,GRAPES全球TLM和ADM增加了相应的切线性和伴随物理过程后,其计算时间不超过无物理过程的1.4倍,最少只有无物理过程的1.1倍,且切线性物理过程下的TLM和ADM的并行加速效率依然很好,这是由于TLM和ADM中没有运行NLM的非线性物理过程(尤其是辐射和边界层CoLM方案), 3个切线性和伴随物理过程都是简化方案,耗时较小。
|
|
表 3 有物理过程的GRAPES全球模式的并行效率(单位:s) Table 3 The parallel efficiency with physical processes (unit:s) |
综上所述,在无物理过程的情况下,纯动力框架的GRAPES全球TLM和ADM相对于NLM模式相当高效。另外,GRPAES切线性物理过程的设计采用2.3节的物理过程轨迹基态计算与扰动计算分离的方法,避免了非线性物理过程带来的低计算效率,并能得到合理的切线性物理过程所需基态,也说明这是一个适合4DVar和奇异向量业务化发展的有效办法。
5 小结本文主要介绍了在GRAPES全球模式2.0版本下发展的TLM和ADM,为面向GRAPES全球4DVar和奇异向量集合预报的业务化发展需求,从效果和计算效率两个方面同时考虑重新设计和开发了GRAPES全球TLM和ADM,相当于TLM和ADM 1.0版本大大提高了并行计算效率,也验证了其合理性和正确性,同时也初步表明增加切线性物理过程能为TLM带来好处,使其朝着业务化发展迈出了一大步,加快了GRAPES全球4DVar的研发进程。当前GRAPES全球TLM和ADM的主要问题和不足:
1)由于半拉格朗日上游点和插值的非线性计算存在着判断(if)开关的不连续性导致的切线性近似的问题,需要从NLM里半拉格朗日的方案上进行更多的研究,采用合理的办法解决该问题。
2) GRAPES全球TLM模式中切线性物理过程仍然很少,且其所采用的物理过程均为极其简单的方案,与目前ECMWF所采用的切线性物理过程差距较大,这方面还需要投入更多的开发力量。
| [1] | Errico R M, VuKicevic T. Sensitivity analysis using an adjoint of the PSU-NCAR mesoscale model. Mon Wea Rev, 1992, 120: 1644–1660. DOI:10.1175/1520-0493(1992)120<1644:SAUAAO>2.0.CO;2 |
| [2] | Errico R M. What is an adjoint model. Bull Amer Meteor Soc, 1997, 78: 2577–2591. DOI:10.1175/1520-0477(1997)078<2577:WIAAM>2.0.CO;2 |
| [3] | Rabier F, Jarvinen H. The ECMWF operational implementation of four-dimensional variational assimilation.Ⅰ:Experimental results with simplified physics. Q J R Met Soc, 2000, 126, (564): 11–43. |
| [4] | Molteni F R B, Palmer T N, Petroloagis T. The ECMWF ensemble prediction system:Methodlogy and validation. Q J R Met Soc, 1996, 119: 269–298. |
| [5] | Cardinal C. Forecast Sensitivity to Observation (FSO) as a Diagnostic Tool. ECMWF Tech Memo, 2009: 26. |
| [6] | 陈德辉, 沈学顺. 新一代数值预报系统GRAPES研究进展. 应用气象学报, 2006, 17, (6): 773–777. |
| [7] | 薛纪善, 陈德辉. 数值预报系统GRAPES的科学设计与应用. 北京: 科学出版社, 2008. |
| [8] | 任迪生, 沈学顺, 薛纪善, 等. GRAPES伴随模式底层数据栈优化. 应用气象学报, 2011, 22, (3): 362–366. |
| [9] | 蒋沁谷, 金之雁. GRAPES全球模式MPI与OpenMP混合并行方案. 应用气象学报, 2014, 25, (5): 581–591. |
| [10] | 刘艳, 薛纪善, 张林, 等. GRAPES全球三维变分同化系统的检验与诊断. 应用气象学报, 2016, 27, (1): 1–15. |
| [11] | 刘永柱, 沈学顺, 李晓莉. 基于总能量模的GRAPES全球模式奇异向量扰动研究. 气象学报, 2013, 71, (3): 517–526. |
| [12] | Janiskova M, Lopez P. Linearized Physics for Data Assimilation at ECMWF. ECMWF Tech Memo, 2012: 26. |
| [13] | Giering R, Kaminski T. Recipes for adjoint code construction. ACM Trans Math Software, 1998, 24, (4): 437–474. DOI:10.1145/293686.293695 |
| [14] | Laurent H.TAPENADE, Automatic Dierentiation by Program Transformation.http://www-sop.inria.fr/tropics.2007. |
| [15] | 宋君强, 伍湘君. GRAPES模式中Helmhothz方程两种求解方法的对比研究. 计算机工程与科学, 2011, 33, (11): 65–70. |
| [16] | 徐国强, 陈德辉. GRAPES物理过程的优化试验及程序结构设计. 科学通报, 2008, (20): 2428–2434. |
| [17] | 张林, 朱宗申. GRAPES模式切线性垂直扩散方案的误差分析和改进. 应用气象学报, 2008, 19, (2): 194–200. |
| [18] | 刘奇俊, 胡志晋, 周秀骥. HLAFS显式云降水方案及其对暴雨和云的模拟(I)云降水显式方案. 应用气象学报, 2003, 14, (增刊): 60–67. |
| [19] | Zou X. Tangent linear and adjoint of "on-off" processes and their feasibility for use in 4-dimensional variational data assimilation. Tells, 1997, 49: 3–31. |
2017, 28 (1): 62-71



