求教，计算水汽通量和水汽通量散度的程序有没有问题

1099221723

program main
implicit none
integer,parameter::nx=33,ny=61,nt=720,nz=8
integer i,j,z,t,irec
real,parameter::Delta=2.5,PI=3.1415926,R=6371e+3
real lon(ny),lat(nx),u(nx,ny,nz,nt),v(nx,ny,nz,nt),q(nx,ny,nz,nt),qu(nx,ny,nz,nt),qv(nx,ny,nz,nt),adq(nx,ny,nz,nt),adqv(nx,ny,nz,nt),div(nx,ny,nz,nt),add(nx,ny,nz,nt)
real fqu(nx,ny),fqv(nx,ny),fadqv(nx,ny),sumu(nx,ny),sumv(nx,ny),sumadqv(nx,ny)
!qu: the vapor flux in the u direction
!qv: the vapor flux in the v direction
!adq: specific humidity advection
!adqv: vapor flux divergence
!div: divergence
!add: the divergence of wind
!经纬度转换成弧度
do i=1,ny
lon(i)=30+(i-1)*Delta
lon(i)=lon(i)*pi/180
end do
do i=1,nx
lat(i)=-20+(i-1)*Delta
lat(i)=lat(i)*pi/180
end do
!读入数据
open(1,file='D:\lunwen\uwnd.grd',form='unformatted',access='direct',recl=nx*ny)
do t=1,nt
do z=1,nz
irec=1
read(1,rec=irec) ((u(i,j,z,t),i=1,nx),j=1,ny)
irec=irec+1
end do
end do
close(1)
open(2,file='D:\lunwen\vwnd.grd',form='unformatted',access='direct',recl=nx*ny)
do t=1,nt
do z=1,nz
irec=1
read(2,rec=irec) ((v(i,j,z,t),i=1,nx),j=1,ny)
irec=irec+1
end do
end do
close(2)
open(3,file='D:\lunwen\shum.grd',form='unformatted',access='direct',recl=nx*ny)
do t=1,nt
do z=1,nz
irec=1
read(3,rec=irec) ((q(i,j,z,t),i=1,nx),j=1,ny)
irec=irec+1
end do
end do
close(3)
!计算各层（8层）的散度
do t=1,nt
do z=1,nz
do j=2,ny-1
do i=2,nx-1
div(i,j,z,t)=(u(i+1,j,z,t)-u(i-1,j,z,t))/(2*R*Delta*PI/180*cos(lat(i)))+(v(i,j+1,z,t)-v(i,j-1,z,t))/(2*R*Delta*PI/180)
enddo
enddo
enddo
enddo
!计算各层（8层）水汽通量以及水汽通量散度
do t=1,nt
do z=1,nz
do j=2,ny-1
do i=2,nx-1
         qv(i,j,z,t)=q(i,j,z,t)*v(i,j,z,t)/9.8
     qu(i,j,z,t)=q(i,j,z,t)*u(i,j,z,t)/9.8
         adq(i,j,z,t)=(u(i,j,z,t)*((q(i+1,j,z,t)-q(i-1,j,z,t))/(2*R*cos(lat(i))*Delta*PI/180))+v(i,j,z,t)*((q(i,j+1,z,t)-q(i,j-1,z,t))/(2*R*Delta*PI/180)))/9.8
         add(i,j,z,t)=(q(i,j,z,t)*div(i,j,z,t))/9.8
         adqv(i,j,z,t)=adq(i,j,z,t)+add(i,j,z,t)
enddo
enddo
enddo
enddo
!输出数据
open(4,file='d:\lunwen\out.grd',form='binary')
do t=1,nt
do z=1,nz
write(4) ((qu(i,j,z,t),i=2,nx-1),j=2,ny-1)
end do
end do
do t=1,nt
do z=1,nz
write(4) ((qv(i,j,z,t),i=2,nx-1),j=2,ny-1)
end do
end do
do t=1,nt
do z=1,nz
write(4) ((adqv(i,j,z,t),i=2,nx-1),j=2,ny-1)
end do
end do
do t=1,nt
do z=1,nz
write(4) ((div(i,j,z,t),i=2,nx-1),j=2,ny-1)
end do
end do
close(4)
end

╯①路霺笑

按理说没错啊

排山倒海

你的4个输出变量div，qu，qv，adqv定义的时候都是：i从1到nx，j从1到ny，但是没有赋初值，虽然后面直接用的是从2到nx-1（ny-1），但是他们的初值还是没有固定，系统容易出错，建议开始定义的时候直接赋初值为0

pcc5201314

看你写的那么多，就知道没错

longlong

{:eb513:}

Mr.Yes

{:eb502:}

tuini

思路很清晰，就是维数设定有点多，大数据不好实现！

windys

我来学习一下！

piloiohe

小白前来学习，还不会怎么写，多谢楼主