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
    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
    !
    !
    Kr = 0.0
    !
    wetup   = 0
    wetdown = 0
    !
    ! 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(i,j)+s%H(i,j)/2>par%eps) then
             wetup(i,j) = 1
          else
             wetup(i,j) = 0
          end if
          if(s%hh(i,j)-s%H(i,j)/2>par%eps) 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
    !
    if (s%H(bl(1,1),1)<=0.01) then
       kr(1) = 0.d0
    else
       iribar = slope/(s%H(bl(1,1)-1,1)*s%k(bl(1,1)-1,1)/(2*par%px))**0.5d0
       alfa   = 0.8
       beta   = 5
       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,Href)

    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
    integer                                                  :: i
    !
    !
    allocate(dist(bl(1,1),1))
    allocate(Htemp(bl(1,1),1))
    !
    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
             Href(bl(1,1)+1-i,1) = Htemp(i,1)
          else
             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((par%ny+1)*(par%nx+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

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

end module reflection_module