为啥输出的结果都是0呢，急急急，请大佬帮忙看看，老师也不知道错哪

sususir

program main
implicit none
integer,parameter::nx=53,ny=29,nz=11,nt=22
real::a=6371000.
integer i,j,iz,it
character*12 timefile(nt)
character*4 levelfile(nz)
integer p(nz)
real::m=17.2693882,n=35.86
real::dx=2.5*2*3.1416/360,dy=2.5*2*3.1416/360
real vort(nx,ny,nz,nt),u(nx,ny,nz,nt),v(nx,ny,nz,nt),lat(ny),div(nx,ny,nz,nt),td(nx,ny,nz,nt),e(nx,ny,nz,nt),q(nx,ny,nz,nt)
real vor_ad(nx,ny,nz,nt),tem(nx,ny,nz,nt),tem_ad(nx,ny,nz,nt),ttd(nx,ny,nz,nt),qu(nx,ny,nz,nt),qv(nx,ny,nz,nt)
data p /1000,925,850,700,500,400,300,250,200,150,100/
data levelfile/'1000','925','850','700','500','400','300','250','200','150','100'/
data timefile/'13052108.000','13052120.000','13052208.000','13052220.000','13052308.000',&
     '13052320.000','13052408.000','13052420.000','13052508.000','13052520.000','13052608.000','13052620.000','13052708.000',&
       '13052720.000','13052808.000','13052820.000','13052908.000','13052920.000','13053008.000','13053020.000','13053108.000','13053120.000'/
!读风速u和v
do iz=1,nz
   do it=1,nt
    open(1,file='D:\micaps/uv/'//trim(levelfile(iz))//'/'//timefile(it))
        do i=1,3
     read(1,*)
    end do
    do j=ny,1,-1
     read(1,*) (u(i,j,iz,it),i=1,nx)
    end do
    do j=ny,1,-1
     read(1,*) (v(i,j,iz,it),i=1,nx)
    end do
        close(1)
   end do
end do
!读温度
do iz=1,nz
   do it=1,nt
    open(3,file='D:\micaps/temper/'//trim(levelfile(iz))//'/'//timefile(it))
        do i=1,4
         read(3,*)
        end do
    do j=ny,1,-1
         read(3,*) (tem(i,j,iz,it),i=1,nx)
        end do
        close(3)
   end do
end do
!读温度露点差
do iz=1,nz
   do it=1,nt
    open(4,file='D:\micaps/t-td/'//trim(levelfile(iz))//'/'//timefile(it))
        do i=1,4
         read(4,*)
        end do
    do j=ny,1,-1
         read(4,*) (ttd(i,j,iz,it),i=1,nx)
        end do
        close(4)
   end do
end do
open(2,file='D:\micaps\all.grd',form='binary')
!计算涡度
do it=1,nt
   do iz=3,7,2
     do j=2,ny-1
           do i=2,nx-1
        vort(i,j,iz,it)=1./(2*a)*((v(i+1,j,iz,it)-v(i-1,j,iz,it))/(dx)/cos(lat(j))-(u(i,j+1,iz,it)&
                -u(i,j-1,iz,it))/dy+2*u(i,j,iz,it)*tan(lat(j)))
           end do
         end do
   end do
end do
!计算散度
do it=1,nt
   do iz=3,7,2
     do j=2,ny-1
           lat(j)=(10+(j-1))*dx
           do i=2,nx-1
                div(i,j,iz,it)=1./(2*a)*((v(i+1,j,iz,it)-v(i-1,j,iz,it))/(dx)/cos(lat(j))-(u(i,j+1,iz,it)&
                -u(i,j-1,iz,it))/dy-2*v(i,j,iz,it)*tan(lat(j)))
           end do
         end do
   end do
end do
!计算温度平流
dO it=1,nt
  dO iz=3,3
    dO j=2,ny-1
         lat(j)=(10+(j-1))*dx
     dO i=2,nx-1
       tem_ad(i,j,iz,it)=1./(2*a)*(-u(i,j,iz,it)*(tem(i+1,j,iz,it)-tem(i-1,j,iz,it))/dx/cos(lat(j))&
           -v(i,j,iz,it)*(tem(i,j+1,iz,it)-tem(i,j-1,iz,it))/dy)
            end do
        end do
  end do
end do
print*,10000000*tem_ad(2,2,3,2)
!计算涡度平流
do it=1,nt
   do iz=5,5
     do j=3,ny-2
           lat(j)=(10+(j-1)*dx)
           do i=3,nx-2
            vor_ad(i,j,iz,it)=1./(2*a)*(-u(i,j,iz,it)*(vort(i+1,j,iz,it)-vort(i-1,j,iz,it))/dx/cos(lat(j))&
                -v(i,j,iz,it)*(vort(i,j+1,iz,it)-vort(i,j-1,iz,it))/dy)
                end do
          end do
        end do
end do
!计算比湿
do it=1,nt
do iz=3,5
  do j=1,ny
   do i=1,nz
    td(i,j,iz,it)=tem(i,j,iz,it)-ttd(i,j,iz,it)
        e(i,j,iz,it)=6.1078*EXP(m*td(i,j,iz,it)/(273.15+td(i,j,iz,it)-n))
        q(i,j,iz,it)=0.622*e(i,j,iz,it)/(iz-0.378*e(i,j,iz,it))
   end do
  end do
end do
end do
!计算水汽通量和水汽通量散度
do it=1,nt
do iz=3,5
   do j=1,ny
    do i=1,nx
         qu(i,j,iz,it)=q(i,j,iz,it)*u(i,j,iz,it)/9.8
     qv(i,j,iz,it)=q(i,j,iz,it)*v(i,j,iz,it)/9.8
        end do
   end do
  end do
end do
!输出结果
open(2,file='D:\micaps\all.grd',form='binary')
  do it=1,nt
    write(2) (((tem_ad(i,j,iz,it),i=2,nx-1),j=2,ny-1),iz=3,3)
  end do
open(21,file='D:\micaps\temper_ad.txt')
open(31,file='D:\micaps\vort.txt')
open(41,file='D:\micaps\div.txt')
open(51,file='D:\micaps\vor_ad.txt')
do it=1,nt
   write(21,*) (((tem_ad(i,j,iz,it),i=2,nx-1),j=2,ny-1),iz=3,3)
  end do
  do it=1,nt
   write(31,*) (((vort(i,j,iz,it),i=2,nx-1),j=2,ny-1),iz=3,7,2)
  end do
  do it=1,nt
   write(41,*) (((div(i,j,iz,it),i=2,nx-1),j=2,ny-1),iz=3,7,2)
  end do
  do it=1,nt
   write(51,*) (((vor_ad(i,j,iz,it),i=3,nx-2),j=3,ny-2),iz=5,5)
  end do
close(2)
close(21)
close(31)
close(41)
close(51)
end

sususir

哈哈哈，我也不知道改了什么，有结果啦