2024年气象统计上机报告

一.资料介绍

现有全球海平面气压场资料，文件名NCEP_slp_30y_Wt.dat，时段：冬季1978~2007年共30年。水平分辨率：7.5*7.5（具体参考NCEP_slp_30y_Wt.ctl文件），格点数：48*24。

要求：

1）利用合成分析方法，分析El Nino和La Nina年（根据实习一中的结果）全球海平面气压异常场之间是否存在显著差异（需检验）？

2）利用相关分析方法，计算热带太平洋Nino3.4区海温指数（实习一中的结果）与全球海平面气压场之间的相关关系。

3）分别计算南方涛动指数SOI与Nino3.4区海温指数各自超前滞后相关系数（自相关）及两者之间的超前滞后相关系数（交叉相关）。注意：这里取滞后时间为1~15年。

这里我们可以根据题2中确定最大/最小值相关系数的中心位置, 计算该位置海平面气压场差值近似定义SOI指数（东太减西太）。

二、程序

第一小题：

program main

implicit none

integer,parameter::nx=48,ny=24,nt=30,n=5 !n为5年的El与La

external average1

integer it,ix,iy

real :: E(nx,ny,5),L(nx,ny,5),El(nx,ny)=0.0,La(nx,ny)=0.0,slp(nx,ny,nt),dif(nx,ny) !E,L为原数据，El，La为距平值

real Es(nx,ny),Ls(nx,ny),t !带s的为均方差

!读取海平面气压

open (1,file='D:\2021tj\NCEP_slp_30y_Wt.dat',form='binary')

do it=1,nt

read(1)((slp(ix,iy,it),ix=1,nx),iy=1,ny)

enddo

close(1)

! 给El 及La赋海平面气压值

do ix=1,nx

do iy=1,ny

E(ix,iy,1)=slp(ix,iy,5)

E(ix,iy,2)=slp(ix,iy,9)

E(ix,iy,3)=slp(ix,iy,14)

E(ix,iy,4)=slp(ix,iy,20)

E(ix,iy,5)=slp(ix,iy,25)

L(ix,iy,1)=slp(ix,iy,7)

L(ix,iy,2)=slp(ix,iy,11)

L(ix,iy,3)=slp(ix,iy,21)

L(ix,iy,4)=slp(ix,iy,22)

L(ix,iy,5)=slp(ix,iy,29)

enddo

!计算出Elnino和Lanina的平均值之后，用差值场求出是否有显著性差异

call average1(E,n,El)

call average1(L,n,La)

do ix=1,nx

do iy=1,ny

dif(ix,iy)=El(ix,iy)-La(ix,iy) !difference为El与La的差值场。

enddo

open(2,file='D:\2021tj\practice2\difference.grd',form='binary')

write(2) ((dif(ix,iy),ix=1,nx),iy=1,ny)

close(2)

call jfc(E,El,nx,ny,n,Es)

call jfc(L,La,nx,ny,n,Ls)

!write(*,*)Ls

call t_test(El,Es,n,La,Ls,n,nx,ny,n,t)

write(*,*)'t=',t

end program

subroutine average1(x,n,y)

implicit none

integer,parameter ::nx=48,ny=24

integer ix,iy,it,n

real:: x(nx,ny,n),y(nx,ny)

y(nx,ny)=0.0

do ix=1,nx

do iy=1,ny

!这五年的El Nino

do it=1,n

if(x(ix,iy,it)<2000.0)then

y(ix,iy)=y(ix,iy)+x(ix,iy,it)

else

n=n-1

endif

enddo

y(ix,iy)=y(ix,iy)/n

enddo

end subroutine

subroutine t_test(xave,xs,n1,yave,ys,n2,nx,ny,nt,t)

implicit none

integer nx,ny,n,nt

real xave(nx,ny),xs(nx,ny),yave(nx,ny),ys(nx,ny),t,s1,s2,s3

integer n1,n2,ix,iy,it

!do it=1,nt

do ix=1,nx

do iy=1,ny

s1=(n1-1)*xs(ix,iy)2

s2=(n2-1)*ys(ix,iy)2

s3=sqrt(1.0/n1+1.0/n2)

t=(xave(ix,iy)-yave(ix,iy))/(sqrt((s1+s2)/(n1+n2-2))*s3)

enddo

! enddo

end subroutine

subroutine jfc(x,ave,nx,ny,nt,s)

implicit none

integer nx,ny,nt,ix,iy,it !nt=n,即只有5年时间

real x(nx,ny,nt),ave(nx,ny),s(nx,ny)

s(nx,ny)=0.0

do ix=1,nx

do iy=1,ny

do it=1,nt

s(ix,iy)=s(ix,iy)+(x(ix,iy,it)-ave(ix,iy))2

enddo

s(ix,iy)=sqrt(s(ix,iy)/nt)

enddo

!s(nx,ny,nt)=sqrt(s(nx,ny,nt)/nt)

end subroutine

第二小题

program main

implicit none

external jvping,jvping_3

integer,parameter::nx=48,ny=24,nt=30

real sst(nt),ssta,sstp(nt)

real slp(nx,ny,nt),slpa(nx,ny),slpp(nx,ny,nt)

real:: k(nx,ny)=0.0,q1(nx,ny)=0.0,q2(nx,ny)=0.0,r(nx,ny)=0.0

integer n,it,ix,iy

open (1,file='D:\2021tj\NCEP_slp_30y_Wt.dat',form='binary')

do it=1,nt

read(1)((slp(ix,iy,it),ix=1,nx),iy=1,ny)

enddo

close(1)

open(3,file='D:\2021tj\practice2\N-original.grd',form='binary')

do it=1,nt

read(3)sst(it)

enddo

close(3)

n=nt

call jvping(sst,n,ssta,sstp)

call jvping_3(slp,slpa,slpp)

do ix=1,nx

do iy=1,ny

k(ix,iy)=0.0

q1(ix,iy)=0.0

q2(ix,iy)=0.0

do it=1,nt

k(ix,iy)=k(ix,iy)+sstp(it)*slpp(ix,iy,it)

q1(ix,iy)=q1(ix,iy)+sstp(it)2

q2(ix,iy)=q2(ix,iy)+slpp(ix,iy,it)2

enddo

r(ix,iy)=k(ix,iy)/sqrt(q1(ix,iy)*q2(ix,iy))

enddo

!write(*,*)r

open(4,file='D:\2021tj\practice2\relative.grd',form='binary')

write(4)((r(ix,iy),ix=1,nx),iy=1,ny)

close(4)

end program

subroutine jvping(x,n,ave,y)

implicit none

integer nt,n

real x(n),y(n),ave

integer it

ave=0.0

do it=1,n

if(x(it)<400)then

ave=ave+x(it)

else

n=n-1

endif

enddo

ave=ave/n

do it=1,n

y(it)=x(it)-ave

enddo

end subroutine

subroutine jvping_3(x,ave,y)

implicit none

integer,parameter::nx=48,ny=24,nt=30

real x(nx,ny,nt),y(nx,ny,nt),ave(nx,ny)

integer ix,iy,it

do ix=1,nx

do iy=1,ny

ave(ix,iy)=0.0

do it=1,nt

ave(ix,iy)=ave(ix,iy)+x(ix,iy,it)/nt

enddo

do ix=1,nx

do iy=1,ny

do it=1,nt

y(ix,iy,it)=x(ix,iy,it)-ave(ix,iy)

enddo

end subroutine

第三小题

program main

implicit none

external jvping,jvping_3

integer,parameter::nx=48,ny=24,nt=30

real sst(nt),ssta,sstp(nt),ssts

real slp(nx,ny,nt),slpa(nx,ny),slpp(nx,ny,nt)

real sstj(15),sstr(15),SOIA,SOI(nt),SOIS,SOIP(nt),SOIJ(15),SOIR(15)

real sxy(15),rxy(15)

integer n,it,ix,iy,j

open (1,file='D:\2021tj\practice2\NCEP_slp_30y_Wt.dat',form='binary')

do it=1,nt

read(1)((slp(ix,iy,it),ix=1,nx),iy=1,ny)

enddo

close(1)

open(3,file='D:\2021tj\practice2\N-original.grd',form='binary')

do it=1,nt

read(3)sst(it)

enddo

close(3)

n=nt

call jvping(sst,n,ssta,sstp)

call jfc(sst,ssta,ssts)

do j=1,15

sstj(j)=0.0

do it=1,nt-j

sstj(j)=sstj(j)+sstp(it)*sstp(it+j)

enddo

sstj(j)=sstj(j)/(nt-j)

sstr(j)=sstj(j)/(ssts2)

enddo

write(*,*)sstr

open(5,file='D:\2021tj\practice2\N_self.grd',form='binary')

write(5)(sstr(j),j=1,15)

close(5)

SOIA=0.0

do it=1,nt

SOI(it)=slp(17,13,it)-slp(32,13,it) !(120E-120W,y==2.5) 的海平面气压值相减得SOI

SOIA=SOIA+SOI(it)/nt

enddo

call jvping(SOI,nt,SOIA,SOIP)

call jfc(SOI,SOIA,SOIS)

do j=1,15

SOIJ(j)=0.0 !1-15年的滞后

do it=1,nt-j

SOIJ(j)=SOIJ(j)+SOIP(it)*SOIP(it+j)

enddo

SOIJ(j)=SOIJ(j)/(nt-j)

SOIR(j)=SOIJ(j)/(SOIS2)

enddo

open(6,file='D:\2021tj\practice2\SLP-self-relative.grd',form='binary')

do j=1,15

write(6)SOIJ(j)

enddo

close(6)

do j=1,15

sxy(j)=0.0

do it=1,nt-j

sxy(j)=sxy(j)+(SOIP(it)*sstp(it+j))

enddo

sxy(j)=sxy(j)/(nt-j)

rxy(j)=sxy(j)/(SOIS*ssts)

ENDDO

open(7,file='D:\2021tj\practice2\cross-r.grd',form='binary')

write(7)(rxy(j),j=1,15)

close(7)

end program

subroutine jvping(x,n,ave,y)

implicit none

integer nt,n

real x(n),y(n),ave

integer it

ave=0.0

do it=1,n

if(x(it)<400)then

ave=ave+x(it)

else

n=n-1

endif

enddo

ave=ave/n

do it=1,n

y(it)=x(it)-ave

enddo

end subroutine

subroutine jvping_3(x,ave,y)

implicit none

integer,parameter::nx=48,ny=24,nt=30

real x(nx,ny,nt),y(nx,ny,nt),ave(nx,ny)

integer ix,iy,it

do ix=1,nx

do iy=1,ny

ave(ix,iy)=0.0

do it=1,nt

ave(ix,iy)=ave(ix,iy)+x(ix,iy,it)/nt

enddo

do ix=1,nx

do iy=1,ny

do it=1,nt

y(ix,iy,it)=x(ix,iy,it)-ave(ix,iy)

enddo

end subroutine

subroutine jfc(x,ave,s)

implicit none

integer,parameter ::n=30

integer it

real x(n),ave,s

s=0.0

do it=1,n

s=s+(x(it)-ave)2

enddo

s=sqrt(s/n)

end subroutine

三、ctl文件与gs文件

第一小题

1.1)difference.ctl

dset D:\2021tj\practice2\difference.grd

undef 1.0E+33

title DIFFERENCE DATA

xdef 48 levels 5 12.5 20 27.5 35 42.5 50 57.5 65 72.5 80 87.5 95 102.5 110 117.5 125 132.5 140 147.5 155 162.5 170 177.5 185 192.5 200 207.5 215 222.5 230 237.5 245 252.5 260 267.5 275 282.5 290 297.5 305 312.5 320 327.5 335 342.5 350 356.25

ydef 24 levels -87.5 -80 -72.5 -65 -57.5 -50 -42.5 -35 -27.5 -20 -12.5 -5 2.5 10 17.5 25 32.5 40 47.5 55 62.5 70 77.5 85

zdef 1 levels 1000

tdef 1 linear 0z01Jan1978 1yr

vars 1

dif 0 99

endvars

1.2)difference.gs

'reinit'

'open D:\2021tj\practice2\difference.ctl'

'set lat -90 90'

'set lon 0 360'

'set lev 1000'

'set t 1'

'set gxout stream'

'd dif'

'draw title difference'

'gxprint D:\2021tj\practice2\difference.png white'

;

第二小题

1.1)correlation.ctl

dset D:\2021tj\practice2\relative.grd

undef 1.0E+33

title relative data

ydef 24 levels -87.5 -80 -72.5 -65 -57.5 -50 -42.5 -35 -27.5 -20 -12.5 -5 2.5 10 17.5 25 32.5 40 47.5 55 62.5 70 77.5 85

zdef 1 levels 1000

tdef 1 linear 0z01Jan1978 1yr

vars 1

r 0 99

endvars

1.2)correlation.gs

'reinit'

'open D:\2021tj\practice2\relative.ctl'

'set lat -87.5 85'

'set lon 5 356.5'

'set lev 1000'

'set t 1'

'd r'

'draw title Correlation between Nino 3.4 & Global SLP'

'gxprint D:\2021tj\practice2\correlation.png white'

;

第三小题

1.1) SLP-self-relative.ctl

dset D:\2021tj\practice2\SLP_self_relative.grd

undef 1.0E+33

title SLP self relative data

xdef 1 linear 5 7.5

ydef 1 linear -87.5 7.5

zdef 1 levels 1000

tdef 15 linear 0z01Jan1978 1yr

vars 1

slpr 0 99

endvars

1.2)SLP-self-relative.gs

'reinit'

'open D:\2021tj\practice2\SLP_self_relative.ctl'

'set x 1'

'set y 1'

'set lev 1000'

'set t 1 15'

'd slpr'

'draw title slp_self_relative'

'gxprint D:\2021tj\practice2\slp_self_relative.ong white'

;

2.1)N_self.ctl

dset D:\2021tj\practice2\N_self.grd

undef 1.0E+33

title relative data

xdef 1 linear 121.88 5.650

ydef 1 linear -18.0950 1.0

zdef 1 levels 1000

tdef 15 linear 0z01Jan1978 1yr

vars 1

nsr 0 99

endvars

2.2)N_self.gs

'reinit'

'open D:\2021tj\practice2\N_self.ctl'

'set x 1'

'set y 1'

'set lev 1000'

'set t 1 15'

'd nsr'

'draw title Nino 3.4 Autocorrelation coefficient'

'gxprint D:\2021tj\practice2\Nino-self-relative.png white'

;

3.1)cross_r.ctl

dset D:\2021tj\practice2\cross_r.grd

undef 1.0E+33

title cross_r data

xdef 1 linear 1 1

ydef 1 linear -87.5 1

zdef 1 levels 1000

tdef 15 linear 0z01Jan1978 1yr

vars 1

cross 0 99

endvars

3.2)cross_r.gs

'reinit'

'open D:\2021tj\practice2\cross_r.ctl'

'set x 1'

'set y 1'

'set lev 1000'

'set t 1 15'

'd cross'

'draw title backward cross between SOI & Nino3.4'

'gxprint D:\2021tj\practice2\cross_r.png white'

;

四、图形及分析

由合成分析可知，两者关系在显著性水平α=0.05的情况下并不显著。

Nino 3.4区海温与全球海平面气压在东半球大致呈正相关，且最大值出现在（0°，160°E附近），在西半球大致呈负相关，且最小值大致出现在（0°，120°W）附近。

知秋君

相关推荐