WASH123D模型-配套程序

时间：大约从2011年8月开始
事件：基于张凡老师最得意的WASH123D模型，进行一系列河道参数试验，为了减少手动工作量同时减少我手动修改可能带来的错误，下面是我当时写的一部分Fortran代码（额，f77格式，基本使用gfortran编译，写得比较乱，在此备份，仅供以后参考或改进）
PS：其实，我一直想不明白，我当时为什么一定要用fortran写呢，唯一的理由可能就是闲的O疼了-_-!，令人欣慰的是，这部分工作为我换来了我的第一篇SCI论文。

! Author: zhb
! Last modified time: 2011.8.20
! Task: 修改指定节点的河道形状(按照梯形拟合)
Program modify
Implicit Real*8(A-H,O-Z)
Integer o1,o2,o3,counts(102,3),num(102,3),Itype,w1,w2,d1,d2,nnode,
$ CS(100,3),TmpI(10),num2(3),cts(3),o4,o5,o6,nLine,o7
Character head(3)*3,buf*30,Ctype(3)*5,Ntype*9,form*20,tc*5,
& onec(130)*1,fname*32,form1*20,form2*20,IC1*2,IC3*1,
$ cbuf(220),fname2*32,ctmp*7,ct1*7,c1*1
Real*8 Area(12,2),Pmeter(12,2),Twidth(12,2),CA(12,2),CP(12,2),
& CT(12,2),Winit,Winit2,TmpD(10),DD,ratio,GA(102,12,2),
& GP(102,12,2),GT(102,12,2)
common /N/degree

o1 = 1
o2 = 2
o3 = 3
o4 = 4
o5 = 5
o6 = 6
o7 = 7
open(unit=o1,file=’ALL-DVAPT.txt’,form=’FORMATTED’,status=’old’,
& access=’sequential’)

Print *,’Please enter the global number of node: ‘
Read *,nnode
Print *,’Please enter degree: ‘
Read *,degree
write(fname,502) nnode,degree
502 Format(‘APT_node’,I2,’_D’,F4.2,’.txt’)
open(unit=o2,file=fname,form=’formatted’,status=
& ‘replace’,access=’sequential’)

Do JK=1,102
IF(JK.le.100) THEN

read(o1,1008) head(1),num(JK,1),counts(JK,1),Ctype(1)
Call crackA(Area)
read(o1,1008) head(2),num(JK,2),counts(JK,2),Ctype(2)
CALL crackA(Pmeter)
read(o1,1008) head(3),num(JK,3),counts(JK,3),Ctype(3)
CALL crackA(Twidth)
ELSE
read(o1,1018) head(1),num(JK,1),counts(JK,1),Ctype(1)
Call crackA(Area)
read(o1,1018) head(2),num(JK,2),counts(JK,2),Ctype(2)
CALL crackA(Pmeter)
read(o1,1018) head(3),num(JK,3),counts(JK,3),Ctype(3)
CALL crackA(Twidth)
ENDIF
Do I=1,12
GA(JK,I,1) = Area(I,1)
GA(JK,I,2) = Area(I,2)
GP(JK,I,1) = Pmeter(I,1)
GP(JK,I,2) = Pmeter(I,2)
GT(JK,I,1) = Twidth(I,1)
GT(JK,I,2) = Twidth(I,2)
ENDDO

ENDDO

!初始化河道编号数组
open(unit=o5,file=’CS.txt’,form=’formatted’,status=’old’,
& access=’sequential’)
DO I=1,100
READ(o5,1009) IC1,IC3,(onec(J),J=1,100)
In = 4
I1 = 1
CALL CrackI(I1,In,onec,TmpI,’INTEGER ‘)
M = TmpI(1)
Do J1=1,3
CS(M,J1) = TmpI(J1+1)
ENDDO
ENDDO

!找到指定节点的河道参数
IA = CS(nnode,1)
IP = CS(nnode,2)
IT = CS(nnode,3)
IndexA = 0
IndexP = 0
IndexT = 0
Icount = 0 !计数器
DO I=1,102
IF(IA.eq.num(I,1)) THEN
IndexA = I
Icount = Icount+1
ENDIF
IF(IP.eq.num(I,2)) THEN
IndexP = I
Icount = Icount+1
ENDIF
IF(IT.eq.num(I,3)) THEN
IndexT = I
Icount = Icount+1
ENDIF
IF(Icount.eq.3) goto 501
ENDDO
501 continue
!提取要修改的参数数组
Do I=1,12
Area(I,1) = GA(IndexA,I,1)
Area(I,2) = GA(IndexA,I,2)
Pmeter(I,1) = GP(IndexA,I,1)
Pmeter(I,2) = GP(IndexA,I,2)
Twidth(I,1) = GT(IndexA,I,1)
Twidth(I,2) = GT(IndexA,I,2)
c print *,’Area:’,Area(I,2)
ENDDO
c read *,II
cts(1) = counts(IndexA,1)
cts(2) = counts(IndexP,2)
cts(3) = counts(IndexT,3)
num2(1) = num(IndexA,1)
num2(2) = num(IndexP,2)
num2(3) = num(IndexT,3)
W = Twidth(3,2)
h = Twidth(3,1)
Winit = 2.0*Area(3,2)/h-w
xie=0.0
IF(Winit.lt.0.000001) THEN
hh=2*Area(3,2)/W
xie=HYPO(hh,(W/2.0))
W=W*degree
xie2=HYPO(hh,(W/2.0))
dd = 2.0*(xie2-xie) !当梯形拟合失效而采用漏斗状拟合时，湿周增加的长度
print *,’Winit 为负值！’
read *,II
ENDIF
!修改顶部河宽
CALL changeT(degree,Winit,Twidth,CT,Winit2)
!修改横断面积
CALL changeA(Area,CA,CT,Winit2)
!修改湿周
CALL changeP(Pmeter,CP,CT,Winit2,Winit,dd)

Call computeD(Area,CA,DD,ratio)
print *,DD,ratio
Call computeD(Pmeter,CP,DD,ratio)
print *,DD,ratio
DO I=1,3
IF(num2(I).lt.1000) THEN
write(o2,1008) head(I),num2(I),cts(I),Ctype(I)
Else
write(o2,1018) head(I),num2(I),cts(I),Ctype(I)
ENDIF
Do J=1,12
tc = Ctype(I)
IF(tc(2:4).eq.’DVA’) THEN
IF(J.le.2) THEN
IF(J.eq.1) write(o2,'(F3.1,1x,F3.1)’) CA(J,1),CA(J,2)
IF(J.eq.2) write(o2,'(F5.3,1x,F3.1)’) CA(J,1),CA(J,2)
ELSE
Call getLength(CA(J,1),w1,4)
Call getLength(CA(J,2),w2,9)
write(form,1012) w1,4,w2,9
write(o2,form) CA(J,1),CA(J,2)
ENDIF
ELSE IF(tc(2:4).eq.’DVP’) THEN
IF(J.le.2) THEN
IF(J.eq.1) write(o2,'(F3.1,1x,F3.1)’) CP(J,1),CP(J,2)
IF(J.eq.2) write(o2,'(F5.3,1x,F3.1)’) CP(J,1),CP(J,2)
ELSE
Call getLength(CP(J,1),w1,4)
Call getLength(CP(J,2),w2,9)
write(form,1012) w1,4,w2,9
write(o2,form) CP(J,1),CP(J,2)
ENDIF
ELSE IF(tc(2:4).eq.’DVT’) THEN
IF(J.le.2) THEN
IF(J.eq.1) write(o2,'(F3.1,1x,E8.1)’) CT(J,1),CT(J,2)
IF(J.eq.2) write(o2,'(F5.3,1x,E8.1)’) CT(J,1),CT(J,2)
ELSE
Call getLength(CT(J,1),w1,4)
Call getLength(CT(J,2),w2,9)
write(form,1012) w1,4,w2,9
write(o2,form) CT(J,1),CT(J,2)
ENDIF
ENDIF
ENDDO
ENDDO

!修改整个1bc文件
write(fname2,1101) nnode,degree
1101 Format(‘n’,I2,’_d’,F4.2,’_LY_1D.1bc’)
write(ct1,'(A3,1X,I3)’) ‘XYS’,num2(1)
open(unit=o6,file=’LY_1D.1bc’,form=’formatted’,
$ status=’old’,access=’sequential’)
Do I=1,4514
read(o6,'(A7)’) cTmp
IF(cTmp.eq.ct1) THEN
nLine = I
goto 1102
ENDIF
ENDDO
1102 continue
close(o6)
open(unit=o6,file=’LY_1D.1bc’,form=’formatted’,
$ status=’old’,access=’sequential’)
open(unit=o7,file=fname2,form=’formatted’,
$ status=’replace’,access=’sequential’)
rewind(o6)
Do I5=1,4512
Do I=1,220
cbuf(I) = ‘ ‘
ENDDO
IF(I5.gt.nLine.and.I5.le.(nLine+38)) goto 1105
IF(nLine.eq.I5) THEN
Do I7=1,39
read(o6,'(A1)’) c1
ENDDO
DO I=1,3
IF(num2(I).lt.1000) THEN
write(o7,1008) head(I),num2(I),cts(I),Ctype(I)
Else
write(o7,1018) head(I),num2(I),cts(I),Ctype(I)
ENDIF
Do J=1,12
tc = Ctype(I)
IF(tc(2:4).eq.’DVA’) THEN
IF(J.le.2) THEN
IF(J.eq.1) write(o7,'(F3.1,1x,F3.1)’) CA(J,1),CA(J,2)
IF(J.eq.2) write(o7,'(F5.3,1x,F3.1)’) CA(J,1),CA(J,2)
ELSE
Call getLength(CA(J,1),w1,4)
Call getLength(CA(J,2),w2,9)
write(form,1012) w1,4,w2,9
write(o7,form) CA(J,1),CA(J,2)
ENDIF
ELSE IF(tc(2:4).eq.’DVP’) THEN
IF(J.le.2) THEN
IF(J.eq.1) write(o7,'(F3.1,1x,F3.1)’) CP(J,1),CP(J,2)
IF(J.eq.2) write(o7,'(F5.3,1x,F3.1)’) CP(J,1),CP(J,2)
ELSE
Call getLength(CP(J,1),w1,4)
Call getLength(CP(J,2),w2,9)
write(form,1012) w1,4,w2,9
write(o7,form) CP(J,1),CP(J,2)
ENDIF
ELSE IF(tc(2:4).eq.’DVT’) THEN
IF(J.le.2) THEN
IF(J.eq.1) write(o7,'(F3.1,1x,E8.1)’) CT(J,1),CT(J,2)
IF(J.eq.2) write(o7,'(F5.3,1x,E8.1)’) CT(J,1),CT(J,2)
ELSE
Call getLength(CT(J,1),w1,4)
Call getLength(CT(J,2),w2,9)
write(form,1012) w1,4,w2,9
write(o7,form) CT(J,1),CT(J,2)
ENDIF
ENDIF
ENDDO
ENDDO
ENDIF
read(o6,'(220A1)’) (cbuf(J5),J5=1,220)
c print *,I5
Do I1=220,1,-1
IF(cbuf(I1).ne.’ ‘) THEN
Ilen = I1
goto 111
ENDIF
ENDDO
111 continue
write(o7,'(220A1)’) (cbuf(I2),I2=1,Ilen)
1105 continue
ENDDO
1001 Format(130A1)
1008 Format(A3,1x,I3,1x,I2,1x,A5)
1009 Format(A2,A1,130A1)
1012 Format(‘(F’,I1,’.’,I1,’,1x’,’,F’,I2,’.’,I1,’)’)
1013 Format(‘(F’,I1,’.’,I1,’)’)
1014 Format(‘(F’,I2,’.’,I1,’)’)
1018 Format(A3,1x,I4,1x,I2,1x,A5)
END
Subroutine crackA(Area)
Implicit Real*8(A-H,O-Z)
Character buf*30,Ntype*9,onec(130)*1
Real*8 Area(12,2),TmpD(10)
Do J=1,12
I1=1
In=2
Ntype = ‘REAL ‘
Read(1,'(130A1)’) (onec(N),N=1,130)
Call CrackD(I1,In,onec,TmpD,Ntype)
Do I2=1,2
Area(J,I2) = TmpD(I2)
ENDDO
ENDDO
END

subroutine CrackI(I1,In,onec,TmpI,NTYPE)
Integer I1,I2,In,TmpI(10)
Character onec(130)*1,Blank*1,Ibuf*130,IFOR*10,NTYPE*9
Blank = ‘ ‘
IF(I1.GT.130) GOTO 299
In1 = In
Do 201 I=1,In1
202 Continue
IF(onec(I1).NE.Blank) goto 215
I1 = I1 + 1
IF(I1.ge.130) THEN
In = I – 1
goto 299
ENDIF
goto 202
215 continue
I2 = I1+1
221 continue
IF(onec(I1).eq.””) THEN
IF(onec(I2).eq.””) goto 225
Else
IF(onec(I2).eq.Blank) goto 225
ENDIF
I2 = I2+1
IF(I2.lt.131) goto 221
I2 = 130
225 continue
IF(onec(I2).eq.””) I1 = I1+1
Length = I2-I1
IF(Length.lt.1) Length = 1
I2 = I2-1

write(Ibuf,'(130A1)’) (onec(J),J=I1,I2)
IF(NTYPE(1:7).eq.’INTEGER’) THEN
WRITE(IFOR,1011) LENGTH
1011 FORMAT(‘(I’,I2,’)’)
Read(Ibuf,IFOR) TmpI(I)
c write(*,*) ‘TmpI’,I,TmpI(I)
c Read *,II
ELSE
goto 299
ENDIF
IF(onec(I2+1).EQ.””) I2=I2+1
I1 = I2+1
201 Continue
299 continue
Return
END
Subroutine CrackD(I1,In,onec,TmpD,NTYPE)
Implicit Real*8(A-H,O-Z)
Integer I1,I2,In
Character onec(130)*1,Blank*1,Ibuf*130,FFOR*10,NTYPE*9
Real*8 TmpD(10)
Blank = ‘ ‘
IF(I1.GT.130) GOTO 399
In1 = In
Do 301 I=1,In1
302 Continue
IF(onec(I1).NE.Blank) goto 315
I1 = I1 + 1
IF(I1.ge.130) THEN
IN = I – 1
goto 399
ENDIF
goto 302
315 continue
I2 = I1+1
321 continue
IF(onec(I1).eq.””) THEN
IF(onec(I2).eq.””) goto 325
Else
IF(onec(I2).eq.Blank) goto 325
ENDIF
I2 = I2+1
IF(I2.lt.131) goto 321
I2 = 130
325 continue
IF(onec(I2).eq.””) I1 = I1+1
Length = I2-I1
IF(Length.lt.1) Length = 1
I2 = I2-1

END
subroutine changeP(Pmeter,CP,CT,Winit2,Winit,dd)
Implicit Real*8(A-H,O-Z)
real*8 Pmeter(12,2),CP(12,2),Bedge,HYPO,CT(12,2),former,Winit2,
$ Winit,dd
common /N/degree
CP(1,1) = Pmeter(1,1)
CP(1,2) = Pmeter(1,2)
CP(2,1) = Pmeter(2,1)
CP(2,2) = Pmeter(2,2)
CP(3,1) = Pmeter(3,1)
former = Pmeter(3,2)
CP(3,2) = former-Winit+Winit2
d = CP(3,2) – Pmeter(3,2)
IF(Winit.lt.0.000001) THEN
d = dd
CP(3,2) = Pmeter(3,2) + d
ENDIF
DO I=4,12
CP(I,1) = Pmeter(I,1)
CP(I,2) = d+Pmeter(I,2)
ENDDO
END

!计算斜边
Real*8 Function HYPO(X,Y)
REAL*8 X,Y
HYPO = SQRT(X**2+Y**2)
END
subroutine changeT(degree,Winit,Twidth,CT,Winit2)
Implicit Real*8(A-H,O-Z)
REAL*8 Twidth(12,2),CT(12,2),Winit

CT(1,1) = Twidth(1,1)
CT(1,2) = Twidth(1,2)
CT(2,1) = Twidth(2,1)
CT(2,2) = Twidth(2,2)
Winit2 = Winit*degree
Do I=3,12
CT(I,1) = Twidth(I,1)
CT(I,2) = Twidth(I,2)*degree
ENDDO
END
Real*8 Function Sladder(X,Y,H)
REAL*8 X,Y,H
Sladder = ((X+Y)*H)/2.0
END
Subroutine getLength(F,w,d)
Real*8 F
character*10 form
character*20 tmp
Integer w,d
write(form,2199) d
2199 Format(‘(F20.’,I1,’)’)
write(tmp,form) F
Do I=1,20
IF(tmp(I:I).ne.’ ‘.and.I.ne.1) THEN
w=20-I+1
goto 2002
ENDIF
ENDDO
2002 continue
END
Subroutine computeD(D1,D2,DD,ratio)
REAL*8 D1(12,2),D2(12,2),tmp,DD,ratio
Integer Ilabel

DD = D1(1,2)-D2(1,2)
Do I=2,12
tmp = DABS(D1(I,2)-D2(I,2))
IF(DD.lt.tmp) Then
DD = tmp
Ilabel = I
ENDIF
ENDDO
ratio = DD/D1(Ilabel,2)
END

!########################################################################
! Author: zhb
! Last modified time: 2011.11.25
! Task: 计算峰值时刻（用户预定义）各节点的水深所对应的河宽
Program look
Implicit Real*8(A-H,O-Z)
Integer o1,o2,o3,counts(3),num(3),Itype,w1,w2,d1,d2,
& Inum,Ia,Ip,It,XYS(102,4),TmpI(10),CS1(100,3),
& Bnum,M,JU(16,2),B1(59,2),B2(25,2),CSN(100)
Character head(3)*3,buf*30,Ctype(3)*5,Ntype*9,form*20,tc*5,
& onec(130)*1,fname*36,IC1*2, IC3*1,C11*6,C21*2
Real*8 Area(12,2),Pmeter(12,2),Twidth(12,2),WA(102,12,2),
& WP(102,12,2),WT(102,12,2),Tmp(100),Ivalue,JUR(16,3),
& Winit,Winit2,TmpD(10),DD,ratio,B1R(59,3),B2R(25,3),
$ XC(100,3),Elem(97,2),stage(100),WSTIME,Wdepth(100)
Common /N1/WA
Common /N2/WP
Common /N3/WT
Common /M1/B1
COmmon /M2/B2
Common /M3/JU
Common /M4/B1R
COmmon /M5/B2R
Common /M6/JUR
Common /N4/Wdepth

o1 = 1
o5 = 5
o2 = 2
o3 = 3
open(unit=o2,file=’LY_1D.1dm’,form=’formatted’,
& status=’old’,access=’sequential’)
read(o2,'(A6)’) C11
read(o2,'(A2)’) C21
read(o2,'(A2)’) C21
read(o2,'(A2)’) C21

!读取单元信息
Do J=1,97
read(o2,1002) IC1,IC3,(onec(I),I=1,130)
I1 = 1
In = 4
Ntype = ‘INTEGER ‘
Call CrackI(I1,In,onec,TmpI,Ntype)
M = TmpI(1)
Elem(M,1) = TmpI(2)
Elem(M,2) = TmpI(3)
ENDDO
!读取节点坐标
Do J=1,100
read(o2,1002) IC1,IC3,(onec(I),I=1,130)
I1 = 1
In = 1
Ntype = ‘INTEGER ‘
Call CrackI(I1,In,onec,TmpI,Ntype)
M = TmpI(1)
In = 3
Ntype = ‘REAL ‘
Call CrackD(I1,In,onec,TmpD,Ntype)
XC(M,1) = TmpD(1)
XC(M,2) = TmpD(2)
XC(M,3) = TmpD(3)
ENDDO

open(unit=o3,file=’fort.20′,FORM=’FORMATTED’,STATUS=’OLD’,ACCESS=
$ ‘SEQUENTIAL’)
!指针到峰值时刻之前
Do I=1,6
read(o3,'(A1)’) IC3
ENDDO
Do I=1,28
Do J=1,101
read(o3,'(A1)’) IC3
ENDDO
ENDDO

!读取峰值时刻，各节点的水位
read(o3,'(A2,I4,D16.8)’) IC1,ISAT,WSTIME
Do I=1,100
read(o3,'(E16.8)’) stage(I)
ENDDO

!计算峰值时刻，各节点的水深，以用于插值
Do I=1,100
Wdepth(I) = stage(I) – XC(I,3)
c print *,I,Wdepth(I)
ENDDO
open(unit=o1,file=’LY_1D.1bc’,form=’FORMATTED’,status=’old’,
& access=’sequential’)

!跳到CS开始之前
Do I=1,117
read(o1,'(A1)’) IC3
ENDDO

!读取CS编号
Do I=118,217
READ(o1,1002) IC1,IC3,(onec(J),J=1,130)
In = 4
I1 = 1
CALL CrackI(I1,In,onec,TmpI,’INTEGER ‘)
M = TmpI(1)
Do 132 J1=2,4
CS1(M,J1-1) = TmpI(J1)
132 continue
ENDDO
Do I=218,474 !到XYS 400之前
read(o1,'(A1)’) IC3
ENDDO
Do I1=1,102
Do I2=1,4
XYS(I1,I2) = -1
ENDDO
ENDDO
!干流,第二项是其在XYS中的编号
Do I=1,16
JU(I,1) = I
JU(I,2) = I
ENDDO

!支流1（42-100）
Do I=1,59
B1(I,1) = I+41
B1(I,2) = I+41
ENDDO

!支流2（17-41）
Do I=1,25
B2(I,1) = I+16
B2(I,2) = I+16
ENDDO

c IL = 0 !汇合点计数器，为了区分出节点100、41和1
Do I0=1,102
read(o1,1008) head(1),(onec(J),J=1,130)
In = 2
I1 = 1
Call CrackI(I1,In,onec,TmpI,’INTEGER ‘)
num(1) = TmpI(1)
Call crackA(Area)
read(o1,1008) head(2),(onec(J),J=1,130)
In = 2
I1 = 1
Call CrackI(I1,In,onec,TmpI,’INTEGER ‘)
num(2) = TmpI(1)
CALL crackA(Pmeter)
read(o1,1008) head(3),(onec(J),J=1,130)
In = 2
I1 = 1
Call CrackI(I1,In,onec,TmpI,’INTEGER ‘)
num(3) = TmpI(1)
CALL crackA(Twidth)
XYS(I0,2) = num(1)
XYS(I0,3) = num(2)
XYS(I0,4) = num(3)
!CSN记录了每个节点所对应的XYS编号
CALL getXYSnum(CSN,num(1),num(2),num(3),CS1,I0)
Do I=1,2
Do J=1,12
WA(I0,J,I) = Area(J,I)
WP(I0,J,I) = Pmeter(J,I)
WT(I0,J,I) = Twidth(J,I)
ENDDO
ENDDO
ENDDO
CSN(41) = CSN(1)
CSN(100) = CSN(1)

Do I=1,16
JU(I,2) = CSN(I)
Bnum = 3
Ivalue=Wdepth(I)
Call inpo(Ivalue,CSN(I),Bnum,I)
ENDDO

Do I=1,59
IF(I.eq.59) THEN
B1(I,2) = JU(1,2)
Bnum = 1
Ivalue=Wdepth(I+41)
Call inpo(Ivalue,JU(1,2),Bnum,I)
goto 911
ENDIF
B1(I,2) = CSN(I+41)
Bnum = 1
Ivalue=Wdepth(I+41)
Call inpo(Ivalue,CSN(I+41),Bnum,I)
911 continue
ENDDO
Do I=1,25
IF(I.eq.25) THEN
B2(I,2) = JU(1,2)
Bnum = 2
Ivalue=Wdepth(I+41)
Call inpo(Ivalue,JU(1,2),Bnum,I)
goto 912
ENDIF
J = CSN(I+16)
B2(I,2) = J
Bnum = 2
Ivalue=Wdepth(I+16)
Call inpo(Ivalue,J,Bnum,I)
912 continue
ENDDO

print *,’Please input the branch number you want to see: 0-2′
read *,Bnum
print *,’Please input the parameter you want to get: ‘
print *,’1:area 2:perimeter 3:topWidth’
read *,II
DO I=1,100
Tmp(I) = -1.0
ENDDO
write(fname,1201) Bnum,II
IF(Bnum.eq.1) THEN
IF(II.eq.1) THEN
DO I=1,59
Tmp(I) = B1R(I,1)
ENDDO
CAll getInfo(59,Tmp,fname)
ELSE IF(II.eq.2) THEN
DO I=1,59
Tmp(I) = B1R(I,2)
ENDDO
CAll getInfo(59,Tmp,fname)
ELSE IF(II.eq.3) THEN
DO I=1,59
Tmp(I) = B1R(I,3)
ENDDO
CAll getInfo(59,Tmp,fname)
ENDIF
ELSE IF(Bnum.eq.2) THEN
IF(II.eq.1) THEN
DO I=1,25
Tmp(I) = B2R(I,1)
ENDDO
CAll getInfo(25,Tmp,fname)
ELSE IF(II.eq.2) THEN
DO I=1,25
Tmp(I) = B2R(I,2)
ENDDO
CAll getInfo(25,Tmp,fname)
ELSE IF(II.eq.3) THEN
DO I=1,25
Tmp(I) = B2R(I,3)
ENDDO
CAll getInfo(25,Tmp,fname)
ENDIF
ELSE IF(Bnum.eq.0) THEN
IF(II.eq.1) THEN
DO I=1,16
Tmp(I) = JUR(I,1)
ENDDO
CAll getInfo(16,Tmp,fname)
ELSE IF(II.eq.2) THEN
DO I=1,16
Tmp(I) = JUR(I,2)
ENDDO
CAll getInfo(16,Tmp,fname)
ELSE IF(II.eq.3) THEN
DO I=1,16
Tmp(I) = JUR(I,3)
ENDDO
CAll getInfo(16,Tmp,fname)
ENDIF
ENDIF

1201 Format(‘B’,I1,’_P’,I1,’.txt’)
1001 Format(130A1)
1002 Format(A2,A1,130A1)
1008 Format(A3,1x,130A1)
1012 Format(‘(F’,I1,’.’,I1,’,1x’,’,F’,I2,’.’,I1,’)’)
END

Subroutine crackA(Area)
Implicit Real*8(A-H,O-Z)
Character buf*30,Ntype*9,onec(130)*1
Real*8 Area(12,2),TmpD(10)
Do J=1,12
I1=1
In=2
Ntype = ‘REAL ‘
Read(1,'(130A1)’) (onec(N),N=1,130)
Call CrackD(I1,In,onec,TmpD,Ntype)
Do I2=1,2
Area(J,I2) = TmpD(I2)
ENDDO
ENDDO
END

Subroutine CrackD(I1,In,onec,TmpD,NTYPE)
Implicit Real*8(A-H,O-Z)
Integer I1,I2,In
Character onec(130)*1,Blank*1,Ibuf*130,FFOR*10,NTYPE*9
Real*8 TmpD(10)
Blank = ‘ ‘
IF(I1.GT.130) GOTO 399
In1 = In
Do 301 I=1,In1
302 Continue
IF(onec(I1).NE.Blank) goto 315
I1 = I1 + 1
IF(I1.ge.130) THEN
IN = I – 1
goto 399
ENDIF
goto 302
315 continue
I2 = I1+1
321 continue
IF(onec(I1).eq.””) THEN
IF(onec(I2).eq.””) goto 325
Else
IF(onec(I2).eq.Blank) goto 325
ENDIF
I2 = I2+1
IF(I2.lt.131) goto 321
I2 = 130
325 continue
IF(onec(I2).eq.””) I1 = I1+1
Length = I2-I1
IF(Length.lt.1) Length = 1
I2 = I2-1

write(Ibuf,'(130A1)’) (onec(j),j=I1,I2)
IF(NTYPE(1:4).eq.’REAL’) THEN
WRITE(FFOR,1031) LENGTH
1031 FORMAT(‘(F’,I2,’.0)’)
READ(Ibuf,FFOR) TmpD(I)
ELSE
goto 399
ENDIF
IF(onec(I2+1).EQ.””) I2=I2+1
I1 = I2+1
301 Continue
399 continue
Return
END
subroutine CrackI(I1,In,onec,TmpI,NTYPE)
Integer I1,I2,In,TmpI(10)
Character onec(130)*1,Blank*1,Ibuf*130,IFOR*10,NTYPE*9
Blank = ‘ ‘
IF(I1.GT.130) GOTO 299
In1 = In
Do 201 I=1,In1
202 Continue
IF(onec(I1).NE.Blank) goto 215
I1 = I1 + 1
IF(I1.ge.130) THEN
In = I – 1
goto 299
ENDIF
goto 202
215 continue
I2 = I1+1
221 continue
IF(onec(I1).eq.””) THEN
IF(onec(I2).eq.””) goto 225
Else
IF(onec(I2).eq.Blank) goto 225
ENDIF
I2 = I2+1
IF(I2.lt.131) goto 221
I2 = 130
225 continue
IF(onec(I2).eq.””) I1 = I1+1
Length = I2-I1
IF(Length.lt.1) Length = 1
I2 = I2-1

Do I=1,100
IF(CS(I,1).eq.Ia.and.CS(I,2).eq.Ip.and.CS(I,3).eq.It) THEN
CSN(I) = I0
ENDIF
ENDDO
901 Continue

END
Subroutine inpo(X,J,Bnum,num)
Integer JU(16,2),B1(59,2),B2(25,2),Bnum,num,J
REAL*8 S1,S2(3,12,2),Y(3),K,Y1,Y2,B,YY,WA(102,12,2),
& WP(102,12,2),WT(102,12,2),X,B1R(59,3),B2R(25,3),JUR(16,3)
Common /N1/WA
Common /N2/WP
Common /N3/WT
Common /M1/B1
COmmon /M2/B2
Common /M3/JU
Common /M4/B1R
COmmon /M5/B2R
Common /M6/JUR
Do I=1,2
Do J1=1,12
S2(1,J1,I) = WA(J,J1,I)
S2(2,J1,I) = WP(J,J1,I)
S2(3,J1,I) = WT(J,J1,I)
ENDDO
ENDDO
S1 = X
Do 120 I1=1,3
IF (S1.gt.S2(I1,12,1)) THEN
print *,’Error: The water stage is larger than top’
goto 999
ENDIF
Do 121 I2=1,11
IF ( S1.ge.S2(I1,I2,1).AND.S1.le.S2(I1,I2+1,1) ) THEN
Y1 = S2(I1,I2,2)
Y2 = S2(I1,I2+1,2)
K = (Y2-Y1)/(S2(I1,I2+1,1)-S2(I1,I2,1))
B = Y1 – K*S2(I1,I2,1)
YY = K*S1+B
Y(I1) = YY
goto 120
ENDIF
121 Continue
120 Continue

IF(Bnum.eq.1) THEN
DO I=1,3
B1R(num,I) = Y(I)
ENDDO
ELSE IF(Bnum.eq.2) THEN
DO I=1,3
B2R(num,I) = Y(I)
ENDDO
ELSE IF(Bnum.eq.3) THEN
DO I=1,3
JUR(num,I) = Y(I)
ENDDO
ENDIF
999 END
Subroutine getInfo(num,B,fname)
REAL*8 B(100),Max,Min,Mean,Total,Tmp(100),ratio,WA(102,12,2),
$ WT(102,12,2),Wdepth(100)
Character*36 fname
Integer MaxID,MinID,JU(16,2),B1(59,2),B2(25,2)
Common /M1/B1
COmmon /M2/B2
Common /M3/JU
Common /N1/WA
Common /N3/WT
Common /N4/Wdepth
Max = B(1)
Min = B(1)
Total = B(1)
Do I=2,num
IF(B(I).lt.0) goto 502
IF(Max.le.B(I)) THEN
Max = B(I)
MaxID = I
ENDIF
IF(Min.ge.B(I)) THEN
Min = B(I)
MinID = I
ENDIF
Total = Total+B(I)
ENDDO
502 continue
Mean = Total/num
open(unit=103,file=’WaterDepth.txt’,form=’formatted’,
&status=’replace’,
& access=’sequential’)
open(unit=101,file=’TopWidth.txt’,form=’formatted’,
&status=’replace’,
& access=’sequential’)
open(unit=102,file=’RiverwayDepth.txt’,form=’formatted’
$ ,status=’replace’,access=’sequential’)
open(unit=100,file=fname,form=’formatted’,status=’replace’,
& access=’sequential’)

c write(100,401) ‘Max: ‘,Max,’ID: ‘,MaxID
c write(100,401) ‘Min: ‘,Min,’ID: ‘,MinID
c write(100,402) ‘Mean: ‘,mean
c write(100,*) ‘ ‘
c write(100,*) ‘Data Part: ‘

Do I=1,num
IF(num.eq.59) THEN
J = I+41
K = B1(I,2)
ELSE IF(num.eq.25) THEN
J = I+16
K = B2(I,2)
ELSE IF(num.eq.16) THEN
J = I
K = JU(I,2)
ENDIF
ratio = B(I)/Max
c write(100,400) I,J,B(I),ratio,WA(K,12,1)
write(100,'(I2,1X,F16.9)’) J,B(I)
write(101,'(F16.9)’) WT(K,12,2)
write(102,'(F6.2)’) WA(K,12,1)
write(103,'(F16.9)’) Wdepth(J)
ENDDO
400 Format(I2,1X,I3,1x,F10.4,2x,F4.2,2x,F6.2)
401 Format(A5,F10.4,3X,A4,I2)
402 Format(A5,F10.4)
END
!########################################################################

! Author: zhb
! Last modified time: 2011.8.20
! Task: 计算模拟结果(stage\velocity\discharge)与基准案例之间的差异——以R方值表征
Program compute
Implicit Real*8(A-H,O-Z)
Character*2 c5,c6
Character*4 c7
Character*7 c1,c2
Character*6 c3,c4
Character*9 c9
Character*80 fname1,fname2
Character*40 c8
integer o11,o12,o21,o22,o31,o32,o41,o42,o51,o52,ISTAT(48)
& ,IST1(48),TL,TL1 !TL–Time Label
REAL*8 VEL(48),DCH(48),VEL1(48),DCH1(48),
& STG(48),STG1(48),TMP,TMP1,
& AMAX,BMAX,ATIME,BTIME,Dvalue,Ratio
Logical fexist

o1 = 1
o2 = 2
o3 = 3
o4 = 4
o5 = 5
o6 = 6
o7 = 7
o8 = 8

c9 = ‘FORMATTED’
!读取原数据，basecase
open(unit=o1,file=’原stage.txt’,FORM=c9,STATUS=’OLD’,ACCESS=
$ ‘SEQUENTIAL’)
open(unit=o2,file=’原velocity.txt’,FORM=c9,STATUS=’OLD’,ACCESS=
$ ‘SEQUENTIAL’)
open(unit=o3,file=’原discharge.txt’,FORM=c9,STATUS=’OLD’,ACCESS=
$ ‘SEQUENTIAL’)

open(unit=o4,file=’stage.txt’,FORM=c9,STATUS=’OLD’,ACCESS=
$ ‘SEQUENTIAL’)
open(unit=o5,file=’velocity.txt’,FORM=c9,STATUS=’OLD’,ACCESS=
$ ‘SEQUENTIAL’)
open(unit=o6,file=’discharge.txt’,FORM=c9,STATUS=’OLD’,ACCESS=
$ ‘SEQUENTIAL’)

Do I=1,48
read(o1,1008) STG(I)
read(o2,1008) VEL(I)
read(o3,1008) DCH(I)
read(o4,1008) STG1(I)
read(o5,1008) VEL1(I)
read(o6,1008) DCH1(I)
ENDDO
open(unit=o7,file=’result_Rsquare.txt’,form=’formatted’
&,status=’replace’,
& access = ‘sequential’)
open(unit=o8,file=’meta.txt’,form=’formatted’,status=’replace’,
& access = ‘sequential’)
R2a = Rsquare(STG1,STG)
R2b = Rsquare(VEL1,VEL)
R2c = Rsquare(DCH1,DCH)
write(o7,1101) R2a
1101 Format(‘水位 R-square: ‘,F11.9)
write(o7,1102) R2b
1102 Format(‘流速 R-square: ‘,F11.9)
write(o7,1103) R2c
1103 Format(‘流量 R-square: ‘,F11.9)

CALL comp(STG,STG1,AMAX,BMAX,TL,TL1)
dvalue = BMAX – AMAX
ratio = dvalue/Amax
write(o7,1013) dvalue,ratio*100
1013 Format(‘峰值水位相差: ‘,F10.4,1x,’比例是: ‘,F10.2,’%’,/)
c write(o7,*) ‘峰值时相标签分别是：’,TL,’ ‘,TL1
write(o8,*) dvalue,ratio
c write(o8,*) ‘ ‘

CALL comp(VEL,VEL1,AMAX,BMAX,TL,TL1)
dvalue = BMAX – AMAX
ratio = dvalue/Amax
write(o7,1014) dvalue,ratio*100
1014 Format(‘峰值速度相差: ‘,F10.4,1x,’比例是: ‘,F10.2,’%’,/)
c write(o7,*) ‘峰值时相标签分别是：’,TL,’ ‘,TL1
write(o8,*) dvalue,ratio
c write(o8,*) ‘ ‘
CALL comp(DCH,DCH1,AMAX,BMAX,TL,TL1)
dvalue = BMAX – AMAX
ratio = dvalue/Amax
write(o7,1015) dvalue,ratio*100
1015 Format(‘峰值流量相差: ‘,F10.4,1x,’比例是: ‘,F10.2,’%’,/)
c write(o7,*) ‘峰值时相标签分别是：’,TL,’ ‘,TL1
write(o8,*) dvalue,ratio
c write(o8,*) ‘ ‘
1001 FORMAT(A7)
1002 FORMAT(A6)
1003 FORMAT(A2,I8)
1004 FORMAT(A4,1X,A40)
1005 FORMAT(A2,I4,D16.8)
1006 FORMAT(2E16.8)
1007 FORMAT(A7,1X,A6)
1008 FORMAT(2E16.8)
1009 Format(A20)
1010 Format(E10.7)
1011 FORMAT(E15.8)

999 END

Real Function Rsquare(Y,YS)
Implicit Real*8(A-H,O-Z)
Real*8 Y(48),YS(48),SSE,SST,meanYS
!YS是原数据
SSE=0.0
SST=0.0
meanYS=0.0
t=0.0
Do I=1,48
t=t+YS(I)
ENDDO
meanYS=t/48.0
Do I=1,48
SSE=SSE+((Y(I)-YS(I)))**2
SST=SST+(YS(I)-meanYS)**2
ENDDO
Rsquare = 1-(SSE/SST)
END
subroutine comp(A,B,AMAX,BMAX,TL,TL1)
Implicit Real*8(A-H,O-Z)
Real*8 A(48),B(48),AMAX,BMAX
Integer TL,TL1
AMAX = A(1)
BMAX = B(1)
J1 = 1
J2 = 1
DO I=2,48
IF(A(I).GT.AMAX) THEN
AMAX=A(I)
J1 = I
ENDIF
IF(B(I).GT.BMAX) THEN
BMAX=B(I)
J2 = I
ENDIF
ENDDO
TL = J1
TL1 = J2
END

!########################################################################
! Author: zhb
! Last modified time: 2011.8.18
! Task: 转换河道参数的格式（分离列）
Program separate
Integer o1,o2,num(3),counts(3)
Character head(3)*3,Ctype(3)*5
Real*8 Ta(12,2),Tp(12,2),Tt(12,2),Ta2(12,2),Tp2(12,2),Tp3(12,2)

o1=1
o2=2
open(unit=o1,file=’source.txt’,form=’formatted’,
$ status=’old’,access=’sequential’)
read(o1,1008) head(1),num(1),counts(1),Ctype(1)
CALL crackA(o1,Ta)
read(o1,1008) head(2),num(2),counts(2),Ctype(2)
CALL crackA(o1,Tp)
read(o1,1008) head(3),num(3),counts(3),Ctype(3)
CALL crackA(o1,Tt)
open(unit=o2,file=’out.txt’,form=’formatted’,
$ status=’replace’,access=’sequential’)
write(o2,'(A2)’) ‘X:’
DO I=1,12
write(o2,'(F16.8)’) Ta(I,1)
ENDDO
write(o2,*) ‘ ‘
write(o2,*) ‘ ‘
write(o2,'(A5)’) ‘Area:’
Do I=1,12
write(o2,'(F16.8)’) Ta(I,2)
ENDDO
write(o2,*) ‘ ‘
write(o2,*) ‘ ‘
write(o2,'(A)’) ‘Wetted Perimeter:’
Do I=1,12
write(o2,'(F16.8)’) Tp(I,2)
ENDDO
write(o2,*) ‘ ‘
write(o2,*) ‘ ‘
write(o2,'(A)’) ‘Top Width:’
Do I=1,12
write(o2,'(F16.8)’) Tt(I,2)
ENDDO

1008 Format(A3,1x,I3,1x,I2,1x,A5)
END
subroutine CrackI(I1,In,onec,TmpI,NTYPE)
Integer I1,I2,In,TmpI(10)
Character onec(130)*1,Blank*1,Ibuf*130,IFOR*10,NTYPE*9
Blank = ‘ ‘
IF(I1.GT.130) GOTO 299
In1 = In
Do 201 I=1,In1
202 Continue
IF(onec(I1).NE.Blank) goto 215
I1 = I1 + 1
IF(I1.ge.130) THEN
In = I – 1
goto 299
ENDIF
goto 202
215 continue
I2 = I1+1
221 continue
IF(onec(I1).eq.””) THEN
IF(onec(I2).eq.””) goto 225
Else
IF(onec(I2).eq.Blank) goto 225
ENDIF
I2 = I2+1
IF(I2.lt.131) goto 221
I2 = 130
225 continue
IF(onec(I2).eq.””) I1 = I1+1
Length = I2-I1
IF(Length.lt.1) Length = 1
I2 = I2-1

!########################################################################
! Author: zhb
! Last modified time: 2011.8.18
! Task: 修改二维参数中的降水曲线
Program modify
Integer o1,o2,RFTYPE(20),Npoint(20),TmpI(10),M,Elem(12209,3),
$ RN(20),RF(12209,2),RFC(12209,2),o3,o4,o5,seeH(20,2),
$ Itest(12209),RNC(16),RTYPEC(16) !RNC是修改后的RN
Character onec(130)*1,IC1*2,IC3*1,NTYPE*9,C1*6,C2*2,RFname*72,
$ XYname*72,post1(9)*9,post2(11)*10,form*24,TmpC*9,
& TmpC2*10,form2*36
Real*8 TmpD(10),RFXYS(20,48,2),XC(6312,3),RPmean,TmpR(48),
$ RFP(20),MaxRF,rate,t1,RFXYSC(20,48,2),Height(12209),X,Y,
$ R1,R2,RFPs(20),Srate,Maxh,Minh,sum,ssum,typeH(15),RFPC(16)
Logical Inquir
Common /H1/Height
Common /R1/RFC
Common /R2/RTYPEC
Common /R3/RFXYS
Common /R4/RFXYSC

o1 = 1
o2 = 2
o3 = 3
o4 = 4
o5 = 5
!读取原降水参数
open(unit=o1,file=’RF-Parameters.txt’,form=’Formatted’,
$ status=’old’,access=’sequential’)
Do J=1,20
read(o1,1111) IC1,IC3,(onec(I),I=1,130)
I1=1
In=2
NTYPE=’INTEGER ‘
Call CrackI(I1,In,onec,TmpI,NTYPE)
RFTYPE(J) = TmpI(1)
Npoint(J) = TmpI(2)
Do I=1,Npoint(J)
Read(o1,'(130A1)’) (onec(K),K=1,130)
I1=1
In=2
NTYPE = ‘REAL ‘
Call CrackD(I1,In,onec,TmpD,NTYPE)
RFXYS(J,I,1) = TmpD(1)
RFXYS(J,I,2) = TmpD(2)
ENDDO
ENDDO
!读取单元和节点坐标
open(unit=o2,file=’mesh_LY.2dm’,form=’formatted’,
& status=’old’,access=’sequential’)

read(o2,1112) C1
read(o2,'(A2)’) C2
read(o2,'(A2)’) C2
read(o2,'(A2)’) C2

Do J=1,12209
read(o2,1111) IC1,IC3,(onec(I),I=1,130)
I1 = 1
In = 5
Ntype = ‘INTEGER ‘
Call CrackI(I1,In,onec,TmpI,Ntype)
M = TmpI(1)
Elem(M,1) = TmpI(2)
Elem(M,2) = TmpI(3)
Elem(M,3) = TmpI(4)
ENDDO
Do J=1,6312
read(o2,1111) IC1,IC3,(onec(I),I=1,130)
I1 = 1
In = 1
Ntype = ‘INTEGER ‘
Call CrackI(I1,In,onec,TmpI,Ntype)
M = TmpI(1)
In = 3
Ntype = ‘REAL ‘
Call CrackD(I1,In,onec,TmpD,Ntype)
XC(M,1) = TmpD(1)
XC(M,2) = TmpD(2)
XC(M,3) = TmpD(3)
ENDDO

!计算各单元的高程
Do I=1,12209
I1 = Elem(I,1)
I2 = Elem(I,2)
I3 = Elem(I,3)
Height(I) = (XC(I1,3)+XC(I2,3)+XC(I3,3))/3.0
ENDDO

!计算各高度段有多少个点
Do I=1,20
seeH(I,1) = I*200
seeH(I,2) = 0
ENDDO
Maxh=Height(1)
Minh=Height(1)
Do I=1,12209
Do J=1,20
IF(Height(I).le.seeH(J,1)) THEN
seeH(J,2) = seeH(J,2)+1
ENDIF
ENDDO
IF(Height(I).ge.Maxh) Maxh=Height(I)
IF(Height(I).le.Minh) Minh=Height(I)
ENDDO

print *,'<=200: ‘,seeH(1,2)
Do I=2,20
I1=seeH(I-1,1)
I2=seeH(I,1)
I3=seeH(I,2)-seeH(I-1,2)
write(*,1122) I1,I2,I3,seeH(I,2)
1122 Format(I4,’-‘,I4,': ‘,I5,2X,I5)
ENDDO
write(*,1120) Maxh
1120 Format(‘Max Height: ‘,F7.2)
write(*,1119) Minh
1119 Format(‘Min Height: ‘,F7.2)
c read *,I
Inquir=.true.
Do 1200 while(Inquir)
Print *,’Do you want to inquire? 1:yes 0:no’
read *,J
IF(J.eq.1) THEN
print *,’p;ease input the height value’
read *,r
Inum=0
Do I=1,12209
IF(Height(I).le.r) Inum = Inum+1
ENDDO
print *,’The number of elements below the height is: ‘,Inum
print *,’ ‘
ELSE IF(J.ne.1) THEN
Inquir = .false.
ENDIF
1200 Continue
!读取并计算原来各降水类型各自包含的单元数
open(unit=o3,file=’RF1.txt’,form=’formatted’,
$ status=’old’,access=’sequential’)
Do I=1,20
RN(I)=0
ENDDO
Do I=1,12209
Itest(I)=0
ENDDO
Do I=1,12208
read(o3,1111) IC1,IC3,(onec(J),J=1,130)
I1=1
In=2
Ntype=’INTEGER ‘
CALL CrackI(I1,In,onec,TmpI,Ntype)
RF(I,1) = TmpI(1) !单元编号
RF(I,2) = TmpI(2) !降水类型
DO J=1,20
IF(RF(I,2).eq.RFTYPE(J)) THEN
RN(J) = RN(J)+1 !各降水类型的单元数
ENDIF
ENDDO
Itest(RF(I,1)) = 1
ENDDO
!找没有赋降水类型的单元号
Isum=0
Do I=1,12209
IF(Itest(I).eq.0) THEN
print *,’Losing element ID: ‘,I
Isum=Isum+1
ENDIF
ENDDO
print *,’Total: ‘,Isum
!计算平均峰值降水
Do J=1,20
Do I=1,48
TmpR(I) = RFXYS(J,I,2)
ENDDO
RFP(J) = MaxRF(TmpR)
print *,J,RFP(J)
RFPs(J) = RFP(J) !保存原降水峰值
ENDDO

Isum=0
Do I=1,20
Isum = Isum+RN(I)
ENDDO
print *,’Isum: ‘,Isum
sum=0.0
Do I=1,20
sum = sum+RN(I)*RFP(I)
ENDDO
RPmean = sum/Isum
print *,’RPmean: ‘,RPmean
c print *,’RN1 RN2 RN3 : ‘,RN1,RN2,RN3
c print *,’RFP1: ‘,RFP1
c print *,’RFP2: ‘,RFP2
c print *,’RFP3: ‘,RFP3
c read *,II
print *,’Please input the additive rate (unit: mm/100m/year)’
read *,rate
Srate = rate
print *,’ ‘
print *,’Study area is divided as: <=200,400,600,800,…,3000,>3000′
print *,’Total rainfall types: 16′

!生成分段高度
Do I=1,15
typeH(I)=I*200.0
ENDDO
!生成新的降水类型编号集合
DO I=1,16
RTYPEC(I) = RFTYPE(I)
ENDDO
!计算修改后的各类型单元数并给RFC赋值
CALL typeCounts(RNC,typeH)
Isum=0
Do I=1,16
print *,’RNC’,I,RNC(I)
Isum=RNC(I)+Isum
ENDDO
print *,’Total RNC: ‘,Isum

!转换rate为适合模型计算的单位
rate = (rate/365.0)/24.0/1000/3600
rate=(200.0/100.0)*rate
!计算修改后的各条曲线的峰值
Isum2=0
Do I=1,15
Isum2=I*RN(16-I)+Isum2
ENDDO
print *,’Isum2: ‘,Isum2
Y=RPmean
X=(Y*Isum+rate*Isum2)/Isum
print *,’X: ‘,X
RFPC(16) = X
Do I=1,15
RFPC(I) = RFPC(16) – (16-I)*rate
ENDDO

DO I=1,16
print *,’RFPC’,I,RFPC(I)
ENDDO
!拟合曲线
Call makeCurve(RFPs,RFPC)

!输出修改后的RF1
open(unit=o4,file=’RF1_MOD.txt’,form=’Formatted’,status=’replace’,
$ access=’sequential’)

IC1=’RF’
IC3=’1′
Do I=1,12209
IF(I.eq.473) goto 1322
I10=IntLength(RFC(I,1))
I20=IntLength(RFC(I,2))
write(form,1323) I10,I20
1323 Format(‘(A2,A1,3X,I’,I1,’,3X,I’,I1,’)’)
write(o4,form) IC1,IC3,RFC(I,1),RFC(I,2)
1322 continue
ENDDO
!输出修改后的降水参数RF_XYS
write(XYname,1123) Srate
1123 Format(‘RFXYS_rate’,ES8.2,’.txt’)
open(unit=o5,file=XYname,form=’formatted’,
$ status=’replace’,access=’sequential’)

IC1=’XY’
IC3=’S’

post1(1) = ‘”Elem._1″‘
I90 = ICHAR(‘1′)
Do I=2,9
TmpC=post1(1)
I91 = I90+I-1
Do J=1,9
IF(J.eq.8) THEN
post1(I)(J:J) = CHAR(I91)
goto 1998
ENDIF
post1(I)(J:J) = TmpC(J:J)
1998 continue
ENDDO
ENDDO
post2(1) = ‘”Elem._10″‘
I90 = ICHAR(‘0′)
Do I=2,10
TmpC2=post2(1)
I91 = I90+I-1
Do J=1,10
IF(J.eq.9) THEN
post2(I)(J:J) = CHAR(I91)
goto 1999
ENDIF
post2(I)(J:J) = TmpC2(J:J)
1999 continue
ENDDO
ENDDO
post2(11) = ‘”Elem._20″‘
Do I=1,16
IF(I.le.9) THEN
K1=IntLength(RFTYPE(I))
K2=IntLength(Npoint(I))
write(form2,1234) K1,K2
1234 Format(‘(A2,A1,1X,I’,I1,’,1X,I’,I1,’,1X,A9)’)
write(o5,form2) IC1,IC3,RFTYPE(I),Npoint(I),post1(I)
Else
K1=IntLength(RFTYPE(I))
K2=IntLength(Npoint(I))
write(form2,1235) K1,K2
1235 Format(‘(A2,A1,1X,I’,I1,’,1X,I’,I1,’,1X,A10)’)
write(o5,form2) IC1,IC3,RFTYPE(I),Npoint(I),post2(I-9)
END IF
Do J=1,48
Id=2
CALL getLength(RFXYSC(I,J,1),Iw,Id)
write(form,1125) Iw,Id
1125 Format(‘(F’,I1,’.’,I1,’,” ”,ES8.2)’)
write(o5,form) RFXYSC(I,J,1),RFXYSC(I,J,2)
ENDDO
ENDDO
Pause
1111 Format(A2,A1,130A1)
1112 Format(A6)
END
subroutine CrackI(I1,In,onec,TmpI,NTYPE)
Integer I1,I2,In,TmpI(10)
Character onec(130)*1,Blank*1,Ibuf*130,IFOR*10,NTYPE*9,Tab*1
Blank = ‘ ‘
Tab=CHAR(9)
IF(I1.GT.130) GOTO 299
In1 = In
Do 201 I=1,In1
202 Continue
IF(onec(I1).NE.Blank.AND.onec(I1).NE.Tab) goto 215
I1 = I1 + 1
IF(I1.ge.130) THEN
In = I – 1
goto 299
ENDIF
goto 202
215 continue
I2 = I1+1
221 continue
IF(onec(I1).eq.””) THEN
IF(onec(I2).eq.””) goto 225
Else
IF(onec(I2).eq.Blank.OR.onec(I2).eq.Tab) goto 225
ENDIF
I2 = I2+1
IF(I2.lt.131) goto 221
I2 = 130
225 continue
IF(onec(I2).eq.””) I1 = I1+1
Length = I2-I1
IF(Length.lt.1) Length = 1
I2 = I2-1

write(Ibuf,'(130A1)’) (onec(J),J=I1,I2)
IF(NTYPE(1:7).eq.’INTEGER’) THEN
WRITE(IFOR,1011) LENGTH
1011 FORMAT(‘(I’,I2,’)’)
Read(Ibuf,IFOR) TmpI(I)
c write(*,*) ‘TmpI’,I,TmpI(I)
c Read *,II
ELSE
goto 299
ENDIF
IF(onec(I2+1).EQ.””) I2=I2+1
I1 = I2+1
201 Continue
299 continue
Return
END
Subroutine CrackD(I1,In,onec,TmpD,NTYPE)
Implicit Real*8(A-H,O-Z)
Integer I1,I2,In
Character onec(130)*1,Blank*1,Ibuf*130,FFOR*10,NTYPE*9,Tab*1
Real*8 TmpD(10)
Blank = ‘ ‘
Tab=CHAR(9)
IF(I1.GT.130) GOTO 399
In1 = In
Do 301 I=1,In1
302 Continue
IF(onec(I1).NE.Blank.AND.onec(I1).NE.Tab) goto 315
I1 = I1 + 1
IF(I1.ge.130) THEN
IN = I – 1
goto 399
ENDIF
goto 302
315 continue
I2 = I1+1
321 continue
IF(onec(I1).eq.””) THEN
IF(onec(I2).eq.””) goto 325
Else
IF(onec(I2).eq.Blank.OR.onec(I2).eq.Tab) goto 325
ENDIF
I2 = I2+1
IF(I2.lt.131) goto 321
I2 = 130
325 continue
IF(onec(I2).eq.””) I1 = I1+1
Length = I2-I1
IF(Length.lt.1) Length = 1
I2 = I2-1

tmp = TmpR(1)
Do I=2,48
IF(tmp.le.TmpR(I)) tmp=TmpR(I)
ENDDO
MaxRF = tmp
END
Subroutine typeCounts(RNC,typeH)
Integer RFC(12209,2),RNC(16),RTYPEC(16)
REAL*8 typeH(15),Height(12209)
Common /H1/Height
Common /R1/RFC
Common /R2/RTYPEC
Do I=1,16
RNC(I) = 0
ENDDO
Do I=1,12209
IF(I.eq.473) goto 1333
IF(Height(I).le.typeH(1)) THEN
RNC(1)=RNC(1)+1
RFC(I,1)=I
RFC(I,2)=RTYPEC(1)
ELSE IF(Height(I).gt.typeH(15)) THEN
RNC(16) = RNC(16)+1
RFC(I,1)=I
RFC(I,2)=RTYPEC(16)
ELSE
Do I1=2,15
IF(Height(I).gt.typeH(I1-1).and.Height(I).le.typeH(I1)) THEN
RNC(I1) = RNC(I1)+1
RFC(I,1)=I
RFC(I,2)=RTYPEC(I1)
ENDIF
ENDDO
ENDIF
1333 continue
ENDDO
END
Subroutine makeCurve(RFPs,RFPC)
REAL*8 RFXYS(20,48,2),RFXYSC(20,48,2),RFPs(20),
$ RFPC(16),r(20),rsum,GR(48),Maxr,rt
Common /R3/RFXYS
Common /R4/RFXYSC

rsum=0.0
DO I=1,48
DO J=1,20
r(J) = RFXYS(J,I,2)/RFPs(J)
rsum=rsum+r(J)
ENDDO
rsum=rsum/20.0
GR(I)=rsum
ENDDO

Maxr=GR(1)
DO I=1,48
IF(Maxr.lt.GR(I)) Maxr=GR(I)
ENDDO
Do I=1,48
GR(I)=GR(I)/Maxr
print *,’GR’,I,GR(I)
ENDDO
Do I=1,48
rt=GR(I)
Do J=1,16
RFXYSC(J,I,1) = RFXYS(J,I,1)
RFXYSC(J,I,2) = rt*RFPC(J)
ENDDO
ENDDO
END
Subroutine getLength(F,w,d)
Real*8 F
character*10 form
character*20 tmp
Integer w,d
write(form,2199) d
2199 Format(‘(F20.’,I1,’)’)
write(tmp,form) F
Do I=1,20
IF(tmp(I:I).ne.’ ‘.and.I.ne.1) THEN
w=20-I+1
goto 2002
ENDIF
ENDDO
2002 continue
END
Integer Function IntLength(IP)
character*20 tmp
Integer width
write(tmp,'(I20)’) IP
Do I=1,20
IF(tmp(I:I).ne.’ ‘.and.I.ne.1) THEN
width=20-I+1
goto 2032
ENDIF
ENDDO
2032 continue
IntLength = width
END

Leave a Reply Cancel reply

Post Navigation