mofangbao 发表于 2012-5-9 09:39
能把你求相关的脚本贴出来吗,我这里也有海温资料,可以试试~
parameter(nx=181,ny=89,nz=1,n=52,undef=32767)
dimension ty(n),tw(nx,ny,nz,n)
real ya,wa(nx,ny,nz),syw(nx,ny,nz),sy,sw(nx,ny,nz),r(nx,ny,nz)
open(1,file='i.txt')
read(1,*) (ty(it),it=1,n)
close(1)
open(2,file='sst.6.grd',form='binary')
read(2) ((((tw(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz),it=1,n)
close(2)
do iz=1,nz
do i=1,nx
do j=1,ny
ya=0.0
wa(i,j,iz)=0.0
do it=1,n
ya=ya+ty(it)
wa(i,j,iz)=wa(i,j,iz)+tw(i,j,iz,it)
end do
ya=ya/real(n)
wa(i,j,iz)=wa(i,j,iz)/real(n)
syw(i,j,iz)=0.0
sy=0.0
sw(i,j,iz)=0.0
do it=1,n
syw(i,j,iz)=syw(i,j,iz)+(ty(it)-ya)*(tw(i,j,iz,it)-wa(i,j,iz))
sy=sy+(ty(it)-ya)*(ty(it)-ya)
sw(i,j,iz)=sw(i,j,iz)+(tw(i,j,iz,it)-wa(i,j,iz))*(tw(i,j,iz,it)-wa(i,j,iz))
end do
syw(i,j,iz)=syw(i,j,iz)/real(n)
sy=sy/real(n)
sw(i,j,iz)=sw(i,j,iz)/real(n)
r(i,j,iz)=0.0
r(i,j,iz)=syw(i,j,iz)/sqrt(sy*sw(i,j,iz))
enddo;enddo;enddo
do it=1,n
do iz=1,nz
do i=1,nx
do j=1,ny
if(tw(i,j,iz,it)==undef)then
r(i,j,iz)=0.0
endif
enddo;enddo;enddo;enddo
|