Fortran 百分位法降水量阈值 所得文件浮点数问题

生气的粉色条纹

这个是子程序，!*-------------------------------------------------------------------------*!*                                                                         *
!*  读取区域站点逐日降水资料，提取连续过程降水日数及降水量，               *
!*  确定各等级百分位阈值，计算满足各阈值条件下各站发生暴雨洪涝的频次。     *         
!*  written by chenxy 2008.5.6  (email:chenxy@cma.gov.cn)                  *
!*  子程序部分                                                             *   
!*-------------------------------------------------------------------------*



                subroutine READAS(filein,ms,cst)
                implicit none
                integer ms
                character*5 cst(1000)
                character*80 filein

                open(22,file=filein,status='old')
                ms=1
100                read(22,*,end=101) cst(ms)
                ms=ms+1
                goto 100
101                close(22)
                ms=ms-1

                return
                end



                subroutine READRD(ms,cst,rd,year1,year2)
                implicit none
                integer ms,year1,year2
                character*5 cst(1000)
                integer rd(1000,year1:year2,366)
                integer iy,i,j,yd
                character*5 cy,filein*80
                logical alive

         rd=-999
         do j=1,ms
                filein=cst(j)//".txt"
                inquire(file=filein,exist=alive)
                if(.not.alive) then
                  print *, filein,' not exist!'
                cycle
                endif
                open(20,file=filein,status="old")
                  do iy=year1,year2
                  yd=365
                  if(mod(iy,400).eq.0.or.(mod(iy,4).eq.0.and.mod(iy,100).ne.0)) yd=366
                  read(20,*) cy
                  read(20,'(30I6)') (rd(j,iy,i),i=1,yd)
                  enddo
                close(20)
          enddo                  

                return
                end


! 统计连续降水日数的过程次数及日数及过程降水量,thr为过程起始开始、开始值
!cc  任意日数的最大连续降水日数rd_con(ns,year1:year2,yd)
!cc  过程累计降水量rs_sum(ns,year1:year2,yd)  
                subroutine RDCONS(ms,rd,year1,year2,rd_con,rd_sum,thr)
                implicit none
                integer ms,year1,year2
                integer rd(1000,year1:year2,366)
                integer iy,yd,i,j,kk,ii
                integer thr

                integer rd_con(1000,year1:year2,366),rd_sum(1000,year1:year2,366)
               
                do iy=year1,year2
                yd=365
                if(mod(iy,400).eq.0.or.(mod(iy,4).eq.0.and.mod(iy,100).ne.0)) yd=366
                do j=1,ms
                do i=1,yd
                rd_con(j,iy,i)=0
                rd_sum(j,iy,i)=0
                  do ii=i,1,-1
                        if(rd(j,iy,ii).ge.thr.and.rd(j,iy,ii).lt.30000) then
                          rd_con(j,iy,i)=rd_con(j,iy,i)+1
                          rd_sum(j,iy,i)=rd_sum(j,iy,i)+rd(j,iy,ii)
                        else
                          exit   
                        endif
                  enddo
                enddo
                enddo  
                enddo

                return
                end


                subroutine floods(fld,ms,year1,year2,rd,rd_con,rd_sum,thresh,daycon,ksy,jsy,per1,per2,FIR)  !FIR 1: thresh<rd<pp1, 0: pp1<rd<pp2, -1:rd>=pp1
                implicit none
                integer ms,year1,year2
                integer rd_con(1000,year1:year2,366),rd_sum(1000,year1:year2,366),rd(1000,year1:year2,366)
                integer rd_con_sum(10000)
                integer thresh,daycon,ksy,jsy,FIR
                integer iy,yd,i,j,kk,ii
                real percentile,per1,per2,pp1,pp2
                real fld(1000)   !统计频次值，一年满足条件的过程次数全部累计

                fld=0
                kk=0
                do iy=ksy,jsy
                yd=365
                if(mod(iy,400).eq.0.or.(mod(iy,4).eq.0.and.mod(iy,100).ne.0)) yd=366
                do j=1,ms
                do i=1,yd-1
                if(daycon.lt.10) then
                  if(rd_con(j,iy,i).eq.daycon.and.rd_con(j,iy,i+1).eq.0) then
                        do ii=i-daycon+1,i
                        if(rd(j,iy,ii).ge.thresh) then
                        kk=kk+1
                        rd_con_sum(kk)=rd_sum(j,iy,i)
                        exit
                        endif                        
                    enddo
                  endif
                else
                  if(rd_con(j,iy,i).ge.daycon.and.rd_con(j,iy,i+1).eq.0) then
                        do ii=i-daycon+1,i
                        if(rd(j,iy,ii).ge.thresh) then
                        kk=kk+1
                        rd_con_sum(kk)=rd_sum(j,iy,i)
                        exit
                        endif                        
                    enddo
                  endif
                endif
                enddo
                enddo
                enddo
               

                if(kk.ge.50) then
                 pp1=percentile(float(rd_con_sum),kk,per1)   !找出百分位为per1%的降水阈值   
                 pp2=percentile(float(rd_con_sum),kk,per2)   !找出百分位为per2%的降水阈值   
                print * ,daycon,per1,'=', pp1/10.0,per2,'=',pp2/10.0

                if(FIR.eq.0)then
                do j=1,ms
                kk=0
                do iy=year1,year2
                yd=365
                if(mod(iy,400).eq.0.or.(mod(iy,4).eq.0.and.mod(iy,100).ne.0)) yd=366
                do i=1,yd-1
                if(daycon.lt.10) then
                  if(rd_con(j,iy,i).eq.daycon.and.rd_con(j,iy,i+1).eq.0) then
                        if(rd_sum(j,iy,i).ge.pp1.and.rd_sum(j,iy,i).lt.pp2) then
                        kk=kk+1
                        endif                        
                  endif
                else
                  if(rd_con(j,iy,i).ge.daycon.and.rd_con(j,iy,i+1).eq.0) then
                        if(rd_sum(j,iy,i).ge.pp1.and.rd_sum(j,iy,i).lt.pp2) then
                        kk=kk+1
                        endif                        
                  endif
                endif
                enddo
                enddo
                fld(j)=float(kk)/(year2-year1+1.0)
                enddo
                elseif(FIR.eq.-1)then
                do j=1,ms
                kk=0
                do iy=year1,year2
                yd=365
                if(mod(iy,400).eq.0.or.(mod(iy,4).eq.0.and.mod(iy,100).ne.0)) yd=366
                do i=1,yd-1
                if(daycon.lt.10) then
                  if(rd_con(j,iy,i).eq.daycon.and.rd_con(j,iy,i+1).eq.0) then
                        if(rd_sum(j,iy,i).ge.pp1) then
                        kk=kk+1
                        endif                        
                  endif
                else
                  if(rd_con(j,iy,i).ge.daycon.and.rd_con(j,iy,i+1).eq.0) then
                        if(rd_sum(j,iy,i).ge.pp1) then
                        kk=kk+1
                        endif                        
                  endif
                endif
                enddo
                enddo
                fld(j)=float(kk)/(year2-year1+1.0)
                enddo

                endif
                endif

                return
                end



                subroutine guiyi(fld1,fld2,fld3,fld4,fld5,ms)
                real fld1(1000),fld2(1000),fld3(1000),fld4(1000),fld5(1000)
                real fmax1,fmax2,fmax3,fmax4,fmax5
                  fmax1=-99999.0
                  fmax2=-99999.0
                  fmax3=-99999.0
                  fmax4=-99999.0
                  fmax5=-99999.0
                do i=1,ms
                   if(fld1(i).gt.-999.0.and.fld1(i).ge.fmax1) fmax1=fld1(i)
                   if(fld2(i).gt.-999.0.and.fld2(i).ge.fmax2) fmax2=fld2(i)
                   if(fld3(i).gt.-999.0.and.fld3(i).ge.fmax3) fmax3=fld3(i)
                   if(fld4(i).gt.-999.0.and.fld4(i).ge.fmax4) fmax4=fld4(i)
                   if(fld5(i).gt.-999.0.and.fld5(i).ge.fmax5) fmax5=fld5(i)
                enddo

                do i=1,ms
           if(fld1(i).gt.-999.0) fld1(i)=fld1(i)/fmax1
           if(fld1(i).gt.-999.0) fld2(i)=fld2(i)/fmax2
           if(fld1(i).gt.-999.0) fld3(i)=fld3(i)/fmax3
           if(fld1(i).gt.-999.0) fld4(i)=fld4(i)/fmax4
           if(fld1(i).gt.-999.0) fld5(i)=fld5(i)/fmax5
                enddo

                return
                end



                subroutine output(filein,ms,cst,year1,year2,rd,fld1,fld2,fld3,fld4,fld5)
                integer year1,year2,ms
                integer rd(1000,year1:year2,366)
                character*5 cst(1000)
                real fld1(1000),fld2(1000),fld3(1000),fld4(1000),fld5(1000)
                character*80 filein

                open(33,file=filein)
                do i=1,ms
                 if(rd(i,year1,1).eq.-999) cycle
                 write(33,"(A5,5F6.2)") cst(i),fld1(i),fld2(i),fld3(i),fld4(i),fld5(i)
                enddo               
                close(33)

                return
                end

!c*******************************
!c       FUNCTION   PERCENTILE  *
!c*******************************

      real function percentile(x, length, per)
      integer length
      real x(length), per
      real xtos(length),bb,cc
      integer nn
       
          MISSING=0.0

      if(per.gt.1.or.per.lt.0) then
        print *, "Function percentile return error: parameter perc"
        stop
      endif

      nn=0
      do i=1, length
        if(x(i).ne.MISSING)then
          nn=nn+1
          xtos(nn)=x(i)
        endif
      enddo
      if(nn.eq.0) then
        percentile=MISSING
      else
        call sort(nn,xtos)
        bb=nn*per+per/3.+1/3.
        cc=real(int(bb))
        if(int(cc).ge.nn) then
          percentile=xtos(nn)
        else
          percentile=xtos(int(cc))+(bb-cc)*(xtos(int(cc)+1)-xtos(int(cc)))
        endif
      endif

          return       
      end

!c*******************************
!c       SUBROUTINE  SORTS      *
!c*******************************

!c---Sorts an array arr(1:n) into ascending numerical order using the Quicksort
!c   algorithm. n is inpu; arr is replace on output by its sorted rearrangement.
!c   Parameters: M is the size of subarrays sorted by straight insertion
!c   and NSTACK is the required auxiliary.
      SUBROUTINE sort(n,arr)
      INTEGER n,M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,temp
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 12 j=l+1,ir
          a=arr(j)
          do 11 i=j-1,1,-1
            if(arr(i).le.a)goto 2
            arr(i+1)=arr(i)
11        continue
          i=0
2         arr(i+1)=a
12      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        temp=arr(k)
        arr(k)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1

          RETURN
      END
不知道是代码的问题还是输入文件的格式问题，但验证了输入文件应该是没问题的，尝试输入几次，最后文件格式都是 51241   NaN  0.00  0.00  0.00  0.00
显示The following floating-point exceptions are signalling: IEEE_INVALID_FLAG
Fortran确实没咋学，但现在也找不到其他资料，求助各位大神