E-P通量的Fortran程序

dovely

先收藏一下~~~~

fancyyrl

谢谢楼主分享

★翼★

我计算ep2的时候，总是出现 Infinity       -Infinity的情况，我读取的数据是grads转化来的二进制，具体程序如下：
      program epflux2
      parameter(nx=810,ny=181,nz=20,nt=1)
      real t(nx,ny,nz,nt),v(nx,ny,nz,nt),avev(ny,nz,nt),avet(ny,nz,nt),
     &avevt(ny,nz,nt),p0(nz),t0(nz,nt),fnn(nz,nt),aveuv(ny,nz,nt)
      data p0/1000.0,975.0,950.0,925.0,900.0,850.0,800.0,750.0,700.0,
     *            650.0,600.0,550.0,500.0,450.0,400.0,350.0,300.0,250.0,
     *        200.0,150.0/
          p0=p0*100
     
      a=6371000        !!!!!地球半径
      omega=7.292e-5   !!!地球旋转角速度
      g=9.80665
      pi=3.1416
      R=287.           !!!气体常数

      TS=240.
      ps=1e5         !!!表面气压


      open (1,file='e:\t.dat',access='direct',form='unformatted',
     &recl=nx*ny)
      m=0
        do kt=1,nt
      do l=1,nz
        m=m+1
      read(1,rec=m,iostat=ferr) ((t(i,j,l,kt),i=1,nx),j=1,ny)
      enddo
      enddo
      close(1)

      open (2,file='e:\v.dat',access='direct',form='unformatted',
     &recl=nx*ny)
      m=0
        do kt=1,nt
      do l=1,nz
        m=m+1
      read(1,rec=m,iostat=ferr) ((v(i,j,l,kt),i=1,nx),j=1,ny)
      enddo
      enddo
      close(2)

      open (3,file='e:\avet.dat',access='direct',form='unformatted',
     &recl=ny)
      m=0
        do kt=1,nt
      do l=1,nz
        m=m+1
      read(3,rec=m,iostat=ferr) (avet(j,l,kt),j=1,ny)
      enddo
      enddo
      close(3)

      open (4,file='e:\avev.dat',access='direct',form='unformatted',
     &recl=ny)
      m=0
        do kt=1,nt
      do l=1,nz
        m=m+1
      read(4,rec=m,iostat=ferr) (avev(j,l,kt),j=1,ny)
      enddo
      enddo
      close(4)

cccccccccccccccccccccccccccccccccccccc
      do n=1,nt
      do k=1,nz
       t0(k,n)=0.0
       do j=1,ny
       do i=1,nx
        t0(k,n)=t0(k,n)+t(i,j,k,n)*(ps/p0(k))**0.286/(nx*ny)
       end do
       end do
      end do
      end do

      do n=1,nt

      do k=2,nz-1
      dx=p0(k+1)-p0(k-1)
      fnn(k,n)=(t0(k+1,n)-t0(k-1,n))/dx
      enddo
      dx=p0(2)-p0(1)
      fnn(1,n)=(t0(2,n)-t0(1,n))/dx
      dx=p0(nz)-p0(nz-1)
      fnn(nz,n)=(t0(nz,n)-t0(nz-1,n))/dx
   
      enddo

c      write(*,*) avev
ccccccccccccccccccccccccccccccccccccccccc
     
      do kt=1,nt
      do l=1,nz
      do j=1,ny
      zz=0.0
      do i=1,nx
      avevt(j,l,kt)=(v(i,j,l,kt)-avev(j,l,kt))*(t(i,j,l,kt)-avet(j,l,kt)
     &)*(ps/p0(l))**0.286+zz
      zz=avevt(j,l,kt)
      enddo
      enddo
      enddo
      enddo
cccccccccccccccccccccccccccccccccccccccc


      open (5,file='e:\ep2.txt',recl=8)

      do kt=1,nt
      do l=1,nz
      do j=1,ny
      f=(2.0*omega*sin((-90+(j-1)*1.0)/180*3.14))!!!科氏参数f
      write(5,*) avevt(j,l,kt)/nx*cos((-90+(j-1)*1.0)/180*3.14)*f*a/
     &fnn(l,kt)
      enddo
      enddo
      enddo
     
      end
具体错误出在哪儿了，请楼主指教

liuchunchun

谢谢分享 非常有用

jolincai

EP1和EP2是什么啊？

wuyao

学习一下，谢谢楼主

聆念LN

你好楼主，我也想要波通量的for程序，可不可以也给我一份，多谢多谢！
596776754@qq.com

signor9

好贵，下下来好好学学

猫猫monty

cold_wq 发表于 2013-9-23 16:39
我有那个For程序，稍后传给你，你留个邮箱地址，与E-P通量公式和主要用途上还是有本质区别的。

您好，现在十分需要三维波活动通量的程序，可以分享吗？邮箱：997538345@qq.com

simba_sy

谢谢 下下来学习下