APO指数年际变化计算

浪子回头

      PROGRAM ebna
implicit real(a-z)
dimension tt(144,73,32),tp(73,32),tr(144,73,32),apo(32)
      dimension ap(32),bz(32),tr1(32),tr2(32)
real fc,a,ave
      open(10,file='d:\xnlw\air2.grd',form='binary',status='old')

      !读数据
do it=1,32
do j=1,73
do i=1,144
read(10)tt(i,j,it)
enddo;enddo;enddo
      write(*,*)'read data ok'
      !纬向平均
do j=1,73
      do it=1,32
      tp(j,it)=0.0
enddo
enddo
      do it=1,32
      do j=1,73
      do i=1,144
      tp(j,it)=tp(j,it)+tt(i,j,it)/144
      enddo;enddo;enddo
print*,tp
      !扰动值
      do it=1,32
do j=1,73
do i=1,144
tr(i,j,it)=tt(i,j,it)-tp(j,it)
enddo;enddo;enddo

       write(*,*)'select region ok111'
       !（所选区域15-50N,180-120W,60E-120E)，区域平均
do it=1,32
tr1(it)=0
enddo
      do it=1,32
do j=43,57
do i=25,49
        tr1(it)=tr1(it)+tr(i,j,it)/375
      enddo
enddo
enddo
do it=1,32
tr2(it)=0
enddo
      do it=1,32
do j=43,57
do i=73,97
        tr2(it)=tr2(it)+tr(i,j,it)/375
      enddo
enddo
enddo
       !区域作差值
      do it=1,32
   apo(it)=tr1(it)-tr2(it)
enddo
print*,apo
      write(*,*)'select region ok222'
!标准化处理
ave=0.0
do it=1,32
   ave=ave+apo(it)/32
enddo
do it=1,32
   ap(it)=apo(it)-ave
enddo
print*,ap
a=0.0
do it=1,32
  a=a+ap(it)**2
enddo
fc=sqrt(a/32)
print*,fc
      write(*,*)'select region ok333'
do it=1,32
bz(it)=ap(it)/fc
enddo
print*,bz
      write(*,*)'select region ok444'
      open(20,file='d:\xnlw\apo.grd',form='binary')
write(20) (apo(it),it=1,32)
open(30,file='d:\xnlw\bz.grd',form='binary')
write(30) (bz(it),it=1,32)
close(20)
close(30)
end

ayzqs

顶，谢谢分享！

言深深

楼主有时间的时候不妨给这些程序加个文件头，说明一下哈

浪子回头

言深深 发表于 2012-12-17 08:23

                               
登录/注册后可看大图

楼主有时间的时候不妨给这些程序加个文件头，说明一下哈

好的，以后会注意的

yh1223

楼主，你是对APO有研究过吗？如果有的话，我想请教你些问题！

mrxiaoan

APO不是指500-200hPa的T吗？那应该有好几层气温吧。但air2.grd里读出来的好像只有一层格点上的多年气温数据，air2.grd是垂直积分后的结果吗？

hzgloria

对啊！楼主是否测试过正确吗？

小蚂蚁

多谢楼主分享