如何做PNA遥相关型（Wallance）

啊哈

选用NCEP位势高度场资料，首先把1966.12-1977.2北半球（set lon 0 360；set lat 0 90）中所有冬季月份（12.1.2月，共45个月）500hPa位势高度资料写入newhgt.grd，然后计算点相关系数（以20°N，160°W为基点A），根据公式 ，点相关系数即计算A点与其余各点的相关系数（对吧？）。把A点这45个月的位势高度看做k变量，其余各点分别看做l变量，分别求出相关系数（Text2.f90），作图（写ctl文件，再写一个gs文件）。但是做出来的图像不对。到底是哪里出了问题，跪求指导！！感激不尽！！

Text2.f90

program dianxiangguan2

integer,parameter::nx=144,ny=37,nt=45,m1=5328

integer i,j,t,n

real tt1(nx,ny,nt),tt2(nx,ny,nt),tt3(nx,ny,nt),tt4(nx,ny,nt),tt5(nx,ny,nt),tt6(nx,ny,nt)

real ss1(nx,ny),ss2(nx,ny),fz(nx,ny),fm(nx,ny),arr(nx,ny)

real sum3(m1),mean2(m1),sum4(m1),sum5(m1),s2(m1)

real::sum1=0,sum2=0

real s1,mean1

open(unit=10,file="c:\Users\user\Desktop\newhgt.grd",form="binary",status="old")

open(unit=11,file="c:\Users\user\Desktop\newnewhgt.grd",form="binary")

read(10) tt1

!基点A（20°N，160°W）对应tt2(81,45,t)

do t=1,nt

         sum1=sum1+tt1(81,45,t)

enddo

mean1=sum1/45

do t=1,nt

         tt2(81,45,t)=tt1(81,45,t)-mean1

         tt3(81,45,t)=tt2(81,45,t)*tt2(81,45,t)

         sum2=sum2+tt3(81,45,t)

enddo

s1=sqrt(sum2/45)

!基点A时间序列的标准差s1，每时刻的距平tt2

ss1(:,:)=s1

do i=1,nx

         do j=1,ny

                   do t=1,nt

                            tt4(i,j,t)=tt2(81,45,t)

                   enddo

         enddo

enddo

n=0

sum3(:)=0

sum4(:)=0

do i=1,nx

         do j=1,ny

                   n=n+1

                   do t=1,nt

                            sum3(n)=sum3(n)+tt1(i,j,t)

                   enddo

                   mean2(n)=sum3(n)/45

                   do t=1,nt

                            tt5(i,j,t)=tt1(i,j,t)-mean2(n)

                            tt6(i,j,t)=tt5(i,j,t)*tt5(i,j,t)

                            sum4(n)=sum4(n)+tt6(i,j,t)

                   enddo

                   s2(n)=sqrt(sum4(n)/45)

         enddo

enddo

n=0

do i=1,nx

         do j=1,ny

                   n=n+1

                   ss2(i,j)=s2(n)

         enddo

enddo

sum5(:)=0

n=0

do i=1,nx

         do j=1,ny

                   n=n+1

                   do t=1,nt

                            sum5(n)=sum5(n)+tt4(i,j,t)*tt5(i,j,t)

                   enddo

                   fz(i,j)=sum5(n)/45

         enddo

enddo

do i=1,nx

         do j=1,ny

                   fm(i,j)=ss1(i,j)*ss2(i,j)

                   arr(i,j)=fz(i,j)/fm(i,j)

         enddo

enddo

write(11) arr

close(10)

close(11)

end

ctl文件：

dset C:\Users\user\Desktop\newnewhgt.grd

title PNA

undef -9.96921e+36

xdef 144 linear 0.0 2.5

ydef 37 linear 0.0 2.5

zdef 1 linear 500

tdef 1 linear Dec1966 1mo

vars 1

hgt 1 99 xiangguanxishu

endvars

Gs文件：

'reinit'

'open c:\Users\user\Desktop\PNA.ctl'

'set lon 0 360'

'set lat 0 90'

'set z 1'

'set mproj nps'

'set mpdset lowres'

'd hgt'

'enable print c:\Users\user\Desktop\PNA.gmf'

'print'

'disable print'

;

mofangbao

哎 这种问题一看就不想去思考  提问不够简洁到位  自己没有把问题理清楚就问 给的程序太长，这样只能打消活跃着回答的积极性。我不是针对楼主你一个人，只是吐槽一下毕业季各种无法让人直视的问题