!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Copyright (C) 2011 UNESCO-IHE, WL|Delft Hydraulics and Delft University !
! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski,            !
! Jaap van Thiel de Vries, Robert McCall                                  !
!                                                                         !
! d.roelvink@unesco-ihe.org                                               !
! UNESCO-IHE Institute for Water Education                                !
! P.O. Box 3015                                                           !
! 2601 DA Delft                                                           !
! The Netherlands                                                         !
!                                                                         !
! This library is free software; you can redistribute it and/or           !
! modify it under the terms of the GNU Lesser General Public              !
! License as published by the Free Software Foundation; either            !
! version 2.1 of the License, or (at your option) any later version.      !
!                                                                         !
! This library is distributed in the hope that it will be useful,         !
! but WITHOUT ANY WARRANTY; without even the implied warranty of          !
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU        !
! Lesser General Public License for more details.                         !
!                                                                         !
! You should have received a copy of the GNU Lesser General Public        !
! License along with this library; if not, write to the Free Software     !
! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307     !
! USA                                                                     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module readbedcomposition_module

use bedcomposition_module
use handlearray_module

!---------Variables declaration---------
implicit none
private
type(bedcomp_data), save, pointer :: bc
logical,            save          :: varthtrlyr
public bc
public varthtrlyr
public readbedcomp
public bedcomp_updmorlyr
public bedcomp_init_dzbed
public bedcomp_init_pbbed
public bedcomp_avalanching_x
!public wave_impact
!---------------------------------------

contains

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

subroutine readbedcomp(s,par)
    !
    use params
    use spaceparams
    use readkey_module
    use filefunctions
    use logging_module
    use interp
    !
    implicit none
    type(parameters)                  :: par
    type(spacepars)                   :: s
    real*8                            :: thtrlyr     ! [m] transport layer thickness in case it is uniform over the domain
    real*8, dimension(:), allocatable :: logsedsig
    integer                           :: istat
    logical                           :: varthtrlyr  ! [T/F] flag: T --> non uniform thtrlyr; F --> uniform thtrlyr
    character(256)                    :: fname1, fname2
    !
    !
    if (trim(par%bedcompfile)/='') then
       call writelog('l','','--------------------------------')
       call writelog('l','','Initializing bed composition')
       call check_file_exist(par%bedcompfile)
       !
       ! readkey_dbl(<file>,<key>,<default>,<min>,<max>)
       ! readkey_int(<file>,<key>,<default>,<min>,<max>)
       ! readkey_name(<file>,<key>)
       !
       istat = 0
       !
       ! Allocate bed compisition structure
       if (istat == 0) allocate (bc         , stat = istat)
       if (istat == 0) allocate (bc%settings, stat = istat)
       if (istat == 0) allocate (bc%state   , stat = istat)
       if (istat == 0) allocate (bc%work    , stat = istat)
       if (istat == 0) allocate (bc%settings%morlyrnum, stat = istat) !  structure containing numerical settings
       if (istat /= 0) then
          call writelog('le','','ERROR: allocation of bed composition structure failed')
          return
       endif
       !
       !
       ! Read bed composition parameters from file
       bc%settings%thunlyr    = readkey_dbl(par%bedcompfile, 'ThUnLyr', 0.d01, 0.d0, 1.d0) !
       bc%settings%thlalyr    = readkey_dbl(par%bedcompfile, 'ThLaLyr', bc%settings%thunlyr , 0.d0, 1.d0) ! [m]  Maximum thickness of Lagrangian layers in case IUnderLyr=2 (default = ThUnLyr)
       bc%settings%theulyr    = readkey_dbl(par%bedcompfile, 'ThEuLyr', bc%settings%thunlyr , 0.d0, 1.d0) ! [m]  Maximum thickness of Eulerian layers in case IUnderLyr=2 (default = ThUnLyr)
       thtrlyr                = readkey_dbl(par%bedcompfile, 'ThTrLyr', 0.d01, 0.d0, 1.d0) ! [m]  thickness of trasport (active) layer; initially it is kept uniform 
       !
       ! TODO: find a way to read a vector of integers instead of reals
       !bc%settings%sedtyp     = readkey_intvec(par%bedcompfile, 'sedtyp',par%ngd,size(bc%settings%sedtyp), 0, 0, 2)  ! sediment type: 0=total/1=noncoh/2=coh
       bc%settings%idiffusion = readkey_int(par%bedcompfile, 'IDiffusion', 0, 0, 10)    ! [-]  Switch for diffusion mechanism (default = 0)
       bc%settings%iporosity  = readkey_int(par%bedcompfile, 'IPorosity', 1, 0, 2)      ! [-]  Switch for porosity mechanism (default = 0)
                                                                                        !      0: based on 'dry bed density'
                                                                                        !      1: linear
                                                                                        !      2: non-linear
       bc%settings%iunderlyr  = readkey_int(par%bedcompfile, 'IUnderLyr', 2, 1, 2)      ! [-]  Flag for underlayer concept 1=well mixed layer, 2=multiple layers
       bc%settings%ttlform    = readkey_int(par%bedcompfile, 'TTLForm', 1, 0, 10)       ! [-]  Flag for transport layer thickness formulation in case iundelyr=2
       bc%settings%keuler     = readkey_int(par%bedcompfile, 'KEuler', 0, 0, 10)        ! [-]  index of first Eulerian (i.e. non-moving) layer
                                                                                        !      2: standard Eulerian, only top layer moves with bed level
                                                                                        !      nlyr: fully Lagrangian (all layers move with bed level)
       bc%settings%nfrac      = readkey_int(par%bedcompfile, 'NFrac', 0, 0, 10)         ! [-]  number of sediment fractions
       bc%settings%mxnulyr    = readkey_int(par%bedcompfile, 'MxNuLyr', 0, 0, 10)       ! [-]  maximum number of underlayers including transport and base layers
       bc%settings%neulyr     = readkey_int(par%bedcompfile, 'NEuLyr', bc%settings%mxnulyr, 0, 10)  ! [-]  Maximum number of Eulerian layers in case IUnderLyr=2 (default = MxNULyr)
       bc%settings%nlalyr     = readkey_int(par%bedcompfile, 'NLaLyr', 0, 0, 10)        ! [-]  Maximum number of Lagrangian layers in case IUnderLyr=2 (default = 0)
       bc%settings%nlyr       = readkey_int(par%bedcompfile, 'NLyr', 0, 0, 10)          ! [-]  number of layers (transport + exchange + under layers)
       bc%settings%ndiff      = readkey_int(par%bedcompfile, 'NDiff', 0, 0, 10)         ! [-]  Number of diffusion coefficients in z-direction (default = 0)  
       bc%settings%nmlb       = readkey_int(par%bedcompfile, 'NMLB', 1, 0, 10)          ! [-]  start index of segments
       bc%settings%nmub       = readkey_int(par%bedcompfile, 'NMUB', 1, 0, 10)          ! [-]  nm end index
       bc%settings%updbaselyr = readkey_int(par%bedcompfile, 'UpdBaseLyr', 1, 1, 4)     ! [-]  Mechanism for computing composition of base layer      
       bc%settings%flufflyr   = readkey_int(par%bedcompfile, 'flufflyr', 0, 0, 1)       ! [-]  Flag to switch to fluf layer concept (1) or not(0) 
       bc%settings%exchlyr    = readkey_logical(par%bedcompfile, 'ExchLyr', .false.)    ! [T/F]  Switch for exchange layer (is not implemented, so default = false)
       varthtrlyr             = readkey_logical(par%bedcompfile, 'varThTrLyr', .false.) ! [T/F]  Switch for non uniform transport layer thickness
       !
       ! Allocate variance of distribution of each sediment class
       allocate(logsedsig(1:par%ngd))
       logsedsig              = readkey_dblvec(par%bedcompfile, 'logsedsig',par%ngd,size(logsedsig),0.d0,0.d0,0.d0) ! [m2] variance of the distribution of each sediment class
       !
       !
       ! read referred files
       !fname1 = readkey_name(par%bedcompfile, 'IniComp')
       !fname2 = readkey_name(par%bedcompfile, 'Diffusion') ![-]  Uniform diffusion coefficient (1 real) or diffusion file <*.ind> with spatially varying values and depth below bed level (1 string)
       ! TODO read these files, if necessary
       !
       ! Numerical settings
       bc%settings%morlyrnum%MinMassShortWarning = 0.d10  ! [m]  minimum erosion thickness for a shortage warning
       bc%settings%morlyrnum%MaxNumShortWarning  = 1      ! [m]  maximum number of shortage warnings remaining
       ! TODO read these values from file
       !
       ! Assign parameters for bed composition to s% and par% variables and vice-versa
       ! TODO: decide if we take all this parmeters from a single file
       par%ngd           = bc%settings%nfrac 
       bc%settings%nlyr  = bc%settings%nlalyr + bc%settings%neulyr + 1
       par%nd            = bc%settings%nlyr 
       par%nd_var        = bc%settings%keuler ! position of the variable layer which is located under the lagrangian layers and the transport layer
       !
       ! Allocate bed composition matrices
       allocate(bc%settings%sedtyp (bc%settings%nfrac                    ) ) ! [-]  sediment type: 0=total/1=noncoh/2=coh
       allocate(bc%settings%kdiff  (bc%settings%nlyr, (s%nx+1) * (s%ny+1)) ) ! [m2/s]  diffusion coefficients for mixing between layers
       allocate(bc%settings%phi    (bc%settings%nfrac                    ) ) ! [-]  D50 diameter expressed on phi scale
       allocate(bc%settings%rhofrac(bc%settings%nfrac                    ) ) ! [kg/m3]  density of fraction (specific density or including pores, at now we consider excluding pores)
       allocate(bc%settings%sigphi (bc%settings%nfrac                    ) ) ! [-]  standard deviation expressed on phi scale
       allocate(bc%settings%thexlyr(bc%settings%nfrac                    ) ) ! [m]  thickness of exchange layer
       allocate(bc%settings%thtrlyr(                  (s%nx+1) * (s%ny+1)) ) ! [m]  Thickness of the transport layer
       allocate(bc%settings%zdiff  (                  (s%nx+1) * (s%ny+1)) ) ! [m]  depth below bed level for which diffusion coefficients are defined
       !
       ! Initialize bed composition matrices
       bc%settings%sedtyp(1) = 2        ! from params.txt or another file ?
       bc%settings%sedtyp(2) = 0
       bc%settings%kdiff     = 0.d0
       bc%settings%phi       = 0.d0     ! phi is automatically calculated in setbedfracprop
       bc%settings%sigphi    = 0.d0     ! sigphi is automatically calculated in setbedfracprop
                                        ! from logsedsig that is assigned independently from the
                                        ! value of D50 and D90.
                                        ! TODO use handle_prob_distr to calculate logsedsig from
                                        ! the values of D50 and D90
       bc%settings%rhofrac   = par%rhos ! from params.txt
                                        ! TODO provide a file with sediment properties ?
       bc%settings%thexlyr   = 0.d0     ! Thickness of exchange layer (not activated)
       bc%settings%thtrlyr   = thtrlyr  ! Thickness of transport layer (chosen uniform)
       bc%settings%zdiff     = 0.d0
       !bc%settings%bfluff0 = 0.d0
       !bc%settings%bfluff1 = 0.d0
       !
       ! Allocation of the bc%state variables
       allocate(bc%state%bodsed  (bc%settings%nfrac,                   (s%nx+1) * (s%ny+1)) ) ! [kg/m2] Array with total sediment per cell per fraction
       allocate(bc%state%dpsed   (                                     (s%nx+1) * (s%ny+1)) ) ! [m] Total depth sediment layer
       allocate(bc%state%msed    (bc%settings%nfrac, bc%settings%nlyr, (s%nx+1) * (s%ny+1)) ) ! [kg/m2] composition of morphological layers: mass of sediment fractions
       allocate(bc%state%preload (                   bc%settings%nlyr, (s%nx+1) * (s%ny+1)) ) ! historical largest load, units : kg
       allocate(bc%state%sedshort(bc%settings%nfrac,                   (s%nx+1) * (s%ny+1)) ) ! sediment shortage in transport layer, units : kg /m2
       allocate(bc%state%svfrac  (                   bc%settings%nlyr, (s%nx+1) * (s%ny+1)) ) ! [-] 1 - porosity coefficient
       allocate(bc%state%thlyr   (                   bc%settings%nlyr, (s%nx+1) * (s%ny+1)) ) ! [m] thickness of morphological layers
       allocate(bc%state%mfluff  (bc%settings%nfrac,                   (s%nx+1) * (s%ny+1)) ) ! [m] thikness of fluff layer (mass ?)
       allocate(bc%state%bfluff0 (bc%settings%nfrac,                   (s%nx+1) * (s%ny+1)) ) ! [kg/m2/s] burial coefficient 1 
       allocate(bc%state%bfluff1 (bc%settings%nfrac,                   (s%nx+1) * (s%ny+1)) ) ! [1/s] burial coefficient 2 
       !
       ! Initialize bc%state variables
       bc%state%bodsed   = 0.d0
       bc%state%dpsed    = 0.d0
       bc%state%msed     = 0.d0
       bc%state%preload  = 0.d0
       bc%state%sedshort = 0.d0
       bc%state%svfrac   = 1-par%por ! porosity is initially chosen uniform everywhere
       bc%state%thlyr    = 0.d0
       bc%state%mfluff   = 0.d0
       bc%state%bfluff0  = 0.d0
       bc%state%bfluff1  = 0.d0
       !
       ! Set sediment fraction properties: we obtain values for phi and sigphi
       call setbedfracprop(bc, bc%settings%sedtyp, par%D50, logsedsig, bc%settings%rhofrac)
       deallocate(logsedsig)
       !
    endif

end subroutine readbedcomp

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

! BAMJ: Subroutine to update the bed composition using the bedcomp approach
! by means of updmorlyr function. Use of allocatable variables instead of 
! pointers because variables change shape

subroutine bedcomp_updmorlyr(s,par,dzg,i,j)
    !
    use params
    use spaceparams
    use message_module
    use bedcomposition_module
    !
    implicit none
    type(spacepars)              :: s
    type(parameters)             :: par
    type(message_stack), pointer :: messages   ! message stack
    real*8, dimension(1:par%ngd) :: dzg        ! [m] bed level variation per sed. class => [kg/m2] as input for updmorlyr  
    real*8, dimension(1)         :: dzb        ! [m] total bed level variation
    real*8                       :: mass_tot
    integer                      :: jd,jg,i,j
    character(message_len)       :: message
    !
    !
    allocate(messages)
    !
    bc%settings%nmlb = (par%nx+1)*(j-1)+i
    bc%settings%nmub = bc%settings%nmlb
    !
    call initstack(messages) ! BAMJ: no clue what it means, copied from example.f90
    !
    ! Input sediment variation in [m] per fraction (excluding pores) is transformed in [kg/m2] per fraction  
    do jg = 1,par%ngd
       dzg(jg) = dzg(jg)*bc%settings%rhofrac(jg) ! from [m] to [kg/m2] (pores are excluded)
    end do
    !
    ! This is the update of the bed composition
    message = 'updating bed composition'
    if (updmorlyr(bc, dzg, dzb, messages) /= 0) call adderror(messages, message)
    !
    ! Assign bc% variables to variables with structure consistent with s% variables
    do jd=1,par%nd
       s%dzbed(i,j,jd) = bc%state%thlyr(jd,bc%settings%nmlb)
       mass_tot        = sum(bc%state%msed(:,jd,bc%settings%nmlb))
       do jg=1,par%ngd
          s%pbbed(i,j,jd,jg) = bc%state%msed(jg,jd,bc%settings%nmlb)/mass_tot
       end do
    end do
    !
    ! Update bed elevation
    s%zb(i,j)          = s%zb(i,j)+dzb(1)
    s%dzbdt(i,j)       = s%dzbdt(i,j)+dzb(1)
    s%sedero(i,j)      = s%sedero(i,j)+dzb(1)
    s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)+dzb(1))
    !
    deallocate(messages)
    !
end subroutine bedcomp_updmorlyr

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

! Subroutine to set s%dzbed in case of bedcomp approach
! At the moment we are using:
! bc%settings%thtrlyr = thickness of top layer
! par%dzg1            = thickness of Lagrangian layer
! par%dzg2            = thickness of Eulerian layer
! par%dzg3            = thickness of Eulerian layers under first transport layer
! NOTE: thlalyr and theulyr represnts the maximum thickness respectively of
! Lagrangian and Eulerian layers: first one can replace par%dzg1
! TODO: decide if these informations have to be provided by bedcomp.txt
!
subroutine bedcomp_init_dzbed(s,par)
    !
    use params
    use spaceparams
    !
    implicit none
    type(spacepars)                                 :: s
    type(parameters)                                :: par
    real*8, dimension(par%nx+1,par%ny+1)            :: thtrlyr_rank2
    real*8, dimension((par%nx+1)*(par%ny+1),par%nd) :: dzbed_rdc_bedcomp ! [m] layer thickness ()
    integer                                         :: cellnumber
    logical, save                                   :: varthtrlyr ! this save is necessary to maintain the value of the variable
    !
    !   
    cellnumber = (par%nx+1)*(par%ny+1)
    bc%settings%nmlb = 1
    bc%settings%nmub = cellnumber
    !
    !allocate(thtrlyr_rank2(1:par%nx+1,1:par%ny+1))
    !allocate(dzbed_rdc_bedcomp(1:cellnumber,1:par%nd))
    !
    ! Initialize transport layer thickness in case it is non uniform
    if (varthtrlyr) then
       call init_thtrlyr(par)      
    end if
    !
    thtrlyr_rank2 = reshape(bc%settings%thtrlyr,shape=(/par%nx+1,par%ny+1/)) ! rank 1-->2
    !
    ! Transport (top) layer thickness
    s%dzbed(:,:,1) = thtrlyr_rank2
    !
    ! Set Lagrangian layers thickness
    if (bc%settings%nlalyr/=0) then
       s%dzbed(:,:,2:bc%settings%nlalyr) = par%dzg1
    end if
    !
    ! Nominal thickness of first (Eulerian) variable layer
    s%dzbed(:,:,bc%settings%keuler) = par%dzg2
    !
    ! Nominal thickness of (Eulerian) underlayers
    s%dzbed(:,:,bc%settings%keuler+1:bc%settings%nlyr) = par%dzg3
    !
    call handlearray(s%dzbed,dzbed_rdc_bedcomp) ! rank 3-->2
    !
    bc%state%thlyr = transpose(dzbed_rdc_bedcomp) ! (nm,jd) --> (jd,nm)
    !
    !deallocate(thtrlyr_rank2)
    !deallocate(dzbed_rdc_bedcomp)
    !
end subroutine bedcomp_init_dzbed

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

! BAMJ: subroutine to initialize nonuniform transport layer thickness
!
subroutine init_thtrlyr(par)
   !
   use params
   use readkey_module
   use logging_module
   use bedcomposition_module
   !
   implicit none
   type(parameters)              :: par
   integer                       :: pos, ier
   real*8, dimension(:), pointer :: thtrlyr
   character(256)                :: fnameTLThick
   !
   !
   fnameTLThick =  readkey_name(par%bedcompfile, 'TLThick') ! it is necessary a file with transport layer thicknesse for each cell
   thtrlyr      => bc%settings%thtrlyr
   !
   open(10,file=fnameTLThick)
   do pos = 1,(par%nx+1)*(par%ny+1)
      read(10,*,iostat=ier)(thtrlyr(pos))
      if (ier .ne. 0) then
         call report_file_read_error(fnameTLThick)
      end if
   end do
   close(10)      
   !
end subroutine init_thtrlyr

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

! BAMJ: subroutine to set bc%state%msed (mass per unit surface area of the
! jg-th sediment class [kg/m2]) in case of bedcomp approach using data 
! provided by s%pbbed (ratio between mass jg-th sediment class and total mass [-])
!
subroutine bedcomp_init_pbbed(s,par)
    !
    use params
    use spaceparams
    !
    implicit none
    type(spacepars)                                         :: s
    type(parameters)                                        :: par
    real*8, dimension((par%nx+1)*(par%ny+1),par%nd,par%ngd) :: pbbed_rdc_bedcomp ! [-] "middle variable" for sed. fractions (nm,jd,jg)
    real*8                                                  :: rhom,invrhom,mass_tot
    integer                                                 :: i,j,jd,jg,pos
    integer                                                 :: cellnumber
    !
    !
    cellnumber = (par%nx+1)*(par%ny+1)
    bc%settings%nmlb = 1
    bc%settings%nmub = cellnumber
    !
    !allocate(pbbed_rdc_bedcomp(1:cellnumber,1:par%nd,1:par%ngd))
    !
    call handlearray(s%pbbed,pbbed_rdc_bedcomp) ! rank 4-->3
    !
    ! Pass from mass fraction [-] to mass [kg/m2]
    do pos=1,cellnumber
       do jd=1,par%nd
          rhom    = 0.0d0 ! averaged density
          invrhom = 0.0d0 ! inverse of average density
          do jg=1,par%ngd
             invrhom = invrhom+pbbed_rdc_bedcomp(pos,jd,jg)/bc%settings%rhofrac(jg)
          end do
          rhom = 1/invrhom
          mass_tot = rhom*bc%state%thlyr(jd,pos)*bc%state%svfrac(jd,pos) ! total mass
          do jg=1,par%ngd
             bc%state%msed(jg,jd,pos) = pbbed_rdc_bedcomp(pos,jd,jg)*mass_tot
          end do
       end do
    end do        
    !
    !deallocate(pbbed_rdc_bedcomp)
    !
end subroutine bedcomp_init_pbbed

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

subroutine bedcomp_avalanching_x(s,par,edg2,edg1,ie,id,j,jdz)
    !
    use params
    use spaceparams
    use bedcomposition_module
    use message_module
    !
    implicit none
    !
    type(spacepars),target       :: s
    type(parameters)             :: par
    type(message_stack), pointer :: messages ! message stack
    real*8, dimension(1,1)       :: dz_bedcomp
    real*8, dimension(par%ngd)   :: edg1,edg2
    real*8                       :: mass_tot
    integer                      :: ie,id,j,jdz
    integer                      :: jg
    character(message_len)       :: message
    !
    !
    allocate(messages)
    !
    ! Pass from [m] to [kg/m2]
    do jg=1,par%ngd
       edg2(jg) = edg2(jg)*bc%settings%rhofrac(jg)*par%dt
       edg1(jg) = edg1(jg)*bc%settings%rhofrac(jg)*par%dt
    end do
    !                     
    ! Update bed in eroding point
    !
    bc%settings%nmlb = (j-1)*(par%nx+1)+ie
    bc%settings%nmub = (j-1)*(par%nx+1)+ie
    !                     
    message = 'updating bed composition'
    if (updmorlyr(bc,-edg2,dz_bedcomp,messages) /= 0) call adderror(messages, message)
    mass_tot            = sum(bc%state%msed(:,jdz,bc%settings%nmub))
    s%pbbed(ie,j,jdz,:) = bc%state%msed(:,jdz,bc%settings%nmub)/mass_tot
    s%dzbed(ie,j,jdz)   = bc%state%thlyr(jdz,bc%settings%nmub)
    !                     
    s%zb(ie,j)          = s%zb(ie,j)+dz_bedcomp(1,1)
    s%dzbdt(ie,j)       = s%dzbdt(ie,j)+dz_bedcomp(1,1)
    s%sedero(ie,j)      = s%sedero(ie,j)+dz_bedcomp(1,1)
    s%structdepth(ie,j) = max(0.d0,s%structdepth(ie,j)+dz_bedcomp(1,1))
    !                    
    ! Update bed in deposition point
    !
    bc%settings%nmlb = (j-1)*(par%nx+1)+id
    bc%settings%nmub = (j-1)*(par%nx+1)+id
    !                     
    message = 'updating bed composition'
    if (updmorlyr(bc,-edg1,dz_bedcomp,messages) /= 0) call adderror(messages, message)
    mass_tot            = sum(bc%state%msed(:,jdz,bc%settings%nmub))
    s%pbbed(id,j,jdz,:) = bc%state%msed(:,jdz,bc%settings%nmub)/mass_tot
    s%dzbed(id,j,jdz)   = bc%state%thlyr(jdz,bc%settings%nmub)
    !                     
    s%zb(id,j)          = s%zb(id,j)+dz_bedcomp(1,1)
    s%dzbdt(id,j)       = s%dzbdt(id,j)+dz_bedcomp(1,1)
    s%sedero(id,j)      = s%sedero(id,j)+dz_bedcomp(1,1)
    s%structdepth(id,j) = max(0.d0,s%structdepth(id,j)+dz_bedcomp(1,1))
    !
    deallocate(messages)
    !
end subroutine bedcomp_avalanching_x

!!==========================================================================!
!!
!! Subroutine to determine wave erosion due to the impact of waves on a steep slope
!! based on Mariotti and Fagherazzi (2010) or Roelvink and Reniers (2012) or new approach
!!
!subroutine wave_impact(s,par,ero_dry,sed_input)
!    !
!    use params
!    use spaceparams
!    !use bedcomposition_module
!    use message_module
!    use handlearray_module
!    use sediment_basics_module
!    use reflection_module
!    !
!    implicit none
!    type(parameters)                    :: par
!    type(spacepars)                     :: s
!    !type(bedcomp_data)                  :: bc
!    real*8, dimension(par%nx+1,par%ny+1,par%ngd), intent(out) :: ero_dry      ! [m3/m2/s] erosion erosion rate from the dry cells
!    real*8, dimension(par%nx+1,par%ny+1,par%ngd), intent(out) :: sed_input    ! [m3/m2/s] flux of material to fluid phase coming from dry cells
!    real*8, dimension(par%ngd)                                :: temp_ero_dry ! [m3/m2/s] total eroded material from dry cells per number of fractions
!    real*8                                                    :: impact_ero   ! [m3/m2/s] eroded mass due to wave impact
!    real*8                                                    :: Pcr          ! [W/m] critical wave power threshold for erosion
!    real*8                                                    :: Pw           ! [W/m] wave enery flux
!    real*8                                                    :: Pr           ! [W/m] roller energy flux
!    real*8, dimension(2,par%ngd)                              :: vfrac        ! [] volume fractions
!    real*8                                                    :: betaMF10     ! [kg/J] "erodibility" parameter from Mariotti and Fagherazzi (2010)
!    real*8                                                    :: betaMal11    ! [m3/J] "erodibility" parameter from Marani et al. (2011)
!    real*8                                                    :: calbeta      ! [-] calibration factor for the erodibility parameter
!    real*8, dimension(2)                                      :: pm           ! [-] fraction of cohesive sediment
!    real*8, dimension(2)                                      :: weight       ! [-] weight function for erosion based on water depth
!    real*8, dimension(par%nx+1,par%ny+1)                      :: hloc
!    real*8, dimension(par%nx+1,par%ny+1)                      :: dzbdxup
!    real*8                                                    :: ee           ! [m3/m/s] erosion rate expressed in meter^2 of sediment (without pores) per second;
!                                                                              ! m^2 are referred as a dx*dz quantity
!    real*8                                                    :: deltaz
!    real*8                                                    :: cell_ero
!    real*8                                                    :: courcheck,deltat
!    integer, dimension(par%nx+1)                              :: cellscarp    ! [-] number of cells composing a scarp
!    integer, dimension(par%nx+1)                              :: toe          ! [-] First cell being part of a scarp
!    integer, dimension(par%nx+1)                              :: top          ! [-] Last cell being part of a scarp
!    integer                                                   :: counter      ! [-] counter of the number of scarps
!    integer                                                   :: myscarp      ! [-] identifier of the scarp at waterline
!    integer                                                   :: subm
!    integer                                                   :: i,j,jg,ii,jj,kk
!    integer, dimension(par%nx+1,par%ny+1)                     :: sl
!    !
!    !                 
!    !allocate(impact_ero(2,bc%settings%nfrac))
!    !
!    hloc           = max(s%hh,par%eps)
!    Pcr            = 0
!    betaMF10       = 1.4E-04
!    betaMal11      = 1.2E-09
!    calbeta        = 10
!    ero_dry(:,:,:) = 0.d0
!    pm(:)          = 0.d0
!    sl             = 0
!    top            = 0
!    toe            = 0
!    cellscarp      = 0
!    dzbdxup        = 0.d0
!    counter        = 0
!    myscarp        = 0
!    !
!    call submergence(s,par,s%wetz,subm)
!    !
!    ! Transformation from [m3/m2/s] to [kg/m2/s]
!    !
!    do jg=1,bc%settings%nfrac
!       s%ero(:,:,jg) = s%ero(:,:,jg)*bc%settings%rhofrac(jg)
!    end do
!    !
!    ! Check if the domain is completely submerged or dry
!    !
!    if (subm==0) then
!       do j=1,par%ny+1
!          do i=1,par%nx-1
!             !
!             ! Identify the cell at the waterline
!             ! TODO: ACCOUNT FOR THE CASE THERE IS A TIDAL POOL "SEPARATING" FLUID PHASE
!             if (s%wetz(i,j)==1 .and. s%wetz(i+1,j)==0) then
!                !
!                ! Calculate wave power at the waterline
!                ! TODO: CORRECT FLUXES
!                Pw = s%E(i,j)*s%cg(i,j)
!                Pr = s%R(i,j)*s%c(i,j)
!                ee = max(0.d0,Pw+Pr-Pcr)*betaMal11*(1-par%por)*calbeta
!                !
!                ! If no erosion occurs, skip the following commands
!                !
!                if (ee>0.d0) then
!                   !
!                   ! Determine the retreat of the scarp
!                   !
!                   if (par%retreat==1) then
!                      !
!                      ! Retreat approach in which two cells at waterline are eroded
!                      ! proportionally to the water depth at the last wet cell
!                      !
!                      call retreat_scheme_1(s,par,hloc,ee,ero_dry,temp_ero_dry,i,j)
!                      ! 
!                      !
!                   else
!                      !
!                      ! Retreat approach in which erosion is proportional to 
!                      ! the domain slope in order to represent the horizontal
!                      ! translation of the bank scarp.
!                      !
!                      call identify_scarp(s,par,cellscarp,toe,top,sl,counter)
!                      !
!                      !
!                      if (sum(cellscarp)==0) then
!                         ! 
!                         ! In case no cell has slope >= scarpsl it is used the approach of retreat 1
!                         ! with no weight based on water depth
!                         !
!                         go to 10 
!                         !call retreat_scheme_1(s,par,hloc,ee,ero_dry,temp_ero_dry,i,j)
!                         ! PERHAPS IT IS NOT A GOOD IDEA TO CALL THE FUNCTION IN CASE OF NO SCARP...
!                         !
!                         !
!                      else
!                         !
!                         do kk=1,counter
!                            if (i>=toe(kk) .and. i<=top(kk)) then
!                                !
!                                ! Retreat scheme is only applied to the scarp intersecting the waterline 
!                                !
!                                call retreat_scheme_2(s,par,ee,ero_dry,temp_ero_dry,toe(kk),top(kk),i,j)
!                                !
!                            end if
!                         end do
!                      end if ! if statement for cellscarp = 0
!                   end if    ! if statement for retreat 1 or 2
!                   !
!                   ! Calculate flux of material eroded from dry cells as input for adv. diff. eq.;
!                   ! material is not added at a single cell but distributed over cells
!                   !
!                   call distr_ero_dry(s,par,temp_ero_dry,sed_input,i,j)
!                   !
!                   !
!                end if       ! if statement for occurrence of erosion               
!             end if          ! if to identify waterline
!          end do             ! do in x direction
!       end do                ! do in y direction
!    end if                   ! if statement for submerged/dry domain
!    !
!    ! Transformation from [kg/m2/s] to [m3/m2/s]
!    !
!10  do jg=1,bc%settings%nfrac
!       s%ero(:,:,jg)  = s%ero(:,:,jg)/bc%settings%rhofrac(jg)
!       ero_dry(:,:,jg) = ero_dry(:,:,jg)/bc%settings%rhofrac(jg)
!    end do
!    !
!end subroutine wave_impact
!
!!==========================================================================!
!
!subroutine retreat_scheme_1(s,par,hloc,ee,ero_dry,temp_ero_dry,i,j)
!    !
!    use params
!    use spaceparams
!    !
!    implicit none
!    type(parameters)                                          :: par
!    type(spacepars)                                           :: s
!    
!    real*8, dimension(par%nx+1,par%ny+1,par%ngd), intent(out) :: ero_dry      ! [m3/m2/s] eroded material from the dry cells
!    real*8, dimension(par%ngd),                   intent(out) :: temp_ero_dry ! [m3/m2/s] total eroded material from dry cells per number of fractions
!    real*8                                                    :: impact_ero   ! [m3/m2/s] erosion rate due to wave impact
!    real*8                                                    :: wave_L       ! [m] wave length
!    real*8                                                    :: wave_L_temp  ! [m] "cumulative" wave length
!    real*8                                                    :: x_temp 
!    real*8, dimension(2)                                      :: weight       ! [-] weight function for erosion based on water depth
!    real*8, dimension(par%nx+1,par%ny+1),         intent(in)  :: hloc
!    real*8,                                       intent(in)  :: ee           ! [m3/m/s] erosion rate expressed in meter^2 of sediment (without pores) per second;
!                                                                              ! m^2 are referred as a dx*dz quantity
!    integer                                                   :: i,j,jg,kk
!    integer                                                   :: ncells
!    integer                                                   :: counter
!    !
!    !
!    ! Determine erosion rate in meter of sediment (no pores) per second
!    !
!    impact_ero = ee/(s%dsz(i,j)+s%dsz(i+1,j))
!    !
!    ! Split erosion rate            
!    !
!    weight(2) = (hloc(i,j))/(max(s%zb(i+1,j)-s%zb(i,j),par%eps))
!    weight(1) = abs(1-weight(2))
!    !
!    if ((weight(1)+weight(2))>1.1 .or. (weight(1)+weight(2))<0.9) then
!        !write(*,*) 'unconsistent impact erosion weight'
!        weight(1) = 0.5
!        weight(2) = 0.5
!    end if
!    !
!    ! Compute erosion rate per cell and store erosion rate
!    ! of the first dry cell to be added in Adv.-Diff.-Eq.
!    !
!    s%ero(i,j,:)   = s%ero(i,j,:)+impact_ero*s%pbbed(i,j,1,:)*weight(1)
!    ero_dry(i,j,:) = impact_ero*s%pbbed(i+1,j,1,:)*weight(2)
!    s%ero(i+1,j,:) = s%ero(i+1,j,:)+ero_dry(i,j,:)
!    do jg=1,par%ngd
!        temp_ero_dry(jg) = ero_dry(i,j,jg)
!    end do
!    !
!    ! In case cell at waterline is lower than previous cell, than material eroded is split
!    ! linearly decreasing over the cells covering a distance equal to a wave length
!    !
!    if (s%zb(i-1,j)>s%zb(i,j)) then
!       wave_L      = 2.d0*par%px/s%k(i,j)
!       wave_L_temp = 0.d0
!       counter     = i
!       do while (wave_L>=wave_L_temp)
!          wave_L_temp = wave_L_temp+s%dsz(counter,j)
!          counter     = counter-1
!          if (counter==0) then
!             exit
!          end if
!       end do
!       ncells       = i-counter
!       s%ero(i,j,:) = s%ero(i,j,:)-impact_ero*s%pbbed(i,j,1,:)*weight(1)
!       x_temp       = s%dsz(i,j)/2
!       do kk=1,ncells
!          s%ero(i+1-kk,j,:) = s%ero(i+1-kk,j,:)+impact_ero*s%pbbed(i+1-kk,j,1,:)*weight(1)/wave_L_temp*2*(1-x_temp/wave_L_temp)
!          !
!          ! Increase x_temp in order to be always at cell center
!          !
!          x_temp = x_temp+s%dsu(i-kk,j)
!          !
!       end do
!    end if
!    !
!    !
!end subroutine retreat_scheme_1
!
!!==========================================================================!
!
!subroutine retreat_scheme_2(s,par,ee,ero_dry,temp_ero_dry,S_toe,S_top,i,j)
!    !
!    use params
!    use spaceparams
!    !
!    implicit none
!    type(parameters)                                          :: par
!    type(spacepars)                                           :: s
!    real*8, dimension(par%nx+1,par%ny+1,par%ngd), intent(out) :: ero_dry      ! [m3/m2/s] eroded material from the dry cells
!    real*8, dimension(par%ngd),                   intent(out) :: temp_ero_dry ! [m3/m2/s] total eroded material from dry cells per fractions
!    real*8                                                    :: impact_ero   ! [m3/m2/s] average erosion rate over the scarp height
!    real*8                                                    :: cell_ero     ! [m3/m2/s] erosion rate per cell
!    real*8                                                    :: deltax       ! [m] length of the scarp (x direction)
!    real*8                                                    :: deltaz       ! [m] heigth of the scarp (z direction)
!    real*8,                                       intent(in)  :: ee           ! [m3/m/s] erosion rate expressed in meter^2 of sediment (without pores) per second;
!                                                                              ! m^2 are referred as a dx*dz quantity
!    real*8                                                    :: courcheck
!    real*8                                                    :: deltat
!    integer,                                      intent(in)  :: S_toe,S_top  ! [-] cell number referred to the toe and top of the scarp at the waterline
!    integer                                                   :: i,j,jg,ii
!    integer                                                   :: cellscarp
!    integer, dimension(par%nx+1,par%ny+1)                     :: sl
!    !
!    !
!    ! Determine scarp length, height and impact_ero which represents a sort of
!    ! celerity of the scarp migrating horizontally
!    ! 
!    deltax     = sum(s%dsz(S_toe:S_top,j))
!    deltaz     = s%zb(S_top,j)-s%zb(S_toe-1,j)
!    impact_ero = ee/deltaZ
!    !
!    ! Initialize total eroded material from dry cells
!    ! 
!    temp_ero_dry = 0.d0
!    !
!    do ii=1,cellscarp
!       do jg=1,par%ngd
!          !
!          ! Erosion rate a single cell
!          !
!          cell_ero = impact_ero*s%pbbed(S_toe+ii-1,j,1,jg)*(s%zb(S_toe+ii-1,j)-s%zb(S_toe+ii-2,j))/s%dsu(S_toe+ii-2,j)
!          !cell_ero = impact_ero/dA(ii)*s%pbbed(toe+ii-1,j,1,jg)*(s%zb(toe+ii-1,j)-s%zb(toe+ii-2,j))/deltaz
!          !
!          courcheck = impact_ero*s%pbbed(S_toe+ii-1,j,1,jg)/s%dsz(S_toe-1+ii,j)*par%dt
!          !courcheck = impact_ero/dA(ii)*s%pbbed(toe+ii-1,j,1,jg)/deltaz*par%dt
!          !deltat    = 1/courcheck*par%dt
!          !
!          s%ero(S_toe+ii-1,j,jg) = s%ero(S_toe+ii-1,j,jg)+cell_ero
!          !
!          ! Account for the material eroded from the dry cells
!          !
!          if (S_toe+ii-1>=i .and. S_toe+ii-1<=S_top) then
!             ero_dry(i,j,jg)  = ero_dry(i,j,jg)+cell_ero
!             temp_ero_dry(jg) = temp_ero_dry(jg)+cell_ero
!          end if
!       end do
!    end do
!
!end subroutine retreat_scheme_2
!
!!==========================================================================!
!
!subroutine retreat_scheme_3
!
!! IDEA: IN CASE OF STRANGE FORM OF THE SCARP OR NEGATIVE SLOPES SHIFT TO RETREAT 2
!
!end subroutine retreat_scheme_3
!
!!==========================================================================!
!
!subroutine identify_scarp(s,par,cellscarp,toe,top,sl,counter)
!    !
!    use params
!    use spaceparams
!    !
!    implicit none
!    type(parameters)                                   :: par
!    type(spacepars)                                    :: s
!    real*8,  dimension(par%nx+1,par%ny+1)              :: dzbdxup
!    integer, dimension(par%nx+1,par%ny+1), intent(out) :: sl        ! [-] flag to identify cells being part of the scarp
!    integer, dimension(par%nx+1),          intent(out) :: cellscarp ! [-] number of cells composing the scarp
!    integer, dimension(par%nx+1),          intent(out) :: toe       ! [-] First cell being part of the scarp
!    integer, dimension(par%nx+1),          intent(out) :: top       ! [-] Last cell being part of the scarp
!    integer,                               intent(out) :: counter   ! [-] identifier of the number of scarps
!    integer                                            :: ii,jj
!    !
!    !
!    ! TODO: CONSIDER TO SET A THRESHOLD FOR SCARP HEIGHT
!    ! Determine local slope of the domain
!    !
!    dzbdxup   = 0.0d0
!    top       = 0
!    toe       = 0
!    cellscarp = 0
!    counter   = 0
!    ! 
!    do jj=1,par%ny+1
!       do ii=2,par%nx
!          !
!          ! The slope of the scarp is determined in an upwind manner
!          !
!          dzbdxup(ii,jj)=(s%zb(ii,jj)-s%zb(ii-1,jj))/s%dsu(ii-1,jj)
!          !
!          ! Identify cells being part of the scarp (slope >= scarpsl)
!          !
!          if (dzbdxup(ii,jj)>=par%scarpsl .and. dzbdxup(ii-1,jj)<par%scarpsl) then
!              !
!              ! At this cell begins a scarp
!              !
!              counter            = counter+1
!              sl(ii,jj)          = 1
!              cellscarp(counter) = cellscarp(counter)+1 ! Number of cells being part of the scarp
!              toe(counter)       = ii
!              !
!          else if (dzbdxup(ii,jj)>=par%scarpsl .and. dzbdxup(ii+1,jj)>=par%scarpsl) then
!              !
!              ! This is a cell within the scarp
!              !
!              sl(ii,jj)          = 1
!              cellscarp(counter) = cellscarp(counter)+1
!              !
!          else if (dzbdxup(ii,jj)>=par%scarpsl .and. dzbdxup(ii+1,jj)<par%scarpsl) then
!              !
!              ! This is the last cell of the scarp
!              !
!              sl(ii,jj)          = 1
!              top(counter)       = ii
!              !
!          else
!              !
!              ! ...I don't know...
!              !
!              sl(ii,jj) = 0
!              !
!          end if
!       end do
!    end do
!    !
!    !
!end subroutine identify_scarp
!
!!==========================================================================!
!
!subroutine distr_ero_dry(s,par,temp_ero_dry,sed_input,i,j)
!    !
!    use params
!    use spaceparams
!    !
!    implicit none
!    type(parameters)                                          :: par
!    type(spacepars)                                           :: s
!    real*8, dimension(par%ngd),                   intent(in)  :: temp_ero_dry ! [m3/m2/s] total eroded material from dry cells per fractions
!    real*8, dimension(par%nx+1,par%ny+1,par%ngd), intent(out) :: sed_input    ! [m3/m2/s] flux of material to fluid phase coming from dry cells
!    real*8                                                    :: wave_L       ! [m] wave length
!    real*8                                                    :: wave_L_temp  ! [m] "cumulative" wave length
!    real*8                                                    :: x_temp
!    integer                                                   :: ncells       ! [-] number of cells covering a wave length
!    integer                                                   :: counter
!    integer                                                   :: i,j,jg,kk
!    !
!    !
!    ! Eroded material from dry cells is distributed linearly decreasing from the waterline up to a 
!    ! distance equal to a wave length ("similar" approach to Nasermoaddeli, 2011)
!    ! TODO: CHECK NASERMOADDELI APPROACH
!    !
!    sed_input = 0.0d0 ! (to be sure)
!    !
!    do jg=1,par%ngd
!       wave_L      = 2.d0*par%px/s%k(i,j)
!       wave_L_temp = 0.d0
!       counter     = i
!       do while (wave_L>=wave_L_temp)
!          wave_L_temp = wave_L_temp+s%dsz(counter,j)
!          counter     = counter-1
!          if (counter==0) then
!             exit
!          end if
!       end do
!       ncells = i-counter
!       x_temp = s%dsz(i,j)/2
!       do kk=1,ncells
!          sed_input(i+1-kk,j,jg) = 2.d0*temp_ero_dry(jg)/wave_L_temp*(1-x_temp/wave_L_temp)
!          !
!          ! Increase x_temp in order to be always at cell center
!          !
!          x_temp = x_temp+s%dsu(i-kk,j)
!          !
!       end do
!    end do
!    
!    !
!    !
!end subroutine distr_ero_dry
!
!!==========================================================================!

end module readbedcomposition_module