run-time error M6201 : MATH 错误，求助大神！！！！

大力

整了好几天，好不容易编好一个程序，到了最后一步，出错！！又找不出来错误，整个人快崩溃了，求助各路大神，帮帮忙！！！具体程序如下：
program main
        implicit none
!   integer,parameter :: imax=1440, jmax=721,kmax=31
        integer,parameter :: imax=500, jmax=200,kmax=31
        real,allocatable:: HGTprs(:,:,:),TMPprs(:,:,:),RHprs(:,:,:),TD(:,:,:),p(:,:,:)
        real,allocatable:: si(:,:),t850(:,:),t500(:,:),pc(:,:),tc(:,:),thse(:,:),rh(:,:),td850(:,:),pp(:,:)
        integer i,j,k
        character*255 filetemp


        allocate(HGTprs(imax,jmax,kmax))
        allocate(TMPprs(imax,jmax,kmax))
        allocate(RHprs(imax,jmax,kmax))
        allocate(TD(imax,jmax,kmax))
        allocate(p(imax,jmax,kmax))       

        allocate(si(imax,jmax))
        allocate(t850(imax,jmax))
        allocate(t500(imax,jmax))
        allocate(pc(imax,jmax))
        allocate(tc(imax,jmax))       
        allocate(thse(imax,jmax))       
        allocate(rh(imax,jmax))
        allocate(td850(imax,jmax))
        allocate(pp(imax,jmax))
               
!        filetemp='D:\data\analysis\20170802\gdas1.fnl0p25.2017080100.f00.grib2'
!        open(21,file=filetemp,form='unformatted',access='direct',recl=imax*jmax*kmax*4,err=1005)
!!        read(21,rec=1,err=1005)(((HGTprs(i,j,k),i=1,imax),j=1,jmax),k=1,kmax)
!        close(21)
!        1005 continue
!!!!!!!!!!温度T       
        filetemp='D:\data\analysis\20170802\gdas1.fnl0p25.2017080100.f00.grib2'
        open(21,file=filetemp,form='unformatted',access='direct',recl=imax*jmax*kmax*4,err=1007)
        read(21,rec=1,err=1007)(((TMPprs(i,j,k),i=1,imax),j=1,jmax),k=1,kmax)
        close(21)
        1007 continue
!!!!!!!!!!相对湿度       
        filetemp='D:\data\analysis\20170802\gdas1.fnl0p25.2017080100.f00.grib2'
        open(21,file=filetemp,form='unformatted',access='direct',recl=imax*jmax*kmax*4,err=1009)
        read(21,rec=1,err=1009)(((RHprs(i,j,k),i=1,imax),j=1,jmax),k=1,kmax)
        close(21)
        1009 continue
       
        do j=1,jmax
        do i=1,imax
        k=6
        t850(i,j)=TMPprs(i,j,k)-273.15
        rh(i,j)=RHprs(i,j,k)
        enddo
        enddo

        do j=1,jmax
        do i=1,imax
        k=13       
        t500(i,j)=TMPprs(i,j,k)-273.15
        enddo
        enddo

        do j=1,jmax
        do i=1,imax
        pp(i,j)=850
        enddo
        enddo

        do j=1,jmax
        do i=1,imax
        if(rh(i,j).gt.0)then
        td850=t850-((14.55+0.114*t850)*(1-0.01*rh)+(2.5+0.007*t850)*(1-0.01*rh)**3+(15.9+0.37*t850)*(1-0.01*rh)**14)
        endif
        enddo
        enddo

        call cal_tc_pc(pp,t850,td850,pc,tc,imax,jmax,kmax,p)
    call cal_thse(pc,tc,thse,imax,jmax)
        call cal_ta_si(t850,pc,tc,thse,t500,si,imax,jmax,kmax,p)

        filetemp='D:\data\analysis\20170802\2017080100.dat'
        open(27,file=filetemp,form='unformatted',access='direct',recl=imax*jmax*kmax*4)
        write(27,rec=1)((si(i,j),i=1,imax),j=1,jmax)
        close(27)


        deallocate(HGTprs)
        deallocate(TMPprs)
        deallocate(RHprs)
        deallocate(TD)
        deallocate(si)
        deallocate(p)       
        deallocate(t850)
        deallocate(t500)
        deallocate(pc)
        deallocate(tc)       
        deallocate(thse)       
        deallocate(rh)
        deallocate(td850)
        deallocate(pp)


end



subroutine cal_ta_si(t850,pc,tc,thse,t500,si,i0,j0,k0,p)                 
          real p(k0)
!***********************************************************************        
!    compute  shawalter index                          by liuhz 2003.11 *      
!     input    t500 =======> temperature at 500 (celsius degree)               
!              t850 =======> temperature at 850 (celsius degree)               
!                pc =======> condesation  pressure (hpa)                        
!                tc =======> condesation  temperature (celsius degree)         
!               thse =======> pseudo-potential temperature at condesation point
!     output     si =======> shawalter index (celsius degree)                  
!***********************************************************************        

      dimension t850(i0,j0),t500(i0,j0),si(i0,j0),tc(i0,j0)                     
      dimension pc(i0,j0),thse(i0,j0)                                          
          real tp0,t2,t1,p_l,thse1,thse2,thse3                                          

      integer kpc                                                               
      pp0=850.                                                                  
        do 200 j=1,j0                                                                  
        do 200 i=1,i0                                                                  
        si(i,j)=9999.                                                                  
      tp0=t850(i,j)                                                            
        if(pc(i,j).gt.9990..or.tp0.lt.-100) goto 200                                   
      kpc=int(pp0/10)-int(pc(i,j)/10)                                          
      p_l=int(pp0/10.)*10.                                                      
        if(kpc.le.0) goto 194                                                         
      do k=1,kpc                                                               
              t1=tp0+(-tp0)/alog(pc(i,j)/pp0)*alog(p_l/pp0)                           
              p_l=p_l-10.                                                              
              if(k.eq.35)then                                                         
                      si(i,j)=t500(i,j)-t1                                                   
                      goto 200                                                               
              endif                                                                    
      enddo                                                                     
194  t1=tc(i,j)                                                               
      p1=pc(i,j)                                                               
        if(t1.lt.-100) goto 200                                                        
      kzero=1                                                                  
!        kzero is a judgement for temperature                                          
10        call cal_thse0_1(p_l,t1,thse1,i0,j0,k0,p)                                    
      call cal_thse0_1(p_l,t1-0.1,thse2,i0,j0,k0,p)                             
        if(thse1.lt.0..or.thse2.lt.0)goto 200                                          
      t2=t1-(thse1-thse(i,j))/(thse1-thse2)*0.1                                 
!        notice difference of thse when tempererture above or beleow zero degree      
      if(t2.lt.0.and.kzero.eq.1)then                                            
      call cal_thse0_1(p_l,t2,thse3,i0,j0,k0,p)                                 
        thse(i,j)=thse3                                                               
      kzero=0                                                                  
      endif                                                                     
      t1=t2                                                                     
      p_l=p_l-10.                                                               
      if(p_l.ge.500.)goto 10                                                   
      si(i,j)=t500(i,j)-t2                                                      
200   continue                                                                  
      return                                                                    
      end                 


        subroutine cal_tc_pc(pp,t,td,pc,tc,i0,j0,k0,p)                           
          real p(k0)
!***********************************************************************        
!    compute  lifting condensation layer and lifting condensation temperature   
!    input  t  =====> lifting temperature (celsius degree)                     
!           td =====>  lifting dew point (celsius degree)                       
!           pp ======>  lifting pressur (hpa)                                   
!    output pc =======> condesation pressure (hpa)                              
!           tc =======> condensation temperature (celsius degree)               
!           thse =====> pseuso equivalent potention temperature (kelvin)        
!***********************************************************************        

      dimension t(i0,j0),td(i0,j0),thse(i0,j0),pp(i0,j0)                        
      dimension tc(i0,j0),pc(i0,j0)                                             
        do 50 j=1,j0                                                                  
        do 50 i=1,i0                                                                  
        tc(i,j)=-120.                                                                  
50        pc(i,j)=9999.                                                               
        do j=1,j0                                                                     
        do i=1,i0                                                                     
        tt=t(i,j)                                                                     
        if(tt.lt.-120.or.td(i,j).lt.-120.or.tt.gt.9990) goto 194                       
      tc(i,j)=tt-(tt-td(i,j))*0.976/(0.976-0.000833*(237.3+td(i,j))**2/(273.15+td(i,j)))                  
      pc(i,j)=pp(i,j)*((273.15+tc(i,j))/(273.15+tt))**3.5                       
!        if(pc(i,j).lt.400)then                                                         
!        tc(i,j)=-120.                                                                  
!      pc(i,j)=9999.                                                            
!        endif                                                                          

194        enddo                                                                     
        enddo                                                                          
!     !计算假相当位温                                                                                       
!      call cal_thse0(pc,tc,thse,i0,j0)                                
      return                                                                                           
      end              


        subroutine cal_thse(pc,tc,thse,i0,j0)                                    
!***********************************************************************        
!    compute  pseuso equivalent potention temperature     by liuhz 2003.12 *   
!     input    tc =======> condesation  temperature (k)                        
!              pc ========>condesation  pressure (hpa)                          
!     output   thse ====> pseudo equivalent potential temperature (k)           
!***********************************************************************        
!      implicit none
       
        integer i0,j0           
    real,dimension(i0,j0) :: tc,thse,pc
       
        integer i,j
        real tt,a,b,e,theta,rmix,paraml                                

      thse=9999.0                                                           
        do j=1,j0                                                                     
        do i=1,i0                                                                     
        tt=tc(i,j)
                                                                            
        if(tt.lt.-100..or.pc(i,j).gt.9990..or.tt.gt.9990.) goto 194                    
                a=7.5*tt/(237.3+tt)                                                      
                e=6.112*10.**a                                                            
                b=1000.0/(pc(i,j)-e)                                                      
                theta=(273.15+tt)*b**0.2856
                rmix=0.622*e/(pc(i,j)-e)
                paramL=2.501*10.**6-2370.*tt                                            
                thse(i,j)=theta*exp(paramL*rmix/(1005.*(273.15+tt)))                                 
194        enddo                                                                     
        enddo                                                                          
    return                                                            
    end

                subroutine cal_thse0_1(pc,tc,thse,i0,j0,k0,p)                                 
          real p(k0)
!***********************************************************************        
!    compute  pseuso equivalent potention temperature   *   
!     input    tc =======> condesation  temperature (k)                        
!              pc ========>condesation  pressure (hpa)                          
!     output   thse ====> pseudo equivalent potential temperature (k)           
!***********************************************************************        

        real pc,tc,thse                                                               

        thse=-100.                                                                     
        tt=tc                                                                          
        if(tt.lt.-100..or.pc.gt.9990.) goto 193                                       
      a=7.5*tt/(237.3+tt)                                                      
      e=6.112*10.**a                                                            
      b=1000.0/(pc-e)                                                           
      thse=(273.15+tt)*b**0.2856*exp((2.501*10.**6-2370.*tt)*0.622*e/(1005.*(273.15+tt)*(pc-e)))                                      
193  continue                                                                  
      return                                                                    
      end

RossilLj

您好，请问您问题解决了吗？我也遇到了相同情况，一直找不到原因