实习一 Fortran文件的读取与处理

实习1a：已知1951-2010年1月蒙古高压强度指数、面积指数、经度指数、纬度指数序列资料分别为p.dat、s.dat、lon.dat和lat.dat。利用Fortran语言编写程序，调用子程序计算这四个指数的气候态（均值）、变率（均方差）和距平值，并将这四个指数的均值和变率写入十进制文件mh1.dat和二进制文件mh1.grd，将这四个指数的距平值写入十进制文件mh2.dat和二进制文件mh2.grd. 文件均保存于sx01文件夹下。

实习1b：已有实习生成的1951-2010年1月蒙古高压强度指数、面积指数、经度指数、纬度指数这四个指数的距平值二进制文件mh2.grd，数据描述文件mh2.ctl 保存于sx01文件夹下。利用Grads下列语句输出这四个指数距平的一维时间分布图，并逐一截屏保存。

1. 编写fortran文件生成数据

program mh
   implicit none
   integer, parameter::ny = 60
   real p(ny), s(ny), lon(ny), lat(ny) !原始数组
   real pa(ny), sa(ny), lona(ny), lata(ny) !距平
   real pav, sav, lonav, latav !均值
   real pd, sd, lond, latd !标准差
   integer i, j, k

!利用open语句打开强度指数数据p.dat
   open (1, file='p.dat')
!利用open语句打开面积指数数据s.dat
   open (2, file='s.dat')
!利用open语句打开经度指数数据lon.dat
   open (3, file='lon.dat')
!利用open语句打开纬度指数数据lat.dat
   open (4, file='lat.dat')

   do i = 1, ny
      read (1, *) p(i)
      read (2, *) s(i)
      read (3, *) lon(i)
      read (4, *) lat(i)
   end do
   close (1); close (2); close (3); close (4)

!调用气候及异常值计算的子程序
   call cha(ny, p, pa, pav, pd)
   call cha(ny, s, sa, sav, sd)
   call cha(ny, lon, lona, lonav, lond)
   call cha(ny, lat, lata, latav, latd)

!将蒙古高压环流指数气候值和标准差写入到mh1.dat!用open语句打开文件mh1.dat
   open (5, file='mh1.dat')
!写入pav，pd
   write (5, *) pav, pd
!写入sav，sd
   write (5, *) sav, sd
!写入lonav，lond
   write (5, *) lonav, lond
!写入latav，latd
   write (5, *) latav, latd
   close (5)

!将蒙古高压环流指数气候值和标准差写入到mh1.grd!用open语句打开文件mh1.grd
   open (6, file='mh1.grd', form='binary')
!写入pav，pd
   write (6) pav, pd
!写入sav，sd
   write (6) sav, sd
!写入lonav，lond
   write (6) lonav, lond
!写入latav，latd
   write (6) latav, latd
   close (6)

!将蒙古高压环流指数距平值写入到mh2.dat，按照pa(60),sa(60),lona(60),lata(60)顺序存放
!用open语句打开文件mh2.dat
   open (7, file='mh2.dat')
   write (7, *) (pa(i), i=1, ny)
   write (7, *) (sa(i), i=1, ny)
   write (7, *) (lona(i), i=1, ny)
   write (7, *) (lata(i), i=1, ny)
   close (7)

!将蒙古高压环流指数距平值写入到mh2.grd，要求利用do循环按照pa(i),sa(i),lona(i),lata(i)顺序存放
!用open语句打开文件mh2.grd
   open (8, file='mh2.grd', form='binary')

   do i = 1, ny
      write (8) pa(i)
      write (8) sa(i)
      write (8) lona(i)
      write (8) lata(i)
   end do
   close (8)
end

subroutine cha(ny, x, xa, xav, xd) Implicit Noneinteger i, ny
real x(ny), xa(ny), xav, xd, sum!依次为原序列、距平、均值、标准差、和
!求原序列的和sum
sum = 0.0
do i = 1, ny
   sum = sum + x(i)
end do
!求平均值xavxav=sum/ny
!求距平xa和标准差xd
xd = 0.0
do i = 1, ny
   xa(i) = x(i) - xav
   xd = xa(i)*xa(i) + xd
end do
xd = sqrt(xd/ny)
return
end

2. 编写gs文件
   mh1.ctl

dset mh1.grd
undef -9.99E+33
title NCEP/NCAR
REANALYSIS PROJECT
xdef 1 linear 60.0 2.5
ydef 1 linear 0.0 2.5
zdef 1 levels 850
tdef 1 linear JAN2002 1mo
vars 8
pav 0 99 p ave
pd 0 99 p d
sav 0 99 s ave
sd 0 99 s d
lonav 0 99 lon ave
lond 0 99 lon d
latav 0 99 lat ave
latd 0 99 lat d
endvars

• 已有实习生成的1951-2010年1月蒙古高压强度指数、面积指数、经度指数、纬度指数这四个指数的距平值二进制文件mh2.grd，数据描述文件mh2.ctl 保存于sx01文件夹下。
• 利用Grads基本语句输出这四个指数距平的一维时间分布图，并逐一截屏保存。

mh2.ctl

dset mh2.grd
undef -9.99E+33
title NCEP/NCAR
REANALYSIS PROJECT
xdef 1 linear 60.0 2.5
ydef 1 linear 0.0 2.5
zdef 1 levels 850
tdef 60 linear jan1951 1yr
vars 4pa 0 99 p anomaly
sa 0 99 s anomaly
lona 0 99 lon anomaly
lata 0 99 lat anomaly
endvars

3. 编写gs文件

'reinit'
'open mh2.ctl'
'set x 1'
'set y 1'
'set z 1'
'set t 1 60'
'd pa'
'c'
'd sa'
'c'
'd lona'
'c'
'd lata'

实习二 数据文件的转换及数据描述文件的建立、基本操作命令

现有sx02文件夹下有十进制月平均风场数据：200hPa纬向风u200.dat、经向风v200.dat；850hPa纬向风u850.dat、经向风v850.dat。时间范围：2002.1-2005.12共48个月，经纬度范围：60-150°E，0-40°N，分辨率为2.5°×2.5°。
（1）编写Fortran程序，将十进制数据资料文件转换成二进制月平均风场文件：uv.grd；
（2）写出二进制文件相应的数据描述文件：uv.ctl；
（3）利用GrADS基本操作命令显示该地区显示2002年7月850hPauv风场和200hPa uv风场，并逐一保存图形。

编程时数组大小：

• X方向格点数：（初始经度-结束经度）/格距+1
  （150-60）/2.5+1=37
• Y方向格点数：（初始纬度-结束纬度）/格距+1
  （40-0）/2.5+1=17
• Z方向层数： u、v为850、 200hPa ，Z为2
• T时次：月平均资料，4年，共48个月

sx02.f90

program convert_to_binary
implicit none

integer i, j, it
integer, parameter::nx=37,ny=17,nt=48
real u200(nx,ny,nt),v200(nx,ny,nt),u850(nx,ny,nt),v850(nx,ny,nt)

  open(1, file='u200.dat')
  open(2, file='v200.dat')
  open(3, file='u850.dat')
  open(4, file='v850.dat')

  do it=1,nt
    read(1,*)((u200(i,j,it),i=1,nx),j=1,ny)
    read(2,*)((v200(i,j,it),i=1,nx),j=1,ny)
    read(3,*)((u850(i,j,it),i=1,nx),j=1,ny)
    read(4,*)((v850(i,j,it),i=1,nx),j=1,ny)
  enddo
  close(1);close(2);close(3);close(4)
  
  open(12, file='uv.grd', form='unformatted', access='stream')
  do it=1,nt
    write(12) ((u850(i,j,it),i=1,nx),j=1,ny)
    write(12) ((u200(i,j,it),i=1,nx),j=1,ny)
    write(12) ((v850(i,j,it),i=1,nx),j=1,ny)
    write(12) ((v200(i,j,it),i=1,nx),j=1,ny)
  enddo
  close(12)
end

uv.ctl

title NCEP/NCAR REANALYSIS PROJECTION
undef -9.99E+33
xdef 37 linear 60.0 2.5
ydef 17 linear 0.0  2.5
zdef 2 levels 850 200
tdef 48 linear JAN2002 1mo
vars 2
u 2 99 u wind(m/s)
v 2 99 v wind(m/s)
endvars

uv-1.gs

'reinit'
'open uv.ctl'
'set grads off'
'set gxout vector'
'set t 1'
'set z 1'
'd u;v'
'draw title uv at 850 in Jan2002'
'printim 1a.png white'

uv-2.gs

'reinit'
'open uv.ctl'
'set grads off'
'set gxout vector'
'set z 1'
'set t 7'
'd u;v'
'draw title uv at 850 in Jul2002'
'printim uv850_1.png white'
'c'
'set grads off'
'set gxout vector'
'set z 2'
'set t 7'
'd u;v'
'draw title uv at 200 in Jul2002'
'printim uv200_1.png white'
'c'


'reinit'
'open uv.ctl'
'set grads off'
'set gxout vector'
'set time jul2002'
'd u;v'
'draw title uv at 850 in Jul2002'
'printim uv850_2.png white'
'c'
'set grads off'
'set gxout vector'
'set lev 200'
'set time jul2002'
'd u;v'
'draw title uv at 200 in Jul2002'
'printim uv200_2.png white'
'c'

实习三 GrADS绘图要素的设置和基础绘图命令

（1）一维单变量绘图
1951-2010年1月蒙古高压强度指数、面积指数、经度指数、纬度指数这四个指数的距平已按照GrADS要求写入二进制文件mh2.grd，其数据描述文件为mh2.ctl，mh2.grd和mh2.ctl保存于sx03文件夹下。请将1951-2010年强度指数距平和面积指数距平绘制为折线，将经度指数距平和纬度指数距平绘制为直方图。

（2）二维单变量绘图
NCEP/NCAR再分析纬向风月平均数据的二进制数据uwnd.mon.mean.nc存于sx03文件夹下。绘制2023年7月120 °E、 0-90N 、1000~100hPa纬向风的纬度-层次剖面图，并以120 °E为标题，将图形保存为u.png。

mh2.gs

'reinit'
*打开 mh2.ctl 文件
'open mh2.ctl'
*设置 x 维数环境
'set x 1'
*设置 y 维数环境
'set y 1'
*设置 z 维数环境
'set z 1'
*设置时间维数环境
'set t 1 60'
***********绘制强度指数距平为折线图
*设置图形出图类型为 line
'set gxout line'
*设置 line 的颜色
'set ccolor 1'
*设置 line 的线形
'set cstyle 1'
*设置 line 的粗细
'set cthick 5'
'set cmark 2'
*显示蒙古高压强度指数距平 pa
'd pa'
***********绘制面积指数距平为折线图
*设置 line 的颜色
'set ccolor 2'
*设置 line 的线形
'set cstyle 3'
*设置 line 的粗细
'set cthick 5'
'set cmark 3'
*显示蒙古高压面积指数距平 sa*注意，若值较小，考虑处理方法
'd sa*10'
*将图形保存为 pasa.png，背景为白色
'printim pasa.png white'
*下面的'c'不可缺少，为什么?
'c'
***********绘制经度指数距平为直方图
*设置图形出图类型为 bar
'set gxout bar'
*设置 bar 的绘制方向
'set barbase 0'
*设置直方条的间隔
'set bargap 50'
*设置 bar 的颜色
'set ccolor 3'
*显示蒙古高压经度指数距平 lona
'd lona'
*将图形保存为 lona.png，背景为白色
'printim lona.png white'
'c'
***********绘制纬度指数距平为直方图
*设置图形出图类型为 bar
'set gxout bar'
*设置 bar 的绘制方向
'set barbase 0'
*设置直方条的间隔
'set bargap 30'
*设置 bar 的颜色
'set ccolor 4'
*显示蒙古高压纬度指数距平 lata
'd lata'
*将图形保存为 lata.png，背景为白色
'printim lata.png white'
;

u.gs

'reinit'
*打开 uwnd.mon.mean.nc 数据
'sdfopen uwnd.mon.mean.nc'
*设置经度为东经 120 度
'set lon 120'
*设置纬度范围为赤道 0 度至北纬 90 度
'set lat 0 90'
*设置等压面从 1000hPa 至 100hPa
'set lev 1000 100'
*设置时间为 2023 年 7 月
'set time jul2023'
*设置 Z 方向取对数尺度
'set zlog on'
*设置出图类型为等值线 contour
'set gxout contour'
*显示纬向风
'd uwnd'
*写标题
'draw title 120`3.`1E'
*存图为 u.png，白色背景
'printim u.png white'
'reinit'
;

实习四 GrADS 变量与函数及描述语言的应用

sx04文件夹下有：
全球月平均降水数据pr_wtr.eatm.mon.mean.nc，1层。
时间从1948年1月开始，水平格距 2.5*2.5

• nc文件打开：sdfopen filename.nc
• nc文件的维数设置、格距、缺测值的查询：
• q ctlinfo

编写gs文件提取1951年至2010年 60年 7月降水数据pre7.grd

'reinit'
'sdfopen pr_wtr.eatm.mon.mean.nc'
'set gxout fwrite'
'set fwrite pre7.grd'
'set x 1 144'
'set y 1 73'
'set z 1'
* 年份循环
iyear=1951
while(iyear<=2010)
'set time jul'iyear''
'd pr_wtr'
iyear=iyear+1
endwhile
'disable fwrite'
'reinit'

* 时次循环
i=43
while(i<=756)
'set t 'i''
'd pr_wtr'
i=i+12
endwhile
'disable fwrite'
'reinit'
;

编写pre7.grd的数据描述文件pre7.ctl；

dset pre7.grd
undef -9.99e+08
xdef 144 linear 0 2.5
ydef 73 linear -90 2.5
zdef 1 linear 0 1
tdef 60 linear jul1951 1yr

vars 1
pre 0 99 precip. at Jul
endvars

3. 计算7月降水的60年平均值，绘制70-140ºE、15-55ºN区域降水的多年平均值，并在我国降水带大值中心（117.5E，32.5N）标注字母W，存为.png图形

'reinit'
'open pre7.ctl'
'set lon 0 360'
'set lat -90 90'
'set z 1'
'set t 1'
'define prclim=ave(pre,t=1,t=60)'
'set grid off'
'set grads off'
'set lon 70 140'
'set lat 15 55'
'set t 1'
'd prclim'
'q w2xy 117.5 32.5'
x1=subwrd(result,3)
y1=subwrd(result,6)
'set string 2 c 8 0'
'set strsiz 0.2'
'draw string 'x1' 'y1' W'
'printim preclim.png white'
'reinit'
;

4）计算1951年至2010年 7月降水距平，绘制1998年7
月全球降水距平场，要求：
A、正距平区填色，画色标；
B、设置等值线线间隔为2；
C、保存为.png图形，白底。
注：1998年为强厄尔尼诺年次年夏天。2000年夏天为
拉尼娜年次年夏天。有兴趣的同学可以再绘制2000年
7月降水距平，与1998年7月降水距平对比看看。

'reinit'
'open pre7.ctl'
'set lon 0 360'
'set lat -90 90'
'set z 1'
'set t 1'
'define prclim=ave(pre,t=1,t=60)'
'set t 1 60'
'define pra=pre-prclim(t=1)'
*----------绘图-----------------
'define_colors'
'set grads off'
'set grid off'
'set lon 0 360'
'set lat -90 90'
'set z 1'
'set time jul1998'
'set gxout shaded'
'set cmin 0'
'set cint 2'
*'set clevs 0 2 4 6 8 10 12 14'
*'set ccols 0 41 42 43 44 45 46 47 48'
'd pra'
'cbarn 1 0 '
'set gxout contour'
'set cint 2'
'd pra'
'draw title Jul1998'
'printim 1998.png white'
'reinit'
;

```glsl
'reinit'
'open pre7.ctl'
'set lon 0 360'
'set lat -90 90'
'set z 1'
'set t 1'
'define prclim=ave(pre,t=1,t=60)'
'set t 1 60'
'define pra=pre-prclim(t=1)'
*----------绘图-----------------
'define_colors'
'set grads off'
'set grid off'
'set lon 0 360'
'set lat -90 90'
'set z 1'
'set time jul1998'
'set gxout shaded'
*'set cmin 0'
*'set cint 2'
'set clevs 0 2 4 6 8 10 12 14'
'set ccols 0 41 42 43 44 45 46 47 48'
'd pra'
'cbarn 1 0 '
'set gxout contour'
'set cint 2'
'd pra'
'draw title Jul1998'
'printim pra.png white'
'reinit'
;

实习五 Nino3.4海温指数与全球降水相关分析

1、所用数据：

• 1951年1月至2013年12月Nino3.4区（热带太平洋西经170度_{西经120度、北纬5度}南纬5度）区域平均海温指数资料Nino34.txt。共有63行数据，每行数据第一个数字为年份，后面12个数字为该年1-12月的海温指数；
• 1951年至2010年 60年7月全球降水数据pre7.grd，x、y方向格点数分别为144和73，水平格距2.5*2.5，为地面变量。相应ctl文件为pre.ctl。

2、实习要求：
1） 用Fortran编写corr.grid.f90文件，计算1951-2010年 1月Nino34海温指数与7月降水的相关系数，计算结果保存于corr.txt和corr.grd 文件中；
2） 编书写corr.grd 的数据描述文件corr.ctl；
3） 编写corr.gs文件，绘制1951-2010年 1月Nino34海温指数与7月降水的相关系数等值线图，通过95%置信度检验的区域填色（60年相关系数95%、99%置信度检验的临界值分别为0.25和0.325）。

program sx05
   integer, parameter::nt = 60, in = 144, jn = 73, nm = 12
   real pre(in, jn, nt), yh(nm, nt), rr(in, jn)
!------------读取 7 月降水场-------------------
   open (40, file='pre7.grd', form='unformatted')
   do it = 1, nt
      read (40) ((pre(ix, iy, it), ix=1, in), iy=1, jn)
   end do
   close (40)
!-------------读取 nino3.4 sst 指数-----------------
   open (2, file='nino34.txt')
!1951 年 1 月-2010 年 12 个月 nino3 区海温指数
   do it = 1, nt
      read (2, *) iyear, (yh(k, it), k=1, nm)
      write (*, *) yh(1, it)
   end do
   close (2)
   do iy = 1, jn
   do ix = 1, in
      call correlation(nt, pre(ix, iy, :), yh(1, :), r)
      rr(ix, iy) = r
   end do
   end do
!---------写出十进制的相关系数文件--------------
   open (3, file='corr.txt')
   do iy = 1, jn
   do ix = 1, in
      write (3, *) rr(ix, iy)
   end do
   end do
   close (3)
!---------写出可用于 Grads 绘图的二进制文件---
   open (4, file='corr.grd', form='unformatted')
   do iy = 1, jn
   do ix = 1, in
      write (4) rr(ix, iy)
   end do
   end do
   close (4)
end
!------求两个一维时间序列的相关系数子程序-
!n 为时间长度，x,y 为两个时间序列，r 为相关系数
subroutine correlation(n, x, y, r) real x(n), y(n)
   ave1 = 0.0; ave2 = 0.0; Var1 = 0.0; Var2 = 0.0
   do i = 1, n
      ave1 = ave1 + x(i)/real(n)
      ave2 = ave2 + y(i)/real(n)
   end do
   do i = 1, n
      Var1 = Var1 + (x(i) - ave1)**2
      Var2 = Var2 + (y(i) - ave2)**2
   end do
   tmp = 0.0
   do i = 1, n
      tmp = tmp + (x(i) - ave1)*(y(i) - ave2)
   end do
   r = tmp/sqrt(Var1*Var2)
end

2、ctl 文件

dset corr.grd
undef -9.99E+33
xdef 144 linear 0 2.5
ydef 73 linear -90 2.5
zdef 1 levels 1000
tdef 1 linear 01jan1951 1mo
vars 1
corr 0 99 q1
endvars

3、gs 文件

'reinit'
'open corr.ctl'
'define_colors'
'set lat -90 90'
'set lon 0 360'
'set t 1'
'set parea 1 10 1 8'
'set xlopts 1 2 0.15'
'set ylopts 1 2 0.15'
'set grads off'
'set grid off'
'set poli on'
'set gxout shaded'
'set clevs -0.325 -0.25 0.25 0.325'
'set ccols 45 43 0 62 64 '
'd corr'
'cbarn'
'set gxout contour'
'd corr'
'printim corr.png white'
'reinit'
;