登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
program dd
parameter(nx=144,ny=73,nz=21,nt=60)
real::h(nx,ny,nt),a,b,c,d,mmx(nx,nt),gs1(nx,nt),gn1(nx,nt),gs2(nx,nt),gn2(nx,nt),gs3(nx,nt),gn3(nx,nt),index1(nx,nt),index2(nx,nt),index3(nx,nt)
integer i,j,k
external smax
!!!读取500hPa的数据
open(20,file='G:\coldwave\blockingindex/500hgt.grd',form='binary')
do k = 1,nt
do j = 1,ny
read(20)(h(i,j,k),i=1,nx)
enddo
enddo
print*,h
do t = 1,nt
do x = 1,nx
!!!计算北指数和南指数,三个扰动值都计算了
d = 0
a = 68 + d
b = 60 + d
c = 52 + d
gs1(x,t) = (h(x,b,t)-h(x,c,t))/20
gn1(x,t) = (h(x,a,t)-h(x,b,t))/20
d = -2
a = 68 + d
b = 60 + d
c = 52 + d
gs2(x,t) = (h(x,b,t)-h(x,c,t))/20
gn2(x,t) = (h(x,a,t)-h(x,b,t))/20
d = 2
a = 68 + d
b = 60 + d
c = 52 + d
gs3(x,t) = (h(x,b,t)-h(x,c,t))/20
gn3(x,t) = (h(x,a,t)-h(x,b,t))/20
!!将符合不同时满足两个条件的三个指数都取为0,这样处理缺测值可以设置为0
if(gs1(x,t) > 0 .and. gn1(x,t) < -10)then
index1(x,t) = gs1(x,t)
else
index1(x,t) = 0
endif
if(gs2(x,t) > 0 .and. gn2(x,t) < -10)then
index2(x,t) = gs2(x,t)
else
index2(x,t) = 0
endif
if(gs3(x,t) > 0 .and. gn3(x,t) < -10)then
index3(x,t) = gs3(x,t)
else
index3(x,t) = 0
endif
!!!!!=====================输出三个指数最大的一个
mmx(x,t) = smax(index1(x,t),index2(x,t),index3(x,t))
enddo !对应经度循环
enddo !对应时间循环
!print*,mmx
open(22,file='G:\coldwave\blockingindex/2018.HGT.blocking.txt',status='replace')
open(44,file='G:\coldwave\blockingindex/2018.HGT.blocking.grd',status='replace',form='binary')
do t=1,nt
do x=1,nx
write(22,*)mmx(x,t)
write(44)mmx(x,t)
enddo
enddo
close(22)
close(44)
end program
!!三个数,返回其中最大值的外部子程序
function smax(a,b,c)
real a,b,c
if(a > b)then
if(a>c)then
e = a
else
e = c
endif
else
if(b>c)then
e = b
else
e = c
endif
endif
smax = e
end function smax
|