module reflection_module

! Module to account for wave reflection in increasing wave shear stress
! at the bank toe

contains

!==========================================================================!

! Check if the domain is completely submerged. In such a case reflection
! is not present (1D case).
!
subroutine submergence(s,par,wete1,subm)

    use params
    use spaceparams
    
    implicit none
    type(spacepars)                             :: s
    type(parameters)                            :: par
    integer                                     :: subm  ! [-] flag to check whether the domain is submerged
    integer, dimension(par%nx+1,par%ny+1)       :: wete1 ! [-] wet and dry points
    integer                                     :: n_tot,i,j
    !
    !
    ! Count submerged cells
    !
    subm  = 0
    n_tot = 0
    do i=1,par%nx+1
       do j=1,par%ny+1
          n_tot = n_tot+wete1(i,j)
       end do
    end do
    !
    ! Assign the flag for the submerged or dry domain
    !
    if (n_tot>=(par%nx-1)*(par%ny+1) .or. n_tot<=2*(par%ny+1)) then
       subm = 1
    else
       subm = 0
    end if

end subroutine submergence

!==========================================================================!

! Subroutine to identify the coordinates of the cells (dry at the waterline)
! (1D and 2D cases). 
!
subroutine boundline(s,par,wete1,bl)
    
    use params
    use spaceparams
    
    implicit none
    type(spacepars)                             :: s
    type(parameters)                            :: par
    integer, dimension(par%nx+1,par%ny+1)       :: wete1 ! [-] wet and dry points
    integer, dimension((par%nx+1)*(par%ny+1),2) :: bl    ! [-] coordinates of waterline, (point number,xcoord,ycoord)
    integer                                     :: i,j,itheta,k,sum
    !
    !
    bl = 0
    k  = 0
    if (par%ny+1==1) then
       !
       ! 1D case
       !
       j = 1
       do i=2,par%nx
          if (wete1(i,j)==1 .and. wete1(i+1,j)==0) then
             k = k+1
             bl(k,1) = i
             bl(k,2) = j
          end if
       end do
    else
       !
       ! 2D case
       !
       do j=2,par%ny
          do i=2,par%nx
             sum = wete1(i-1,j)+wete1(i+1,j)+wete1(i,j+1)+wete1(i,j-1)
             if (sum<4 .and. sum>0) then
                k = k+1
                bl(k,1) = i
                bl(k,2) = j
             end if
          end do
       end do
    end if
    
end subroutine boundline  

!==========================================================================!

! Determine the relection coefficient Kr based on the local bathymetry (1D case)
!
subroutine ref_par_1D(s,par,bl,Kr)
    
    use params
    use spaceparams
    
    implicit none
    type(spacepars)                                          :: s
    type(parameters)                                         :: par
    real*8,  dimension((par%ny+1)*(par%nx+1)),   intent(out) :: Kr          ! [-] reflection coefficient
    real*8                                                   :: alfa,beta   ! [-] parameters
    real*8                                                   :: slope       ! [-] tan(beta)
    real*8                                                   :: iribar      ! [-] local iribarrean number
    real*8                                                   :: wave_L      ! [m] local wave length
    real*8                                                   :: wave_L_temp ! [m] local wave length
    real*8                                                   :: myH         ! [m] local max wave height
    real*8                                                   :: eqH         ! [m] local equivalent wave height
    real*8                                                   :: eqk         ! [m] local equivalent wave number
    integer, dimension(par%nx+1,par%ny+1)                    :: wetup       ! [-] flag wet/dry for h+H/2
    integer, dimension(par%nx+1,par%ny+1)                    :: wetdown     ! [-] flag wet/dry for h-H/2
    real*8,  dimension(:), allocatable                       :: xx
    real*8,  dimension(:), allocatable                       :: yy
    integer, dimension((par%ny+1)*(par%nx+1),2), intent(in)  :: bl 
    integer                                                  :: i,j,nn,nl,counter
    !
    !
    Kr = 0.0
    !
    wetup   = 0
    wetdown = 0
    !
    ! Identify the maximum local wave height
    !
    wave_L      = 2.d0*par%px/s%k(bl(1,1)-1,bl(1,2))
    wave_L_temp = 0.0d0
    counter     = 0
    myH         = 0.0d0
    do while (wave_L>=wave_L_temp)
       if (bl(1,1)-counter<=1) then
          exit
       end if
       if (s%H(bl(1,1)-counter,bl(1,2))>myH) then
          myH = s%H(bl(1,1)-counter,bl(1,2))
       end if
       wave_L_temp = wave_L_temp+s%dsz(bl(1,1)-counter,bl(1,2))
       counter     = counter+1
    end do
    !
    ! Identify the cells used to determine the local slope of the domain.
    ! Cells included between wave crest and trough are selected.
    !
    do j=1,par%ny+1
       do i=1,par%nx+1
          if(s%hh(bl(1,1),bl(1,2))+s%zb(bl(1,1),bl(1,2))+myH/2>s%zb(i,j)) then
             wetup(i,j) = 1
          else
             wetup(i,j) = 0
          end if
          if(s%hh(bl(1,1),bl(1,2))+s%zb(bl(1,1),bl(1,2))-myH/2>s%zb(i,j)) then
             wetdown(i,j) = 1
          else
             wetdown(i,j) = 0
          end if
       end do
    end do
    !
    nn = sum(wetup)-sum(wetdown)+1
    if (nn<=4) then
       nn = 4
    end if
    !
    allocate(xx(nn))
    allocate(yy(nn))
    !
    ! Take the selected points above and below the waterline
    !
    do i=1,nn
       if (bl(1,1)-2+i>=par%nx+1) then
          xx(i) = s%x(par%nx+1,1)
          yy(i) = s%zb(par%nx+1,1)
       else
          xx(i) = s%x(bl(1,1)-2+i,1)
          yy(i) = s%zb(bl(1,1)-2+i,1)
       end if
    end do
    !
    ! Find the mean slope
    !
    call linreg_1D(xx,yy,slope,nn)
    !
    ! Find Iribarrean number and reflection coefficient.
    ! In this case H and k are calculated before the last wet cell
    !
    eqH = 0.d0
    eqk = 0.d0
    if (s%H(bl(1,1),1)<=0.01) then
       kr(1) = 0.d0
    else
       eqH = sum(s%H(bl(1,1)-counter:bl(1,1)-1,bl(1,2)))/counter
       eqk = sum(s%k(bl(1,1)-counter:bl(1,1)-1,bl(1,2)))/counter 
       !iribar = slope/(s%H(bl(1,1)-1,1)*s%k(bl(1,1)-1,1)/(2*par%px))**0.5d0
       iribar = slope/(eqH*eqk/(2*par%px))**0.5d0
       alfa   = 0.35
       beta   = 15
       Kr(1)  = alfa*iribar**2/(iribar**2+beta)
    end if
    !
    deallocate(xx)
    deallocate(yy)
    !
end subroutine ref_par_1D

!==========================================================================!

! Determine the relection coefficient Kr based on the local bathymetry (2D case)
!
subroutine ref_par_2D(s,par,Kr)
    
    use params
    use spaceparams
    
    implicit none
    type(spacepars)     :: s
    type(parameters)    :: par
    real*8, intent(out) :: Kr        ! [-] reflection coefficient
    real*8              :: alfa,beta ! [-] parameters
    
    !
    ! TO BE IMPLEMENTED
    !
    
end subroutine ref_par_2D

!==========================================================================!

! Subroutine to determine the value of the domain slope through a least
! square linear regression
!
subroutine linreg_1D(xx,yy,slope,nn)

   implicit none
   real*8, dimension(nn) :: xx,yy
   real*8, intent(out)   :: slope
   real*8                :: sumx,sumx2,sumxy,sumy,sumy2
   integer               :: i,nn
   !
   !
   sumx  = 0.0
   sumx2 = 0.0
   sumxy = 0.0
   sumy  = 0.0
   sumy2 = 0.0
   !
   do i=1,nn
      sumx  = sumx+xx(i)
      sumy  = sumy+yy(i)
      sumx2 = sumx2+xx(i)*xx(i)
      sumy2 = sumy2+yy(i)*yy(i)
      sumxy = sumxy+xx(i)*yy(i)
   end do
   !
   slope = (nn*sumxy-sumx*sumy)/(nn*sumx2-sumx**2)
   !
end subroutine linreg_1D

!==========================================================================!

! Subroutine to determine the reflected height distribution for the 2D case
!
subroutine distr_refl(s,par,bl,Kr)

    use params
    use spaceparams
    
    implicit none 
    
    type(parameters)                            :: par
    type(spacepars)                             :: s
    real*8,  dimension((par%ny+1)*(par%nx+1))   :: Kr
    integer, dimension((par%ny+1)*(par%nx+1),2) :: bl ! [-] coordinates of waterline
    integer                                     :: i,j,itheta,Nref
    !
    !
    if (par%ny+1==1) then
       i    = bl(1,1)
       j    = bl(1,2)
       Nref = floor(par%px/s%dtheta)
       do itheta = 1,s%ntheta
       
       ! TO BE COMPLETED
       
       ! TODO :change the difference has not to be done with respect to thetamean
          if (s%theta(itheta)<s%thetamean(i,j)+90*par%px/180 .and. s%theta(itheta)>s%thetamean(i,j)) then
             s%ee(i,j,Nref-itheta+1) = s%ee(i,j,itheta)*Kr(1)**2
          elseif  (s%theta(itheta)>s%thetamean(i,j)+270*par%px/180 .and. s%theta(itheta)<s%thetamean(i,j)+2*par%px) then
             s%ee(i,j,3*Nref-itheta+1) = s%ee(i,j,itheta)*Kr(1)**2
          end if
       end do
    end if

end subroutine distr_refl

!==========================================================================!

! Subroutine to calculate the reflected wave height
!
subroutine reflected_H(s,par,Kr,bl)

    use params
    use spaceparams
    
    implicit none 
    
    type(parameters)                                         :: par
    type(spacepars)                                          :: s
    real*8 , dimension((par%ny+1)*(par%nx+1))  , intent(in)  :: Kr    ! [-] reflection coefficient
    !real*8 , dimension(par%nx+1,par%ny+1)      , intent(out) :: Href  ! [m] reflected height component
    integer, dimension((par%ny+1)*(par%nx+1),2), intent(in)  :: bl    ! [-] coordinates of waterline
    real*8 , dimension(:,:), allocatable                     :: dist
    real*8 , dimension(:,:), allocatable                     :: Htemp
    !real*8                                                   :: k_avg
    integer                                                  :: i
    !
    !
    allocate(dist(bl(1,1),1))
    allocate(Htemp(bl(1,1),1))
    !
    !k_avg = sum(s%k*s%wetz)/sum(s%wetz)
    if (par%ny==0) then 
       !
       ! 1D case
       !
       do i=1,bl(1,1)
          !
          ! Create a vector with increrasing distance offshore from the waterline
          !
          dist(i,1)  = s%sdist(bl(1,1),1)-s%sdist(bl(1,1)+1-i,1)
          !
          !Htemp(i,1) = s%H(bl(1,1)+1-i,1)*Kr(1)*abs(cos(s%k(i,1)*dist(i,1)))*exp(-s%k(i,1)*dist(i,1)/(3.d0*par%px))
          ! Removed the cosine dependency
          !
          ! Value of reflected wave height from to waterline toward offshore direction
          ! 
          Htemp(i,1) = s%H(bl(1,1)+1-i,1)*Kr(1)*exp(-s%k(i,1)*dist(i,1)/(3.d0*par%px))
          !
       end do
       do i=1,par%nx+1
          !
          ! Assign the value Htemp to the vector containing the reflected wave height component
          ! all over the domain (of course it is zero shoreward the waterline)
          !                              
          if (i<=bl(1,1)) then
             s%Href(bl(1,1)+1-i,1) = Htemp(i,1)
          else
             s%Href(i,1) = 0.d0
          end if
       end do
       !   
    else
       !
       ! 2D case
       ! TO BE IMPLEMENTED
       !
    end if
    !   
    deallocate(dist)
    deallocate(Htemp)
    !
end subroutine reflected_H

!==========================================================================!

function reduced_lastwetcell(s,par,bl,myvar)
    !
    ! Subroutine to reduce the value of the variable myvar at last wet cell.
    ! Such value is weighted averaged on a wave length backward. In case the
    ! value of the variable at last wet cell is lower than the average, then
    ! its value is maintained
    !
    ! TODO: insert a criterion on bed slope to do it...
    !
    use params
    use spaceparams
    !
    implicit none 
    !
    !
    type(parameters)                                         :: par
    type(spacepars)                                          :: s
    real*8                                                   :: reduced_lastwetcell  ! [general] value of myvar to be substituted at the last wet cell
    real*8                                                   :: wave_L        ! [m] wave length
    real*8                                                   :: wave_L_temp   ! [m] "cumulative" wave length
    real*8,  dimension(par%nx+1,par%ny+1),       intent(in)  :: myvar         ! [general] variable to be reduced
    integer, dimension(1,2), intent(in)                      :: bl            ! [-] coordinates of waterline
    integer                                                  :: ncells        ! number of cells on which operate average
    integer                                                  :: counter
    !
    !
    wave_L      = 0.0d0
    wave_L_temp = 0.0d0
    ncells      = 0
    wave_L      = 2.d0*par%px/s%k(bl(1,1),bl(1,2))
    counter     = bl(1,1)
    do while (wave_L>=wave_L_temp)
       wave_L_temp = wave_L_temp+s%dsz(counter,bl(1,2))
       counter     = counter-1
       if (counter==0) then
          exit
       end if
    end do
    ncells              = bl(1,1)-counter
    reduced_lastwetcell = sum(myvar(bl(1,1)-ncells:bl(1,1)-1,bl(1,2)))/ncells
    if (myvar(bl(1,1),bl(1,2))<=reduced_lastwetcell) then
       reduced_lastwetcell = myvar(bl(1,1),bl(1,2))
    end if
    !
    !
end function reduced_lastwetcell

!==========================================================================!

subroutine check_wave_steepness(s,par)
    !
    ! Subroutine to check the stepèness of waves...
    !
    use params
    use spaceparams
    !
    implicit none 
    !
    !
    type(parameters)                     :: par
    type(spacepars)                      :: s
    real*8, dimension(par%nx+1,par%ny+1) :: wave_L
    real*8, dimension(par%nx+1,par%ny+1) :: H_corr
    real*8, dimension(par%nx+1,par%ny+1) :: urms_corr
    real*8, dimension(par%nx+1,par%ny+1) :: check_LT
    real*8, dimension(par%nx+1,par%ny+1) :: check_H
    integer                              :: i,j
    !
    ! Check wave steepness
    !
    wave_L    = 0.d0
    H_corr    = 0.d0
    urms_corr = 0.d0
    check_LT  = 0.d0
    check_H   = 0.d0
    do j=1,par%ny+1
       do i=1,par%nx+1
          wave_L(i,j) = 2.d0*par%px/s%k(i,j)
          check_H(i,j)  = s%H(i,j)/s%hh(i,j)
          check_LT(i,j) = s%H(i,j)*wave_L(i,j)**2/s%hh(i,j)**3
          if (check_LT(i,j)>26) then
             write(*,*) 'no linear theory'
          end if
          if (s%H(i,j)/wave_L(i,j)>0.142d0*tanh(s%k(i,j)*s%hh(i,j))) then
             !write(*,*) 'too steep wave'
             H_corr(i,j) = wave_L(i,j)*0.142d0*tanh(s%k(i,j)*s%hh(i,j))
             !s%H(i,j) = wave_L(i,j)*0.142d0*tanh(s%k(i,j)*s%hh(i,j))
             urms_corr(i,j) = par%px*H_corr(i,j)/par%Trep/(sqrt(2.d0)*sinh(s%k(i,j)*max(s%hh(i,j),par%delta*H_corr(i,j))))
          end if
       end do
    end do
    !
end subroutine check_wave_steepness

!==========================================================================!

subroutine increase_urms(par,s,bl,Kr,local_wetz)
   !
    use params
    use spaceparams
    !
    implicit none 
    !
    !
    type(parameters)                                        :: par
    type(spacepars)                                         :: s
    real*8 , dimension(:,:), allocatable                    :: dist
    real*8 , dimension(par%nx+1,par%ny+1)                   :: myurms1
    real*8 , dimension(par%nx+1,par%ny+1)                   :: myurms2
    real*8 ,                                     intent(in) :: Kr
    real*8                                                  :: k_avg
    integer, dimension(par%nx+1,par%ny+1)                   :: local_wetz
    integer, dimension((par%ny+1)*(par%nx+1),2), intent(in) :: bl    ! [-] coordinates of waterline
    integer                                                 :: i,j
    !
    allocate(dist(bl(1,1),1))
    k_avg   = sum(s%k*local_wetz)/sum(local_wetz)
    myurms1 = 0.d0
    myurms2 = 0.d0
    if (par%ny==0) then 
       !
       ! 1D case
       j = 1
       do i=1,bl(1,1)
          !
          ! Create a vector with increrasing distance offshore from the waterline
          dist(i,1)  = s%sdist(bl(1,1),1)-s%sdist(bl(1,1)+1-i,1)
       end do
       do i=1,bl(1,1)
          !
          ! Determine increased root mean squared velocity
          ! Three options:
          !     a) the |sin(ks)| function is applied to the overall value of urms
          !     b) the |sin(ks)| function is applied only to the reflected value of urms
          !     c) the |sin(ks)| function is not applied
          !     d) linear combination of exp(-ks) and exp(-ks)*|sin(ks)|
          !s%urms(i,j) = par%px * ( s%H(i,j)+s%Href(i,j) ) * abs(sin(s%k(i,j)*dist(bl(1,1)-i+1,j))) &
          !              /par%Trep/(sqrt(2.d0)*sinh(s%k(i,j)*max(s%hh(i,j),par%delta*s%H(i,j))))
          s%urms(i,j) = par%px * ( s%H(i,j) + s%Href(i,j)*abs(sin(s%k(i,j)*dist(bl(1,1)-i+1,j))) ) &
                        /par%Trep/(sqrt(2.d0)*sinh(s%k(i,j)*max(s%hh(i,j),par%delta*s%H(i,j))))
          !s%urms(i,j) = par%px*(s%H(i,j)+s%Href(i,j))/par%Trep/(sqrt(2.d0)*sinh(s%k(i,j)*max(s%hh(i,j),par%delta*s%H(i,j))))
          !
          !myurms1(i,j) = par%px*(s%H(i,j)+s%Href(i,j))/par%Trep/(sqrt(2.d0)*sinh(s%k(i,j)*max(s%hh(i,j),par%delta*s%H(i,j))))
          !myurms2(i,j) = par%px*(s%H(i,j)+s%Href(i,j))*abs(sin(s%k(i,j)*dist(bl(1,1)-i+1,j)))/par%Trep/(sqrt(2.d0)*sinh(s%k(i,j)*max(s%hh(i,j),par%delta*s%H(i,j))))
          !myurms1(i,j) = par%px*(s%H(i,j)+s%Href(i,j))/par%Trep/(sqrt(2.d0)*sinh(s%k(i,j)*max(s%hh(i,j),par%delta*s%H(i,j))))
          !myurms2(i,j) = par%px*(s%H(i,j)+s%Href(i,j))*abs(sin(k_avg*dist(bl(1,1)-i+1,j)))/par%Trep/(sqrt(2.d0)*sinh(s%k(i,j)*max(s%hh(i,j),par%delta*s%H(i,j))))
          !s%urms(i,j)  = myurms1(i,1)*(1-Kr) + myurms2(i,j)*Kr                       
       end do
    end if
    
    !
    deallocate(dist)
    !
    !
end subroutine increase_urms

!==========================================================================!

end module reflection_module