DO k=1,NZ ! level
DO l=1,NT ! time
open(1,file='d:\micaps/height/'//trim(levelfile(k))//'/'//timefile(l))
do ii=1,4
read(1,*)
enddo
do j=ny,1,-1
read(1,*)(h(i,j,k,l),i=1,NX)
enddo
enddo
enddo
close(1)
open(2,file='d:/micaps/5/height.grd',form='binary')
dO l=1,nt
write(2)(((h(i,j,k,l),i=1,nx),j=1,ny),k=1,nZ)
enddo
close(2)
open(2,file='d:/micaps/5/height.grd',form='binary')
do l=1,nt
read(2) (((h(i,j,k,l),i=1,nx),j=1,ny),k=1,nz)
enddo
h=h*10
!边界上的liu值
liu=0
do j=1,ny
F=Omega*2*sin(lat(j))
enddo
do l=1,nt
do k=1,nz
do i=1,nx
! j=ny上的liu
liu(i,1,k,l)=h(i,1,k,l)/(Omega*2*sin(lat(1)))
liu(i,ny,k,l)=h(i,ny,k,l)/(Omega*2*sin(lat(ny)))
enddo
enddo
enddo
do l=1,nt
do k=1,nz
do j=1,ny
! i=1.and.i=nx上的liu
liu(1,j,k,l)=h(1,j,k,l)/F
liu(nx,j,k,l)=h(nx,j,k,l)/F
enddo
enddo
enddo
!-------------------计算流函数--------------
n=0
RR=0.0
RRmax=1.e-5
!迭代次数
do l=1,nt
do k=1,nz
do while(RRmax(1,1,k,l)>1.e-7)
do j=1+1,ny-1
do i=1+1,nx-1
RR(i,j,k,l)=(liu(i+1,j,k,l)+liu(i-1,j,k,l))/((R*delta*pi/180*cos(lat(j)))**2)+(liu(i,j+1,k,l)+liu(i,j-1,k,l))/((R*delta*pi/180)**2)-(2.0/(R*delta*pi/180*cos(lat(j)))**2+2.0/(a*delta*pi/180)**2)*liu(i,j,k,l)-voro(i,j,k,l)
enddo
enddo
do j=1+2,ny-2
do i=1+2,nx-2
liu(i,j,k,l)=liu(i,j,k,l)+RR(i,j,k,l)/(2./((R*delta*pi/180*cos(lat(j)))**2)+2./((R*delta*pi/180)**2))
enddo
enddo
RRmax=abs(RR)
do j=1+2,ny-2
do i=1+2,nx-2
if(RRmax(1,1,k,l)<=RRmax(i,j,k,l))then
RRmax(1,1,k,l)=abs(RRmax(i,j,k,l))
endif
enddo
enddo
n=n+1
enddo
enddo
enddo
open(22,file='d:/micaps/5/liu.grd',form='binary')
dO l=1,nt
write(22)(((liu(i,j,k,l),i=1,nx),j=1,ny),k=1,nZ)
enddo
open(23,file='d:/micaps/5/woro.txt')
do j=1,ny
write(23,*)(liu(i,j,5,20),i=1,nx)
enddo
print*,n