请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 6696|回复: 0

fortran计算出的流函数只有边界值有值

[复制链接]

新浪微博达人勋

发表于 2017-12-8 14:48:01 | 显示全部楼层 |阅读模式
1金钱
我用第二边界条件,和利布曼法计算流函数,但是我流函数计算出来,只有边界值有值,中间值为0,想了很久了,实在想不出来,求大家帮帮忙,谢谢大家
流函数的边条件

     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



    end

密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表