! *********************************************************
! **      CALCULATION FOR TIDAL RIVER                    **
! **                                                     **
! **                             '18.7.18                **
! **                                                     **
! **                            Coded by Koji ASAI       **
! **                            Yamaguchi Univ.          **
! *********************************************************

module comm
 implicit none

 real(8) :: dx   ! grid spacing(m)
 real(8) :: dt   ! time interval(sec)
 real(8) :: q0   ! unit width discharge(m2/s)
 real(8) :: s0   ! bed slop(-)
 real(8) :: h0   ! initial water depth(m)
 real(8) :: grv  ! gravity acceleration(m/s2)
 real(8) :: n    ! Manning's roughness parameter
 real(8) :: amp  ! tidal amplitude(m)
 real(8) :: t0   ! tidal period(sec)
 real(8) :: pi   ! the circular constant
 real(8) :: time ! time(sec)

 integer :: imax ! the number of grid
 integer :: per  ! the number o time step in one tide
 integer :: ends ! the number of end time step
 integer :: ebb  ! output in ebb tide
 integer :: low  ! output in low tide
 integer :: fld  ! output in flood tide
 integer :: hig  ! output in high tide
 integer :: k    ! time step
 integer :: i    ! grid number

 real(8), dimension(:), allocatable :: r1, r2, c, h, u, q, z, x
 real(8), dimension(:), allocatable :: r1old, r2old

end module comm


! --------------------------------------------------------------------
program main
 use comm
 implicit none

 call iniset

 do k=1, ends
   time=time+dt
   if(mod(k,100)==0) write(6,*) k, time
   call getuc
   call next
 end do

 write(6,*) '  CALCULATION HAS FINISHED SUCCESFULLY! '
 stop

end program main


! --------------------------------------------------------------------
subroutine iniset
 use comm
 implicit none

 dt=1.0d0
 dx=40.0d0
 n =0.03d0
 q0=0.005d0
 s0=0.0010d0
 h0=2.0d0
 amp=0.50d0
 t0=43200.0d0
 pi=4.0d0*atan(1.0d0)
 grv=9.8d0
 time=0.0d0
 imax=64

 per=int(t0/dt +0.001)
 ends=per *6
 hig=ends-per+per/4
 ebb=ends-per+per/4*2
 low=ends-per+per/4*3
 fld=ends


 allocate (h(0:imax+1))        ! water depth
 allocate (u(0:imax+1))        ! velocity
 allocate (q(0:imax+1))        ! discharge
 allocate (c(0:imax+1))        ! celerity
 allocate (z(0:imax+1))        ! bed elevation
 allocate (x(0:imax+1))        ! x-coordinate
 allocate (r1(0:imax+1))       ! u+c
 allocate (r2(0:imax+1))       ! u-c
 allocate (r1old(0:imax+1))    ! u+c at previous step
 allocate (r2old(0:imax+1))    ! u-c at previous step

! ----- clear -----
 do i=0,imax+1
  h(i)=0.0d0
  u(i)=0.0d0
  q(i)=0.0d0
  c(i)=0.0d0
  z(i)=0.0d0
  x(i)=0.0d0
  r1(i)=0.0d0
  r2(i)=0.0d0
  r1old(i)=0.0d0
  r2old(i)=0.0d0
 end do

! ----- x-coordinate -----
 x(1)=0.0d0
 do i=2,imax
  x(i)=x(i-1) +dx
 end do

! ----- bed elevation -----
 z(1)=0.0d0     ! assumed that this position is at a reference level (up stream side).
 do i=2,imax
  z(i)=z(i-1) -s0*dx
 end do

! ----- initial condition -----
 do i=1, imax
  h(i)=h0
  q(i)=q0
  u(i)=q0/h0
  c(i)=dsqrt(grv*h(i))
  r1old(i)=u(i)+2.0d0*c(i)
  r2old(i)=u(i)-2.0d0*c(i)
 end do

end subroutine iniset


! --------------------------------------------------------------------
subroutine getuc
 use comm
 implicit none
 real(8) :: r1foot, r2foot, alpha, beta, hold

 do i=2, imax-1
  alpha=(u(i)+c(i))*dt/dx
   if(dabs(alpha) >1.0d0) then
     write(6,*) k,i, ' alpha',alpha, u(i), c(i)
     stop
   end if
   if(alpha >=0.0)  r1foot=r1old(i)*(1.0d0-dabs(alpha)) +r1old(i-1)*dabs(alpha)
   if(alpha <0.0)   r1foot=r1old(i)*(1.0d0-dabs(alpha)) +r1old(i+1)*dabs(alpha)
  r1(i)=r1foot +grv*s0*dt -grv*n*n*dabs(u(i))*u(i)/(h(i)**(4./3.))*dt

  beta=(u(i)-c(i))*dt/dx
   if(dabs(beta) >1.0d0) then
     write(6,*) k,i,' beta=',beta, u(i),c(i)
     stop
   end if
   if(beta >0.0)  r2foot=r2old(i)*(1.0d0-dabs(beta)) +r2old(i-1)*dabs(beta)
   if(beta <=0.0) r2foot=r2old(i)*(1.0d0-dabs(beta)) +r2old(i+1)*dabs(beta)
  r2(i)=r2foot +grv*s0*dt -grv*n*n*u(i)*dabs(u(i))/(h(i)**(4./3.))*dt
 end do

! ------ u, c, h, q at new time step except both boundaries -----
 do i=2,imax-1
  u(i)=0.50d0*(r1(i)+r2(i))
  c(i)=0.25d0*(r1(i)-r2(i))
  h(i)=c(i)*c(i)/grv
  q(i)=u(i)*h(i)
 end do

! ----- boundary condition at upstream (i=1) -----
 hold=h(1)
 q(1)=q0
 h(1)=hold -dt/dx*(q(2)-q(1))
 u(1)=q(1)/h(1)
 c(1)=dsqrt(grv*h(1))

! ----- boundary condition at downstream (i=imax) -----
 hold=h(imax)
 h(imax)=h0 +amp*dsin(2*pi/t0*time)
 q(imax)=q(imax-1) -dx/dt*(h(imax)-hold)
 u(imax)=q(imax)/h(imax)
 c(imax)=dsqrt(grv*h(imax))

end subroutine getuc


! --------------------------------------------------------------------
subroutine next
 use comm
 implicit none

 do i=1,imax
  r1old(i)=u(i)+2.0d0*c(i)
  r2old(i)=u(i)-2.0d0*c(i)
 end do


 if(k == hig ) then
  open(16, file='hig_h.txt')
  open(26, file='hig_u.txt')
  do i=1,imax
   write(16,*) x(i), z(i), z(i)+h(i)
   write(26,*) x(i), u(i)
  end do
  close(16)
  close(26)
 end if

 if(k == ebb ) then
  open(16, file='ebb_h.txt')
  open(26, file='ebb_u.txt')
  do i=1,imax
   write(16,*) x(i), z(i), z(i)+h(i)
   write(26,*) x(i), u(i)
  end do
  close(16)
  close(26)
 end if

 if(k == low ) then
  open(16, file='low_h.txt')
  open(26, file='low_u.txt')
  do i=1,imax
   write(16,*) x(i), z(i), z(i)+h(i)
   write(26,*) x(i), u(i)
  end do
  close(16)
  close(26)
 end if

 if(k == fld ) then
  open(16, file='fld_h.txt')
  open(26, file='fld_u.txt')
  do i=1,imax
   write(16,*) x(i), z(i), z(i)+h(i)
   write(26,*) x(i), u(i)
  end do
  close(16)
  close(26)
 end if

end subroutine next