- 积分
- 377
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2022-12-14
- 最后登录
- 1970-1-1
![未绑定新浪微博用户 新浪微博达人勋](source/plugin/sina_login/img/gray.png)
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
这个是子程序,!*-------------------------------------------------------------------------*!* *
!* 读取区域站点逐日降水资料,提取连续过程降水日数及降水量, *
!* 确定各等级百分位阈值,计算满足各阈值条件下各站发生暴雨洪涝的频次。 *
!* 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确实没咋学,但现在也找不到其他资料,求助各位大神
|
|