module gvariables implicit none real(8),parameter::pi=3.14159265358979323846D0 real(8),parameter::g=9.80665D0 !Gravity acceleration integer,parameter::NUMS=10 !Limit number of section definition points of Surge Tank integer,parameter::NUMQ=20 !Limit number of discharge definition points character(len=256)::Title !Title integer::ICT !Calculation case (1: normal, 2: AFC) real(8)::AFCA !Discharge of half amplitude of AFC(m3/s) real(8)::AFCT !Period of AFC(s) !Definition of times real(8)::TMAX !End time for calculation real(8)::dt !Time increment for calculation real(8)::DTWR !Time increment for print out !Definition of properties of orifice (port) real(8)::PAA !Section area of port real(8)::PCI !Discharge coefficient of port for in-flow real(8)::PCO !Discharge coefficient of port for out-flow !Definition of propertoes of pressure tunnel & reservior real(8)::RWL !Water level at reservoir (EL) real(8)::TNL !Length of pressure tunnel real(8)::TNA !Section area of presure tunnel real(8)::TNC !Head loss coefficient of pressure tunnel !Definition of properties of shaft integer::NST !Number of sections real(8)::SEL(NUMS+2) !Level definition of sections of shaft (EL) real(8)::SAA(NUMS+2) !Section area of specified level in the shaft character(len=50)::SLB(NUMS+1) !Label for section (text: for drawing) !Difinition of discharge integer::NQT !Number of discharge input points real(8)::QTQ(NUMQ+1) !Discharge at specified time point real(8)::QTI(NUMQ+1) !Time at specified time point !Others real(8)::AFCS !Starting time of interception during AFC operation real(8)::QI !Initial discharge character(len=50)::fnameR,fnameW !Input & output file name end module gvariables program f90_SURGE use gvariables implicit none integer::iret call getarg(1,fnameR) call getarg(2,fnameW) call DINP() !Data input call PRE() !Pre-treatment call CALC(iret) !Execution of calculation print "('End of calculation')" if(iret==0)print "(i0,': Succsessfull termination')",iret if(iret==1)print "(i0,': Water Level is greater than Top of the shaft')",iret if(iret==2)print "(i0,': Water Level is less than Bottom of the shaft')",iret stop contains subroutine DINP() integer::i open(11,file=fnameR,status='old') read(11,'(a)') Title read(11,*) ICT,AFCA,AFCT read(11,*) TMAX,dt,DTWR read(11,*) PAA,PCI,PCO read(11,*) RWL,TNL,TNA,TNC read(11,*) NST do i=1,NST read(11,*) SAA(i),SEL(i),SLB(i) end do read(11,*) NQT do i=1,NQT read(11,*) QTQ(i),QTI(i) end do close(11) end subroutine DINP subroutine PRE() integer::i,j real(8)::wa,wb character(len=50)::wc !QI : Initial discharge !AFCS : Starting time of interception (normal condition: AFCS=-1.0) !In normal condition, starting time of interception is set on QTI(1). !In AFC operation, starting time of interception is set on QTI(2). QI=QTQ(1) select case(ICT) case(1) AFCS=-1.0D0 !Normal case(2) AFCS=QTI(2) !AFC end select !To verify the order of section levels of shaft (from bottom to top) do i=1,NST do j=i+1,NST if(SEL(i)>=SEL(j))then !To sort section numbers in small order of elevation wa=SEL(j) wb=SAA(j) wc=SLB(j) SEL(j)=SEL(i) SAA(j)=SAA(i) SLB(j)=SLB(i) SEL(i)=wa SAA(i)=wb SLB(i)=wc end if end do end do SAA(NST+1)=SAA(NST)!Dummy for judgement of top level (section area) SEL(NST+1)=SEL(NST)!Dummy for judgement of top level (EL) end subroutine PRE subroutine CALC(iret) integer:: j,ic,np,npe,npn,npw,iw real(8) dto,qn,vi,zi,zo,tmi,tmx,zmi,zmx,dts,fa,dv,dz,zw,vw,dti,dtw integer,intent(out)::iret iret=0 open(12, file=fnameW, status='replace') !To print input data call WOUT1() np=int(DTWR/dt+0.00001D0) !Step number of print out dto=dt !Time increment qn=QI !Initial discharge vi=qn/TNA !Initial velocity of tunnel discharge zi=TNC*vi*vi !Initial head loss due to initial velocity if(vi<0.0D0)zi=-TNC*vi*vi !Definition of sign of initial head loss (direction of flow) zo=RWL-zi !Initial water level in Surge Tank !To print initial values (time=0) dts=0.0D0 zw=zo vw=vi zi=zi qn=qn call WOUT2(dts,zw,vw,qn) !To search the section number with initial water level do j=NST,1,-1 if(SEL(j)NST)iw=NST if(zw=np)then npw=0 zw=RWL-zi vw=vi call WOUT2(dts,zw,vw,qn) end if end do iret=iret close(12) end subroutine CALC subroutine RUNGE(qo,vo,zo,fa,dtt,dv,dz,dts) real(8),intent(inout)::qo real(8),intent(out)::dv,dz real(8),intent(in)::vo,zo,fa,dtt,dts real(8)::wrk,dz0,dv0,cf,dtw,v1,z1,dv1,dz1,dv2,dz2,dv3,dz3 wrk=0.0D0 dz0=dtt*FZ(vo,qo,fa) cf=PCI*PAA !Water level rising (in-flow from port) if(dz0>0.0D0)cf=PCO*PAA !Water level drop (out-flow from port) wrk=FK(vo,qo,cf) dv0=dtt*FV(zo,vo,wrk) v1=vo+0.5D0*dv0 z1=zo+0.5D0*dz0 dtw=dts-dtt*0.5D0 qo=QNCAL(dtw) wrk=FK(v1,qo,cf) dv1=dtt*FV(z1,v1,wrk) dz1=dtt*FZ(v1,qo,fa) v1=vo+0.5D0*dv1 z1=zo+0.5D0*dz1 wrk=FK(v1,qo,cf) dv2=dtt*FV(z1,v1,wrk) dz2=dtt*FZ(v1,qo,fa) v1=vo+dv2 z1=zo+dz2 dtw=dts qo=QNCAL(dtw) wrk=FK(v1,qo,cf) dv3=dtt*FV(z1,v1,wrk) dz3=dtt*FZ(v1,qo,fa) dv=(dv0+2.0D0*dv1+2.0D0*dv2+dv3)/6.0D0 dz=(dz0+2.0D0*dz1+2.0D0*dz2+dz3)/6.0D0 end subroutine RUNGE real(8) function FZ(v,q,fa) ! Calculation of dz real(8),intent(in)::v,q,fa FZ=(q-TNA*v)/fa end function FZ real(8) function FK(v,q,cf) ! Calculation of k real(8),intent(in)::v,q,cf FK=abs((TNA*v-q)/cf)*(TNA*v-q)/cf/2.0D0/g end function FK real(8) function FV(z,v,wrk) ! Calculation of dv real(8),intent(in)::z,v,wrk FV=(z-TNC*abs(v)*v-wrk)/TNL*g end function FV real(8) function QNCAL(tt) ! Calculation of discharge by interpolation method real(8),intent(in)::tt real(8)::qn integer::i if(tt>=AFCS)then do i=2,NQT if(QTI(i-1)