GrADS任意方向垂直剖面.docx

GrADS任意方向垂直剖面.docx

《GrADS任意方向垂直剖面.docx》由会员分享，可在线阅读，更多相关《GrADS任意方向垂直剖面.docx（11页珍藏版）》请在冰豆网上搜索。

GrADS任意方向垂直剖面

GrADS下任意方向垂直剖面的实现

缪明余清平

（解放军理工大学气象学院）

廖德敏

（93966部队气象室）

引言

GrADS（GridAnalysisandDisplaySystem）是当前气象学中应用最为广泛的图形图像处理软件之一，它功能强大，使用方便而深受用户的欢迎。

利用它可直接实现图像的动画显示，在三维数值模拟研究中，在固定模式高度的情况下，可较为方便地显示某一模式预报量在该层次上的水平分布；在固定经（纬）度的情况下，还可方便地显示某一模式预报量沿该经（纬）度的经（纬）向垂直剖面分布。

但在实际问题的三维数值模拟研究中，用户关心的不仅仅是模式预报量沿某一经纬向的垂直剖面分布，更多的时候用户关心的是沿某一方向的垂直剖面分布，但在GrADS下如何实现这一功能，这可能是困扰一般的GrADS用户的问题。

其实，GrADS除了强大的图形图象的分析处理能力外，还有较强的数据处理能力，我们可利用GrADS的交互性和数据处理的特点，利用GrADS脚本语言的编程技术，按GrADS的数据格式，通过提取所要分析的模式预报量沿某一方向上的值（ascII码形式），然后再将提取ascII码值转换成GrADS支持的二进制格式，并为转化后的数据编写相应的数据描述文件（ctl文件），在GrADS下打开它就可实现某一模式预报量沿该方向的垂直剖面分布的分析。

一、实现的思想和步骤

在介绍实现的思想和gs编程（GrADSScriptprograming）之前先说明GrADS下某些基本的常识和某些基本命令的含义和功能。

1．GrADS基本常识和基本命令的简介

（1）进入GrADS后，在“ga->”后每执行一个GrADS命令，在文本窗口都会返回1行或2行的的文本代码，返回的代码即为GrADS的默认变量“result”值，每执行一次GrADS命令，原来的“result”值将被刷新成刚执行的GrADS命令返回码。

在gs编程中，可通过代码“sayresult”看到默认变量“result”对应的返回文本代码。

（2）命令“qbpos”：

功能：

确定鼠标在图形窗口中所点击处的位置（页面坐标）

返回码：

执行完该命令后，有一行形如“position=x0y010”的文本返回码，其中x0、y0即为用户在图形窗口中点击处对应的页面坐标。

注意：

执行该命令后，GrADS将等待用户的响应，直到用户用鼠标点击图形窗口后才可继续执行其它的命令或返回到“ga->”状态下。

（3）在GrADS下成功打开一个.ctl描述文件对应的数据文件后，在某一水平显示环境下（即t、z一定的情况下）一旦执行了“df”（f为.ctl文件中的某个变量名），图形窗口任何一点就在页面坐标、经纬坐标和格点坐标三者间建立起一一对应的关系，这种对应关系只有在一定的缩放比例环境下才有意义，打开.ctl描述文件对应的数据文件后，执行命令“df”的一个目的就是要在页面坐标与经纬坐标间建立一定的缩放比例。

执行“clear”命令或“reset”命令后，页面坐标与经纬坐标间建立一定的缩放比例自动终止。

（4）命令“qxy2wx0y0”：

功能：

实现页面坐标（x0，y0）与对应的经纬坐标（lon0，lat0）之间的转换；

返回码：

在页面坐标与经纬坐标间建立一定的缩放比例后，执行该命令，在文本窗口中将有一行形如“Lon=lon0Lat=lat0”文本返回码，其中lon0、lat0为与页面坐标（x0，y0）对应的经纬度。

注意：

只有在页面坐标与经纬坐标间建立起一定的缩放比例的前提下，才会有“Lon=lon0Lat=lat0”的返回码，也才能实现页面坐标（x0，y0）与对应的经纬坐标（lon0，lat0）之间的转换。

（5）几个GrADS脚本编程语言的内部函数（intrinsicfunction）

sublin（string,nth_line）:

功能：

从字符串“string”中提取第“nth_line”行子字符串。

注意：

“nth_line”为正整数。

举例：

假设string=“Demonstratehowtouse

theintrinsicfunction！

”

则：

sublin（string,2）=“theintrinsicfunction！

”

subwrd（string,nth_word）：

功能：

从字符串“string”中提取第“nth_word”个“字（word）”。

注意：

“nth_word”为正整数；在GrADS下，用空格符来作为“字（word）”与“字（word）”间的分隔符，因此“字（word）”是指任何两相邻两空格间的非空隔字符或非空隔字符串。

举例：

如：

string的值同上，则：

subwrd（string,6）=“intrinsic”；

如：

string=“Demonstratehowtouse

theintrinsic-function！

”，即“intrinsic”与“function”用连字符“-”相连后，则：

subwrd（string,6）=“intrinsic-function！

”

再如：

“Demonstratehowtouse

theintrinsicfunction！

”，即“intrinsic”与“function”间多加入几个空隔后，则：

subwrd（string,6）=“intrinsic”。

用户通过这两个内部函数可将执行完GrADS命令后返回的文本代码中有用的部分提取出来以便进一步的处理。

2．实现思想

任意模式预报量f可看成是空间和时间的函数f（x,y,z,t），在GrADS下，视模式预报量f为x,y,z,t的一个四维变量，当固定时间（t=t0）和垂直方向的层次（z=z0），通过命令“displayf”的可得到模式预报量f在t=t0时刻，z=z0模式高度上的二维水平分布；当固定所有维，即x=x0,y=y0,z=z0和t=t0一定的情况下，“displayf”得到的是模式预报量f在时间t=t0时刻、在空间某一点M（x0,y0,z0）的值，用户可利用GrADS语言的编程技术编写一个简单的gs程序将所关心的模式预报量f在固定时空的情况下按照GrADS的数据格式依次地沿某一方向、自下而上地将该模式预报量的值完全提取出来，并将它以文件的形式保存，然后通过一个简单的Fortran程序将提取出的模式预报量f的ascII码形式的值转换成二进制，并为转化后的数据编写相应的数据描述文件，在GrADS下打开该描述文件，即可显示该模式预报量f沿该方向的垂直剖面分布。

3．实现的步骤

（1）确定要分析的垂直剖面的方向，即确定垂直剖面的剖线AB。

（2）按一定的精度（水平分辨率）确定剖线上点A0、A1、A2…An（A0和An分别与A和B重合）的位置选取的水平分辨率最好等于或略大于格点资料的水平分辨率。

（3）按先从A0到An、然后从下到上的顺序提取某一时刻模式预报量f在直线AB上各点（A0、A1、A2…An）在各个高度上的值，并将提取的值以文件的形式保存以便处理。

（4）编写Fortran程序将提取的值转化成二进制的的数据格式。

（5）为转化后的GrADS二进制数据编写相应的.ctl描述文件。

（6）在GrADS打开.ctl描述文件，对它进行正常操作，完成物理量f沿剖线AB的垂直剖面分布的分析。

二、实例分析与应用

1．可靠性检验

（1）利用直接法与利用本法的经纬向的垂直解剖面分布的比较

如图1（a）和图（b）均为某一时刻三维模拟的预报量u沿同一纬度的纬向垂直剖面分布，其中图（a）是直接利用GrADS自身支持的沿某一经纬向的垂直解剖面分布显示功能得到的结果，图（b）是利用本文所介绍方案的一个显示结果，从中可看出，利用本文所介绍的方案显示的结果与直接利用GrADS自身支持的沿某一经纬向的垂直解剖面分布显示功能得到的结果是较为一致的。

从图2（a）和（b）中也可得到类似的结论。

可见，本文所介绍方案是可靠的。

（a）（b）

图1模拟u分量的纬向垂直剖面分布

（a）直接法（b）间接法

（a）（b）

图2模拟u分量的经向垂直剖面分布

（a）直接法（b）间接法

2．实际应用

如图3，该图是一次飑线过程三维数值模拟得到的某一水平截面的垂直速度场分布，图中呈西南-东北走向带状分布的垂直上升运动即为与飑线相联系的强对流带，在飑线的垂直结构分析中，分析研究某些模式预报量沿横越飑线方向（与强对流带垂直）的垂直分布是十分必要的，利用本文所介绍的方案很容易地实现模式预报量沿横越飑线方向的垂直剖面分布的分析。

按照上述所介绍的方法步骤，即可显示模式预报量沿横越飑线方向的垂直剖面分布。

图1：

垂直剖面剖线示意图，其中AB为所选取的剖线.

如图4，该就是按文中介绍方法的到的沿直线AB方向的相

对系统（飑线）u、垂直速度w和u-w流场垂直分布。

（a）水平速度沿直线AB分量（m/s）

（b）垂直速度W（m/s）

（c）U-W流场

图4模拟结果沿直线AB的垂直剖面分布

三结束语

本文介绍了在GrADS下如何实现任意方向的垂直剖面分析，通过本法与GrADS自身支持的沿经纬向的垂直剖面的分析比较，表明该法在分析沿某一方向的垂直剖面分析中是可信的。

结合笔者三维模拟飑线中遇到的一个实际问题，利用该法较好地分析了相对飑线的水平速度u分量、垂直速度w和垂直环流结构，发现该法在实际应用中是可行的。

参考文献：

[1]BrianDoty，GrADS用户手册。

附录：

数据提取gs源程序

functionmain（）

'openXXX\xxx.ctl'//打开数据描述文件（“XXX”为路径名，“xxx.ctl”文件名）

'dvarnam'//显示描述文件中某一物理量“varnam”以便更好确定剖线位置和建立//图形窗口中页面坐标与经纬坐标的缩放比例环境（scalingenvironment）

******************************确定剖线位置***********************************

say'Clickgraphicwindowstospecifythefirstpointcrosssection_line!

'//文本提示，在文//本窗口看到“Clickgraphicwindowstospecifythefirstpointcrosssection_line!

”//的提示后，鼠标点击图形窗口以确定剖线的起点位置A

'qbpos'

x0=subwrd（result,3）

y0=subwrd（result,4）//图形窗口上鼠标点击处（剖线起点A）对应的页面坐标（x0,y0）

say'Clickgraphicwindowstospecifythesecondpointcross_section_line!

'

'qbpos'

x1=subwrd（result,3）

y1=subwrd（result,4）//剖线终点B对应的页面坐标（x1,y1）

'drawline'x0''y0''x1''y1''//在图形窗口中画出剖线，即连接点A（x0,y0）到点

//B（x1,y1）的线段注意此处的单引号

*****求解点A（x0,y0）和点B（x1,y1）对应的经纬坐标A（lon0,lat0）和B（lon1,lat1）*****

'qxy2w'x0''y0''//将页面坐标（x0,y0）转换为对应的经纬坐标（lon0,lat0）

lon0=subwrd（result,3）

lat0=subwrd（result,6）

'qxy2w'x1''y1''//将页面坐标（x0,y0）转换为对应的经纬坐标（lon1,lat1）

lon1=subwrd（result,3）

lat1=subwrd（result,6）

**********求解点A（lon0,lat0）到点B（lon1,lat1）的实际距离（L）*******************************************************

Lx=（lon1-lon0）*DisPerLon//DisPerLon为单位经度对应的实际距离（112km）

Ly=（lat1-lat0）*DisPerLat//DisPerLat为单位纬度对应的实际距离（223km）

//实际应用中必须给出DisPerLon和DisPerLat的值

'dpow（'Lx',2）+pow（'Ly',2）'函数的用法

L=subwrd（result,4）

'dpow（'L',1/2）'

LAB=subwrd（result,4）//LAB为点A（lon0,lat0）到点B（lon1,lat1）的实际距离

*********求取相邻各等距点间的经度增量（inclon）和纬度增量（inclat）*********

decimal=L/horRes//horRes为直线上相邻各点间的距离（水平分辨率）

//实际应用中必须给出horRes的值，给出的值略大于模式水平分辨//率为宜

decimal_part=substr（decimal,3,20）

n=decimal-decimal_part

inclat=（lat1-lat0）/n

inclon=（lon1-lon0）/n

***计算直线AB上各等分点A0、A1、A2…An（A0和An分别与A和B重合）所在经纬度****

i=1

while（i<=Num_of_point）

lon.i=lon0+（i-1）*inclon

lat.i=lat0+（i-1）*inclat

i=i+1

endwhile

**********按GrADS数据格式提取t=t0时刻模式预报量“varnam”**********

*********在直线AB上各点（A0、A1、A2…An）在各个高度上的值**********

'settt0'

nz=1

while（nz<=NZ）//NZ为模式的垂直分层大小

i=1

'setz'nz''

while（i<=n）

'setlon'lon.i''

'setlat'lat.i''

'dvarnam'

dummy=sublin（result,2）

data=subwrd（dummy,4）

write（'XXX\gad.dat',data,‘append’）//将提取的数据保存到文件“gad.dat”中（“XXX”为路径名）//

i=i+1

endwhile

nz=nz+1

endwhile

close（'XXX\varnam.dat'）

return