谁知道计算这个计算湿位涡的程序错误在哪？

良辰

http://www.06climate.com/view/1462.html，在资料站下载的，运行时出错，哪位前辈用过，希望指点出来错误之处！万分感激！

abd

楼主都不贴出是什么错误，还让人家自己下载、运行、找错，修改？呵呵，帮楼主顶了！O(∩_∩)O~

良辰

abd 发表于 2012-9-27 21:31

                               
登录/注册后可看大图

楼主都不贴出是什么错误，还让人家自己下载、运行、找错，修改？呵呵，帮楼主顶了！O(∩_∩)O~

我的意思是用过这个程序的人肯定是知道错误之处的！那好吧，我贴程序和错误！
程序如下：
C     ***************************************************
C     *    CALCULATE PERTENCIAL VORTICITY (MPV,PV)      *
C     ***************************************************
c$large
       PROGRAM MAIN
       PARAMETER(n=31,m=26,L=21,NT=7)
       DIMENSION T(N,M,L,NT),U(N,M,L,NT),V(N,M,L,NT),VEL(N,M,L,NT),
     *       RH(N,M,L,NT),ALF(m),IP(L),tp(n,m,l,nt),d(m),TD(N,M,L,NT),
     *       QE(N,M,L,NT),SM1(N,M,L,NT),SM2(N,M,L,NT),S(N,M,L,NT),ff(m)
  real DY
c       CHARACTER*50 DS1,MS,DSS
c       CHARACTER*4 DDHH
       DATA IP/1000,975,950,925,900,850,800,750,700,650,600,550,500,
     *     450,400,350,300,250,200,150,100/
       WRITE(*,*) 'pause'
       pause
       OPEN(1,FILE='F:\taifeng\blw_data\TMPprs.dat',form='binary')
      READ(1) ((((t(i,j,k,it),i=1,n),j=1,m),k=1,l),it=1,nt)
  close (1)   
OPEN(2,FILE='F:\taifeng\blw_data\uwind.dat',form='binary')
      READ(2) ((((u(i,j,k,it),i=1,n),j=1,m),k=1,l),it=1,nt)
  close (2)   
OPEN(3,FILE='F:\taifeng\blw_data\vwind.dat',form='binary')
      READ(3) ((((v(i,j,k,it),i=1,n),j=1,m),k=1,l),it=1,nt)
  close (3)   

OPEN(4,FILE='F:\taifeng\blw_data\RHprs.dat',form='binary')
      READ(4) ((((rh(i,j,k,it),i=1,n),j=1,m),k=1,l),it=1,nt)
  close (4)   
      F1=0.0001454        ! f1=2*omege
      do j=1,m
alf(j)=(j+14)*1.0
  FF(j)=F1*SIN(3.14159*ALF(J)/180.0)    !   ff地转参数  
D(j)=6378.2*1000.*cos(alf(j)*3.14159/180)*1.0*2.0*3.14159/360.
enddo
       DY=6378.2*1000.0*1.0*2.0*3.14159/360.
c      write(*,*)  (ff(j),j=1,21)
       write(*,*)'2.    CALCULATE  AND STORE QSE,TP,QE...........'
       CALL  JQSE(T,TD,RH,IP,tp,QE,N,M,L,NT)
      
c write(*,*) ((qe(i,j,1),i=1,n),j=1,m)
     
CC------------------------------------------------------CC
       wRITE(*,*) 'pause2'
       pause
       write(*,*)'3.    CALCULATE POTENTIAL VORTICITY ...........'
  !  计算位涡S,湿位涡SM
       CALL  JMPV(TP,QE,U,V,IP,VEL,SM1,SM2,S,d,dy,ff,N,M,L,NT)
c      open (321,file='d:\0601grd\sm1.dat')
c      write(321,700) ((((sm1(i,j,k,it),i=1,n),j=1,m),k=1,l),it=1,nt)
c 700  format(1x,31f8.3)
c      close(321)  
CC------------------------------------------------------CC
       wRITE(*,*) 'pause3'
       pause
       write(*,*)'4.    SAVE MOIST POTENTIAL VORTICITY(MPV)......'
c write(*,*) (((sm(i,j,k),i=1,n),j=1,m),k=1,3)
OPEN(20,FILE='F:\taifeng\blw_data\MVP1.dat',form='binary')
      WRITE(20) ((((sm1(i,j,k,it),i=1,n),j=1,m),k=1,l),it=1,nt)
      CLOSE(20)
     
OPEN(30,FILE='F:\taifeng\blw_data\mpv2.DAT',form='binary')
      WRITE(30) ((((sm2(i,j,k,it),i=1,n),j=1,m),k=1,l),it=1,nt)
      CLOSE(30)
OPEN(40,FILE='F:\taifeng\blw_data\pv.DAT',form='binary')
      WRITE(40) ((((s(i,j,k,it),i=1,n),j=1,m),k=1,l),it=1,nt)            
      CLOSE(40)
      end
      
    !主程序结束  
      
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
       SUBROUTINE JQSE(T,TD,RH,IP,tp,QE,N,M,L,NT)
C      PARAMETER(N=21,M=29,L=7)
       DIMENSION T(N,M,L,NT),TD(N,M,L,NT),IP(L),tp(n,m,l,nt),
     *           QE(N,M,L,NT),RH(N,M,L,NT)
  DO 40 IT=1,NT
       DO 40 K=1,L
       DO 40 I=1,N
       DO 40 J=1,M
      IF(T(I,J,K,IT).gt.-15.0) THEN
c 水面处理
BA=17.2693882
      BB=35.86
else
c 冰面处理
      BA=21.3745584
      BB=7.66
endif
c由相对湿度rh用迭代法求露点温度
Es=6.11*exp(BA*T(I,J,K,IT)/(273.16-BB+T(I,J,K,IT)))
         E=RH(I,J,K,IT)*ES/1.
122        IF(E.Lt.Es) then
     T(i,j,k,it)=T(i,j,k,it)-0.05
          Es=6.11*exp(BA*T(I,J,K,IT)/(273.16-BB+T(I,J,K,IT)))
     goto 122
     else     
          TD(I,J,K,IT)=t(i,j,k,it)
     endif
c 以下注释掉的和上面的部分作用相同，得到的结果一样
c       IF(T(I,J,K).gt.-15.0) THEN
c         BA=17.2693882
c         BB=35.86
c         Es=6.11*exp(BA*T(I,J,K)/(273.16-BB+T(I,J,K)))
c         E=RH(I,J,K)*ES/100.
c 122        IF(E.Lt.Es) then
c     T(i,j,k)=T(i,j,k)-0.05
c          Es=6.11*exp(BA*T(I,J,K)/(273.16-BB+T(I,J,K)))
c     goto 122
c     else     
c          TD(I,J,K)=t(i,j,k)
c     endif
c       ELSE
c         BA=21.3745584
c         BB=7.66
c         Es=6.11*exp(BA*T(I,J,K)/(273.16-BB+T(I,J,K)))
c         E=RH(I,J,K)*ES/100.
c 222      IF(E.Lt.Es) thenc
c     T(i,j,k)=T(i,j,k)-0.05
c          Es=6.11*exp(BA*T(I,J,K)/(273.16-BB+T(I,J,K)))
c     goto 222
c     else     
c          TD(I,J,K)=t(i,j,k)
c     endif
c       ENDIF
40     CONTINUE
       DO 920 IT=1,NT
       DO 920 K=1,L
       DO 920 I=1,N
       DO 920 J=1,M
       IF (TD(I,J,K,IT).GT.-15.) THEN
       QP1=6.11*EXP(17.26*TD(I,J,K,IT)/(TD(I,J,K,IT)+237.3)) !273.16-35.86=237.3
       Q=0.622*QP1/(iP(K)-0.378*QP1)               !QP1为水汽压
                                               !Q为比湿
       ELSE
       QP1=6.11*EXP(21.87*TD(I,J,K,IT)/(TD(I,J,K,IT)+265.5))  !273.16-7.66=265.5
       Q=0.622*QP1/(iP(K)-0.378*QP1)
       ENDIF
       TT=273.16+T(I,J,K,IT)
  TTD=273.16+TD(i,j,k,it)
       TP(i,j,k,it)=TT*((1000.0/iP(K))**0.286)              !tp是位温
TL=0.622*(597.3-0.566*t(i,j,k,it))*TTD
     */(0.622*(597.3-0.566*t(i,j,k,it))+0.24*TD(i,j,k,it)*log(TT/TTD))
c       TL=338.72-0.24*TT+1.24*Td(i,j,k,it)                !tl是凝结高度温度
       QE(I,J,K,IT)=TP(I,J,K,IT)*EXP((597.3-0.566*t(i,j,k,it))
     *  *Q/(0.2403*TL))                               ! qe是相当位温
                                                         
c      QSE(I,J,K)=TT*EXP(0.28586*ALOG(1000./iP(K))+2500.*Q/TL)
920    CONTINUE
       RETURN
       END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
       SUBROUTINE JMPV(TP,QE,U,V,IP,VEL,SM1,SM2,S,d,dy,ff,N,M,L,NT)
       DIMENSION U(N,M,L,NT),V(N,M,L,NT),QE(N,M,L,NT),TP(N,M,L,NT),
     *           IP(L),SM1(N,M,L,NT),SM2(N,M,L,NT),S(N,M,L,NT),
     *           d(m),ff(m),VEL(M,N,L,NT)
  real DY
       DO 168 IT=1,NT
       DO 168 K=1,L
       IF(K.EQ.1) THEN
         KK1=K+1
         KK0=K
         DP=IP(KK1)-IP(KK0)       ! 前差
       ELSE IF(K.EQ.L) THEN
         KK1=K
         KK0=K-1
         DP=IP(KK1)-IP(KK0)       ! 后差
       ELSE
         KK1=K+1
         KK0=K-1
         DP=(IP(KK1)-IP(KK0))*2   ! 中央差
       ENDIF
       DO 168 I=2,n-1
       DO 168 J=2,m-1
c       F=F1*SIN(3.14159*ALF(I)/180.)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!   请完成.............
  VEL(I,J,K,IT)=(V(I+1,J,K,IT)-V(I-1,J,K,IT))/(2*D(J))
     &         -(U(I,J+1,K,IT)-U(I,J-1,K,IT))/(2*DY)      ! VEL涡度
      
       S(I,J,K,IT)=-9.8*(FF(J)+VEL(I,J,K,IT))*(TP(I,J,KK1,IT)
     &   -TP(I,J,KK0,IT))/DP +9.8*((V(I,J,KK1,IT)-V(I,J,KK0,IT))/DP)*
     &    ((TP(I+1,J,K,IT)-TP(I-1,J,K,IT))/(2*D(J)))-9.8*((U(I,J,KK1,IT)
     &    -U(I,J,KK0,IT))/DP)*((TP(I,J+1,K,IT)-TP(I,J-1,K,IT))/(2*DY))
  SM1(I,J,K,IT)=-9.8*(FF(J)+VEL(I,J,K,IT))*
     &             ((QE(I,J,KK1,IT)-QE(I,J,KK0,IT))/DP)
  SM2(I,J,K,IT)=9.8*(((V(I,J,KK1,IT)-V(I,J,KK0,IT))/DP)*
     &    (QE(I+1,J,K,IT)-QE(I-1,J,K,IT))/(2*D(J))-((U(I,J,KK1,IT)
     &    -U(I,J,KK0,IT))/DP)*(QE(I,J+1,K,IT)-QE(I,J-1,K,IT))/(2*DY))

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
168    CONTINUE
cccccc计算边界值     
       DO 169 IT=1,NT
       DO 169 K=1,L
       IF(K.EQ.1) THEN
         KK1=K+1
         KK0=K
         DP=IP(KK1)-IP(KK0)       ! 前差
       ELSE IF(K.EQ.L) THEN
         KK1=K
         KK0=K-1
         DP=IP(KK1)-IP(KK0)       ! 后差
       ELSE
         KK1=K+1
         KK0=K-1
         DP=(IP(KK1)-IP(KK0))*2   ! 中央差
       ENDIF
c       F=F1*SIN(3.14159*ALF(I)/180.)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!   请完成.............
  VEL(1,1,K,IT)=0.0    ! VEL涡度
  S(1,1,K,IT)=0.0
  SM1(1,1,K,IT)=0.0
  SM2(1,1,K,IT)=0.0
  VEL(1,m,K,IT)=0.0          ! VEL涡度
  S(1,m,K,IT)=0.0
  SM1(1,m,K,IT)=0.0
  SM2(1,m,K,IT)=0.0
       VEL(n,1,K,IT)=0.0          ! VEL涡度
  S(n,1,K,IT)=0.0
  SM1(n,1,K,IT)=0.0
  SM2(n,1,K,IT)=0.0
       VEL(n,m,K,IT)=0.0          ! VEL涡度
  S(n,m,K,IT)=0.0
  SM1(n,m,K,IT)=0.0
  SM2(n,m,K,IT)=0.0


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
169    CONTINUE
       RETURN
       END
个人感觉是  VEL(I,J,K,IT)=(V(I+1,J,K,IT)-V(I-1,J,K,IT))/(2*D(J))
     &         -(U(I,J+1,K,IT)-U(I,J-1,K,IT))/(2*DY)      ! VEL涡度是这有问题，维数越界了吧，I定义了，这里出现了I+1等等。请用过的人或大神指点！
错误如图：

mofangbao

楼主可以在程序里多加print*的提示，看看是运行到哪里开始出错的，然后针对出错的部分仔细分析，看提示是数组越界，但是实际上可能是别的原因引起的

abd

我也一般用管理员说的方法，每一步用print输出看看数据，一步一步确认，楼主加油！

lgg00828

良辰 发表于 2012-9-28 00:11
我的意思是用过这个程序的人肯定是知道错误之处的！那好吧，我贴程序和错误！
程序如下：
C     **** ...

我也遇到这个问题了，m取值在20以下就没问题，超过20就出问题了，说是exp函数的问题，我也不懂，不知道楼主解决了吗

north23

看看啊           

读书一摞摞

同被湿位涡折磨