数值winds组合

数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-1-目录目录.
1绪论.
51、什么是数值预报和数值模拟.
52、两者有何异同.
53、课程目的、要求、内容.
54、参考书目.
6第一章大气运动方程组.
7第一节基本方程组.
71、运动方程.
72、连续方程.
83、状态方程.
84、热力学方程.
85、水汽方程.
96、闭合方程组.
9第二节球坐标系方程组.
101、球坐标系.
102、局地直角坐标系.
11第三节铅直坐标变换.
111、z坐标系.
112、P坐标系123、δ坐标系.
124、θ坐标系.
135、混合坐标系.
136、pzδ、、坐标系的比较.
14第四节地图投影坐标系.
151、为何采用地图投影坐标系.
152、常用的地图投影.
153、地图投影坐标系中的大气运动方程组.
16第二章大气运动的整体性质与能量约束.
19第一节大尺度大气运动的特性.
191、准地转特性.
192、准静力平衡.
203、准水平特征.
214、准无辐散.
22第二节小尺度运动的特性(包率内斯克Boussinesq近似)221、特性.
222、Boussinesq近似22目目录录-2-第三节大气中的几个守恒量.
231、总质量守恒.
232、绝热大气位温与位温平方守恒.
233、绝热无摩擦大气总能量守恒.
244、绝热无摩擦大气位涡和位涡拟能守恒.
24第三章大气中的波动.
25第一节大气中的混合波解.
251、大气波动基本方程.
252、波解分析.
26第二节各种波动特征.
261、声波.
262、重力内波.
263、惯性波.
264、惯性重力内波.
275、Rossby波.
27第四章数值计算法.
29第一节水平空间差分网格及其差分计算格式.
29第二节时间差分.
31第三节各种常用的数值时间积分格式.
321、取两个时间层的情况.
332、取三个时间层的格式.
35第四节差分近似的相容性、精确性、收敛性和稳定性.
361、差分近似的相容性.
362、差分近似的精确性.
373、差分近似的收敛性.
374、各种时间积分格式的稳定性.
38第五节起步问题.
441、多步法起步.
442、用稳定的两个时间层格式起步.
45第六节常用差分格式的系统性误差.
451、波速偏慢问题.
452、差分方程的计算解(虚假的计算波)45第七节非线性计算不稳定.
46第八节平滑和滤波.
471、5点平滑公式.
472、9点平滑公式.
48第九节水平边界条件.
481、固定边界和刚体边界.
482、海绵边界(可比较稳定)493、外推边界.
494、网格嵌套.
50第十节垂直方向的差分.
501、垂直分层.
50数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-3-2、垂直差分.
51第五章正压原始方程模式.
53第一节历史回顾.
53第二节过滤模式.
54第三节正压原始方程组.
55第四节正压原始方程模式的物理属性.
55第六章斜压原始方程模式.
57第一节模式方程组及其积分关系.
571、垂直区间约束关系.
572、质量守恒关系.
573、变量的个别变化及其通量形式.
574、水平气压梯度力的垂直积分.
585、总能量守恒.
58第二节变量配置与垂直差分格式.
581、垂直分层与变量配置.
582、垂直差分格式的构造.
58第三节斜压原始方程模式计算方案.
59第七章初值处理.
61第一节客观分析.
611、水平插值.
612、垂直插值.
65第二节资料初始化.
661、静处理方法.
662、动处理方法.
683、变分处理方法.
68第三节四维同化.
681、间歇资料同化.
682、连续资料同化.
69第八章模式中物理过程的处理方法.
71第一节辐射过程:711、太阳辐射:712、地气系统的红外热辐射(长波辐射)79第二节湍流和湍流扩散作用.
83第三节边界层处理.
85第四节地形作用.
881、地形的热力作用.
882、地形动力作用.
883、模式中对地形的处理.
88第五节凝结加热.
89附录一Fortran在大气科学数值模拟中的运用.
91一、FORTRAN基础911、字符集.
912、程序书写格式.
91目目录录-4-3、变量.
924、变量名称的取名原则.
935、程序结构.
936、输入/输出947、选择结构.
958、循环结构.
969、数组的声明与使用.
9610、文件的基本读写操作.
9711、子程序SUBROUTINE9712、函数FUNCTION98二、差分法Fortran实例.
98附录二Grads的应用.
991、基本使用操作.
992、Fortran写Grads数据文件1023、Matlab写Grads数据文件.
103实习一差分方法.
105实习二时间积分格式.
107实习三水平内插.
113实习四USTC9层模式115数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-5-绪论1、什么是数值预报和数值模拟过去讲数值预报是指数值天气预报,即用电子计算机来求解大气动力学方程,从而制作天气预报的方法;近年来发展出气候数值预报和环境数值预报.

数值模拟主要是采用数值方法来模拟大气的运动和变化规律.
大气运动非常复杂,描述它的方程组也非常复杂,目前还没有解析解,只能求它的数值解.
计算机是必备条件.
2、两者有何异同相同点:采用的基本方程、计算方法相同,得到的结果有不少是一样的不同点:1、两者的目的和任务不一样数值模拟:研究、实验先行,为预报打基础.
如:⑴设计模式,对各种计算方法(如各种差分格式、谱模式等,包括计算Jacobi的方法就好几种)进行实验,找出最佳.
⑵对各种物理因子,首先设计计算方案,然后实验,看其对预报是否重要.
举例:O3、SST…搞好了,提供给预报用.
数值预报:已有成型的模式和方案,每天把初值送入,计算结果,分析预报2、两者使用的资料不一样:预报用当时实况,模拟则不一定3、成功的数值模式就可转去应用,就变成数值预报模式3、课程目的、要求、内容目的:学会使用计算机作数值预报、数值模拟的原理和方法;熟悉Fortran、Matlab、GrADS等软件在大气科学研究中的应用内容:针对本科,本课程主要讲解数值天气预报.
绪绪论论-6-要求:(1)认真听讲,作笔记(2)期末考试+平时实习(一份实习报告)4、参考书目《数值天气预报基础》周毅等编著气象出版社《数字天气预报》沈桐立等编著气象出版社数值天气预报的基础知识(大气运动基本方程组)数值计算方法(差分法求解偏微分方程组)数值天气预报模式(正压、斜压原始方程模式)初始条件和边界条件(初值、侧边界条件)模式物理过程的处理方法数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-7-第一章大气运动方程组数值天气预报:第一节基本方程组1、运动方程惯性参考系,又名静止参考系、绝对参考系牛顿第二运动定律:单位质量空气所受合外力aaiidVfdt=∑,其中aadVdtt=+地球表面(旋转参考系):3aVVr=+Ω*其中Ω为地球旋转角速度,517.
29210sΩ*=,3V为相对速度即:绝对速度=相对速度+因地球自转而产生的牵连速度对于任意矢量A有adAdAAdtdt=+Ω*,(注:令Ar=可推得上式)对大气,令aAV=,有332()aadVdVVrdtdt=+Ω*+Ω*Ω*其中各项依次为绝对加速度,相对加速度,科氏加速度,向心加速度初始条件、边界条件数值积分基本方程组未来时刻气象要素的空间分布动力热力动量守恒定律(牛2)质量守恒定律气体实验定律热力学第一定律水汽守恒定律……第第一一章章大大气气运运动动方方程程组组-8-33333112()2adVPVrgFPgVFdtρρ=Ω*Ω*Ω*Ω*+各项依次为气压梯度力、科氏力、惯性离心力、地球引力、摩擦力,其中科氏力和惯性离心力合称视视力以区别于真实力,惯性离心力和地球引力合称重力g,3aGMgrr=,气象上通常取1g9.
8ms=)2、连续方程(质量守恒定律的数学表达式)0ddtρδτ=(xyzδτδδδ=有限物质体积元)推得:330dVdtρρ+=或330Vtρρ+=其中33Vρ为质量散度,速度散度33uvwVxyz=++,表示密度变化由幅散幅合引起.
3、状态方程理想气体状态方程:pRTρ=(R为比气体常数)干空气:11287RJkgK=湿空气:(10.
61)wRqR=+,其中比湿wqρρ=引进虚温:(10.
61)vTqT=+湿空气状态方程:wvpRTRTρρ==4、热力学方程理想气体热力学第一定律:vdTdCpQdtdtα=+数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-9-其中质量定容热容11717vCJkgK=,比容1αρ=,各项为:单位质量理想气体的内能变化率,理想气体压缩功率,非绝热加热率(水汽相变潜热,辐射湍流对热量的输送等)11,1004pvpRTCCRJkgKα∵质量定压热容pdTdpCQdtdtα∴=,用位温表示lnpdCTQdtθ=对于干绝热过程Q=0,则有lnln0pCdTRdp=定义位温00()pRCpTpθ=,其中001000phPa=则有位温守恒:lnlnln0ppCdCdTRdpθ==对于非绝热过程,定义比熵lnpsCθ=,则热力学方程为:dsQdtT=5、水汽方程类似于连续方程(但考虑源汇):33wwdVsdtρρ每单位体积中提供的水汽质量)将比湿wqρρ=代入得:331qsVqStρ+==(水汽源,单位时间内为单位质量空气提供的水汽量)其中33Vq为水汽输送量6、闭合方程组基本方程组:运动方程33312dVPgVFdtρΩ*+连续方程330dVdtρρ+=第第一一章章大大气气运运动动方方程程组组-10-热力学方程pdTdpCQdtdtα=水汽方程1dqSdt=状态方程pRTρ=1FQS、、已知→3TVpρ、、、、q闭合模式大气为绝热,无摩擦,干燥大气10FQS=、、→3TVpρ、、、闭合第二节球坐标系方程组1、球坐标系原点O:地心λ经度φ(θ)纬度(余纬)r径向长度cos,,ddrurvrwdtdt===取薄层近似,即rza=+(z为海拔高度,6371akm=),zara球坐标系中基本方程组取薄层近似后的形式为:21tancos1tan1(tan)0coszpdupuvfvFdtaadvpufuFdtaadwpgFdtzduvwvdtaazadTdpCQdtdtpRTλρλρρρρλαρ=+==其中,2sinf=Ω,cosduvwdttaazλ+++虽然简化,但仍精确方程组仍十分复杂,在极点为奇点,解析分析和数值计算难度大数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-11-2、局地直角坐标系球坐标系的简化形式:保持了球面坐标系的框架,忽略了球面曲率的影响111()0xyzpdupfvFdtxdvpfuFdtydwpgFdtzduvwdtxyzdTdpCQdtdtpRTρρρρραρ=++=+=+==此方程组常用于尺度不太大的问题,如台风、积云数字模拟,但由于这两问题旋转和轴对称性较强,所以常常将它们转换到组坐标系来设计模式第三节铅直坐标变换Z坐标系便于显示大气运动演变图象,但在数值天气预报中有时并不方便.

坐标变换:P坐标系,θ坐标系,δ坐标系变换推导见P12~151、z坐标系(见前节球坐标系和局地直角坐标系)边界条件:上边界:,0zwρ→∞=下边界:0,0,sssszwzzwVz====平坦地:起伏地:球坐标系局地直角坐标大气运动基本方程组(海拔高度z坐标系)第第一一章章大大气气运运动动方方程程组组-12-2、P坐标系只有在静力平衡时才使用pgzρ=2221011pppTqdVgzfkVFdtVpdTdQCFdtdtzpgpRTdqScFdtωωρρρ+==+====+其中ppdVdttpω=++边界条件:上边界:0,0pω→→下边界:,ssssspppVptω==+平坦地:0sω=起伏地:ssssssspVpwgwzωρ3、δ坐标系又名地形坐标系令TsTppppδ=,取0Tp=,则sppδ=数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-13-(ln)()()01ssssTppqdVgzRTpfkVFdtppVptzRTgdTRTdpdQFdtCpdtCdtdqcFdtpRTδδδδδδρ==++=+=其中1()ssdVdttdpddtpdtδδδδδωδ=++==边界条件:上边界0,0δδ→→,下边界1,0sδδ→→4、θ坐标系(等熵坐标系,略,见P21~22)5、混合坐标系P-δ混合(上层P坐标,下层δ坐标)msmppppδ=或1mmsmmppppδδδ=+要使mp处从上面算得mp的w值与从下面算上去的w值相等,保持大气运动连续性.
边界条件:0,0TTppw===第第一一章章大大气气运运动动方方程程组组-14-,mmppww==0,0mmδδ==spp=1,0ssδδδ===6、pzδ、、坐标系的比较δ坐标系:优点:1、地形处理方便,在地形面上0oδ=2、连续方程变成地面气压倾向方程3、可以做谱展开缺点:1、气压梯度变成两大项之小差,不易算准,设计模式时要小心2、只能描述静力平衡运动z坐标系:优点:1、是严格的正交坐标系2、下边界不随时间变化3、等高面水平面4、可描述准静力平衡和非静力平衡运动缺点:1、与山地相截,计算面出现空洞,在其上计算梯度,选取边界条件困难2、不能做谱展开p坐标系:优点:1、等压面接近水平2、气压梯度力形式简单3、连续方程形式简单,便于计算4、等压面上的Tq、等量是观测量,直接得到,便于计算(降水等)缺点:1、下边界面不是坐标面,p随时空变化,地形起伏区等压面也与地形相截形成空洞,且空洞范围随时间变化,使山脉边缘的格点有时会增减,故计算和编程困难2、不能做谱展开3、不能描述非静力平衡运动数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-15-第四节地图投影坐标系1、为何采用地图投影坐标系1、天气工作需要(天气图的应用)2、经纬度网格由于经线向极地汇集,使差分计算不便,甚至引起误差(如散度计算)2、常用的地图投影主要是正形投影(即投影后微元的形状不变,两条微分线段所夹的角度投影不变,投影前描写方程的正交曲线坐标投影后不变)如果原来地球上的距离为,eedxdy,投影到地图上的距离为,dxdy的话,则有:eedxmdxdymdy==m称为地图投影坐标放大系数(或地图放大系数)①极射赤面投影(极射赤平投影)示意图见P27fig1.
3ccoscoscossinxmaymaλλ==01sin1sinm+=+常用于北半球图,南半球图②Lambert投影(等角圆锥投影)示意图见P27fig1.
3b常用于中纬度1tan()cos221tan()sin22xAllyAllπλπλ==001tan()tansin222()cossintan2llAmaπθθθθ==第第一一章章大大气气运运动动方方程程组组-16-01212909045θ==或,可以等于,当1230,60==时,0.
7156l=单位经度圆锥面所张的平面角(弧度)00tan2lDAθ=,0D为标准纬度012(或)在图面上的径向距离(即极点到30或60的距离)11423.
5Akm③麦卡托Mercator投影(一种圆柱投影)示意图见P27fig1.
3a,fig1.
5常用于低纬天气如季风、台风等的研究地图面为圆柱面与地球赤道相切,光源放在地心.
地图面展开后为一矩形,经纬度成矩形分布1sinlncos1cosxayamλ=+==放大倍速3、地图投影坐标系中的大气运动方程组⑴梯度、散度和涡度的表达式3mimjkxyzφφφφ=++23()()uvwVmxmymz=++2311ijkmmVmxyzuvwmm*=2222322()mxyzφφφ=++其中数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-17-11,,dxdydzuvwmdtmdtdt===⑵方程组22211[()]11[()]1()0xyzTpdumdpvfmvuFdtdxxmymdvmdpufmvuFdtdyxmymdwdpgFdtdzuvwmtxmymzdTRTddQFdtCdtdtdmumvwdttxyzρρρρρρρρρ=++++==+=+++其中由于极射赤面投影投影和麦卡托投影的m比较简单,可代入运动方程的[]中,得到:极射赤面投影的运动方程:20201[]1sin1[]1sinxydumdpxvyuvfFdtdxadvmdpxvyuufFdtdyaρρ++麦克托投影的第一、二运动方程:(2)sincos(2)sincosxydumdpuvFdtdxadvmdpuuFdtdyaρρ=+Ω++=Ω++xyFF、等用地图投影坐标系中的表达式数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-19-第二章大气运动的整体性质与能量约束设计数值模式,出于描述物理过程的需要和实际计算的可能,难免要对完全的大气动力学和热力学方程进行某些简化,并在求解时采取某些近似方法.
如此做法可能使得模式大气与实际大气在物理性质上不一致,使大气原有的能量关系遭到破坏,造成虚假的能量的源和汇,从而给计算结果带来严重误差.
研究大气的整体性质和能量关系,可以更好地认识模式大气的特点,考察它是否合理地描述了大气运动的性质.
一个与大气原有的物理性质保持协调一致的数值模式,有利于减少误差和增加计算的稳定性.
第一节大尺度大气运动的特性1、准地转特性地转关系:==yPfuxPfvggρρ11(Z坐标系)或==yzgfuxzgfvgg(P坐标系,z为等压面高度)在自由大气中,yxFF,等可略,局地直角坐标系下的运动方程可写成:+====ggfufuyPfudtdvfvfvxPfvdtduρρ11说明dtdu和dtdv是由于风场u、v和气压场不满足地转平衡而造成的.
fv与xPρ1整体属性能量关系大气数值模拟设计约束第第二二章章大大气气运运动动的的整整体体性性质质与与能能量量约约束束-20-(或fu与yPρ1)之差值称为地转偏差.
在大尺度运动中地转偏差较小,一般dtdu比fv和xPρ1小一个量级.
因此说大尺度大气运动是准地转平衡的.
而且大气有自己调节的能力:当出现地转偏差时,气压场和风场会进行调整,而达到一个新的平衡状态.
这就是地转适应过程.

准地转性的应用大尺度运动的准地转性对我们讨论问题提供了很大的方便,在不少情况下,我们可以用地转风代替实际风(如地转涡度,Zfgyuxvggg2==ζ代替实际涡度;初值处理时,用地转风近似求初始风场).
但要注意,不能用地转风散度代替实际风散度.
因为实际散度在天气变化中非常重要,但地转风散度几乎为0,只有地球曲率项:ctgavafvyfvyfxzfgfyxzgyxzfgyxzfgyvxuggggg=Ω=Ω==++=+1cos2cos21222在中纬地区,地转风散度约为1.
6*10-7/s准地转:很接近地转,但又不完全是地转.
有些量可以用地转公式来算,而有些量(如散度)就不能,需通过其他途径来计算.
2、准静力平衡从第三运动方程,在不考虑摩擦力对垂直运动影响时:01=++zPgdtdwρ如果把气压场看成有一个静态的气压P和一个扰动的气压P′,即PPP′+=.
则数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-21-有:011=′+++zPzPgdtdwρρ根据实测资料分析,第2、3项最大,且近于平衡,即:01=+zPgρ(在大气静止时,垂直气压梯度力与向下的重力平衡).
也就是说zPdtdw′ρ1~(气压场的扰动引起垂直运动.
对于大尺度,这个量很小,因此把它略去).
于是,gzPρ=.
这就是静力平衡方程.
1.
只有在静力平衡情况下,才可以采用P坐标和σ坐标;2.
做静力平衡近似,并不等于认为没有垂直运动或垂直速度不变.
只是它在第三运动方程中不出现,w可以通过连续方程或热通量方程把它求出来.
所以叫着准静力平衡.
3、准水平特征大尺度运动中空气速度的垂直分量与水平分量之比非常小,空气运动差不多是平行于地面的.
例如:一个典型的水平尺度L~103km的天气系统,其垂直尺度H通常受到对流层厚度D的限制,即H≤D,而D~10km量级.
空气块穿过这个典型天气系统,其典型水平速度和垂直速度分别为U和W,气块穿过水平距离L所需时间为:T≈L/U.
同时气块还移动了垂直距离H,所以T≈H/W≤D/W.
于是有L/U≤D/W,则W/U≤D/L~10-2.
W/U是空气运动轨迹线的斜率,可见空气运动基本上是水平的.
空气的垂直运动虽小,但它很重要.
只有W不等于0,才会发生各种重要的天气现象.
因此,我们说大尺度大气运动接近水平,但又不完全水平,所以说是准水平,它与静力平衡是一致的.
第第二二章章大大气气运运动动的的整整体体性性质质与与能能量量约约束束-22-4、准无辐散大尺度运动的水平散度很小,约为10-6/s,比涡度小一个量级以上.
对水平运动方程进行散度运算可得散度方程:0)()(2)(22=Φ+++++yuxvfuyfxvyuyvxupvypuxDdtdDωω式中:VD=2——水平散度在大尺度运动中,散度方程中前三项比后三项小两个量级以上,所以基本不可能从散度方程来求散度及其变化.
若在方程中只保留最大的项和比它们小一个量级的项,可得到平衡方程:0)()(22=Φ+yuxvfuyfxvyuyvxu平衡方程中再也不出现带有散度和垂直速度的项,因此平衡方程表明:大尺度运动不但是准水平的,而且是准无辐散的.
用平衡方程描述的运动称为准水平无辐散运动.
平衡方程是散度方程的一级近似,它的精度比地转关系高一个量级,因此用平衡方程描述的准水平无辐散运动比准地转运动更接近实际大气的运动.

特别是在低纬地区,地转关系不太好,平衡方程更为常用;在数值预报中有时用平衡方程来进行初值处理.
第二节小尺度运动的特性(包率内斯克Boussinesq近似)1、特性准地转、准水平、准无辐散、准静力平衡都不满足2、Boussinesq近似对于小尺度运动,特别是热对流运动(如积云问题),气块受热产生密度扰动.
此密度扰动即使微小,对垂直运动也会有大的影响,因为它可以产生较大的浮力,因此很重要.
表示方法:设′+=′+=ρρρ)()(zPzPP数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-23-(1)水平运动方程中可以略去密度扰动的影响:+=++=yxFyPfudtdvFxPfvdtduρρ11(2)在垂直方向上,P和ρ满足静力平衡:gzPρ=.
但dtdw不可略,第三运动方程变为:gzPdtdwρρ′′=~dtdw是由于扰动P′和ρ′引起(3)连续方程不可简化,保留可压的形式.
这种近似在研究对流问题(如积云、雷暴、龙卷、局地暴雨等问题)时常用.

第三节大气中的几个守恒量大气中有几个守恒量,在设计模式和进行数值模拟过程中需要注意保持它们的守恒.
1、总质量守恒由连续方程可推得0=+∫∫∞∞sszzVdzdztρρ(具体推导见p.
35~p.
36),其垂直边界条件为:∞→→==zwzzzVwssss0ρ对于全球大气,速度矢量在边界面外法线方向上分量为0,于是对上式全球大气积分,有:0=∫∫∫∞szAdAdztρ这表明,全球大气的总质量是守恒的.
这与大气下边界地形的起伏状态无关.
2、绝热大气位温与位温平方守恒有绝热大气的热力学方程(0=dtdθ)和连续方程可得:第第二二章章大大气气运运动动的的整整体体性性质质与与能能量量约约束束-24-0)()(33=+Vtρθρθ对上式进行全球大气积分,利用垂直边界条件,可得:0=∫∫∫∞szAdAdztρθ这表明,在绝热条件下全球大气的位温是守恒的.
同样,通过02=dtdθ可以推得:02=∫∫∫∞szAdAdztρθ,以此类推.
对任意幂次的以θ为函数的代数多项式,其全球积分也是守恒的.
3、绝热无摩擦大气总能量守恒通过运动方程、连续方程和热力学方程可得总能量(动能、位能和内能的总和)方程的通量形式(具体推导见p.
36~p.
38):QFpVVVTcKTcKtvvρρρρ+=+Φ+++Φ++3333333)()()(对上式进行全球大气积分,利用垂直边界条件,并设T、3V有界,可得:∫∫∫∫∫∫∞∞+=Φ++sszAzAvdAdzQFVdAdzTcKt)()(33ρρ上式说明,全球大气的动能、内能和位能之和的变化是由于非绝热加热过程和摩擦消耗所造成的,在绝热无摩擦条件下,全球大气总能量是守恒的.
4、绝热无摩擦大气位涡和位涡拟能守恒可推得0)(=+∫∫∫∞szAdAdzzftθζ(具体推导见p.
38)0])[(2=+∫∫∫∞szAdAdzzftθζ上式表明,在绝热无摩擦条件下大气的位涡和位涡拟能(位势涡度的平方)是守恒的.
数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-25-第三章大气中的波动大气动力学方程组描写了大气中的各类波动,包含了多种尺度运动,其物理因子很复杂.
主要有:声波、重力波、惯性波、涡旋慢波(长波、罗斯贝波)和超长波.
其中,重力波又分为重力外波和重力内波,重力波与惯性波相结合,成为重力惯性波.
声波——由大气的可压缩性产生;重力波——由于地球重力场产生;惯性波——地球旋转效应产生;Rossby波——考虑到地球的球面效应,在科氏参数的变化条件下产生;声波、重力外波都是快波,它们对天气变化特别是大尺度天气变化作用不大(重力内波与雷暴、龙卷有关),但对数值计算有影响,容易引起计算不稳定,所以常常要把它们滤掉.
第一节大气中的混合波解1、大气波动基本方程采用小扰动法,对绝热无摩擦的大气动力学基本方程组(p.
51)进行线性化,并假设基本流场为静止的(p.
523.
1),可推得描写大气波动的基本方程(具体推导见p.
52~p.
54):01222212212222222221=++++++VtczcggNztftxtNtssλλλβλ其中0RTcsκ=为绝热声速,vpcc=κ,zgN=02lnθ,N为浮力振荡频率,即Brunt-Vaisala频率.
1λ和2λ为示踪系数,可取1或者0.
01=λ表示取静力平衡近似,02=λ表示取滞弹性近似.
不同尺度的大气运动各类波动联系第第三三章章大大气气中中的的波波动动-26-2、波解分析设波动解代入大气波动基本方程得到ω的5次方程,即频率方程:022422132521=++RsBRsANωωωωωωωλωωωλλ(其中各量参见p.
54中的3.
23).
分别可得到对于高频波的解2,1ω(式3.
26声波频率)、解4,3ω(式3.
27惯性重力内波频率),以及对于低频波的解5ω(式3.
28Rossby波频率).
因此,一般大气系统的方程组中包含声波、惯性重力波和Rossby波.
第二节各种波动特征1、声波声波是由于大气的可压缩性引起的.
当空气受声源扰动压缩时,会引起密度和气压扰动的局地变化,并使四周空气产生膨胀压缩的交替变化,这种过程的传播即形成声波.
声波是一种纵波.
滤波:(1)假定大气不可压(或三维无辐散,02=λ);(2)假定静力平衡(01=λ),滤垂直方向声波(垂直方向两力平衡,没有外力引起空气上下压缩)2、重力内波重力波是在重力作用下引起的波动,它是一种横波.
在大气内部密度不连续面或风不连续面(如锋面)形成,或在稳定层结下气块受扰动在垂直方向上振荡形成的波动,都与重力有关.
3、惯性波大气中的扰动在科氏力作用下会产生惯性,这种振荡在空间的传播上就形成惯性波,它是一种横波.
数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-27-4、惯性重力内波在地球大气中,地球旋转和大气层结的共同作用,使得惯性波与重力内波混合在一起,形成惯性重力内波.
5、Rossby波从天气图上可以看到,在北半球中高纬上空有绕极的西风带,西风带中约有3~5个波,这些波动具有重要的天气和气候意义.
Rossby首先从理论上解释了这种波动的成因,指出决定此波产生的重要因子是地转参数f随纬度的变化.
这种波动被称为Rossby波,又称为大气长波或行星波.
数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-29-第四章数值计算法差分法(FiniteDifferenceMethod)数值天气预报和数值模拟问题主要解偏微分方程的问题.
它的解和它所要求的初值,都不是自变量的连续函数,直接用它们做预报会遇到困难.
原因是:1)观测站的分布是离散的,不可能给出连续的初值.
2)方程组为复杂的非线性方程组,目前还未有一个普适的解析求解方法.
因此,为便于求解方程,首先要将方程离散化.
离散化方法有三:(1)差分法;(2)谱方法;(3)有限元法,即把偏微分方程变成相应的泛函极小问题来求解.
目前,差分法最成熟,用的最多.
谱方法也已广泛应用.
有限元法用得较少.

差分法有如下优点:1)用差分代替微分,比较直观2)在球面或地图面上取网格进行差分运算,比较容易3)差分法既可用于全球问题,也可用于局地问题.
对垂直坐标系没有特别的要求.
第一节水平空间差分网格及其差分计算格式(TheGridPointMethod&FiniteDifferenceSchemes)在水平面上,通常把计算区域取成矩形或扇形(按地球经纬度).
另外还有三角形或菱形,但编程复杂,少用.
第第四四章章数数值值计计算算法法-30-如图,令NjyjyMixix,.
.
.
2,1,.
.
.
2,1=Δ==Δ=划分成网格一般如取矩形计算区,则取dyx=Δ=Δ-格距方形网格如取经纬度计算区,则cos,=Δ=Δdxdy经纬度网格由此,每个位置yx,都可用ji,来表示.
ji,为整数的点称为网格点,计算区内的点又叫内点,边界上的点叫边界点.
函数),(yxf在ji,点上的值用jif,来表示.
(f代表.
.
.
,,,,Tpwvu)水平空间差分有三种表示方法(FiniteDifferenceSchemes):X方向为:jixjijijifdffxxf,,,1,1)(1++→Δ≈向前差(ForwardDifference)jixjijijifdffxxf,,1,,1)(1→Δ≈向后差(BackwardDifference)jixjijijifdffxxf,,1,1,21)(21→Δ≈+中央差(CenteredDifference)同理,在Y方向为:jiyjijijifdffyyf,,1,,1)(1++→Δ≈向前差(ForwardDifference)jiyjijijifdffyyf,1,,,1)(1→Δ≈向后差(BackwardDifference)jiyjijijifdffyyf,1,1,,21)(21→Δ≈+中央差(CenteredDifference)用差商来代替微商必然会带来一定的误差,称为截断误差(TruncationError),它反应了差商的精度.
之所以会产生截断误差,是由于我们在构造差分格式时,从函数的Taylor展开式中所取的次数有限.
通常在所研究的气象问题中,实际使用的网格距xΔ与方程所描述运动的特征尺度相比往往是一个很小的量,因此,对于函数),(txf可展开成Taylor级数:数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-31-)1(.
.
.
!
3)(!
2)(),(),(333222+Δ+Δ+Δ+=Δ+xxfxxfxxftxftxxf)2(.
.
.
!
3)(!
2)(),(),(333222+ΔΔ+Δ=Δxxfxxfxxftxftxxf由(1)、(2)两式可得)(.
.
.
.
!
2),(),(122xOxffxxfxtxftxxfxfiiΔ+Δ=ΔΔΔ+=+ττ)(.
.
.
.
!
2),(),(122xOxffxxfxtxxftxfxfiiΔ+Δ=Δ+ΔΔ=ττ这就是向前差和向后差格式,若R为它们的截断误差,由上两式可见,其截断误差为xΔ的量级.
将(1)和(2)相减,可得:))((2.
.
.
.
!
3)(2),(),(211233xOxffxxfxtxxftxxfxfiiΔ+Δ=Δ+ΔΔΔ+=+ττ这就是中央差格式,其截断误差R的量级为2)(xΔ将(1)、(2)相加,可得22xf的中央差格式为:))(()(2221122xOxfffxfiiiΔ+Δ+=+τττ其截断误差也是2)(xΔ的量级所以中央差比向前差和向后差精确些,因此一般水平差分通常采用中央差格式.

第二节时间差分(FiniteTimeDifference)对于tf也用差商来表示,如取三个时间层1,,1+τττ,则有:第第四四章章数数值值计计算算法法-32-)(11τττffttfΔ≈+时间向前差)(11Δ≈τττffttf时间后差)(2111+Δ≈τττffttf时间中央差这三种表示法截断误差也不一样,分别为)(tOΔ,)(tOΔ和))((2tOΔ,因此中央差的精度高些.
第三节各种常用的数值时间积分格式对于大气动力学方程组,可写成:fvxpzuwyuvxuutu+++=ρ1)(.
.
.
=tv.
.
.
=tT左边随时间变化,右边随空间变化,对于已知的空间分布,可以先用空间差分把等式右边的项算出来,然后等式左边写成上面的时间差分形式,就可计算另一时刻的要素值.
如上面u的方程,如果等式右边算出的值记做),,,(tzyxF,则第一运动方程可简写为:),,,(tzyxFtu=计算tu有三种差分方法,计算F也有三种空间差分方法,把它们分别组合在一起,有九种解方程的方法,但因时间后差对我们做预报没什么意义,所以一般不用,于是就有以下6种:时间向前差,空间向前差*时间向前差,空间向后差*时间向前差,空间中央差时间中央差,空间向前差时间中央差,空间向后差*时间中央差,空间中央差其中又以*的三个用得较多,用时间差分格式和空间差分格式分别组成的解偏微数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-33-分方程的格式,叫做时间积分格式.
时间积分格式又分显式和隐式两种,显式格式是空间差分用当时时刻的资料(要素值)来算,而隐式格式的空间差分部分用下一时刻(预报时刻)的要素值来算.
1、取两个时间层的情况(Two-Level-Schemes):1,+ττ对于时间步长为tΔ,τ和1+τ记积分步数,则∫Δ+Δ+++Δ+=+=ttFFtudttzyxFuu)1(11)(),,,(τττττττβα这里1=+βαA)向前差格式:(显式)——用当前时刻值算F,再算下一时刻时间值.
即取0,1==βα,则tFuuΔ+=+τττ1,截断误差为:)(tOΔ向前差格式又可分两种,一种的τF用向后差,另一种τF用中央差,分别叫做迎风格式和欧拉(前差)格式.
如对前面的运动方程,其差分格式分别为:1、迎风格式:(时间前差,空间后差).
.
.
)()(1,,,,1,,,1,ΔΔ=+ττττττττjijijijijijijijiuudtvuudtuuu这表明u的变化只受到其上游u值的影响,所以叫迎风格式.
2、欧拉向前差格式:(时间前差,空间中央差,但用后一时刻值算).
.
.
)(2)(21,1,,,1,1,,1,ΔΔ=+++ττττττττjijijijijijijijiuudtvuudtuuu第第四四章章数数值值计计算算法法-34-B)向后差格式:(隐式)时间向前差,空间中央差,但用后一时刻值算取1,0==βα,即tFuuΔ+=++11τττ截断误差为)(tOΔ由于在τ时刻,1+τF步还不能直接算,因为111,,+++τττpvu等值还没有,只能先把差分式子写出,如:.
.
.
)(21,11,11,,1,Δ=+++++τττττjijijijijiuudtuuu然后把1+τu移到等式左边,用迭代法来算.
因此麻烦,但以后要讲,它比向前差格式稳定,所以有时选用它.
C)梯形格式:(隐式)时间向前差,空间中央差(但用当时值和下一步值各1/2来算)取21==βα,即)(211+++Δ+=ττττFFtuu,截断误差:))((2tOΔ此格式计算比较麻烦,但精确度较高D)欧拉后差格式:(Matsuno)格式是两步显式格式考虑到隐式向后差格式稳定性好,但计算不方便,故设计此格式.
第一步:先用向前差格式求出一个带"*"的量.
.
.
)(2)(21,1,,,1,1,,*,ΔΔ=++τττττττjijijijijijijijiuudtvuudtuuu第二步:将带"*"的量用于求空间差分,求出下一步1+τu.
.
.
)(2)(2*1,*1,*,*,1*,1*,,1,ΔΔ=+++jijijijijijijijiuudtvuudtuuuττ写得简单点:ττtFuuΔ+=*)]1(,[**+Δ=τtuFF数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-35-*1tFuuΔ+=+ττ这是个常用格式.
E)欧拉——梯形格式(又叫Heun's格式)ττtFuuΔ+=*)(2*1FFtuu+Δ+=+τττ2、取三个时间层的格式(Three-Level-Schemes):取1,,1+τττ步为例,对于),,,(tzyxFtu=的积分:)(),,,(111)1()1(11+Δ+Δ+++Δ+=+=∫ττττττττγβαFFFtudttzyxFuutt其中2=++γβα,截断误差为))((2tOΔ有以下三种:A)中央差格式(又叫蛙跳格式),时间中央差,空间中央差.
是个显式格式.
取2,0,0===βγα,用当前时刻要素值算FτττtFuuΔ+=+211例,第1运动方程:.
.
.
)()(1,1,,,1,1,1,1,ΔΔ=+++ττττττττjijijijijijijijiuudtvuudtuuu第第四四章章数数值值计计算算法法-36-B)Adams后差格式(显式)0,23,21===γβα)3(2111++Δ+=ττττFFtuuC)Milue-Simpson格式(隐式))2(21111++++Δ+=τττττFFFtuu其中蛙跳格式用得最多.
第四节差分近似的相容性、精确性、收敛性和稳定性当我们把大气动力学方程组中的微分运算用差分近似来代替时,需要考虑差分近似的相容性、精确性、收敛性和稳定性:1、差分近似的相容性差分近似的相容性即是差分系统在某种程度上是否近似于微分系统,或者说它们是否相协调.
若是相容的话,则在小的时间步长和空间步长趋于零的极限条件下,差分系统应等同于微分系统.
在上述讨论中,若取时间向前差,空间向前或向后差格式,则其截断误差的量级为:)()(xOtORΔ+Δ=在时间向前差、空间中央差的格式中截断误差为:)()(2xOtORΔ+Δ=在蛙跳格式中:)()(22xOtORΔ+Δ=当tΔ和xΔ都→0时,这几个差分方程的截断误差R都趋近于零,所以对于我们讨论的微分方程而言,这几种差分格式都是相容的.
数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-37-2、差分近似的精确性差分近似的精确性问题,即是讨论差分系统相对于微分系统的精确程度.
差分近似的精确性取决于差分近似时的截断误差和计算时的舍入误差.
为了减少截断误差,需要在计算条件允许的情况下,尽可能地减小tΔ和xΔ、yΔ的大小,以使截断误差尽可能地小.
对于舍入误差,现在计算机一般不严重,但要注意有时算术运算会出一些问题,编程时要注意.
(如很小的数与很大的数相除后再来一个大数,最好先乘后除).
虽然计算机的舍入误差不大,但因为差分方程求解是一步一步进行.
如果误差不断的积累,也会增长很厉害,而造成计算不稳定.

3、差分近似的收敛性考虑微分方程初值问题的真解),(txifΔΔτ和与之相容的差分方程的解τif之差),(txiffiΔΔττ在求解区域中的极值)),((,txiffMaxiiΔΔτττ,当tΔ和xΔ→0时,它也趋近于0,则称这个差分格式是收敛的,或者称差分方程的解收敛于微分方程的解.
差分近似的相容性不能保证其收敛性.
下面用一维平流方程来讨论此问题.
一维平流方程为:0=+xuctu其通解是:)(),(ctxtxu=右图中,常数=ctx的直线是这个方程的特征线(即线上的点满足方程).
若用时间向前差,空间中央差、向前差或向后差,三种差分格式写出一维平流方程的差分形式分别为:02111=Δ+Δ++xuuctuuiiiiττττ(A)011=Δ+Δ++xuuctuuiiiiττττ(B)常数=ctxxt第第四四章章数数值值计计算算法法-38-011=Δ+Δ+xuuctuuiiiiττττ(C)由差分方程(A)可知,图中τ,i点的u值τiu依赖于11111,,+τττiiiuuu,而这三点的值又由222122122,,,,++τττττiiiiiuuuuu的值算出,……,如此倒推,τiu的值决定于三角形内所有交点的值,这三角形区域称为数值解的依存域.
τiu的值与这个依存域外的点无关系,但是一维平流方程的真解在τ,i的值即),(txiuΔΔτ只与初值P点有关.
P点是常数=ctx的特征线与x轴的交点.
如果P点落在三角形依存域外,则无论tΔ和xΔ取得如何小,τiu都不可能收敛于真解.
特征线的斜率为c1,只有xtcΔΔ>1时,特征值才会在数值解的依存域内,所以差分近似的收敛性的必要条件是xtcxtcΔ≤Δ≤ΔΔ或1同样可证明差分格式(B)和(C)的收敛性必要条件为:1≤ΔΔxtc上述收敛性条件称为柯朗条件.
4、各种时间积分格式的稳定性为了简便起见,用一维线性平流方程来讨论各时间积分格式的线性计算稳定性.
一维线性平流方程为0=+xuctu,c为常数(实数)设)(ktxieu+=α,α是波数,k是圆频率.
代入方程得:ckα=(实数),所以此方程解是稳定的.
c就是波速.
常数=ctxxtPτ,iτ11+iii12112+++τiiiii数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-39-(一)几个显式格式的稳定性1、欧拉向前差格式(即时间向前差,空间中央差,F用当前时刻值计算)一维线性平流方程的差分近似为:0)(2)(1111=Δ+Δ++ττττnnnnuuxcuut注:因后面要用i表示虚数,此处用n表示空间位置令上式基本解为:)(xntkinAeuΔ+Δ=αττ这里α是波数(相当于Lπ2),k是圆频率,irikkk+=.
当0≠ik时,有:)(xntkitknrieAeuΔ+ΔΔ=ατττ除了是一个波动解外,τnu还以e指数变化,增长率为tkieΔτ.
如果0xtcc且时才稳定,数值预报中c不一定都大于0,所以通常不用这格式.
而是在预报员计算平流项时(如算T平流,水汽平流等),才用这种格式,且把计算点放在风的下游,以满足0>c,因而叫做迎风格式.
3、欧拉后差(松野)格式(两步显式格式,或叫迭代格式)一维平流方程用欧拉后差格式写出为:第一步:0)(2)(1111*=Δ+Δ++ττττnnnnuuxcuut第二步:0)(2)(11*11*11=Δ+Δ+++ττττnnnnuuxcuut两式中把1*+τ的项消去,则有:)2(4)(2)(12222111τττττττ++++ΔΔ=Δ+Δnnnnnnnuuuxtcuuxcuut也是相当于在欧拉向前差格式中人为地加了一个平滑项,只是它不是用1+n点和1n点来平滑,而是用2+n点和2n点来平滑,它对高频波衰减很有意义,故常用.
用上面地类似方法可以证明它也是条件稳定的,其稳定条件也为1≤ΔΔxtc.
4、蛙跃格式(时间和空间皆用中央差)用这种格式,一维平流方程变成:0)(2)(211111=Δ+Δ++ττττnnnnuuxcuut同样令其基本解为:)(xntkinAeuΔ+Δ=αττ代入得:0)(2)(21)1()1()1()1(=Δ+ΔΔΔ+ΔΔΔΔ+xnixnitikxnitiktikeeAexceAeAetααταττ第第四四章章数数值值计计算算法法-42-0)()(=ΔΔ+ΔΔΔΔΔΔΔΔxixixnitikxnitiktiktikeeeAextceeeAeααατατ因为tiktiktiktikeeeeΔΔΔΔ=12和xixixxixeexixiΔ=Δ+ΔΔ+Δ=ΔΔαααααααsin2sincossincos所以,0)sin2(12=ΔΔΔ+ΔΔxiextcetiktkiα或写成:01)sin2()(2=ΔΔΔ+ΔΔtiktikexxtcieα设λ=ΔΔxtc,得:xxietikΔ±Δ=Δαλαλ22sin1sin由此,差分方程有两个解,而微分方程只有一个解,哪个真解,哪个假解,如何处理设irikkk+=,xxieetktikirΔ±Δ=ΔΔαλαλ22sin1sin或:xxietkitktkrriΔ±Δ=Δ+ΔΔαλαλ22sin1sin)sin(cos下面分两种情况来讨论:(1)xΔ),取kmx500~300:Δ.
对于次天气尺度(km1000≤),取kmx200~100:Δ.
对于小尺度,xΔ取得更小.
(2)提高格式的精度,多用中央差格式.
2、差分方程的计算解(虚假的计算波)在讨论蛙跃格式的稳定性中,我们看到一维线形平流方程用蛙跃格式时,其差分近式可得到(A)式,即:xxietikΔ±Δ=Δαλαλ22sin1sin,其中xtcΔΔ=λ也就是说,差分方程有两个解,而微分方程只有一个解.
为讨论方便,设:tikeBΔ=,xΔ=Ωαλsin,上式变成:21Ω±Ω=iB讨论稳定的情况,即1≤ΔΔxtc或211Ω+Ω=iB221ΩΩ=iB第第四四章章数数值值计计算算法法-46-考虑到)sin(arcsin)sin(sinsin1Ω=Ω=Δ=Ωxαλ)cos(sin)(sinsin111122Ω=Ω=ΩβieiB=ΩΩ=)sin(sin)cos(sin111)(11112)sinsin()sincos()cos(sin)sin(sinβπππ+=Ω++Ω+=ΩΩ=ieiiB其中)sinarcsin(sin1xΔ=Ω=αλβ由tikeBΔ=11,可写成β=Δtk1tikeBΔ=22βπ+=Δtk2于是,差分方程解为:xniiineEeMeuΔ++=ατβπβττ)()(其中M和E为待定常数,但它们有一定关系,可由初值来定.
在初始时刻0=τ时,xnixnineEMAeuΔΔ+==αα)(0,所以AEM=+又因1=πie,所以有:)()()1()(βτατβτατ+ΔΔ+=xnixninEeeEAu即差分方程的解包含两个波,第一个波振幅为)(EA,相速为tΔαβ.
当1→ΔΔxtc时,xxxtcΔ→ΔΔΔ=ααβ)sinarcsin(,相速ctxt→ΔΔ→Δαααβ,所以第一个波的解与微分方程的波解相同,称为物理解.
第二个波振幅为E,相速为tΔαβ,而且每个时间步长其位相改变一个π(即每一步位相相反一次),出现时间振荡,振荡周期为)2(tΔ,这不是真的解,叫做计算解,或计算波,是由于用了时间中央差,其差分方程变成了二阶差分方程,于是有了增根.
实际工作中要把这个增根去掉,采用时间平滑的办法,格式为:)2(2111τττττuuuuu++=+第七节非线性计算不稳定非线性计算不稳定出现在非线性差分方程中,虽然tΔ满足线性方程计算稳定性要求,但还可以出现非线性计算不稳定.
数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-47-形成原因:由于波与波的非线性作用引起的波的混淆现象造成.
简例:一个波11()ixteακ,另一波22()ixteακ,经过非线性项xuu运算后,会出来一个新波:1212ixteαακκ++新波的波数2121ααααα或>+=,2Lπα=,因此新波的波长小于1波的波长,也小于2波的波长.
但由于网格只能分辨波长为2d的波,短波无法分辨,而且会被误认为是可分辨的波.
如右图,细实线的波被误认为是细虚线的波,于是能量加到粗实线的波里,积分步数越多,加得越多,造成计算不稳,虚假能量增长.
第八节平滑和滤波为了消除由于中央差引起的计算波和非线性计算不稳定产生的虚假的能量增长,通常用平滑的办法来滤波.
由于造成的误差多是波长较短的扰动,天气意义不大,但影响运算,所以把它们滤掉.
平滑方法:1、5点平滑公式,,1,1,,1,1,(4)4ijijijijijijijsFffffff++s为平滑级数,如s=1,则用四周点的平均值代替中间点的值;s=0.
5,四周点平均值取1/2,中间点取1/2.
第第四四章章数数值值计计算算法法-48-2、9点平滑公式,,1,11,1,1,1,11,1!
,!
,1,(8)8ijijijijijijijijijijijsfffffffffff一般认为9点平滑太强,多用5点平滑.
对模式积分一段时间,平滑一次.
第九节水平边界条件在做有限区域的数值模拟时,要有侧边界条件.
例如,用空间中央差做一时间步长后,除边界外,内部点的量都变了,但边界点未计算到,所以在算下一步时要相应给出侧边界的值.
对于原始方程模式而言,侧边界的给法如下:1、固定边界和刚体边界若模拟计算区是矩形,边界和yx,轴平行,则有下列几种:1.
在DyLx,0;,0==的边界上,0uvzttt===即边界值固定不变.
2.
在常数=y的边界上000uyvzy===在常数=x的边界上数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-49-000vxuzx===3.
在常数=y的边界上000uuzugtxxvzfugy++==+=在常数=x的边界上000vvzvgtyyuzfvgx++==+=2.
和3.
为刚体边界条件,优点:简单方便;缺点:易不稳定(举例:5层模式).
2、海绵边界(可比较稳定)边界上取固定边界,边界向内两三圈设过渡带,过渡带内边界用模拟值,在二者之间逐渐过渡,如:(1)1(1)AaAaAτττ++=+对于边界点"1",a=1;"4"点,a=0;第"2"和第"3"点a可以线性变化,如第"2"点a=2/3,第"3"点a=1/3.
3、外推边界(1)111jjuuττ++=(2)1()()jjuuttττ=1234第第四四章章数数值值计计算算法法-50-(3)2jjuuττ=(4)122jjjuuuτττ=(5)12()2()()jjjuuutttτττ=(6)11()()jjjjuCuutdτττ=(7)121()(34)2jjjjjuCuuutdττττ=+(8)11111122(1)jjjjjCtCtuuuddτττττ+ΔΔ=+这8种和固定值边界做过试验.
固定值最差,做约40步以后偏离真解很多,且振荡,第4—8种最好,第1种也不错,第2种稍差,第3种不好,后来发散很厉害.
4、网格嵌套由粗网格模式向细网格模式提供边值,主要问题:边界值要处理得好.
方法:1.
用内插法:把粗网格上的值插到细网格上,方法简单,但会产生边界值与内部值不协调.
2.
用"1.
"法后,在边界加过渡区,类似海绵边界处理.
这个方法不仅可以处理边界问题,在污染物输送模式中可用此法为计算区内提供,,,wvu降水等场.
第十节垂直方向的差分1、垂直分层设计模式时,为了能考虑到大气在垂直方向上不均匀的特性(与大气斜压性有关).
要把大气分成若干层,以便计算z或p项,通常分层办法为:先将大数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-51-气分为K层,在每层中加一个计算面,所以有2K+1个面,包括大气层顶和地面.
如下图:2、垂直差分将,,,,uvTqφ等变量放在奇数面上,,wω等放在偶数面上如第1层与第3层之间垂直差分,以P坐标为例:(1)连续方程中ppΔ02~ωωω(2)个别变化项中的puω等一般把运动方程与连续方程一起,将运动方程写成通量形式,所以此项写成)(1~)(22443uuppuωωωΔ因在第2、4层无u,所以4u、2u(或4v、2v等)用5u、3u和3u、1u插值来求.
插值方法有各种,以后在适当章节中讲.
01234地00=ω00=wρTPP=或11111,,,,φqTvu2ω2w33333,,,,φqTvu4ω4wsωsw数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-53-第五章正压原始方程模式正压大气和斜压大气正压大气是指大气的密度ρ的分布只依赖于气压P的大气即ρ≡ρ(P).
这是一种理想化的大气,比较简单,但对于研究大气中某些物理过程是很有用的.

在这种理想化的大气模式中,等压面,等密度面和等温面都是重合的.
在等压面上没有等温线,水平温度梯度为零.
因此,地转风随高度的变化为零即不存在热成风.
在这种模式中,等压面之间厚度处处相等.
因而往往可以用某层等压面的运动状态代表整层大气运动状态(常有500百帕等压面),这种模式称为正压模式.
斜压大气则是指密度ρ的分布依赖于气压和温度的大气,即ρ=ρ(T,P),实际大气都是斜压大气.
在斜压大气中,等压面与等密度面和等温面彼此相交,形成许多"力管",等压面上等温线越密集表示大气的斜压性越强.
在这种大气中,地转风随高度有变化,即存在热成风.
在天气图上凡是锋区,急流区等温线密集的区域,都表明斜压性较强.
第一节历史回顾挪威锋面气旋学派创始人V.
Bjerkness首先提出数值天气预报的思想;英国数学家Richardson实践了V.
Bjerkness的思想,试验失败,主要因为缺少可用的初始资料(尤其是探空资料),方案存在严重问题;无人再试图进行数值预报,大气动力学理论得到了发展;基于高空探测资料的分析,Rossby提出了著名的长波理论,认为大尺度环流的演变主要由大气长波所控制,而声波、重力波的作用并不重要;Charney发表了著名的"大气运动的尺度"论文,利用尺度分析的方法,证明了把静力平衡、地转平衡的假设系统地引进动力方程组,便能"滤去"声波和重力波而使方程组简化;Charney等利用滤声波和重力波的正压无辐散模式成功预报;瑞典军事气象局使用正压模式第一次实现了业务数值预报过滤模式转向原始方程模式,经过30多年的发展,数值预报得到了长足的发展实际的复杂大气模式大气热力状态和运动状态理想化、简单化大气模式或预报模式控制模式大气运动的方程组20世纪初1916~1918之后20多年20世纪30年代1948年1950年1954年.
.
.
.
.
.
第第五五章章正正压压原原始始方方程程模模式式-54-第二节过滤模式结合尺度分析,对涡度方程和散度方程进行简化,可以得到滤去声波、重力波的过滤模式——平衡模式和准地转模式.
准地转正压模式方程为正压准地转无辐散涡度方程:0)()(=++=+ΨfVtfdtdζζζ其中Ψ*=ΨkV旋转风f1=Ψ流函数Ψ=2ζ模式中不包含垂直运动,不包含上下层之间的相互作用和热力因子影响.
因此,准地转正压模式的预报能力是很有限的,适合用于制作500hPa等压面上的形势预报.
而且模式中只包含具有天气意义的大尺度波动——Rossby波,声波和重力波被过滤了.
这一简单模式尚能描述中纬度大尺度运动准水平、准静力平衡、准地转平衡的特征,能够反映大气运动的某些整体性质,及控制大范围环流演变的主要波动——Rossby波.
求准地转正压模式数值解应考虑的问题:(1)客观分析——把实际观测的资料插到网格点上(2)初值化——对已插到网格点上的初值进行处理,减少由于观测误差或客观分析误差导致的虚假小扰动.
譬如采用5点平滑.
(3)侧边界条件——对于有限区域预报,必须给出侧边界条件,譬如采用固定边界条件.
(4)空间网格距和时间步长——为保证计算稳定,要根据线性计算稳定性条件确定网格距和时间步长.
准地转正压过滤模式的具体计算方案和步骤参见p.
111.
因为滤去了声波和重力波的过滤模式不能描写大气中的地转适应过程,而实略去平衡方程中的非线性项非线性平衡模式线性平衡模式部分考虑了地转偏差的作用故可称为非地转过滤模式绝对涡度平流相对涡度的局地变化数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-55-际大气运动总是处于地转平衡不断破坏,又不断向地转平衡调整之中,在这个过程中往往伴随着天气剧烈变化,这正是需要预报的.
第三节正压原始方程组+++===+==*+pyvxutdtddtdpVpRTpVfkdtdVωθω0ln0原始方程模式取了准静力平衡近似,滤去了声波(Lamb波除外),保留了大气长波和重力惯性波.
原始方程模式同时描写了大气运动的准地转演变过程和地转适应,这是与过滤模式最基本的区别.
原始方程模式的物理性能比过滤模式更完善,所能描述的大气动力过程更接近实际过程.
对于具有自由表面的均质不可压缩流体,从原始方程组可推导得到浅水方程组(参见p.
113):=++=+++=++0)()(vyuxtyfuyvvxvutvxfvyuvxuutuss它可用于大气中500hPa等压面.
在给定初始条件),(),,(),,(,0000yxyxvvyxuut====和侧边界条件下即可求解.
第四节正压原始方程模式的物理属性正压模式假定流体是均匀不可压缩的,所以模式中不会存在声波、重力内波解,但在自由面上会产生重力外波,由于考虑了f的作用,出现的应是重力惯性外波.
除此以外,因为模式中还考虑了f随纬度的变化,因而模式中还会产生第第五五章章正正压压原原始始方方程程模模式式-56-Rossby波.
原始方程模式的时间步长比准地转模式小得多,因而模式的计算量大大增加.
模式大气整体性质(积分性质):性质1:质量守恒性质2:平流过程中全球大气动量守恒性质3:全球大气总能量守恒性质4:全球大气位涡和位涡拟能守恒性质5:不计地球自转作用,无辐散大气涡度拟能守恒和水平尺度守恒(平均波数守恒)模式大气的这些整体性质,不但对我们了解模式大气动力学特征有意义,对于制作数值预报也是重要的.
具体制作数值预报时,总是要将连续空间离散化,进行差分近似计算,如果差分近似计算(差分格式)破坏了这些整体性质,就相当于引入了虚假的"源"与"汇",预报结果就难以让人置信.
所以,这些整体性质就为我们设计差分格式提供了约束条件.
数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-57-第六章斜压原始方程模式第一节模式方程组及其积分关系通过σ坐标系中通量形式的模式方程组(p.
141之6.
1~6.
6)可以推导出斜压大气整体性质的一些积分关系,其中主要包括垂直区间的约束关系、质量守恒关系、保守量的守恒性、水平气压梯度力的垂直积分关系以及总能量守恒等.

1、垂直区间约束关系根据上下边界0=σ和1=σ,垂直方向上的积分满足:110=∫σd2、质量守恒关系0=sPt在静力平衡条件下,地面气压变化取决于空气柱质量的变化,地面平均气压不随时间变化表示全球大气质量守恒.
3、变量的个别变化及其通量形式对于保守量B(空气微团的动量、位温或比湿),有∫∫=AdABdt010σπ.
说明物理量B对全球大气是守恒的.
正压原始方程模式把大气视为一层流体只能描述二维运动斜压原始方程模式实际大气三维空间运动,必须考虑温度平流、垂直运动以及各种物理量在垂直方向上的变化,需要采用多层模式第第六六章章斜斜压压原原始始方方程程模模式式-58-保守量的任一函数F(B)也是保守量.
可推得∫∫=AdAFdt010σπ,说明保守量B的任一函数F(B),如果其积分存在,则对全球大气也是守恒的.
4、水平气压梯度力的垂直积分地表平坦时,水平气压梯度力的垂直积分沿水平面任意闭合曲线积分等于零.
(特殊积分条件,参见p.
144)5、总能量守恒∫∫=++AspdAdTcKt0)(10σπ表明在绝热无摩擦情况下全球大气总能量守恒.
以上讨论了连续大气的几个重要积分关系,这些都是设计差分格式时所必须满足的约束条件.
第二节变量配置与垂直差分格式1、垂直分层与变量配置见第四章第十节2、垂直差分格式的构造连续方程个别变化项气压梯度力动能产生项热力学方程静力方程数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-59-第三节斜压原始方程模式计算方案(略)数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-61-第七章初值处理对大气动力方程组求解需要初值,因此,初值处理是数值模拟和数值预报的一个重要的部分,必不可少的部分.
工作量也相当大,通常一个数值模式往往还有初值处理的配套程序.
一般初值处理包括两个部分:(1)是将不规则分布的测站资料处理成格点资料和模式层的资料,这即是客观分析.
(2)是把格点上的资料进行初值化处理,使其能适用于模式,不因初值的误差而影响计算过程和计算结果.
第一节客观分析1、水平插值如有水平面(或等压面)上各测站资料(或与本模式格点不同的格点资料),插值到模式的格点上.
如图,格点g周围有若干点(1)水平内插法:一个最简单又常用的插值方法.
选靠近g的N各测站资料∑∑===NiiNiiigwwff11其中iw为权重.
譬如211iirw+=,(Ni,.
.
.
,2,1=)表明各点的f值对gf的贡献大小,如果0=ir,则贡献为100%.
如r大则贡献小.
如果得到的不是测站资料,而是格点资料,但其格点与模式格点不一样,也可用类似方法来做(如网格嵌套为其他模式提供资料).
对于数值预报,初值要求更严些,常常要做些订正.
如算出高度场后要用风来订正(地转风原理).
g第第七七章章初初值值处处理理-62-如图,风订正(实线)与不订正不大一样.
对短期数值预报很重要,但对长期数值模拟试验则不是很重要.
用这个方法插值时一定要注意:周围各点与格点g一定要在同一高度或同一等压面上.
如果不在同一等压面上,则一定要先求出周围各点在与格点同一气压处的值再进行插值.
所以通常是先在等压面上插值到各格点,再垂直插值到各模式层.
要是反过来,先垂直插值再水平插值会很麻烦的.
如σ面上格点,不好用σ面上其他点的值来插(因为各点上P不一样)要先在P1面和P2面上把观测值插到A1、B1、C1上和A2、B2、C2上,再垂直插值.
(2)多项式法:假设某气象场(如位势高度Z)可用一个x和y的多项式来描述,如:MjijiyxaZjijiij≤+≥=∑,0,,如3=M时为三次曲面,2=M时为二次曲面.
例:2=M时,20201112201000yayaxyaxaxaaZ+++++=只要求出ija,从各格点的x、y,就可以求出格点的Z值.
为了要定出系数ija,对于2=M时,至少要有6点的观测值,但有时多于6点,则更好.
可以部分地去掉观测值的随机误差.
如果有一组(N个测点)的观测值obNobobZZZ,.
.
.
,,21,同时根据这N个点的位置(x,y)可以用上式写出:20201112201000220220122112220210002210210111112120110001NNNNNNNyayayxaxaxaaZyayayxaxaxaaZyayayxaxaxaaZ+++++=+++++=+++++=用最小二乘法就可以求出ija.
即令:A1B1C1P1A2B2C2P2580579581581582ABCσ数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-63-∑==NnnobnZZE12)(达极小求0=ijaE时的ija.
20≤+≤ji由此可有ija的6个方程,可求出ija.
如果用上风的资料,计算的Z更好,此时可用上地转风近似.
∑∑==++=NnnobnnobnVNnnobnvvuuwZZE122212])()[()(其中obnu和obnv为实测的.
nu和nv用地转风公式.
2Vw是风的观测值对高度测值的权重(单位为秒).
由nZ来算:如2=M时,)2()2(101120010211ayaxafgxZfgvayaxafgyZfgunnnnnnnn++==++==同样求E最小时的ija,即可.
实际工作中,这N个测站要选在要计算的格点附近,不要离得太远.
但是,有时候,在缺少测站的地方,在格点附近测站不到6个(2=M时),无法确定ija的值,这时可以这样做:①先做测站多的地区的格点;②把计算好的格点值当成测站值,加上另一些测站值来计算与它们相邻的格点的值.
此法优点:考虑了地转风关系和网格点附近的水平方向的连续性,数学上比较严格.

缺点:没有考虑时间连续性;在记录少的地区可能发生较大误差,产生计算不稳定;运算量比较大.
(3)逐步订正法:为了克服多项式在记录稀少地区应用的困难,设计了逐步订正法.
此法的思第第七七章章初初值值处处理理-64-路是,先作一个预备场,再利用预备场和观测值的差,进行订正.
预备场可以是12小时前的预备值,或气候平均值,或最近一次的分析值,或这几种值的加权平均.
以高度场为例,订正方法如下:①如果某测站Si只有高度的观测值obSiZ,附近格点预备场的高度值为GZ.
从GZ可以用前面的内插值法求出Si点处预备场的值GSiZ,于是,观测值obSiZ和GSiZ有一差值,用这差值来订正格点值GZ.
订正值为)(GSobSiZiiiZZwC=,其中iw为权重因子.
如果在格点G周围选N个测站,则总订正值为∑∑===NiiNiZiGwCC11.
订正后的GZ为GGGCZZ+=~②如果测站S中有高度观测值,又有测风值,则订正值为:])([GSobSobSobSiZViiiiiZyuxvgfkZwCΔΔ+=其中k表示实测风相对于高度的重要性,可由平均的uug/来决定.
取某一范围的半径R,取此范围内的测站来订正.
若此范围内有ZN个点只有测高资料,ZVN个点有高度和风资料,则总的订正值为:∑∑∑∑====++=ZVZZVZNiiNiiNiZViNiZiGwwCCC1111iw有时用下式计算:RrwRrrRrRwii>==≤(0一般取为20°或30°纬度)②211212121fffffff>=+时,平衡方程是个椭圆型方程,可以解出Ψ.
再由,uvyxΨΨ==可求出初值,uv解平衡方程比较麻烦,要先把它写成差分形式,然后用迭代法求Ψ.
(3)地转关系与平衡方程混合使用由于低纬地区风场比气压场清楚(等压线稀疏,不好分析),所以可用风场反算Z场,具体做法:①用实测风求Ψ先用实测风求出散度D和涡度,再由速度势与D的关系和流函数与涡第第七七章章初初值值处处理理-68-度的关系求出χ和ψ2Dχ=,边界上0χ=,解泊松方程求出χ后,由nVnψχ=+求出边界上的值再由2ψ=求出ψ,再求出风Vψ②由ψ用平衡方程反算高度场Z③加权平均求高度在低纬用Vψ,在中高纬用地转风gV,它们之间用Vψ与gV的加权平均.
2、动处理方法(研究生阶段讲)3、变分处理方法(研究生阶段讲)第三节四维同化以上用客观分析处理的资料基本上是定时的常规观测资料.
自从气象卫星升空(第一颗是1960年Tiros号——泰尔斯号)以后,不定时的非常规资料大量增加.
如果能利用起来,对弥补洋面和沙漠地区观测资料的不足是很有用的.
但有一个最大的问题,对于极轨卫星来说,它只能观测到其轨道下方(星下点)附近的带状区域,即某时刻只能观测某地区.
随着地球自转,虽然可以得到全球资料,但是不同地区的资料观测时间不一样.
所谓四维同化就是把不同时刻(t)、不同地区(x,y)、不同高度(Z)、不同性质(常规或非常规)的气象观测资料源源不断地输入计算机,通过一定的程序,把它们协调起来,溶合成常规的定时资料,为数值预报提供初值(或更新预报值).
进行四维同化主要分为间歇资料同化和连续资料同化:1、间歇资料同化可用于提供初值和在数值预报中间某一时间间隔引入卫星观测资料,它包含数数值值预预报报与与数数值值模模拟拟方方法法——周任君博士-69-三部分工作:①采用时间插值,得到指定时刻的值(卫星观测值)②空间插值,把卫星资料插值到格点上(用前述客观分析方法)③调整模式预报值,向观测值逼近,且不激发气象噪音2、连续资料同化即在预报过程中不断引入卫星观测资料,因此必须使模式对新资料的"冲击"进行调整,否则很容易产生"噪音".
目前预报中一般不采用这种方法.
但有人在发展一种伴随模式来进行连续同化,计算量很大,将来可能很有发展前途.