1. 引言
在计算科学领域，自动微分(automatic differentiation，简称为AD)也称为算法微分或计算微分，是一种用数值方法获得由计算机程序描述的函数导数的方法[1] 。AD通过链式求导法则用解析方法获得计算机精度的函数一阶或高阶导数，因而避免了差分近似数值离散时带来的截断误差。AD的作用对象为已分解为初等函数和基本运算的程序代码，因而避免了符号微分(计算机代数)求导运算时可能出现无效代码及难于转换为单个表达式的问题。利用AD实现源程序到源程序的自动转换，可极大提高代码转换效率，以前利用手工方法可能需要几年才能完成的工作利用AD技术可在几个月甚至几天内完成。AD有两种基本运算方式，即向前运算和逆运算。在向前运算方式下，导数运算通过链式法则向前传递，一阶导数矢量称为切向量。在逆运算方式下，对中间变量的导数运算通过链式法则反向传递，一阶导数矢量称为梯度。向前运算方式称为切线性方式，逆运算方式称为伴随方式。这两种方法都可用来分析输出对输入的敏感性，但切线性方式适于输出变量多、输入变量少情形，而伴随方式适于输出变量少、输入变量多情形。与有限差分等计算敏感性的常规方法相比，AD具有计算时间少和计算精度高的优点。这是因为，伴随方法的计算量与敏感性所针对的输入变量个数无关，AD计算导数没有截断误差。由于AD两种基本运算方式互为逆运算，AD软件一般可同时给出程序代码的切线性运算程序和伴随运算程序。在地球科学(如大气科学)领域，数值求解预报方程的程序一般称为数值模式，其相应的切线性运算程序和伴随运算程序分别称为切线性模式和伴随模式。当前，AD工具已有多种，如在地球科学领域应用较多的ADIFOR [2] 、Tapenade [3] 、TAF [4] 及OpenAD/F [5] 等。其中，OpenAD/F是开源软件。使用开源软件的优点不仅在于成本低廉，而且在使用过程中可对其进行发展，使其更适用于解决特定问题，因而具有更大的灵活性。
伴随模式是分析输出对输入敏感性的有力工具，可应用于最优化问题、稳定性分析和参数估计等研究[6] 。在过去20多年里，伴随方法已经在大气科学[7] -[10] 、计算流体力学[11] 、核物理[12] 及工程设计优化[13] 等许多领域得到应用。伴随通过求取微分运算方法获得，而微分运算分在连续条件下进行和离散条件下进行两种情况。相应地，获得伴随方法有两种，分别称之为连续伴随和离散伴随。连续伴随先对连续形式方程求微分获得连续形式的伴随方程，然后再离散化。而离散伴随是在离散形式方程的基础上求微分。如前所述，利用AD技术可以极高的效率实现离散伴随。
数值天气预报是一种大规模科学计算问题，涉及的动力理论、数值计算等方面的问题较复杂。因此，为突出物理意义及增加所研究问题的针对性，在理论分析时常采用简化但保留关键动力学过程的方程组开展研究工作。正压大气原始方程组是在静力平衡条件下把大气流体视为有一自由面的均匀不可压缩流体所得到的简化结果。正压原始方程组虽然形式简单，但解中包含非线性效应以及与斜压原始方程组解中相似的波动结构，因而在地球物理流体力学研究以及数值算法的试验中得到广泛的应用[14] 。
本文将在开源自动微分工具辅助下实现正压大气原始方程模式的伴随模式，并在此基础上开展台风强度对初值和地形敏感性试验，并与采用非开源自动微分工具实现的伴随模式试验结果[15] 进行对比。文章第2部分给出正压大气原始方程模式伴随模式的实现，第3部分给出数值试验及结果，最后(第4部分)为论文总结。
2. 伴随模式实现及试验设计
2.1. 自动微分工具OpenAD
OpenAD(www.mcs.anl.gov/OpenAD/)是美国多家大学和研究机构合作发展出的自动微分产品，分Fortran语言(OpenAD/F)和c/c++语言(ADIC)两个软件包，二者基于相同的框架体系。开发OpenAD的基本设想是：此工具主要用于解决大气和海洋大规模科学计算数值模式资料同化中的梯度计算问题，使用者具有程序调试及编写经验。因而特别强调其应具有使用灵活性、采用开源组件及模块化特点，而不强调其对语言特征的支持能力。所以，使用OpenAD时，一般需对数值计算程序进行修改，以避免使用不常用的语言特征，从而满足OpenAD的要求。即不能象使用某些商业软件那样将其简单地视为黑箱，这正是许多研究者所希望的。本文使用的是OpenAD/F [5] ，其由前端(针对Fortran语言)、抽象层、代码分析和AD转换引擎四个模块组成。前端模块利用Open64 (mfef90)和whirl对Fortran程序进行分析并转换成特定内部格式(也称为whirl格式)，也可将特定内部格式转换成Fortran程序(whirl2f)，所以前端既是入口也是出口。抽象层模块的主要功能是定义一种扩展置标语言抽象接口格式(简称为XAIF)，从而使不同的自动微分算法可用一种中性语言形式表达，此模块利用工具whirl2xaif(将特定内部格式转换成中性语言形式)和xaif2whirl(将中性语言形式转换成特定内部格式)与前端发生联系。代码分析模块采用OpenAnalysis实现别名分析、DU(Define-Use)/UD(Use-Define)链及输入输出分析的功能，OpenAnalysis询问前端并将分析结果返回前端。AD转换引擎模块使用xaifBooster实现数值内核的自动微分变换，AD转换是在中性语言形式下进行的。
2.2. 伴随模式实现
正压原始方程组及其离散方案参见文献15。数值模式采用Fortran77语言编写，OpenAD/F适用Fortran90标准程序。如前所述，OpenAD/F不特别强调其对语言特征的支持能力，要求程序只使用Fortran90常用语言特征。因此，首先需对原源程序进行改写。改写之处主要有：将common方式数据传递方式改为module方式；将隐式数据类型定义方式该为显式方式，尤其明确浮点数据类型使用(一般采用real*8 (64位浮点数)，如果数据输入输出需使用real*4 (32位浮点数)，也要明确定义。否则，后面利用OpenAD/F进行代码转换时数据类型会出现问题)；调整时间积分实现方式，将1重循环积分计算改写为2重循环(例如，如果总共积分1200步，可改写为在内重循环中积分40步，在外重循环中循环30次方式。两个循环运算均编写为独立的子程序(为便于说明，分别称其为外循环子程序和内循环子程序)，即在外循环子程序中循环调用30次内循环子程序，在内循环子程序中积分40步。这样做的目的是，方便在伴随模式实现时采用2层检查点(checkpointing)方式，以节省存储空间。如不这样做，可能会因为存储数组过大超过计算机许可范围而使伴随模式无法运行)；去除goto等Fortran90标准不建议使用的语句。
完成前面数值模式改写之后，在程序中添加注释指令(Fortran编译器视其为注释，而OpenAD/F视其为指令)，主要分为两类，一类为在实现时间积分的子程序(即调用外循环子程序的子程序)中添加：
c$openad INDEPENDENT(z0)
c$openad INDEPENDENT(zo)
c$openad DEPENDENT(fc)
其中，前2条语句声明模式初值(位势高度z0)和地形高度(zo)为独立变量，第3条语句声明代价函数(fc)为依赖变量。另一类为在每个子程序定义语句前添加模板注释指令，有以下四种：
c$openad XXX Template ad_template.joint_split_oif.f
c$openad XXX Template ad_template.joint_split_iif.f
c$openad XXX Template ad_template_timing.joint.f
c$openad XXX Template OADrts/ad_template.split.f
这里前2条语句分别声明外循环子程序和内循环子程序模板，第3句声明与外循环子程序处于同等地位但没有循环积分运算的子程序模板，第4句声明在循环内部参与时间积分运算的子程序模板。使用者可根据需要对这些模板进行修改或给出新的模板。然后利用OpenAD/F系统中的python程序对子程序进行预处理，预处理的目的是对源程序施加检查并将其转换成标准形式，以便于后续处理。OpenAD/F前端将预处理后的程序转为内部格式，OpenAnalysis对此结果进行分析，利用OpenAD Fortran工具包whirl2xaif工具将分析结果写成中间格式文件。在此基础上，利用AD转换引擎实现伴随变换，结果仍被写为中间格式。之后，利用xaif2whirl将中间格式转换为内部格式，再利用whirl2f反向分析得到的内部格式文件，生成Fortran语言文件。最后，利用python程序对子程序进行后处理，得到最终的Fortran子程序文件。
根据OpenAD/F数据类型定义方式，在主程序中调用OpenAD/F数据类型定义模块并对主程序数据类型重新定义，以使得与生成的Fortran程序数据类型定义方式一致。此外，还要在主程序中添加一些运行伴随模式所需要的初始化操作语句及运行方式选择语句等。之后，对主程序及生成的子程序进行编译、链接，即可生成可执行文件。
在借助OpenAD/F生成的数值计算程序中，既包含原来的数值预报模式(经过了OpenAD/F改写，但数值计算实现没有变化。一般称其为正演模式)，还包括与数值预报模式相对应的伴随模式。所以适当设置模式参数，就可利用借助OpenAD/F生成的数值预报模式及其伴随模式开展数值模拟试验。
2.3. 数值模拟试验设计
模式网格距取为50 km，计算区域为Lambert地图投影平面上以北纬40度、东经110度为中心的131 × 151网格范围。模式中定义台风强度(用台风中心位势高度表征)为代价函数。以2006年5月17日20时MACC(Monitoring Atmospheric Composition and Climate)数据集500百帕等压面资料为初值，针对0601号台风“珍珠”个例开展数值模拟试验。此试验正演模式、网格配置及天气个例选取与文献15相同，但伴随模式不同。文献15 伴随模式是在非开源自动微分工具Tapenade辅助下实现的。此处选取相同天气个例及正演模式的原因是可以利用文献15结果对本文结果进行验证。
3. 数值模拟试验结果
台风“珍珠”是2006年登陆我国的第一号台风(标号为0601)。其形成于5月10日，在15日发展为强台风。在500百帕等压面上，5月17日20时涡旋中心处于北纬22.1度、东经116.6度位置，12小时后涡旋中心移至北纬24.8度、东经117.5度处，强度明显减弱[15] 。与MACC数据集结果相比，模拟的12小时后的涡旋中心强度偏强，移动速度略偏慢(图略，或参见文献[15] 图2(b))。
图1为积分伴随模式12小时得到的代价函数对初值和地形的梯度。可以看出，2006年5月18日08时登陆台风中心位势高度对初始时刻(即2006年5月17日20时)台风中心位势高度的梯度为正值，而对初始时刻台风中心西南和东北方向的外围区域位势高度的梯度为负值(图1(a))。这说明：初始时刻台风强度越强则模拟的12小时后登陆台风中心强度越强；但初始时刻台风中心外围不同方位环流对登陆台风中心强度影响是有差异的，西南和东北方位环流越强则模拟的登陆台风中心强度越强，而台风中心北部和南部环流越强则模拟的登陆台风中心强度越弱。登陆台风中心位势高度对登陆区域及台湾岛地形高度的梯度为正值，但对陆地上113˚E~115˚E区域地形高度的梯度为负值(图1(b))。这意味着，台湾岛和登陆处地形越高则模拟的登陆台风中心强度越弱，但台风涡旋影响范围其它区域地形可有与此相异作用。这些结果与文献15完全一致，二者差异(见图2)数值上比梯度值计算结果小2个数量级以上(比较图2与图1)。出现差异的主要原因是两个模式浮点数据类型的精度不一样，本文中为双精度而文献15中为单精度。
4. 小结与结论
自动微分工具已成为辅助伴随模式实现的有力手段，开展相关技术应用研究具有重要实用价值。
(a) (b)
Figure 1. Gradients of the cost function to the initial value and the topography. (a) to the initial value, the contour interval is 0.005; (b) to the topography, the contour interval is 1.0e−6 gpm/m
图1. 代价函数对初值和地形的梯度；(a)为对初值的梯度，等值线间隔为0.005；(b)为对地形的梯度，等值线间隔为1.0e−6 gpm/m
(a) (b)
Figure 2. Differences of results from the two adjoint models. (a) gradient to the initial value, the contour interval is 5.0e−8; (b) gradient to the topography, the contour interval is 3.0e−9 gpm/m
图2. 两个伴随模式结果差异；(a) 为对初值的梯度，等值线间隔为5.0e−8；(b) 为对地形的梯度, 等值线间隔为3.0e−9 gpm/m
OpenAD/F因使用开源组件而具有使用灵活、模块化及可发展等特点，更受到具有编程基础且强调针对性使用并发展(而不是简单地接受一个黑箱)的科技工作者的欢迎。论文在该开源自动微分工具辅助下开展了伴随模式实现技术的探索工作，所用数值模式为大气科学领域常用来讨论理论问题的正压原始方程模式。利用所得伴随模式进行了关于登陆台风强度对初值和地形敏感性方面的分析工作，所得结论与采用一个非开源自动微分工具得到的伴随模式分析结果一致。这表明，此次伴随模式的实现是成功的，这为利用OpenAD/F实现更复杂数值模式伴随模式工作打下了基础。需要说明的是，与利用Tapenade实现伴随模式工作相比，利用OpenAD/F实现伴随模式时需对原模式进行较多修改，且模式运行速度慢。运行速度慢与模式采用双精度数据类型有关，但更多地可能是优化工作有待深入。
基金项目
本文受国家自然科学基金项目(41276190)和国家公益行业专项(GYHY201506011)资助。