任意方向垂直剖面

PIZI

function main()

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

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

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

say ' Click graphic windows to specify the first point cross section_line! '  // 文本提示，在文// 本窗口看到“Click graphic windows to specify the first point cross section_line!”// 的提示后，鼠标点击图形窗口以确定剖线的起点位置A

' q bpos '

x0=subwrd (result,3)

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

say ' Click graphic windows to specify the second point cross_section_line! '

' q bpos '

x1=subwrd (result,3)

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

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

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

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

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

lon0=subwrd (result,3)

lat0=subwrd (result,6)

' q  xy2w  '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的值

' d pow ('Lx',2)+pow ('Ly',2)'函数的用法

L=subwrd (result,4)

'd pow ('L',1/2)'

L=subwrd (result,4)  // L为点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上各等分点A、A、AA(A和A分别与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上各点（A、A、AA）在各个高度上的值   **********

'set t t0'

nz=1

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

i=1

'set z 'nz' '

while (i<=n)

'set lon 'lon.i' '

'set lat 'lat.i' '

'd varnam'

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