!======================================================================
! CROCO is a branch of ROMS developped at IRD, INRIA, 
! Ifremer, CNRS and Univ. Toulouse III  in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org/
!
! These Sediment Dynamics routines are based on the Community Sediment
! Transport Model (CSTM) as implemented in ROMS by J. Warner (USGS)
!
! Adapted from John Warner's code to CROCO by R. Benshila, CNRS 2015
! -- modifs:
!   P. Marchesiello, Dec 2015
!                    Jun 2020: Van der A (2013) bedload
!   P. Marchesiello & G. Morvan, Sep 2021: 5th order interpolation
!     of bedload fluxes; DUNE test case in common with MUSTANG model
!
!=====================================================================
!
#include "cppdefs.h"
#ifdef SEDIMENT

# ifdef BEDLOAD /* for bedload scheme */
#  define bnew nnew
# else
#  define bnew nstp
# endif

# ifdef BEDLOAD
#  ifndef BEDLOAD_VANDERA
#   define BSTRESS_UPWIND
#  endif
#  ifdef ANA_DUNE
#   undef  BSTRESS_UPWIND
#   undef  BEDLOAD_SLUMP
#   undef  BEDLOAD_SHALLOW_LIMIT
#  endif
#  ifdef BEDLOAD_VANDERA /* wave half-cycle concept */
#   define BEDLOAD_VANDERA_ZEROCURR
#   ifndef BEDLOAD_VANDERA_ZEROCURR
#    define BEDLOAD_VANDERA_FD_MADSEN
#    undef  BEDLOAD_VANDERA_FD_SANTOSS
#   endif
#   define BEDLOAD_VANDERA_ASYM_LIMITS
#   undef  BEDLOAD_VANDERA_SURFACE_WAVE
#   undef  BEDLOAD_VANDERA_WAVE_AVG_STRESS
#   define BEDLOAD_SLUMP
#   define BEDLOAD_SHALLOW_LIMIT
#  endif
# endif
!
# undef  LINEAR_CONTINUATION /* for settling scheme */
# undef  NEUMANN
!
# define MORPH_FAC *morph_fac

      subroutine sediment (tile)
!
!=======================================================================
!                                                                      !
!  This routine it is the main driver for the sediment-transport       !
!  model. Currently, it includes calls to the following routines:      !
!                                                                      !
!  * Vertical settling of sediment in the water column.                !
!  * Erosive and depositional flux interactions of sediment            !
!    between water column and the bed.                                 !
!  * Transport of multiple grain sizes.                                !
!  * Bed layer stratigraphy.                                           !
!  * Bed morphology.                                                   !
!  * Bedload based on Meyer Peter Mueller.                             !
!  * Bedload based on Soulsby combined waves + currents                !
!    (p166 Soulsby 1997)                                               !
!  * Bedload slope term options: Nemeth et al, 2006, Coastal           !
!    Engineering, v 53, p 265-275; Lesser et al, 2004, Coastal         !
!    Engineering, v 51, p 883-915.                                     !
!                                                                      !
!  * Seawater/sediment vertical level distribution:                    !
!                                                                      !
!         W-level  RHO-level                                           !
!                                                                      !
!            N     _________                                           !
!                 |         |                                          !
!                 |    N    |                                          !
!          N-1    |_________|  S                                       !
!                 |         |  E                                       !
!                 |   N-1   |  A                                       !
!            2    |_________|  W                                       !
!                 |         |  A                                       !
!                 |    2    |  T                                       !
!            1    |_________|  E                                       !
!                 |         |  R                                       !
!                 |    1    |                                          !
!            0    |_________|_____ bathymetry                          !
!                 |/////////|                                          !
!                 |    1    |                                          !
!            1    |_________|  S                                       !
!                 |         |  E                                       !
!                 |    2    |  D                                       !
!            2    |_________|  I                                       !
!                 |         |  M                                       !
!                 |  NLAY-1 |  E                                       !
!        Nbed-1   |_________|  N                                       !
!                 |         |  T                                       !
!                 |  NLAY   |                                          !
!         Nbed    |_________|                                          !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G.    !
!    Arango, 2008:  Development of a three-dimensional,  regional,     !
!    coupled wave, current, and sediment-transport model, Computers    !
!    & Geosciences, 34, 1284-1306.                                     !
!                                                                      !
!=======================================================================

      implicit none
# include "param.h"
# include "ocean3d.h"
# include "scalars.h"
# include "sediment.h"
# include "private_scratch.h"
      integer tile, trd, omp_get_thread_num, itrc
# include "compute_tile_bounds.h"
!
      trd=omp_get_thread_num()
# ifdef BEDLOAD
      call sed_bedload_tile (Istr,Iend,Jstr,Jend,
     &                       A2d(1,1,trd), A2d(1,2,trd),
     &                       A2d(1,3,trd), A2d(1,4,trd),
     &                       A2d(1,5,trd),
#  ifdef BEDLOAD_VANDERA
     &                       A2d(1,6,trd), A2d(1,7,trd),
     &                       A2d(1,8,trd)
#  else
     &                       A2d(1,6,trd), A2d(1,7,trd)
#  endif
     &                      )
C$OMP BARRIER
# endif /* BEDLOAD */
# ifdef SUSPLOAD
#  ifdef SED_FLOCS
      call sed_flocmod_tile (Istr,Iend,Jstr,Jend)
#  endif
C$OMP BARRIER
      call sed_settling_tile (Istr,Iend,Jstr,Jend,
     &                        A2d(1,1,trd), A2d(1,2,trd),
     &                        A2d(1,3,trd), A2d(1,4,trd),
     &                        A2d(1,5,trd), A2d(1,6,trd),
     &                        A2d(1,7,trd), A2d(1,8,trd),
     &                        A2d(1,9,trd), A2d(1,10,trd),
     &                        B2d(1,trd) )
C$OMP BARRIER
      call sed_fluxes_tile (Istr,Iend,Jstr,Jend,A2d(1,1,trd))
C$OMP BARRIER
# endif /* SUSPLOAD */

# if defined COHESIVE_BED || defined MIXED_BED
      call sed_bed_cohesive_tile (Istr,Iend,Jstr,Jend,
     &                   A2d(1,1,trd), A2d(1,2,trd) )
# else
      call sed_bed_tile (Istr,Iend,Jstr,Jend,
     &                   A2d(1,1,trd), A2d(1,2,trd) ,
     &                   A2d(1,3,trd), A2d(1,4,trd) )
# endif
C$OMP BARRIER
!
      call sed_surf_tile (Istr,Iend,Jstr,Jend)
C$OMP BARRIER
# if defined TRACERS
      do itrc=1,NT
#  if defined EW_PERIODIC || defined NS_PERIODIC || defined MPI
#   ifdef THREE_GHOST_POINTS_TS
        call  exchange_r3d_3pts_tile (Istr,Iend,Jstr,Jend,
     &          t(START_2D_ARRAY,1,nnew,itrc))
#   else
        call  exchange_r3d_tile (Istr,Iend,Jstr,Jend,
     &          t(START_2D_ARRAY,1,nnew,itrc))
#   endif
#  endif
      enddo
# endif

!
      return
      end

# ifdef BEDLOAD
!======================================================================
!                                                                      !
!  This routine computes sediment bedload transport using the Meyer-   !
!  Peter and Muller (1948) formulation for unidirectional flow,        !
!  or SANTOSS formulation (Van der A 2013), that accounts for waves by !
!  treating phase-lag and wave asymetry effects in half-cycles         !
!  (from Warner et al. implementation in COAWST)                       !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Van der A, D. A., Ribberink, J. S., van der Werf, J. J.,            !
!    O’Donoghue, T., Buijsrogge, R. H., and Kranenburg, W. M., 2013:   !
!    Practical sand transport formula for non-breaking waves and       !
!    currents. Coastal Engineering, 76, 26-42                          !
!                                                                      !
!  Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G.    !
!    Arango, 2008:  Development of a three-dimensional,  regional,     !
!    coupled wave, current, and sediment-transport model, Computers    !
!    & Geosciences, 34, 1284-1306.                                     !
!                                                                      !
!  Meyer-Peter, E. and R. Muller, 1948: Formulas for bedload transport !
!    In: Report on the 2nd Meeting International Association Hydraulic !
!    Research, Stockholm, Sweden, pp 39-64.                            !
!                                                                      !
!  Wu W. & Lin Q., 2014: Nonuniform sediment transport under non-      !
!    breaking waves and currents. Coastal Engineering, 90, pp 1-11.    !
!                                                                      !
!======================================================================
!
      subroutine sed_bedload_tile (Istr, Iend, Jstr, Jend,
     &                             tau_w,FX,FE,FX_r,FE_r
#  ifdef BEDLOAD_VANDERA
     &                            ,tau_m,phic,phicw
#  else
     &                            ,angleu,anglev
#  endif
     &                            )
!
      implicit none
      integer Istr,Iend,Jstr,Jend
      real  FX(PRIVATE_2D_SCRATCH_ARRAY)
     &     ,FE(PRIVATE_2D_SCRATCH_ARRAY)
     &     ,FX_r(PRIVATE_2D_SCRATCH_ARRAY)
     &     ,FE_r(PRIVATE_2D_SCRATCH_ARRAY)
     &     ,tau_w(PRIVATE_2D_SCRATCH_ARRAY)
#  ifdef BEDLOAD_VANDERA
     &     ,tau_m(PRIVATE_2D_SCRATCH_ARRAY)
     &     ,phic(PRIVATE_2D_SCRATCH_ARRAY)
     &     ,phicw(PRIVATE_2D_SCRATCH_ARRAY)
#  else
     &     ,angleu(PRIVATE_2D_SCRATCH_ARRAY)
     &     ,anglev(PRIVATE_2D_SCRATCH_ARRAY)
#  endif
!
      integer i, j, k, ised
      real eps, Dstp
      real cff, cff1, cff2, cff3, cff4, cff5,cff6
      real a_slopex, a_slopey, sed_angle
      real bedld, bedld_mass, dzdx, dzdy, dzdxdy
      real smgd, smgdr, osmgd, Umag
      real Ua, Ra, phi, Clim
      real kdmax,Kbh,Kbh2,Kdh,Fwave
      real rhs_bed, bed_change
      real twopi, otwopi, sqrt2
      integer imin,imax,jmin,jmax
      parameter (eps = 1.e-20)
      parameter (kdmax = 100.)
      real K1, K2, K3, K4, K5, K6
      parameter (K1=0.6666666666, K2=0.3555555555,
     &           K3=0.1608465608, K4=0.0632098765,
     &           K5=0.0217540484, K6=0.0065407983)
#  if defined BEDLOAD_WENO5
      real :: flux5_weno
#   define FLUX5 flux5_weno
#  elif defined BEDLOAD_UP5
      real :: flux5
#   define FLUX5 flux5
#  endif
#  ifdef BEDLOAD_VANDERA
      real wavecycle, alphac, alphaw
      real ksd_wbl, ustrc_wbl, thck_wbl, udelta_wbl
      real phic_sgwbl, fd_wbl
      real Hs, Td, depth, k_wn, c_w
      real d50, d50_mix, d90, rhos
      real urms, uhat, ahat, phi_curwave
      real theta, fd, ksw, eta, alpha, tau_wRe
      real ursell_no, RR_asymwave, beta_asymwave
      real r, Su, Au, Sk, Ak, T_cu, T_tu
      real ucrest_r, utrough_r, T_crest, T_trough
      real umax, umin, uhat_c, uhat_t, mag_uc, mag_ut
      real dsf_c, dsf_t, theta_c, theta_t
      real theta_cx, theta_cy, theta_tx, theta_ty
      real mag_theta_c, mag_theta_t, mag_bstrc
      real om_cc, om_tt, om_ct, om_tc
      real bedld_cx, bedld_cy, bedld_tx, bedld_ty
      real bedld_x, bedld_y
      real kh_calc    ! functions
      real mu_calc, ksw_calc, fw_calc, fwi_calc
      real dsf_calc, w_s_calc, w_sc_calc, theta_cr_calc
#  endif
#  ifdef BEDLOAD_SLUMP
      real sedslope_crit, slopefac, slopefac_local
#  endif
#  include "param.h"
#  include "scalars.h"
#  include "grid.h"
#  include "ocean3d.h"
#  include "ocean2d.h"
#  include "sediment.h"
#  include "forces.h"
#  ifdef WKB_WWAVE
#   include "wkb_wwave.h"
#  endif
#  include "bbl.h"
#  include "mixing.h"
#  include "compute_auxiliary_bounds.h"
!
!-----------------------------------------------------------------------
! Set boundary Indices
!-----------------------------------------------------------------------
!
#  ifdef EW_PERIODIC
#   define I_EXT_RANGE Istr-1,Iend+1
#  else
      if (WESTERN_EDGE) then
        imin=Istr
      else
        imin=Istr-1
      endif
      if (EASTERN_EDGE) then
        imax=Iend
      else
        imax=Iend+1
      endif
#   define I_EXT_RANGE imin,imax
#  endif
#  ifdef NS_PERIODIC
#   define J_EXT_RANGE Jstr-1,Jend+1
#  else
      if (SOUTHERN_EDGE) then
        jmin=Jstr
      else
        jmin=Jstr-1
      endif
      if (NORTHERN_EDGE) then
        jmax=Jend
      else
        jmax=Jend+1
      endif
#   define J_EXT_RANGE jmin,jmax
#  endif

#  if defined BEDLOAD_WENO5 || defined BEDLOAD_UP5
#   define FX_r bedload_FX
#   define FE_r bedload_FE
#  endif
!
!-----------------------------------------------------------------------
! Set parameters
!-----------------------------------------------------------------------
!
      twopi=2.*pi
      otwopi=1./twopi
      sqrt2=SQRT(2.)
!
#  ifdef BEDLOAD_VANDERA
      alphac=1.0   ! bedload scale fac for curr contrib
      alphaw=1.0   ! bedload scale fac for wave contrib
#  endif
!
!-----------------------------------------------------------------------
! Compute bottom stress components at rho points
!    tau_m: form stress for flow
!    tau_w: skin stress for sediments
!-----------------------------------------------------------------------
!
      do j=J_EXT_RANGE
        do i=I_EXT_RANGE
#  if defined BEDLOAD_VANDERA && !defined BEDLOAD_VANDERA_ZEROCURR
          tau_m(i,j)=0.5*SQRT( (bustr(i,j)+bustr(i+1,j))**2
     &                        +(bvstr(i,j)+bvstr(i,j+1))**2 )
#  endif
          tau_w(i,j)=SQRT( bustrw(i,j)*bustrw(i,j)
     &                   + bvstrw(i,j)*bvstrw(i,j) )
#  if !defined BEDLOAD_VANDERA && defined BSTRESS_UPWIND
#   define tau_wX angleu
#   define tau_wE anglev
          cff1=0.5*(1.0+SIGN(1.0,bustr(i+1,j)))
          cff2=0.5*(1.0-SIGN(1.0,bustr(i+1,j)))
          cff3=0.5*(1.0+SIGN(1.0,bustr(i  ,j)))
          cff4=0.5*(1.0-SIGN(1.0,bustr(i  ,j)))
          tau_wX(i,j)=cff3*(cff1*bustr(i,j)+
     &                cff2*0.5*(bustr(i,j)+bustr(i+1,j)))+
     &                cff4*(cff2*bustr(i+1,j)+
     &                cff1*0.5*(bustr(i,j)+bustr(i+1,j)))
          cff1=0.5*(1.0+SIGN(1.0,bvstr(i,j+1)))
          cff2=0.5*(1.0-SIGN(1.0,bvstr(i,j+1)))
          cff3=0.5*(1.0+SIGN(1.0,bvstr(i,j)))
          cff4=0.5*(1.0-SIGN(1.0,bvstr(i,j)))
          tau_wE(i,j)=cff3*(cff1*bvstr(i,j)+
     &                cff2*0.5*(bvstr(i,j)+bvstr(i,j+1)))+
     &                cff4*(cff2*bvstr(i,j+1)+
     &                cff1*0.5*(bvstr(i,j)+bvstr(i,j+1)))
#  endif
        enddo
      enddo
!
!-----------------------------------------------------------------------
! Compute angles
!-----------------------------------------------------------------------
!
! Compute some constant bed slope parameters.
!
      sed_angle=DTAN(33.*pi/180.)
!
!  Compute angle between currents and waves (radians).
!
      do j=J_EXT_RANGE
       do i=I_EXT_RANGE
#  ifdef BEDLOAD_VANDERA
!
! Compute angle between currents and waves
! .. from current direction toward wave vector.
!
          cff1=0.5*(bustr(i,j)+bustr(i+1,j))
          cff2=0.5*(bvstr(i,j)+bvstr(i,j+1))
          if (cff1.eq.0.) then
            phic(i,j)=0.5*pi*SIGN(1.,cff2)
          else
            phic(i,j)=ATAN2(cff2,cff1)
          endif
          phicw(i,j)=Dwave(i,j)-phic(i,j)  ! Dwave assumed in CROCO
                                           ! grid convention
!
! Convert phicw to angle between waves and current
!  .. from wave direction towards current vector
!
          phicw(i,j)=twopi-phicw(i,j)
!
#  elif !defined BSTRESS_UPWIND
          cff1=0.5*(bustr(i,j)+bustr(i+1,j))
          cff2=0.5*(bvstr(i,j)+bvstr(i,j+1))
          Umag=SQRT(cff1*cff1+cff2*cff2)+eps
          angleu(i,j)=cff1/Umag
          anglev(i,j)=cff2/Umag
#  endif /* BEDLOAD_VANDERA */
        enddo
      enddo
!
      do ised=1,NST   ! <============== sed loop
!
!-----------------------------------------------------------------------
! Compute magnitude of bedload at rho points
!
!  bedld has dimensions of kg m-1 s-1.
!  FX_r and FE_r have dimensions of kg.
!-----------------------------------------------------------------------
!
#  ifdef BEDLOAD_MARIEU
        smgdr=Srho(ised)
#  else
        smgd=(Srho(ised)/rho0-1.)*g*Sd(ised)
        osmgd=1./smgd
        smgdr=SQRT(smgd)*Sd(ised)*Srho(ised)
#  endif
!
#  ifdef BEDLOAD_VANDERA
        rhos=Srho(ised)          ! (kg/m3)
        d50=Sd(ised)             ! (m)
        d90=1.3*d50              ! (m)
#  endif
!
        do j=J_EXT_RANGE     ! <===== i,j loops
          do i=I_EXT_RANGE
!
#  ifdef BEDLOAD_VANDERA
!
! -----------------  Van der A formulation  ---------------------
!
! VA-2013 equation 1 is solved:
! The transport formula is based on the wave “half-cycle” concept:
! The wave-averaged net sand transport rate in the wave boundary
! layer is described as the difference between sand transported
! during the positive “crest” and the negative “trough” half-cycles.
! Phase-lag and wave asymetry effects are taken into account.
!
! Wave & sediment parameters
!
            depth=z_w(i,j,N)+h(i,j)             ! (m)
            Hs=2.*Awave(i,j)                    ! (m)
            Td=MAX(Pwave(i,j),1.)               ! (s)
            phi_curwave=phicw(i,j)              ! angle cw
            k_wn=kh_calc(Td,depth)/depth        ! Wave number
            c_w=twopi/(k_wn*Td)                 ! Wave speed

            if(NST>1) then
              d50_mix=Sdens(i,j)
            else
             d50_mix=d50
            endif
!
! Compute bed wave orbital velocity (m/s) & excursion amplitude
!
            urms=twopi/Td*Awave(i,j)/(SINH(k_wn*depth)+eps)
#   ifdef MASKING
     &                                          *rmask(i,j)
#   endif
#   ifdef WET_DRY
     &                                      *rmask_wet(i,j)
#   endif
            uhat=urms*sqrt2
            ahat=uhat*Td*otwopi
!
! Magnitude of stress for computing current velocity
! at the wave boundary layer
!
#   ifndef BEDLOAD_VANDERA_ZEROCURR
            mag_bstrc=tau_m(i,j)
#    ifdef MASKING
     &                             *rmask(i,j)
#    endif
#    ifdef WET_DRY
     &                         *rmask_wet(i,j)
#    endif
#   else
            mag_bstrc=0.
#   endif
!
! Compute near bottom current at wave boundary layer (wbl)
! Use wbl thickness based on Madsen to get near bottom current velocity.
! Make sure that wbl is under total depth and greater than apparent
! roughness.
!
#   ifdef BEDLOAD_VANDERA_ZEROCURR
            udelta_wbl=0.
            ksd_wbl=0.
            thck_wbl=0.
            ustrc_wbl=0.
#   else
#    ifdef BBL
            ksd_wbl=30.*Zbapp(i,j)
#    else
            ksd_wbl=30.*Zob(i,j)
#    endif
            thck_wbl=0.09*ksd_wbl*(ahat/ksd_wbl)**0.82
            !thck_wbl=0.1
            ustrc_wbl=sqrt(mag_bstrc)
            cff1=MAX(MIN(0.9*depth, thck_wbl), 1.1*ksd_wbl)
            cff2=LOG(30.*cff1/ksd_wbl)
            udelta_wbl=(ustrc_wbl/vonKar)*cff2
#   endif
!
! Ruessink et al. provides equations for calculating skewness parameters
! Uses Malarkey and Davies equations to get "r" and "phi"
! common to both crest and trough.
!
            call skewness_params(Hs, Td, depth, r, phi, ursell_no)
!
! Abreu et al. use skewness params (r, phi) to get representative critical
! orbital velocity for crest and trough cycles.
!
            call abreu_points(r, phi, uhat, Td,
     &                        T_crest, T_trough,
     &                        T_cu, T_tu, umax, umin,
     &                        RR_asymwave, beta_asymwave)
!
! Representative crest half cycle water particle velocity
! as well as full cycle orbital velocity and excursion.
! VA2013 Equation 10, 11.
!
            uhat_c=umax
            uhat_t=umin
            ucrest_r=0.5*sqrt2*uhat_c
            utrough_r=0.5*sqrt2*uhat_t
!
! Full wave cycle
!
            call full_wave_cycle_stress_factors(d50, d90, osmgd,
     &                                           Td, c_w, depth,
     &                                              phi_curwave,
     &                                  RR_asymwave, uhat, ahat,
     &                                    umax, umin, mag_bstrc,
     &                            T_crest, T_trough, T_cu, T_tu,
     &                                       ksd_wbl, ustrc_wbl,
     &                             thck_wbl, udelta_wbl, fd_wbl,
     &                                 alpha, eta, ksw, tau_wRe )
!
! Bed shear stress (Shields parameter) for Crest half cycle
! alpha VA2013 Eqn. 19
!
            call half_wave_cycle_stress_factors( T_cu, T_crest,
     &                                               ahat, ksw,
     &                                           fd_wbl, alpha,
     &                                          alphac, alphaw,
     &                                              d50, osmgd,
     &               ucrest_r, uhat_c, udelta_wbl, phi_curwave,
     &                                         tau_wRe,  dsf_c,
     &                         theta_cx, theta_cy, mag_theta_c )
!
! Compute sediment load entrained during each crest half cycle
!
            wavecycle=1.
            call sandload_vandera( wavecycle,
     &                             Hs, Td,  depth, RR_asymwave,
     &                                 d50, d50_mix, rhos, c_w,
     &                                              eta, dsf_c,
     &                                   T_crest, T_cu, uhat_c,
     &                               mag_theta_c, om_cc, om_ct )
!
! Bed shear stress (Shields parameter) for Trough half cycle
! alpha VA2013 Eqn. 19
!
            call half_wave_cycle_stress_factors( T_tu, T_trough,
     &                                                ahat, ksw,
     &                                            fd_wbl, alpha,
     &                                           alphac, alphaw,
     &                                               d50, osmgd,
     &               utrough_r, uhat_t, udelta_wbl, phi_curwave,
     &                                           tau_wRe, dsf_t,
     &                          theta_tx, theta_ty, mag_theta_t )
!
! Compute sediment load entrained during each trough half cycle
!
            wavecycle=-1.
            call sandload_vandera( wavecycle,
     &                              Hs, Td,  depth, RR_asymwave,
     &                                  d50, d50_mix, rhos, c_w,
     &                                               eta, dsf_t,
     &                                   T_trough, T_tu, uhat_t,
     &                                mag_theta_t, om_tt, om_tc )
!
! Compute sediment load entrained during crest and trough half cycles
!
            cff1=MAX(0.5*T_crest/T_cu, 0.)
            cff2=sqrt(mag_theta_c)*(om_cc+cff1*om_tc)
            cff3=(theta_cx/mag_theta_c)
            bedld_cx=cff2*cff3
            cff3=(theta_cy/mag_theta_c)
            bedld_cy=cff2*cff3
!
            cff1=MAX(0.5*T_trough/T_tu, 0.)
            cff2=sqrt(mag_theta_t)*(om_tt+cff1*om_ct)
            cff3=(theta_tx/mag_theta_t)
            bedld_tx=cff2*cff3
            cff3=(theta_ty/mag_theta_t)
            bedld_ty=cff2*cff3
!
! Total bedload (kg m-1 sec-1). Use velocity-load equation 1 of VA2013
!
            bedld_x=smgdr*( bedld_cx*T_crest+
     &                      bedld_tx*T_trough )/Td
            bedld_y=smgdr*( bedld_cy*T_crest+
     &                      bedld_ty*T_trough )/Td
!
! Partition bedload into xi and eta directions, still at rho points
! (FX_r and FE_r have dimensions of kg).
! --> project on CROCO grid (bedld is wave-aligned in Van Der A)
!
            FX_r(i,j)=(bedld_x*COS(Dwave(i,j)) -
     &                 bedld_y*SIN(Dwave(i,j)))*on_r(i,j)*dt
            FE_r(i,j)=(bedld_x*SIN(Dwave(i,j)) +
     &                 bedld_y*COS(Dwave(i,j)))*om_r(i,j)*dt
!
! -----------------  END Van der A formulation  ------------------
!
#  else
!
!                    Use simpler formulations
!
! -------------  Meyer-Peter Muller or Wu & Lin  -----------------
!
#   ifdef BSTRESS_UPWIND
!
! Use partitions of stress from upwind direction, still at rho points.
!
#    ifdef BEDLOAD_MPM
            bedld=8.*(MAX((ABS(tau_wX(i,j))*osmgd-0.047),0.)**1.5)
#    elif defined BEDLOAD_WULIN
            cff=ABS(tau_wX(i,j))/tau_ce_2d(i,j,ised)-1.
            bedld=0.0053*(MAX(cff,0.)**2.2)
#    elif defined BEDLOAD_MARIEU
            bedld=(0.001*ubar(i,j,knew)**3)
#    endif
     &                                 *smgdr*SIGN(1.,tau_wX(i,j))
            FX_r(i,j)=bedld*on_r(i,j)*dt
#    ifdef BEDLOAD_MPM
            bedld=8.*(MAX((ABS(tau_wE(i,j))*osmgd-0.047),0.)**1.5)
#    elif defined BEDLOAD_WULIN
            cff=ABS(tau_wE(i,j))/tau_ce_2d(i,j,ised)-1
            bedld=0.0053*(MAX(cff,0.)**2.2)
#    elif defined BEDLOAD_MARIEU
            bedld=0.
#    endif
     &                                 *smgdr*SIGN(1.,tau_wE(i,j))
            FE_r(i,j)=bedld*om_r(i,j)*dt
#    undef tau_wX
#    undef tau_wE
#   else
!
! Use mean stress tau_w at rho points
!
#    ifdef BEDLOAD_MPM
            bedld=8.*(MAX((tau_w(i,j)*osmgd-0.047),0.)**1.5)*smgdr
#    elif defined BEDLOAD_WULIN
            cff=tau_w(i,j)/tau_ce_2d(i,j,ised)-1.
            bedld=0.0053*(MAX(cff,0.)**2.2)*smgdr
#    elif defined BEDLOAD_MARIEU
            bedld=(0.001*ubar(i,j,knew)**3)*smgdr
#    endif
!
! Partition bedld into xi and eta directions, still at rho points.
! (FX_r and FE_r have dimensions of kg).
!
            FX_r(i,j)=angleu(i,j)*bedld*on_r(i,j)*dt
            FE_r(i,j)=anglev(i,j)*bedld*om_r(i,j)*dt
#   endif /* BSTRESS_UPWIND */

#  endif /* BEDLOAD_VANDERA */
!
!-----------------------------------------------------------------------
! Correct for along-direction slope.
! Limit slope to 0.9*sed angle.
!-----------------------------------------------------------------------
!
#  if defined SLOPE_NEMETH || defined SLOPE_LESSER
            cff1=0.5*(1.+SIGN(1.,FX_r(i,j)))
            cff2=0.5*(1.-SIGN(1.,FX_r(i,j)))
            cff3=0.5*(1.+SIGN(1.,FE_r(i,j)))
            cff4=0.5*(1.-SIGN(1.,FE_r(i,j)))
#  endif
#  ifdef SLOPE_NEMETH
            dzdx=(h(i+1,j)-h(i  ,j))/om_u(i+1,j)*cff1+
     &           (h(i-1,j)-h(i  ,j))/om_u(i  ,j)*cff2
            dzdy=(h(i,j+1)-h(i,j  ))/on_v(i,j+1)*cff3+
     &           (h(i,j-1)-h(i,j  ))/on_v(i  ,j)*cff4

            a_slopex=1.7*dzdx
            a_slopey=1.7*dzdy
!
! Add contribution of bed slope to bed load transport fluxes.
!
            FX_r(i,j)=FX_r(i,j)*(1.+a_slopex)
            FE_r(i,j)=FE_r(i,j)*(1.+a_slopey)
!
#  elif defined SLOPE_LESSER
            dzdx=MIN(((h(i+1,j)-h(i  ,j))/om_u(i+1,j)*cff1+
     &                (h(i  ,j)-h(i-1,j))/om_u(i  ,j)*cff2),0.52)*
     &                SIGN(1.,FX_r(i,j))
            dzdy=MIN(((h(i,j+1)-h(i,j  ))/on_v(i,j+1)*cff3+
     &                (h(i,j  )-h(i,j-1))/on_v(i  ,j)*cff4),0.52)*
     &                SIGN(1.,FE_r(i,j))
            cff=ATAN(dzdx)
            a_slopex=sed_angle/(COS(cff)*(sed_angle-dzdx))
            cff=ATAN(dzdy)
            a_slopey=sed_angle/(COS(cff)*(sed_angle-dzdy))
!
! Add contribution of bed slope to bed load transport fluxes.
!
            FX_r(i,j)=FX_r(i,j)*a_slopex
            FE_r(i,j)=FE_r(i,j)*a_slopey
!
#  endif /* SLOPE_NEMETH || SLOPE_LESSER */
!
!-----------------------------------------------------------------------
! Apply morphology factor.
!-----------------------------------------------------------------------
!
            FX_r(i,j)=FX_r(i,j) MORPH_FAC
            FE_r(i,j)=FE_r(i,j) MORPH_FAC
!
!-----------------------------------------------------------------------
! Apply bedload transport rate coefficient.
! Also limit bedload to the fraction of each sediment class.
!-----------------------------------------------------------------------
!
            FX_r(i,j)=FX_r(i,j)*bedload_coeff*bed_frac(i,j,1,ised)
            FE_r(i,j)=FE_r(i,j)*bedload_coeff*bed_frac(i,j,1,ised)
!
!-----------------------------------------------------------------------
! Limit bed load to available bed mass.
!-----------------------------------------------------------------------
!
            bedld_mass=ABS(FX_r(i,j))+ABS(FE_r(i,j))+eps
            FX_r(i,j)=MIN(ABS(FX_r(i,j)),
     &                    bed_mass(i,j,1,nstp,ised)*
     &                    om_r(i,j)*on_r(i,j)*ABS(FX_r(i,j))/
     &                    bedld_mass)*
     &                    SIGN(1.,FX_r(i,j))
            FE_r(i,j)=MIN(ABS(FE_r(i,j)),
     &                    bed_mass(i,j,1,nstp,ised)*
     &                    om_r(i,j)*on_r(i,j)*ABS(FE_r(i,j))/
     &                    bedld_mass)*
     &                    SIGN(1.,FE_r(i,j))
#  ifdef MASKING
            FX_r(i,j)=FX_r(i,j)*rmask(i,j)
            FE_r(i,j)=FE_r(i,j)*rmask(i,j)
#  endif
#  ifdef WET_DRY
            FX_r(i,j)=FX_r(i,j)*rmask_wet(i,j)
            FE_r(i,j)=FE_r(i,j)*rmask_wet(i,j)
#  endif
!
          enddo ! I_EXT_RANGE
        enddo  ! J_EXT_RANGE  <===== i,j loops
!
!-----------------------------------------------------------------------
! Apply boundary conditions
!-----------------------------------------------------------------------
!
#  ifndef EW_PERIODIC
        if (WESTERN_EDGE) then
          do j=J_EXT_RANGE
            FX_r(Istr-1,j)=FX_r(Istr,j)
            FE_r(Istr-1,j)=FE_r(Istr,j)
          enddo
        endif
        if (EASTERN_EDGE) then
          do j=J_EXT_RANGE
            FX_r(Iend+1,j)=FX_r(Iend,j)
            FE_r(Iend+1,j)=FE_r(Iend,j)
          enddo
        endif
#  endif
#  ifndef NS_PERIODIC
        if (SOUTHERN_EDGE) then
          do i=I_EXT_RANGE
            FX_r(i,Jstr-1)=FX_r(i,Jstr)
            FE_r(i,Jstr-1)=FE_r(i,Jstr)
          enddo
        endif
        if (NORTHERN_EDGE) then
          do i=I_EXT_RANGE
            FX_r(i,Jend+1)=FX_r(i,Jend)
            FE_r(i,Jend+1)=FE_r(i,Jend)
          enddo
        endif
#  endif

#  if !defined EW_PERIODIC && !defined NS_PERIODIC
        if ((SOUTHERN_EDGE).and.(WESTERN_EDGE)) then
          FX_r(Istr-1,Jstr-1)=0.5*(FX_r(Istr  ,Jstr-1)+
     &                             FX_r(Istr-1,Jstr  ))
          FE_r(Istr-1,Jstr-1)=0.5*(FE_r(Istr  ,Jstr-1)+
     &                             FE_r(Istr-1,Jstr  ))
        endif
        if ((SOUTHERN_EDGE).and.(EASTERN_EDGE)) then
          FX_r(Iend+1,Jstr-1)=0.5*(FX_r(Iend  ,Jstr-1)+
     &                             FX_r(Iend+1,Jstr  ))
          FE_r(Iend+1,Jstr-1)=0.5*(FE_r(Iend  ,Jstr-1)+
     &                             FE_r(Iend+1,Jstr  ))
        endif
        if ((NORTHERN_EDGE).and.(WESTERN_EDGE)) then
          FX_r(Istr-1,Jend+1)=0.5*(FX_r(Istr-1,Jend  )+
     &                             FX_r(Istr  ,Jend+1))
          FE_r(Istr-1,Jend+1)=0.5*(FE_r(Istr-1,Jend  )+
     &                             FE_r(Istr  ,Jend+1))
        endif
        if ((NORTHERN_EDGE).and.(EASTERN_EDGE)) then
          FX_r(Iend+1,Jend+1)=0.5*(FX_r(Iend+1,Jend  )+
     &                             FX_r(Iend  ,Jend+1))
          FE_r(Iend+1,Jend+1)=0.5*(FE_r(Iend+1,Jend  )+
     &                             FE_r(Iend  ,Jend+1))
        endif
#  endif

#  undef I_EXT_RANGE
#  undef J_EXT_RANGE
!
!-----------------------------------------------------------------------
! Compute face fluxes at u and v points before taking divergence.
!
!  --> Extended stencil in case of high-order interpolation.
!-----------------------------------------------------------------------
!
#  if defined BEDLOAD_WENO5 || defined BEDLOAD_UP5
C$OMP BARRIER
#   if defined BD_EXCHANGE
        if (IstrU.le.Iend) then
          call exchange_r2d_3pts_tile (Istr,Iend,Jstr,Jend,FX_r)
        endif
        if (JstrV.le.Jend) then
          call exchange_r2d_3pts_tile (Istr,Iend,Jstr,Jend,FE_r)
        endif
#   endif

#   ifdef EW_PERIODIC
        imin=1
        imax=LOCALLM+1
#   else
#    ifdef MPI
        if (WESTERN_EDGE) then
          imin=3
        else
          imin=1
        endif
        if (EASTERN_EDGE) then
          imax=Lmmpi-1
        else
          imax=Lmmpi+1
        endif
#    else
        imin=3
        imax=Lm-1
#    endif
#   endif
#   ifdef NS_PERIODIC
        jmin=1
        jmax=LOCALMM+1
#   else
#    ifdef MPI
        if (SOUTHERN_EDGE) then
          jmin=3
        else
          jmin=1
        endif
        if (NORTHERN_EDGE) then
          jmax=Mmmpi-1
        else
          jmax=Mmmpi+1
        endif
#    else
        jmin=3
        jmax=Mm-1
#    endif
#   endif
#  endif /* BEDLOAD_WENO5 || BEDLOAD_UP5 */
!
        do j=Jstr,Jend
          do i=Istr,Iend+1
#  if defined BEDLOAD_WENO5 || defined BEDLOAD_UP5
            cff=SIGN(1.,FX_r(i,j))
            cff1=0.5*(1.+SIGN(1.,FX_r(i,j)))
            cff2=0.5*(1.-SIGN(1.,FX_r(i,j)))
            IF ( (( i.lt.imin ) .or. ( i.gt.imax ))
     &                          .and.( imin.lt.imax ) ) THEN
               FX(i,j)=0.5*(1.+SIGN(1.,FX_r(i-1,j)))*
     &              (cff1*FX_r(i-1,j)+
     &               cff2*0.5*(FX_r(i-1,j)+FX_r(i,j)))+
     &                 0.5*(1.-SIGN(1.,FX_r(i-1,j)))*
     &              (cff2*FX_r(i  ,j)+
     &               cff1*0.5*(FX_r(i-1,j)+FX_r(i,j)))
            ELSE
               FX(i,j)=FLUX5(FX_r(i-3,j),FX_r(i-2,j),
     &                       FX_r(i-1,j),FX_r(i  ,j),
     &                       FX_r(i+1,j),FX_r(i+2,j),cff)
            ENDIF
#  elif defined BEDLOAD_UP1
            cff1=0.5*(1.+SIGN(1.,FX_r(i,j)))
            cff2=0.5*(1.-SIGN(1.,FX_r(i,j)))
            FX(i,j)=0.5*(1.+SIGN(1.,FX_r(i-1,j)))*
     &              (cff1*FX_r(i-1,j)+
     &               cff2*0.5*(FX_r(i-1,j)+FX_r(i,j)))+
     &              0.5*(1.-SIGN(1.,FX_r(i-1,j)))*
     &              (cff2*FX_r(i  ,j)+
     &               cff1*0.5*(FX_r(i-1,j)+FX_r(i,j)))
#  else
            FX(i,j)=0.5*(FX_r(i-1,j)+FX_r(i,j))
#  endif
!
#  ifdef SLOPE_KIRWAN
            dzdx=(h(i,j)-h(i-1,j))/om_u(i,j)
            a_slopex=(MAX(0.,abs(dzdx)-0.01)
     &               *SIGN(1.,dzdx)*0.1)
     &               *om_r(i,j)*dt
            FX(i,j)=FX(i,j)+a_slopex MORPH_FAC
#  endif /* SLOPE_KIRWAN */
#  ifdef MASKING
            FX(i,j)=FX(i,j)*umask(i,j)
#  endif
#  ifdef WET_DRY
            FX(i,j)=FX(i,j)*umask_wet(i,j)
#  endif
          enddo
        enddo
        do j=Jstr,Jend+1
          do i=Istr,Iend
#  if defined BEDLOAD_WENO5 || defined BEDLOAD_UP5
            cff=SIGN(1.,FE_r(i,j))
            cff1=0.5*(1.+SIGN(1.,FE_r(i,j)))
            cff2=0.5*(1.-SIGN(1.,FE_r(i,j)))
            IF ( (( j.lt.jmin ) .or. ( j.gt.jmax ))
     &                          .and.( jmin.lt.jmax ) ) THEN
               FE(i,j)=0.5*(1.+SIGN(1.,FE_r(i,j-1)))*
     &              (cff1*FE_r(i,j-1)+
     &               cff2*0.5*(FE_r(i,j-1)+FE_r(i,j)))+
     &                 0.5*(1.-SIGN(1.,FE_r(i,j-1)))*
     &              (cff2*FE_r(i  ,j)+
     &               cff1*0.5*(FE_r(i,j-1)+FE_r(i,j)))
            ELSE
               FE(i,j)=FLUX5(
     &                      FE_r(i,j-3),FE_r(i,j-2),
     &                      FE_r(i,j-1),FE_r(i  ,j),
     &                      FE_r(i,j+1),FE_r(i,j+2),cff)
            ENDIF
#  elif defined BEDLOAD_UP1
            cff1=0.5*(1.+SIGN(1.,FE_r(i,j)))
            cff2=0.5*(1.-SIGN(1.,FE_r(i,j)))
            FE(i,j)=0.5*(1.+SIGN(1.,FE_r(i,j-1)))*
     &              (cff1*FE_r(i,j-1)+
     &               cff2*0.5*(FE_r(i,j-1)+FE_r(i,j)))+
     &              0.5*(1.-SIGN(1.,FE_r(i,j-1)))*
     &              (cff2*FE_r(i  ,j)+
     &               cff1*0.5*(FE_r(i,j-1)+FE_r(i,j)))
#  else
            FE(i,j)=0.5*(FE_r(i,j-1)+FX_r(i,j))
#  endif
#  ifdef SLOPE_KIRWAN
            dzdy=(h(i,j)-h(i,j-1))/on_v(i,j)
            a_slopey=(MAX(0.,abs(dzdy)-0.01)
     &               *SIGN(1.,dzdy)*0.1)
     &               *on_r(i,j)*dt
            FE(i,j)=FE(i,j)+a_slopey MORPH_FAC
#  endif /* SLOPE_KIRWAN */
#  ifdef MASKING
            FE(i,j)=FE(i,j)*vmask(i,j)
#  endif
#  ifdef WET_DRY
            FE(i,j)=FE(i,j)*vmask_wet(i,j)
#  endif
          enddo
        enddo

#  if defined BEDLOAD_WENO5 || defined BEDLOAD_UP5
#   undef FX_r
#   undef FE_r
#  endif
!
#  ifdef BEDLOAD_SLUMP
!-----------------------------------------------------------------------
!  Sed slump computation to allow slumping (from COAWST).
!     sedslope is the critical slope to allow slumping.
!     slopefac is the scale factor for sediment movement.
!-----------------------------------------------------------------------
!
        sedslope_crit=0.1  ! Critical bed slope for slumping
        slopefac=0.1       ! Bedload bed slumping factor
!
!  U-direction slumping
!
        do j=Jstr,Jend
          do i=Istr,Iend+1
            dzdx=(h(i,j)-h(i-1,j))/om_u(i,j)
            dzdy=(h(i,j)-h(i,j-1))/on_v(i,j)
            dzdxdy=sqrt(dzdx**2.+dzdy**2.)
            cff=Srho(ised)*(1.-bed_poros(i,j,1))
            slopefac_local=cff*dt*slopefac MORPH_FAC
            cff=ABS(dzdxdy)-sedslope_crit
            cff1=(0.5+SIGN(0.5,cff))
            cff=0.5*cff*cff1/(pm(i,j)*pn(i,j))
            cff2=cff*slopefac_local*SIGN(1.,dzdx)
            FX(i,j)=FX(i,j)+cff2
          end do
        end do
!
!  V-direction slumping
!
        do j=Jstr,Jend+1
          do i=Istr,Iend
            dzdy=(h(i,j)-h(i,j-1))/on_v(i,j) ! positive is downhill
            dzdx=(h(i,j)-h(i-1,j))/om_u(i,j)
            dzdxdy=sqrt(dzdx**2.+dzdy**2.)
            cff=Srho(ised)*(1.-bed_poros(i,j,1))
            slopefac_local=cff*dt*slopefac MORPH_FAC
            cff=ABS(dzdxdy)-sedslope_crit
            cff1=(0.5+SIGN(0.5,cff))
            cff=0.5*cff*cff1/(pm(i,j)*pn(i,j))
            cff2=cff*slopefac_local*SIGN(1.,dzdy)
            FE(i,j)=FE(i,j)+cff2
          end do
        end do
#  endif /* BEDLOAD_SLUMP */
!
!-----------------------------------------------------------------------
! Limit fluxes to prevent bottom from breaking thru water surface.
!-----------------------------------------------------------------------
!
        do j=Jstr,Jend
          do i=Istr,Iend
!
! Total thickness available.
!
#  ifdef WET_DRY
            Dstp=MAX(z_w(i,j,N)-z_w(i,j,0)-1.5*Dcrit(i,j),0.)
#  else
            Dstp=MAX(z_w(i,j,N)-z_w(i,j,0)-0.01,0.)  !!!
#  endif
!
! Thickness change that wants to occur.
!
            rhs_bed=(FX(i+1,j)-FX(i,j)+
     &               FE(i,j+1)-FE(i,j))*pm(i,j)*pn(i,j)
            bed_change=rhs_bed/(Srho(ised)*(1.-bed_poros(i,j,1)))
!
! Limit that change to be less than available (if accretional)
!
            cff=MAX(bed_change-Dstp,0.)
            cff1=cff/ABS(bed_change+eps)
            cff1=cff1*0.5*(1.-SIGN(1.,bed_change))
            FX(i+1,j  )=FX(i+1,j  )*(1.-cff1)
            FX(i  ,j  )=FX(i  ,j  )*(1.-cff1)
            FE(i  ,j+1)=FE(i  ,j+1)*(1.-cff1)
            FE(i  ,j  )=FE(i  ,j  )*(1.-cff1)
          end do
        end do
!
!-----------------------------------------------------------------------
! Apply closed boundary conditions on (FX,FE) if any
!-----------------------------------------------------------------------
!
#  ifndef EW_PERIODIC
#   ifndef OBC_WEST
        if (WESTERN_EDGE) then
          do j=Jstr-1,Jend+1
            FX(Istr,j)=0.
          enddo
        endif
#   endif
#   ifndef OBC_EAST
        if (EASTERN_EDGE) then
          do j=Jstr-1,Jend+1
            FX(Iend+1,j)=0.
          enddo
        endif
#   endif
#  endif
#  ifndef NS_PERIODIC
#   ifndef OBC_SOUTH
        if (SOUTHERN_EDGE) then
          do i=Istr-1,Iend+1
            FE(i,Jstr)=0.
          enddo
        endif
#   endif
#   ifndef OBC_NORTH
        if (NORTHERN_EDGE) then
          do i=Istr-1,Iend+1
            FE(i,Jend+1)=0.
          enddo
        endif
#   endif
#  endif
!
!-----------------------------------------------------------------------
!  Determine flux divergence and evaluate change in bed properties.
!-----------------------------------------------------------------------
!
        do j=Jstr,Jend
          do i=Istr,Iend
            cff=(FX(i+1,j)-FX(i,j)+
     &           FE(i,j+1)-FE(i,j))*pm(i,j)*pn(i,j)
#  ifdef SPONGE_SED
     &                             *sed_sponge(i,j)
#  endif

#  ifdef BEDLOAD_SHALLOW_LIMIT
            cff=cff/max(1.,(0.5/(z_w(i,j,N)-z_w(i,j,0)
#   ifdef WET_DRY
     &                           -Dcrit(i,j)+eps
#   endif
     &            )**2))
     &
#  endif

            bed_mass(i,j,1,nnew,ised)=MAX(bed_mass(i,j,1,nstp,ised)-
     &                                    cff,0.)
#  if !defined SUSPLOAD
            if (NLAY.gt.1) then
              do k=2,NLAY
                bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k,nstp,ised)
              enddo
            endif
#  endif
            bed_thick(i,j,1)=MAX(bed_thick(i,j,1)-
     &                           cff/(Srho(ised)*
     &                           (1.-bed_poros(i,j,1))), 0.)
          enddo
        enddo
!
!-----------------------------------------------------------------------
!  Output bedload transports [kg/m2/s].
!-----------------------------------------------------------------------
!
        cff=1./dt
        do j=Jstr,Jend
          do i=Istr,Iend
            bedldu(i,j,ised)=(FX(i,j)-FX(i+1,j))*pm(i,j)*pn(i,j)*cff
          enddo
        enddo
        do j=Jstr,Jend
          do i=Istr,Iend
            bedldv(i,j,ised)=(FE(i,j)-FE(i,j+1))*pm(i,j)*pn(i,j)*cff
          enddo
        enddo
!
      enddo  ! ised =========
!
!-----------------------------------------------------------------------
!  Update mean surface properties.
!-----------------------------------------------------------------------
!
      do j=Jstr,Jend
        do i=Istr,Iend
          cff3=0.0
          do ised=1,NST
            cff3=cff3+bed_mass(i,j,1,nnew,ised)
          end do
          if (cff3.eq.0.0) then
            cff3=eps
          end if
          do ised=1,NST
            bed_frac(i,j,1,ised)=bed_mass(i,j,1,nnew,ised)/cff3
          end do
!
          cff1=1.
          cff2=1.
          cff3=1.
          cff4=1.
          do ised=1,NST
            cff1=cff1*tau_ce_2d(i,j,ised)**bed_frac(i,j,1,ised)
            cff2=cff2*Sd(ised)**bed_frac(i,j,1,ised)
            cff3=cff3*(wsed(ised)+eps)**bed_frac(i,j,1,ised)
            cff4=cff4*Srho(ised)**bed_frac(i,j,1,ised)
          enddo
          taucb(i,j)=cff1
          Ssize(i,j)=cff2
          w_set(i,j)=cff3
          Sdens(i,j)=MAX(cff4,1050.)
        enddo
      enddo
!
! Set boundary conditions
!
#  ifndef EW_PERIODIC
      if (EASTERN_EDGE) then
        do j=Jstr,Jend
           taucb(Iend+1,j)=taucb(Iend,j)
           Ssize(Iend+1,j)=Ssize(Iend,j)
           w_set(Iend+1,j)=w_set(Iend,j)
           Sdens(Iend+1,j)=Sdens(Iend,j)
        enddo
      endif
      if (WESTERN_EDGE) then
        do j=Jstr,Jend
           taucb(Istr-1,j)=taucb(Istr,j)
           Ssize(Istr-1,j)=Ssize(Istr,j)
           w_set(Istr-1,j)=w_set(Istr,j)
           Sdens(Istr-1,j)=Sdens(Istr,j)
        enddo
      endif
#  endif
#  ifndef NS_PERIODIC
      if (NORTHERN_EDGE) then
        do i=Istr,Iend
          taucb  (i,Jend+1)=taucb(i,Jend)
          Ssize  (i,Jend+1)=Ssize(i,Jend)
          w_set  (i,Jend+1)=w_set(i,Jend)
          Sdens  (i,Jend+1)=Sdens(i,Jend)
        enddo
      endif
      if (SOUTHERN_EDGE) then
        do i=Istr,Iend
          taucb  (i,Jstr-1)=taucb(i,Jstr)
          Ssize  (i,Jstr-1)=Ssize(i,Jstr)
          w_set  (i,Jstr-1)=w_set(i,Jstr)
          Sdens  (i,Jstr-1)=Sdens(i,Jstr)
        enddo
      endif
#  endif
#  if !defined EW_PERIODIC && !defined NS_PERIODIC
      if (SOUTHERN_EDGE.and.WESTERN_EDGE) then
        taucb(Istr-1,Jstr-1)=0.5*(taucb(Istr,Jstr-1)+taucb(Istr-1,Jstr))
        Ssize(Istr-1,Jstr-1)=0.5*(Ssize(Istr,Jstr-1)+Ssize(Istr-1,Jstr))
        w_set(Istr-1,Jstr-1)=0.5*(w_set(Istr,Jstr-1)+w_set(Istr-1,Jstr))
        Sdens(Istr-1,Jstr-1)=0.5*(Sdens(Istr,Jstr-1)+Sdens(Istr-1,Jstr))
      endif
      if (SOUTHERN_EDGE.and.EASTERN_EDGE) then
        taucb(Iend+1,Jstr-1)=0.5*(taucb(Iend,Jstr-1)+taucb(Iend+1,Jstr))
        Ssize(Iend+1,Jstr-1)=0.5*(Ssize(Iend,Jstr-1)+Ssize(Iend+1,Jstr))
        w_set(Iend+1,Jstr-1)=0.5*(w_set(Iend,Jstr-1)+w_set(Iend+1,Jstr))
        Sdens(Iend+1,Jstr-1)=0.5*(Sdens(Iend,Jstr-1)+Sdens(Iend+1,Jstr))
      endif
      if (NORTHERN_EDGE.and.WESTERN_EDGE) then
        taucb(Istr-1,Jend+1)=0.5*(taucb(Istr,Jend+1)+taucb(Istr-1,Jend))
        Ssize(Istr-1,Jend+1)=0.5*(Ssize(Istr,Jend+1)+Ssize(Istr-1,Jend))
        w_set(Istr-1,Jend+1)=0.5*(w_set(Istr,Jend+1)+w_set(Istr-1,Jend))
        Sdens(Istr-1,Jend+1)=0.5*(Sdens(Istr,Jend+1)+Sdens(Istr-1,Jend))
      endif
      if (NORTHERN_EDGE.and.EASTERN_EDGE) then
        taucb(Iend+1,Jend+1)=0.5*(taucb(Iend,Jend+1)+taucb(Iend+1,Jend))
        Ssize(Iend+1,Jend+1)=0.5*(Ssize(Iend,Jend+1)+Ssize(Iend+1,Jend))
        w_set(Iend+1,Jend+1)=0.5*(w_set(Iend,Jend+1)+w_set(Iend+1,Jend))
        Sdens(Iend+1,Jend+1)=0.5*(Sdens(Iend,Jend+1)+Sdens(Iend+1,Jend))
      endif
#  endif

#  if defined BD_EXCHANGE
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 taucb(START_2D_ARRAY))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 Ssize(START_2D_ARRAY))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 w_set(START_2D_ARRAY))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 Sdens(START_2D_ARRAY))
#  endif
!
!-----------------------------------------------------------------------
! Set boundary conditions
!-----------------------------------------------------------------------
!
#  ifndef EW_PERIODIC
      if (EASTERN_EDGE) then
        do j=Jstr,Jend
           do k=1,NLAY
             do ised=1,NST
                bed_mass(Iend+1,j,k,nnew,ised)=
     &              bed_mass(Iend,j,k,nnew,ised)
                bed_frac(Iend+1,j,k,ised)=
     &              bed_frac(Iend,j,k,ised)
             enddo
           enddo
           bed_thick(Iend+1,j,1)=bed_thick(Iend,j,1)
         enddo
      endif
      if (WESTERN_EDGE) then
        do j=Jstr,Jend
          do k=1,NLAY
            do ised=1,NST
                bed_mass(Istr-1,j,k,nnew,ised)=
     &              bed_mass(Istr,j,k,nnew,ised)
                bed_frac(Istr-1,j,k,ised)=
     &              bed_frac(Istr,j,k,ised)
             enddo
           enddo
          bed_thick(Istr-1,j,1)=bed_thick(Istr,j,1)
        enddo
      endif
#  endif
#  ifndef NS_PERIODIC
      if (NORTHERN_EDGE) then
        do i=Istr,Iend
          do k=1,NLAY
            do ised=1,NST
               bed_mass(i,Jend+1,k,nnew,ised)=
     &              bed_mass(i,Jend,k,nnew,ised)
               bed_frac(i,Jend+1,k,ised)=
     &              bed_frac(i,Jend,k,ised)
            enddo
          enddo
          bed_thick(i,Jend+1,1)=bed_thick(i,Jend,1)
        enddo
      endif
      if (SOUTHERN_EDGE) then
        do i=Istr,Iend
          do k=1,NLAY
            do ised=1,NST
               bed_mass(i,Jstr-1,k,nnew,ised)=
     &              bed_mass(i,Jstr,k,nnew,ised)
               bed_frac(i,Jstr-1,k,ised)=
     &              bed_frac(i,Jstr,k,ised)
            enddo
          enddo
          bed_thick(i,Jstr-1,1)=bed_thick(i,Jstr,1)
        enddo
      endif
#  endif
#  if !defined EW_PERIODIC && !defined NS_PERIODIC
      if (SOUTHERN_EDGE.and.WESTERN_EDGE) then
        bed_mass(Istr-1,Jstr-1,:,nnew,:)=
     &  0.5*(bed_mass(Istr,Jstr-1,:,nnew,:)
     &      +bed_mass(Istr-1,Jstr,:,nnew,:))
        bed_frac(Istr-1,Jstr-1,:,:)=
     &  0.5*(bed_frac(Istr,Jstr-1,:,:)
     &      +bed_frac(Istr-1,Jstr,:,:))
        bed_thick(Istr-1,Jstr-1,1)=
     &  0.5*(bed_thick(Istr,Jstr-1,1)
     &      +bed_thick(Istr-1,Jstr,1))
      endif
      if (SOUTHERN_EDGE.and.EASTERN_EDGE) then
        bed_mass(Iend+1,Jstr-1,:,nnew,:)=
     &  0.5*(bed_mass(Iend,Jstr-1,:,nnew,:)
     &      +bed_mass(Iend+1,Jstr,:,nnew,:))
        bed_frac(Iend+1,Jstr-1,:,:)=
     &  0.5*(bed_frac(Iend,Jstr-1,:,:)
     &      +bed_frac(Iend+1,Jstr,:,:))
        bed_thick(Iend+1,Jstr-1,1)=
     &  0.5*(bed_thick(Iend,Jstr-1,1)
     &      +bed_thick(Iend+1,Jstr,1))
      endif
      if (NORTHERN_EDGE.and.WESTERN_EDGE) then
        bed_mass(Istr-1,Jend+1,:,nnew,:)=
     &  0.5*(bed_mass(Istr,Jend+1,:,nnew,:)
     &      +bed_mass(Istr-1,Jend,:,nnew,:))
        bed_frac(Istr-1,Jend+1,:,:)=
     &  0.5*(bed_frac(Istr,Jend+1,:,:)
     &      +bed_frac(Istr-1,Jend,:,:))
        bed_thick(Istr-1,Jend+1,1)=
     &  0.5*(bed_thick(Istr,Jend+1,1)
     &      +bed_thick(Istr-1,Jend,1))
      endif
      if (NORTHERN_EDGE.and.EASTERN_EDGE) then
        bed_mass(Iend+1,Jend+1,:,nnew,:)=
     &   0.5*(bed_mass(Iend,Jend+1,:,nnew,:)
     &       +bed_mass(Iend+1,Jend,:,nnew,:))
        bed_frac(Iend+1,Jend+1,:,:)=
     &   0.5*(bed_Frac(Iend,Jend+1,:,:)
     &       +bed_Frac(Iend+1,Jend,:,:))
        bed_thick(Iend+1,Jend+1,1)=
     &   0.5*(bed_thick(Iend,Jend+1,1)
     &       +bed_thick(Iend+1,Jend,1))
      endif
#  endif

#  if defined BD_EXCHANGE
      do ised=1,NST
        do k=1,NLAY
          call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_mass(START_2D_ARRAY,k,nnew,ised))
          call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_frac(START_2D_ARRAY,k,ised))
        enddo
        call exchange_u2d_tile (Istr,Iend,Jstr,Jend,
     &                 bedldu(START_2D_ARRAY,ised))
        call exchange_v2d_tile (Istr,Iend,Jstr,Jend,
     &                 bedldv(START_2D_ARRAY,ised))
      enddo
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_thick(START_2D_ARRAY,1))
#  endif

      return
      end
# endif /* BEDLOAD */

# ifdef SUSPLOAD
!======================================================================
!                                                                      !
!  This routine computes the vertical settling (sinking) of suspended  !
!  sediment via a semi-Lagrangian advective flux algorithm. It uses a  !
!  parabolic,  vertical reconstructuion of the suspended  sediment in  !
!  the water column with PPT/WENO constraints to avoid oscillations.   !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Colella, P. and P. Woodward, 1984: The piecewise parabolic method   !
!    (PPM) for gas-dynamical simulations, J. Comp. Phys., 54, 174-201. !
!                                                                      !
!  Liu, X.D., S. Osher, and T. Chan, 1994: Weighted essentially        !
!    nonoscillatory shemes, J. Comp. Phys., 115, 200-212.              !
!                                                                      !
!  Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G.    !
!    Arango, 2008:  Development of a three-dimensional,  regional,     !
!    coupled wave, current, and sediment-transport model, Computers    !
!    & Geosciences, 34, 1284-1306.                                     !
!                                                                      !
!=======================================================================
!
      subroutine sed_settling_tile (Istr,Iend,Jstr,Jend,
     &                              Hz_inv,Hz_inv2,Hz_inv3,qc,
     &                              FC,qR,qL,WR,WL,tau_w,ksource)
!
      implicit none
#  include "param.h"
#  include "ocean3d.h"
#  include "sediment.h"
#  include "bbl.h"
#  include "scalars.h"
#  include "grid.h"
#  include "mixing.h"
#  ifndef BBL
#   include "forces.h"
#  endif
#  ifdef DIAGNOSTICS_TS
#   include "diagnostics.h"
#  endif
!
      integer Istr,Iend,Jstr,Jend
      real Hz_inv(PRIVATE_1D_SCRATCH_ARRAY,N),
     &    Hz_inv2(PRIVATE_1D_SCRATCH_ARRAY,N),
     &    Hz_inv3(PRIVATE_1D_SCRATCH_ARRAY,N),
     &         qc(PRIVATE_1D_SCRATCH_ARRAY,N),
     &         FC(PRIVATE_1D_SCRATCH_ARRAY,0:N),
     &         qR(PRIVATE_1D_SCRATCH_ARRAY,N),
     &         qL(PRIVATE_1D_SCRATCH_ARRAY,N),
     &         WR(PRIVATE_1D_SCRATCH_ARRAY,N),
     &         WL(PRIVATE_1D_SCRATCH_ARRAY,N)
      real tau_w(PRIVATE_2D_SCRATCH_ARRAY)
      integer  ksource(PRIVATE_1D_SCRATCH_ARRAY,N),
     &         ised,indx,i,j,k,ks
      real dltR, dltL,cff,cffR,cffL,cu
!
!-----------------------------------------------------------------------
!  Compute bottom skin-friction stress due to combined maximum wave
!  and current interaction
!-----------------------------------------------------------------------
!
      do j=Jstr,Jend
        do i=Istr,Iend
          tau_w(i,j)=sqrt(bustrw(i,j)**2
     &                   +bvstrw(i,j)**2)
#  ifdef WET_DRY
          tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
#  endif
        enddo
      enddo
      do j=Jstr,Jend                        ! Auxiliary step: save
        do k=1,N                            ! inverses of grid box
          do i=Istr,Iend                    ! heights to avoid
            Hz_inv(i,k)=1./Hz(i,j,k)        ! repeated divisions
          enddo
        enddo
        do k=1,N-1
          do i=Istr,Iend
            Hz_inv2(i,k)=1./(Hz(i,j,k)+Hz(i,j,k+1))
          enddo
        enddo
        do k=2,N-1
          do i=Istr,Iend
            Hz_inv3(i,k)=1./(Hz(i,j,k-1)+Hz(i,j,k)+Hz(i,j,k+1))
          enddo
        enddo

!
! Vertical sinking of suspended particles:   Copy concentration of
!========= ======= == ========= ==========   suspended sediment into
! scratch array "qc" (q-central, restrict it to be positive) which
! are hereafter interpreted as a set of grid-box averaged values for
! concentration. then reconstruct vertical profile of "qc" in terms
! of a set of parabolic segments within each grid box; and, finally,
! compute semi-Lagrangian flux due to sinking.
!
        do ised=1,NST
          indx=itrc_sed+ised-1
          do k=1,N
            do i=Istr,Iend
              qc(i,k)=t(i,j,k,nnew,indx)
            enddo
          enddo

          do k=N-1,1,-1
            do i=Istr,Iend
              FC(i,k)=(qc(i,k+1)-qc(i,k))*Hz_inv2(i,k)
            enddo
          enddo
          do k=2,N-1
            do i=Istr,Iend
              dltR=Hz(i,j,k)*FC(i,k)
              dltL=Hz(i,j,k)*FC(i,k-1)
              cff=Hz(i,j,k-1)+2.*Hz(i,j,k)+Hz(i,j,k+1)
              cffR=cff*FC(i,k)
              cffL=cff*FC(i,k-1)            ! Apply PPM monotonicity
                                            ! constraint to prevent
              if (dltR*dltL .le. 0.) then   ! oscillation within the
                dltR=0.                     ! grid box
                dltL=0.
              elseif (abs(dltR) .gt. abs(cffL)) then
                dltR=cffL
              elseif (abs(dltL) .gt. abs(cffR)) then
                dltL=cffR
              endif                         ! Compute right and left
                                            ! side values qR,qL of
              cff=(dltR-dltL)*Hz_inv3(i,k)  ! parabolic segments
              dltR=dltR-cff*Hz(i,j,k+1)     ! within grid box Hz(k)
              dltL=dltL+cff*Hz(i,j,k-1)     ! (WR,WL are measures of
              qR(i,k)=qc(i,k)+dltR          ! quadratic variations).
              qL(i,k)=qc(i,k)-dltL
              WR(i,k)=( 2.*dltR-dltL )**2   ! NOTE: Although each
              WL(i,k)=( dltR-2.*dltL )**2   ! parabolic segment is
            enddo                           ! monotone within its
          enddo          !--> discard FC    ! grid box, monotonicity
                                            ! of the whole profile is
          cff=1.0E-14                       ! not guaranteed, because
          do k=2,N-2                        ! qL(k+1)-qR(k) may still
            do i=Istr,Iend                  ! have different sign
              dltL=max(WL(i,k),   cff)      ! than qc(k+1)-qc(k)...
              dltR=max(WR(i,k+1), cff)
              qR(i,k)=(dltR*qR(i,k)+dltL*qL(i,k+1))/(dltR+dltL)
              qL(i,k+1)=qR(i,k)
            enddo                           !   ...this possibility
          enddo      !--> discard WR,WL     ! is excluded, after qL
                                            ! and qR are reconciled
          do i=Istr,Iend                    ! using WENO procedure.
            FC(i,N)=0.     !<-- no-flux BC
#  if defined LINEAR_CONTINUATION
            qL(i,N)=qR(i,N-1)
            qR(i,N)=2.*qc(i,N)-qL(i,N)
#  elif defined NEUMANN
            qL(i,N)=qR(i,N-1)
            qR(i,N)=1.5*qc(i,N)-0.5*qL(i,N)
#  else
            qR(i,N)=qc(i,N)                 ! Strictly monotone
            qL(i,N)=qc(i,N)                 ! version as the default:
            qR(i,N-1)=qc(i,N)               ! distributions at top...
#  endif
#  if defined LINEAR_CONTINUATION
            qR(i,1)=qL(i,2)
            qL(i,1)=2.*qc(i,1)-qR(i,1)
#  elif defined NEUMANN
            qR(i,1)=qL(i,2)
            qL(i,1)=1.5*qc(i,1)-0.5*qR(i,1)
#  else
            qL(i,2)=qc(i,1)                 !     ...and bottom grid
            qR(i,1)=qc(i,1)                 ! boxes re assumed to be
            qL(i,1)=qc(i,1)                 ! piecewise constant.
#  endif
          enddo

          do k=1,N                          ! Since the reconciled
            do i=Istr,Iend                  ! interfacial values may
              dltR=qR(i,k)-qc(i,k)          ! cause non-monotonic
              dltL=qc(i,k)-qL(i,k)          ! behavior of parabolic
              cffR=2.*dltR                  ! segments inside grid
              cffL=2.*dltL                  ! box apply monotonicity
                                            ! constraint again.
              if (dltR*dltL .lt. 0.) then
                dltR=0.
                dltL=0.
              elseif (abs(dltR) .gt. abs(cffL)) then
                dltR=cffL
              elseif (abs(dltL) .gt. abs(cffR)) then
                dltL=cffR
              endif

              qR(i,k)=qc(i,k)+dltR
              qL(i,k)=qc(i,k)-dltL
            enddo
          enddo     !--> discard everything, except qR,qL
!
! After this moment reconstruction is considered complete. The next
! stage is to compute vertical advective fluxes FC. It is expected
! that sinking may occurs relatively fast, the algorithm is designed
! to be free of CFL criterion, which is achieved by allowing
! integration bounds for semi-Lagrangian advective flux to use as
! many grid boxes in upstream direction as necessary.
!
          cff=dt*abs(Wsed(ised))          ! In the two code segments
          do k=1,N                        ! WL is z-coordinate of the
            do i=Istr,Iend                ! departure point for grid
              FC(i,k-1)=0.                ! box interface z_w with
              WL(i,k)=z_w(i,j,k-1)+cff    ! the same indices;
#  ifdef WET_DRY
     &                *rmask_wet(i,j)
#  endif
              WR(i,k)=Hz(i,j,k)*qc(i,k)   ! FC is finite volume flux;
              ksource(i,k)=k              ! ksource(:,k) is index of
            enddo                         ! vertical grid box which
          enddo                           ! contains the departure
          do k=1,N                        ! point (restricted by N);
            do ks=k,N-1                   ! During the search: also
              do i=Istr,Iend
                if (WL(i,k) .gt. z_w(i,j,ks)) then
                  ksource(i,k)=ks+1
                  FC(i,k-1)=FC(i,k-1)+WR(i,ks)
                endif
              enddo                       ! add in content of whole
            enddo                         ! grid boxes participating
          enddo       !--> discard WR     ! in FC.

          do k=1,N                        ! Finalize computation of
            do i=Istr,Iend                ! flux: add fractional part
              ks=ksource(i,k)
              cu=min(1.,(WL(i,k)-z_w(i,j,ks-1))*Hz_inv(i,ks))
              FC(i,k-1)=FC(i,k-1) + Hz(i,j,ks)*cu*( qL(i,ks)
     &                             +cu*( 0.5*(qR(i,ks)-qL(i,ks))
     &             -(1.5-cu)*(qR(i,ks)+qL(i,ks)-2.*qc(i,ks))  ))
            enddo
          enddo
          cff=1/tau_cd(ised)
          do i=Istr,Iend
            !FC(i,0)=FC(i,0)*min(1.,REAL(int(cff*tau_w(i,j))))
            do k=1,N,+1
               t(i,j,k,nnew,indx)=qc(i,k)
     &                   +(FC(i,k)-FC(i,k-1))*Hz_inv(i,k)
#  ifdef SPONGE_SED
     &                                       *sed_sponge(i,j)
#  endif
            enddo
            settling_flux(i,j,ised)=FC(i,0)
#  ifdef SPONGE_SED
     &                                       *sed_sponge(i,j)
#  endif
          enddo
!
!  Diagnostics for settling: use TForc array
!
#  ifdef DIAGNOSTICS_TS
          do k=1,N
            do i=Istr,Iend
              TForc(i,j,k,indx)=(FC(i,k)-FC(i,k-1))/dt
#   ifdef MASKING
     &                                    * rmask(i,j)
#   endif
            enddo
          enddo
#  endif /* DIAGNOSTICS_TS */

        enddo   ! ised loop
      enddo     ! j loop

!
!  Update Trate diagnostics
!
#  ifdef DIAGNOSTICS_TS
      do ised=1,NST
        indx=itrc_sed+ised-1
        do k=1,N
          do j=Jstr,Jend
            do i=Istr,Iend
              Trate(i,j,k,indx)=(Hz(i,j,k)*t(i,j,k,nnew,indx)
     &                          -Hz_bak(i,j,k)*t(i,j,k,nstp,indx))
     &                                                         /dt
#   ifdef MASKING
     &                                                 *rmask(i,j)
#   endif
            enddo
          enddo
        enddo
      enddo
#  endif /* DIAGNOSTICS_TS */
!
!  Exchange for output purpose
!
#  if defined BD_EXCHANGE
      do ised=1,NST
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                          settling_flux(START_2D_ARRAY,ised))
      enddo
#  endif

      return
      end
!
!=======================================================================
!                                                                      !
!  This computes sediment erosion flux and resuspended load.           !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G.    !
!    Arango, 2008:  Development of a three-dimensional,  regional,     !
!    coupled wave, current, and sediment-transport model, Computers    !
!    & Geosciences, 34, 1284-1306.                                     !
!                                                                      !
!=======================================================================
!
      subroutine sed_fluxes_tile (Istr,Iend,Jstr,Jend, tau_w)
!
      implicit none
#  include "param.h"
#  include "grid.h"
#  include "bbl.h"
#  include "forces.h"
#  include "sediment.h"
#  include "ocean3d.h"
#  include "scalars.h"
#  include "mixing.h"
#  ifdef DIAGNOSTICS_TS
#   include "diagnostics.h"
#  endif
      integer Istr,Iend,Jstr,Jend
      integer i,j,ised,ised0
      integer indx
      integer imin,jmin,imax,jmax
      real cff,cff1,cff2,cff3,cff4,cff5,tau_water,Zr
      real tau_w(PRIVATE_2D_SCRATCH_ARRAY)
#  ifdef BED_ARMOR
      integer ised0
      real ph,pe
#  endif
!
!-----------------------------------------------------------------------
!  Compute bottom skin-friction stress due to combined maximum wave
!  and current interaction
!-----------------------------------------------------------------------
!
      do j=Jstr,Jend
        do i=Istr,Iend
          tau_w(i,j)=sqrt(bustrw(i,j)**2
     &                   +bvstrw(i,j)**2)
#  ifdef WET_DRY
          tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
#  endif
        enddo
      enddo
!
!-----------------------------------------------------------------------
!  Compute erosion flux and resuspended load
!-----------------------------------------------------------------------
!
      do j=jstr,jend
        do ised=1,NST
          do i=istr,iend
            indx=itrc_sed+ised-1
!
!  Calculate critical shear stress in Pa
!  => change for COHESIVE bed here
!

#  if defined COHESIVE_BED
            cff = rho0/bed_bctr(i,j,1)
#  elif defined MIXED_BED
            cff = MAX(bot_dprp(i,j)*bed_bctr(i,j,1)/rho0+
     &            (1.-bot_dprp(i,j))*tau_ce(ised),
     &            tau_ce(ised))
            cff=1./cff
#  else
            cff=1./tau_ce_2d(i,j,ised)
#  endif
!
!  Compute erosion, ero_flux (kg/m2).
!
            cff1=(1.0-bed_poros(i,j,1))*bed_frac(i,j,1,ised)
            cff2=dt*Erate(ised)*cff1
            cff3=Srho(ised)*cff1
            cff4=bed_mass(i,j,1,bnew,ised)
            cff5=settling_flux(i,j,ised)
#  if defined SED_TAU_CD_CONST || defined SED_TAU_CD_LIN
             if (tau_w(i,j).GT.tau_cd(ised) ) then
                cff5=0.0
             endif
#  endif
#  if defined SED_TAU_CD_LIN
             if (tau_w(i,j).LE.tau_cd(ised).AND.
     &                       tau_cd(ised).GT.0.0) then
                cff5=(1.0-tau_w(i,j)/tau_cd(ised))*cff5
             endif
#  endif
            ero_flux(i,j,ised)=
     &            MIN(MAX(0.,cff2*(cff*tau_w(i,j)-1.)),
     &                MIN(cff3*bot_thick(i,j) MORPH_FAC ,cff4)+
     &                cff5)
            ero_flux(i,j,ised)=ero_flux(i,j,ised)
#  ifdef SPONGE_SED
     &                             *sed_sponge(i,j)
#  endif
!
!  Update global tracer variables (m Tunits for nnew indx, Tuints for 3)
!  for erosive flux.
!
            t(i,j,1,nnew,indx)=t(i,j,1,nnew,indx)+ero_flux(i,j,ised)
     &                                                    /Hz(i,j,1)
#  if defined SED_TAU_CD_CONST || defined SED_TAU_CD_LIN
     &                       + (settling_flux(i,j,ised)-cff5)/Hz(i,j,1)
            settling_flux(i,j,ised)=cff5
#  endif
          enddo  ! i loop

!
!  Update erosion/settling diagnostics (bottom level only)
!
#  ifdef DIAGNOSTICS_TS
          do i=Istr,Iend
            TForc(i,j,1,indx) = TForc(i,j,1,indx) +
     &                          ero_flux(i,j,ised)/dt
#   ifdef MASKING
     &                                   * rmask(i,j)
#   endif
          enddo
#  endif /* DIAGNOSTICS_TS */

        enddo    ! ised loop
      enddo      ! j loop

!
!  Update Trate diagnostics (bottom level only)
!
#  ifdef DIAGNOSTICS_TS
      do ised=1,NST
        indx=itrc_sed+ised-1
        do j=Jstr,Jend
          do i=Istr,Iend
            Trate(i,j,1,indx)=(Hz(i,j,1)*t(i,j,1,nnew,indx)
     &                        -Hz_bak(i,j,1)*t(i,j,1,nstp,indx))
     &                                                       /dt
#   ifdef MASKING
     &                                               *rmask(i,j)
#   endif
          enddo
        enddo
      enddo
#  endif /* DIAGNOSTICS_TS */
!
!  Exchange for output purpose
!
#  if defined BD_EXCHANGE
      do ised=1,NST
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                          ero_flux(START_2D_ARRAY,ised))
      enddo
#   if defined SED_TAU_CD_CONST || defined SED_TAU_CD_LIN
      do ised=1,NST
        call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                          settling_flux(START_2D_ARRAY,ised))
      enddo
#   endif
#  endif

      return
      end
# endif /* SUSPLOAD */
!
!=======================================================================
!                                                                      !
!  This routine computes sediment bed layer stratigraphy.              !
!                                                                      !
!  Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G.    !
!    Arango, 2008:  Development of a three-dimensional,  regional,     !
!    coupled wave, current, and sediment-transport model, Computers    !
!    & Geosciences, 34, 1284-1306.                                     !
!                                                                      !
!=======================================================================
!
# if !defined COHESIVE_BED && !defined MIXED_BED

      subroutine sed_bed_tile (Istr,Iend,Jstr,Jend,
     &                         tau_w,dep_mass,Awrk1,Awrk2)
!
      implicit none
#  include "param.h"
#  include "bbl.h"
#  include "forces.h"
#  include "grid.h"
#  include "ocean3d.h"
#  include "scalars.h"
#  include "sediment.h"
      integer Istr,Iend,Jstr,Jend
      integer i,j,k,ks,ised,Ksed, NLAYm1
      real tau_w(PRIVATE_2D_SCRATCH_ARRAY)
      real dep_mass(PRIVATE_1D_SCRATCH_ARRAY,NST)
      real Awrk1(PRIVATE_2D_SCRATCH_ARRAY)
      real Awrk2(PRIVATE_2D_SCRATCH_ARRAY)
      real eps,cff,cff1,cff2,cff3
      real thck_avail,thck_to_add,newlayer_thick
      integer imin,imax,jmin,jmax
      parameter ( eps=1.D-20, newlayer_thick=0.001 )
!
!-----------------------------------------------------------------------
!  Compute bottom skin-friction stress due to combined maximum wave
!  and current interaction
!-----------------------------------------------------------------------
!
      do j=Jstr,Jend
        do i=Istr,Iend
          tau_w(i,j)=sqrt(bustrw(i,j)**2
     &                   +bvstrw(i,j)**2)
!
!-----------------------------------------------------------------------
!  Calculate active layer thickness
!     --> Harris and Wiberg 1997: za = k1 (tau-tc)rho + k2 Sd50
!-----------------------------------------------------------------------
!
          bot_thick(i,j)=MAX(newlayer_thick,
     &                       0.007*(tau_w(i,j)-taucb(i,j))*rho0)+
     &                       6.*Ssize(i,j)
        enddo
      enddo
!-----------------------------------------------------------------------
!  Update bed properties according to ero_flux and dep_flux.
!-----------------------------------------------------------------------
!
!
      NLAYm1=max(NLAY-1,1)  ! for NLAY=1 case
!
      do j=Jstr,Jend
!
        do ised=1,NST
          do i=Istr,Iend
            dep_mass(i,ised)=0.0

#  ifdef SUSPLOAD
!  Suspended erosion/deposition. Apply morphology factor.
!
!  The deposition and resuspension of sediment on the bottom "bed"
!  is due to precipitation flux, already computed, and the
!  resuspension (erosion, hence called ero_flux). The resuspension is
!  applied to the bottom-most grid box value qc(:,1) so the total mass
!  is conserved. Restrict "ero_flux" so that "bed" cannot go negative
!  after both fluxes are applied.
!
            ero_flux(i,j,ised)=ero_flux(i,j,ised) MORPH_FAC
            settling_flux(i,j,ised)=settling_flux(i,j,ised) MORPH_FAC
#  endif
!
! Net deposited sediment, including bedload
!
            cff=
#  ifdef SUSPLOAD
     &          settling_flux(i,j,ised)-ero_flux(i,j,ised)
#  endif
#  ifdef BEDLOAD
     &          +bed_mass(i,j,1,nnew,ised)-bed_mass(i,j,1,nstp,ised)
#  endif

            if ( cff .gt. 0. ) then
!
!  If first time step of deposit, then store deposit material in
!  temporary array, dep_mass.
!
!  Use initial resolution Bthk(1) to trigger new layer formation
!  (prevents stacking all layers near the surface - PM).
!
              if ( time .gt. (bed_age(i,j,1)+1.1*dt) .and.
!     &             bed_thick(i,j,1) .gt. newlayer_thick ) then
     &             bed_thick(i,j,1) .gt. Bthk(1) ) then
                dep_mass(i,ised)=cff
              endif
              !bed_age(i,j,1)=time   !!!???
            endif
!
!  Update bed mass arrays.
!
            bed_mass(i,j,1,nnew,ised)=MAX(bed_mass(i,j,1,nstp,ised)+
     &                                    cff, 0.)
            do k=2,NLAY
              bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k,nstp,ised)
            enddo
          enddo   ! i loop
        enddo     ! ised loop
!
!  If first time step of deposit, create new layer and combine bottom
!  two bed layers.
!
        do i=Istr,Iend
          cff=0.
!
!  Determine if deposition ocurred here.
!
          if (NLAY.gt.1) then
!
            do ised=1,NST
              cff=cff+dep_mass(i,ised)
            enddo
!
            if (cff.gt.0.) then
!
!  Combine bottom layers.
!
              bed_poros(i,j,NLAY)=0.5*(bed_poros(i,j,NLAYm1)+
     &                                 bed_poros(i,j,NLAY))
              bed_age  (i,j,NLAY)=0.5*(bed_age  (i,j,NLAYm1)+
     &                                 bed_age  (i,j,NLAY))
              do ised=1,NST
                bed_mass(i,j,NLAY,nnew,ised)=
     &                             bed_mass(i,j,NLAYm1,nnew,ised)+
     &                             bed_mass(i,j,NLAY  ,nnew,ised)
              enddo
!
!  Push layers down.
!
              do k=NLAY-1,2,-1
                bed_poros(i,j,k)=bed_poros(i,j,k-1)
                bed_age  (i,j,k)=bed_age  (i,j,k-1)
                do ised =1,NST
                  bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k-1,nnew,ised)
                enddo
              enddo
!
!  Set new top layer parameters.
!
              k=min(2,NLAY) ! for compilation if NLAY=1
              do ised=1,NST
                bed_mass(i,j,k,nnew,ised)=MAX(bed_mass(i,j,k,nnew,ised)-
     &                                        dep_mass(i,ised),0.)
                bed_mass(i,j,1,nnew,ised)=dep_mass(i,ised)
              enddo
              bed_age(i,j,1)=time  !!! update here
!
            endif ! new deposition cff>0
          endif  ! NLAY=1
!
!  Recalculate thickness and fractions for all layers.
!
          do k=1,NLAY
            cff3=0.
            do ised=1,NST
              cff3=cff3+bed_mass(i,j,k,nnew,ised)
            enddo
            if (cff3.eq.0.) then
              cff3=eps
            endif
            bed_thick(i,j,k)=0.
            do ised=1,NST
              bed_frac(i,j,k,ised)=bed_mass(i,j,k,nnew,ised)/cff3
              bed_thick(i,j,k)=MAX(bed_thick(i,j,k)+
     &                         bed_mass(i,j,k,nnew,ised)/
     &                         (Srho(ised)*(1.-bed_poros(i,j,k))),0.)
            enddo  ! ised loop
          enddo    ! k loop
        enddo      ! i loop
      enddo        ! j loop
!
!  Ensure top bed layer thickness is greater or equal than active layer
!  thickness. If need to add sed to top layer, then entrain from lower
!  levels. Create new layers at bottom to maintain Nlay
!
      do j=Jstr,Jend
        do i=Istr,Iend
          if (bot_thick(i,j).gt.bed_thick(i,j,1)) then
            if (NLAY.eq.1) then
              bot_thick(i,j)=bed_thick(i,j,1)
            else
              thck_to_add=bot_thick(i,j)-bed_thick(i,j,1)
              thck_avail=0.
              Ksed=1       ! initialize
              do k=2,NLAY
                if (thck_avail.lt.thck_to_add) then
                  thck_avail=thck_avail+bed_thick(i,j,k)
                  Ksed=k
                endif
              enddo
!
!  Catch here if there was not enough bed material.
!
              if (thck_avail.lt.thck_to_add) then
                bot_thick(i,j)=bed_thick(i,j,1)+thck_avail
                thck_to_add=thck_avail
              endif
!
!  Update bed mass of top layer and fractional layer.
!
              cff2=MAX(thck_avail-thck_to_add,0.)/
     &             MAX(bed_thick(i,j,Ksed),eps)
              do ised=1,NST
                cff1=0.
                do k=1,Ksed
                  cff1=cff1+bed_mass(i,j,k,nnew,ised)
                enddo
                cff3=cff2*bed_mass(i,j,Ksed,nnew,ised)
                bed_mass(i,j,1   ,nnew,ised)=cff1-cff3
                bed_mass(i,j,Ksed,nnew,ised)=cff3
              enddo
!
!  Update thickness of fractional layer ksource_sed
!
              bed_thick(i,j,Ksed)=MAX(thck_avail-thck_to_add,0.)
!
!  Upate bed fraction of top layer.
!
              cff3=0.
              do ised=1,NST
                cff3=cff3+bed_mass(i,j,1,nnew,ised)
              enddo
              if (cff3.eq.0.) then
                cff3=eps
              endif
              do ised=1,NST
                bed_frac(i,j,1,ised)=bed_mass(i,j,1,nnew,ised)/cff3
              enddo
!
!  Upate bed thickness of top layer.
!
              bed_thick(i,j,1)=bot_thick(i,j)
!
!  Pull all layers closer to the surface.
!
              ks=Ksed-2  ! ks is the number of new layers
              do k=Ksed,NLAY
                bed_thick(i,j,k-ks)=bed_thick(i,j,k)
                bed_poros(i,j,k-ks)=bed_poros(i,j,k)
                bed_age  (i,j,k-ks)=bed_age  (i,j,k)
                do ised=1,NST
                  bed_frac(i,j,k-ks,ised)=bed_frac(i,j,k,ised)
                  bed_mass(i,j,k-ks,nnew,ised)=bed_mass(i,j,k,nnew,ised)
                enddo
              enddo
!
!  Add new layers onto the bottom. Split what was in the bottom layer to
!  fill these new empty cells.
!
              !ks=0
              cff=1./REAL(ks+1)
              do k=NLAY,NLAY-ks,-1
                bed_thick(i,j,k)=bed_thick(i,j,NLAY-ks)*cff
                bed_age  (i,j,k)=bed_age  (i,j,NLAY-ks)
                do ised=1,NST
                  bed_frac(i,j,k,ised)=bed_frac(i,j,NLAY-ks,ised)
                  bed_mass(i,j,k,nnew,ised)=
     &                             bed_mass(i,j,NLAY-ks,nnew,ised)*cff
                enddo   ! ised loop
              enddo     ! k loop
            endif       ! NLAY > 1
          endif         ! increase top bed layer
        enddo           ! i loop
      enddo             ! j loop
!
! Set boundary conditions
!
#  ifndef EW_PERIODIC
      if (EASTERN_EDGE) then
        do j=Jstr,Jend
          do k=1,NLAY
            do ised=1,NST
            bed_mass(Iend+1,j,k,nnew,ised)=
     &                       bed_mass(Iend,j,k,nnew,ised)
            bed_frac(Iend+1,j,k,ised)=
     &                        bed_frac(Iend,j,k,ised)
             enddo
             bed_thick(Iend+1,j,k)=bed_thick(Iend,j,k)
             bed_poros(Iend+1,j,k)=bed_poros(Iend,j,k)
             bed_age(Iend+1,j,k)=bed_age(Iend,j,k)
           enddo
           bot_thick(Iend+1,j)=bot_thick(Iend,j)
         enddo
        endif
      if (WESTERN_EDGE) then
        do j=Jstr,Jend
          do k=1,NLAY
            do ised=1,NST
            bed_mass(Istr-1,j,k,nnew,ised)=
     &                       bed_mass(Istr,j,k,nnew,ised)
            bed_frac(Istr-1,j,k,ised)=
     &                        bed_frac(Istr,j,k,ised)
             enddo
             bed_thick(Istr-1,j,k)=bed_thick(Istr,j,k)
             bed_poros(Istr-1,j,k)=bed_poros(Istr,j,k)
             bed_age(Istr-1,j,k)=bed_age(Istr,j,k)
           enddo
           bot_thick(Istr-1,j)=bot_thick(Istr,j)
        enddo
      endif
#  endif
#  ifndef NS_PERIODIC
      if (NORTHERN_EDGE) then
        do i=Istr,Iend
          do k=1,NLAY
            do ised=1,NST
              bed_mass (i,Jend+1,k,nnew,ised)=
     &                            bed_mass(i,Jend,k,nnew,ised)
              bed_frac (i,Jend+1,k,ised)=
     &                            bed_frac(i,Jend,k,ised)
            enddo
            bed_thick(i,Jend+1,k)=bed_thick(i,Jend,k)
            bed_poros(i,Jend+1,k)=bed_poros(i,Jend,k)
            bed_age  (i,Jend+1,k)=bed_age(i,Jend,k)
          enddo
          bot_thick  (i,Jend+1)=bot_thick(i,Jend)
        enddo
      endif
      if (SOUTHERN_EDGE) then
        do i=Istr,Iend
          do k=1,NLAY
            do ised=1,NST
              bed_mass (i,Jstr-1,k,nnew,ised)=
     &                            bed_mass(i,Jstr,k,nnew,ised)
              bed_frac (i,Jstr-1,k,ised)=
     &                            bed_frac(i,Jstr,k,ised)
            enddo
            bed_thick(i,Jstr-1,k)=bed_thick(i,Jstr,k)
            bed_poros(i,Jstr-1,k)=bed_poros(i,Jstr,k)
            bed_age  (i,Jstr-1,k)=bed_age(i,Jend,k)
          enddo
          bot_thick  (i,Jstr-1)=bot_thick(i,Jstr)
        enddo
      endif
#  endif
#  if !defined EW_PERIODIC && !defined NS_PERIODIC
      if (SOUTHERN_EDGE.and.WESTERN_EDGE) then
        bed_mass(Istr-1,Jstr-1,:,nnew,:)=
     &  0.5*(bed_mass(Istr,Jstr-1,:,nnew,:)
     &      +bed_mass(Istr-1,Jstr,:,nnew,:))
        bed_frac(Istr-1,Jstr-1,:,:)=
     &  0.5*(bed_frac(Istr,Jstr-1,:,:)
     &      +bed_frac(Istr-1,Jstr,:,:))
        bed_thick(Istr-1,Jstr-1,:)=
     &  0.5*(bed_thick(Istr,Jstr-1,:)
     &      +bed_thick(Istr-1,Jstr,:))
        bed_poros(Istr-1,Jstr-1,:)=
     &  0.5*(bed_poros(Istr,Jstr-1,:)
     &      +bed_poros(Istr-1,Jstr,:))
        bed_age(Istr-1,Jstr-1,:)=
     &  0.5*(bed_age(Istr,Jstr-1,:)
     &      +bed_age(Istr-1,Jstr,:))
      endif
      if (SOUTHERN_EDGE.and.EASTERN_EDGE) then
        bed_mass(Iend+1,Jstr-1,:,nnew,:)=
     &  0.5*(bed_mass(Iend,Jstr-1,:,nnew,:)
     &      +bed_mass(Iend+1,Jstr,:,nnew,:))
        bed_frac(Iend+1,Jstr-1,:,:)=
     &  0.5*(bed_frac(Iend,Jstr-1,:,:)
     &      +bed_frac(Iend+1,Jstr,:,:))
        bed_thick(Iend+1,Jstr-1,:)=
     &  0.5*(bed_thick(Iend,Jstr-1,:)
     &      +bed_thick(Iend+1,Jstr,:))
        bed_poros(Iend+1,Jstr-1,:)=
     &  0.5*(bed_poros(Iend,Jstr-1,:)
     &      +bed_poros(Iend+1,Jstr,:))
        bed_age(Iend+1,Jstr-1,:)=
     &  0.5*(bed_age(Iend,Jstr-1,:)
     &      +bed_age(Iend+1,Jstr,:))
      endif
      if (NORTHERN_EDGE.and.WESTERN_EDGE) then
        bed_mass(Istr-1,Jend+1,:,nnew,:)=
     &  0.5*(bed_mass(Istr,Jend+1,:,nnew,:)
     &      +bed_mass(Istr-1,Jend,:,nnew,:))
        bed_frac(Istr-1,Jend+1,:,:)=
     &  0.5*(bed_frac(Istr,Jend+1,:,:)
     &      +bed_frac(Istr-1,Jend,:,:))
        bed_thick(Istr-1,Jend+1,:)=
     &  0.5*(bed_thick(Istr,Jend+1,:)
     &      +bed_thick(Istr-1,Jend,:))
        bed_poros(Istr-1,Jend+1,:)=
     &  0.5*(bed_poros(Istr,Jend+1,:)
     &      +bed_poros(Istr-1,Jend,:))
        bed_age(Istr-1,Jend+1,:)=
     &  0.5*(bed_age(Istr,Jend+1,:)
     &      +bed_age(Istr-1,Jend,:))
      endif
      if (NORTHERN_EDGE.and.EASTERN_EDGE) then
        bed_mass(Iend+1,Jend+1,:,nnew,:)=
     &   0.5*(bed_mass(Iend,Jend+1,:,nnew,:)
     &       +bed_mass(Iend+1,Jend,:,nnew,:))
        bed_frac(Iend+1,Jend+1,:,:)=
     &   0.5*(bed_frac(Iend,Jend+1,:,:)
     &       +bed_frac(Iend+1,Jend,:,:))
        bed_thick(Iend+1,Jend+1,:)=
     &   0.5*(bed_thick(Iend,Jend+1,:)
     &       +bed_thick(Iend+1,Jend,:))
        bed_poros(Iend+1,Jend+1,:)=
     &   0.5*(bed_poros(Iend,Jend+1,:)
     &       +bed_poros(Iend+1,Jend,:))
        bed_age(Iend+1,Jend+1,:)=
     &   0.5*(bed_age(Iend,Jend+1,:)
     &       +bed_age(Iend+1,Jend,:))
      endif
#  endif

#  if defined BD_EXCHANGE
      do ised=1,NST
        do k=1,NLAY
          call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_frac(START_2D_ARRAY,k,ised))
          call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_mass(START_2D_ARRAY,k,nnew,ised))
        enddo
      enddo
      do k=1,NLAY
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_poros(START_2D_ARRAY,k))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_thick(START_2D_ARRAY,k))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_age(START_2D_ARRAY,k))
      enddo
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bot_thick(START_2D_ARRAY))
#  endif
!
!-----------------------------------------------------------------------
!  Store total bed thickness
!-----------------------------------------------------------------------
!
#  ifdef MORPHODYN
      do j=Jstr,Jend
        do i=Istr,Iend
          bed_thick_tot(i,j,nnew)=0.
          do k=1,NLAY
            bed_thick_tot(i,j,nnew)=bed_thick_tot(i,j,nnew)+
     &                              bed_thick(i,j,k)
          enddo
        enddo
      enddo
!
! Set boundary conditions
!
#   ifndef EW_PERIODIC
      if (EASTERN_EDGE) then
        do j=Jstr,Jend
          bed_thick_tot(Iend+1,j,nnew)=bed_thick_tot(Iend,j,nnew)
        enddo
      endif
      if (WESTERN_EDGE) then
        do j=Jstr,Jend
          bed_thick_tot(Istr-1,j,nnew)=bed_thick_tot(Istr,j,nnew)
        enddo
      endif
#   endif
#   ifndef NS_PERIODIC
      if (NORTHERN_EDGE) then
        do i=Istr,Iend
          bed_thick_tot(i,Jend+1,nnew)=bed_thick_tot(i,Jend,nnew)
        enddo
      endif
      if (SOUTHERN_EDGE) then
        do i=Istr,Iend
          bed_thick_tot(i,Jstr-1,nnew)=bed_thick_tot(i,Jstr,nnew)
        enddo
      endif
#   endif
#   if !defined EW_PERIODIC && !defined NS_PERIODIC
      if (SOUTHERN_EDGE.and.WESTERN_EDGE) then
        bed_thick_tot(Istr-1,Jstr-1,nnew)=
     &  0.5*(bed_thick_tot(Istr,Jstr-1,nnew)
     &      +bed_thick_tot(Istr-1,Jstr,nnew))
      endif
      if (SOUTHERN_EDGE.and.EASTERN_EDGE) then
        bed_thick_tot(Iend+1,Jstr-1,nnew)=
     &  0.5*(bed_thick_tot(Iend,Jstr-1,nnew)
     &      +bed_thick_tot(Iend+1,Jstr,nnew))
      endif
      if (NORTHERN_EDGE.and.WESTERN_EDGE) then
        bed_thick_tot(Istr-1,Jend+1,nnew)=
     &  0.5*(bed_thick_tot(Istr,Jend+1,nnew)
     &      +bed_thick_tot(Istr-1,Jend,nnew))
      endif
      if (NORTHERN_EDGE.and.EASTERN_EDGE) then
        bed_thick_tot(Iend+1,Jend+1,nnew)=
     &   0.5*(bed_thick_tot(Iend,Jend+1,nnew)
     &       +bed_thick_tot(Iend+1,Jend,nnew))
      endif
#   endif
!
#   if defined BD_EXCHANGE
       call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                         bed_thick_tot(START_2D_ARRAY,nnew))
#   endif
!
!-----------------------------------------------------------------------
!  Store bed evolution
!-----------------------------------------------------------------------
!
      do j=Jstr-1,Jend+1
        do i=Istr-1,Iend+1
          dh(i,j)=bed_thick_tot(i,j,nstp)-bed_thick_tot(i,j,nnew)
        enddo
      enddo
!
#  endif /* MORPHODYN */

      return
      end

# else
**********************************************************************
      subroutine sed_bed_cohesive_tile (
     &                         Istr,Iend,Jstr,Jend,tau_w,dep_mass)
!***********************************************************************
!
      implicit none
#  include "param.h"
#  include "bbl.h"
#  include "forces.h"
#  include "grid.h"
#  include "ocean3d.h"
#  include "scalars.h"
#  include "sediment.h"
      integer Istr,Iend,Jstr,Jend
      real tau_w(PRIVATE_2D_SCRATCH_ARRAY)
      real dep_mass(PRIVATE_1D_SCRATCH_ARRAY,NST)
      integer Ksed, i, ised, j, k, ks, nnn, indx
      real cff, cff1, cff2, cff3,cff4,eps
      real thck_avail, thck_to_add
      real nlysm(NST)   !!!!???(OpenMP)
      integer imin,imax,jmin,jmax
      real add,newlayer_thick

!     parameter (eps = 1.0E-14,newlayer_thick=0.0005) !BEDTOY
      parameter (eps = 1.0E-14,newlayer_thick=0.001)

#  if defined MIXED_BED
      real pcoh
#  endif
#  if defined COHESIVE_BED || defined MIXED_BED
      real alpha, bzactv, frt, tcb_top, tcb_bot
      real tcr(NLAY)  !!!!???(OpenMP)
#  endif
#  if defined COHESIVE_BED || defined MIXED_BED
      real bmz(NLAY)
#  endif
#  if defined SED_FLOCS && defined SED_DEFLOC
      real masseq(NMUD)
#  endif

! KLUDGE alert
! minlayer_thick = 0.0005
      minlayer_thick = newlayer_thick
!
!-----------------------------------------------------------------------
! Compute sediment bed layer stratigraphy.
!-----------------------------------------------------------------------
!
! why MPM ???
#  if defined BEDLOAD_MPM || defined SUSPLOAD
      do j=Jstr-1,Jend+1
        do i=Istr-1,Iend+1
          tau_w(i,j)=SQRT(bustrw(i,j)**2
     &                   +bvstrw(i,j)**2)
#   ifdef WET_DRY
          tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
#   endif
        enddo
      enddo
#  endif
!
!-----------------------------------------------------------------------
!  Update bed properties according to ero_flux and dep_flux.
!-----------------------------------------------------------------------
!
#  ifdef SUSPLOAD
      J_LOOP : do j=Jstr,Jend
!
!  The deposition and resuspension of sediment on the bottom "bed"
!  is due to precipitation flux FC(:,0), already computed, and the
!  resuspension (erosion, hence called ero_flux). The resuspension is
!  applied to the bottom-most grid box value qc(:,1) so the total mass
!  is conserved. Restrict "ero_flux" so that "bed" cannot go negative
!  after both fluxes are applied.
!
      do i=Istr,Iend
         SED_LOOP: do ised=1,NST
            dep_mass(i,ised)=0.0
#   ifdef SED_MORPH
            ! Apply morphology factor.
             ero_flux(i,j,ised)=ero_flux(i,j,ised)   MORPH_FAC
             settling_flux(i,j,ised)=settling_flux(i,j,ised)  MORPH_FAC
#   endif

            !  Update bed mass arrays.
            bed_mass(i,j,1,nnew,ised)=MAX(bed_mass(i,j,1,bnew,ised)-
     &                                   (ero_flux(i,j,ised)-
     &                                    settling_flux(i,j,ised)), 0.)
            do k=2,NLAY
               bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k,nstp,ised)
            enddo
         enddo SED_LOOP

         cff3=0.0
         do ised=1,NST
            cff3=cff3+bed_mass(i,j,1,nnew,ised)
         enddo

         if (cff3.eq.0.0) then
            cff3=eps
         endif
         bed_thick(i,j,1)=0.0
         do ised=1,NST
            bed_frac(i,j,1,ised)=bed_mass(i,j,1,nnew,ised)/cff3
            bed_thick(i,j,1)=MAX(bed_thick(i,j,1)+
     &        bed_mass(i,j,1,nnew,ised)/
     &        (Srho(ised)*
     &        (1.0-bed_poros(i,j,1))),0.0)
         enddo
      enddo

      enddo J_LOOP
#  endif /* SUSPLOAD section */
!
!-----------------------------------------------------------------------
!  At this point, all deposition or erosion is complete, and
!  has been added/subtracted to top layer. Thickness has NOT been corrected.
!-----------------------------------------------------------------------
!
!!!!! taucb ????
      J_LOOP_CB : do j=Jstr,Jend
        I_LOOP_CB: do i=Istr,Iend
!       Calculate active layer thickness, bot_thick(i,j).
!       (trunk version allows this to be zero...this has minimum of 6*D50)
          bot_thick(i,j)=MAX(0.0,
     &                          0.007*
     &                          (tau_w(i,j)-taucb(i,j))*rho0)+
     &                          6.0*Ssize(i,j)
#  if defined COHESIVE_BED
          bot_thick(i,j)=MAX(0.0,
     &                          0.007*
     &                          (tau_w(i,j)*rho0-bed_bctr(i,j,1)))+
     &                          6.0*Ssize(i,j)
#  endif
#  if defined MIXED_BED
          cff1=             MAX(0.0,
     &                          0.007*
     &                          (tau_w(i,j)*rho0-bed_bctr(i,j,1)))+
     &                          6.0*Ssize(i,j)
          cff2=             MAX(0.0,
     &                          0.007*
     &                          (tau_w(i,j)-taucb(i,j))*rho0)+
     &                          6.0*Ssize(i,j)
          bot_thick(i,j)=MAX(cff1,cff2)
#  endif
#  ifdef MORPHODYN
           ! Apply morphology factor.
           bot_thick(i,j)=MAX( bot_thick(i,j)*morph_fac,
     &          bot_thick(i,j))
#  endif
#  if defined COHESIVE_BED || defined MIXED_BED
!  Find first layer with tc > tb
!  Remember, the critical stresses apply to the TOP of each layer
          Ksed = 1
          bzactv = 0.0
          frt = 0.0 ! CRS
          cff1 = rho0*tau_w(i,j)
          tcb_top = bed_bctr(i,j,1)
          tcb_bot = bed_bctr(i,j,2)
#   if defined MIXED_BED
!         Calc. tau crit for bottom of next layer
!         Update cohesive property of seds in top layer
          cff3 = 0.0
          cff4 = 1.0
          do ised=imud1,imud2
             cff3=cff3+bed_frac(i,j,1,ised)
             cff4=cff4*tau_ce(ised)**bed_frac(i,j,1,ised)
          enddo
          do ised=isand1,isand2
             cff4=cff4*tau_ce(ised)**bed_frac(i,j,1,ised)
          enddo
          pcoh=min(max((cff3-transN)/(transC-transN)
     &      ,0.0),1.0)
          tcb_top = pcoh*bed_bctr(i,j,1)+(pcoh-1.0)*cff4*rho0   !BF
#   endif
          if(cff1 .GT. tcb_top)then
!            Calculate tcb_temp for next layer
             tcb_bot = bed_bctr(i,j,Ksed+1)
#   if defined MIXED_BED
!            Recalculate cohesive fraction and mean grain tau crit
!            Note that Ksed is used for grain props, and Ksed+1 for
!            bed tau crit at bottom of layer Ksed
             cff3 = 0.0
             cff4 = 1.0
             do ised=imud1,imud2
                cff3=cff3+bed_frac(i,j,Ksed,ised)
                cff4=cff4*tau_ce(ised)**bed_frac(i,j,Ksed,ised)
             enddo
             do ised=isand1,isand2
                cff4=cff4*tau_ce(ised)**bed_frac(i,j,Ksed,ised)
             enddo
!            Calculate cohesive behavior and blended tau crit
             pcoh=min(max((cff3-transN)/(transC-transN),
     &            0.0),1.0)
             tcb_bot =pcoh*bed_bctr(i,j,Ksed+1)+(pcoh-1.0)*cff4*rho0  !BF
#   endif
             do while ( (Ksed.LE.(NLAY-1)) .AND.
     &            (cff1 .GT. tcb_bot))
!ALA   Instead of adding entire 2nd layer, add just what is needed
!ALA!               Add entire layer
!ALA                bzactv = bzactv + bed_thick(i,j,Ksed)
!       Add thickness equivalent to difference between cff1 and tcb_bot
                if (cff1.GT.bed_bctr(i,j,Ksed+1)) then
                   bzactv = bzactv + bed_thick(i,j,Ksed)
                   tcb_top = tcb_bot
                   tcb_bot = bed_bctr(i,j,Ksed+1)
                else
                   bzactv = bzactv + MIN(bed_thick(i,j,Ksed),
     &                      (MAX(0.0,0.007*(cff1-tcb_bot))))
                   tcb_top = cff1
                   tcb_bot = bed_bctr(i,j,Ksed+1)
                endif
#   if defined MIXED_BED
!               Recalculate cohesive fraction and mean grain tau crit
                cff3 = 0.0
                cff4 = 1.0
                do ised=imud1,imud2
                   cff3=cff3+bed_frac(i,j,Ksed,ised)
                   cff4=cff4*tau_ce(ised)**bed_frac(i,j,Ksed,ised)
                enddo
                do ised=isand1,isand2
                   cff4=cff4*tau_ce(ised)**bed_frac(i,j,Ksed,ised)
                enddo
!               Calculate cohesive behavior and blended tau crit
                pcoh=min(max((cff3-transN)/(transC-transN),
     &               0.0),1.0)
                tcb_bot = pcoh*bed_bctr(i,j,Ksed+1)+(pcoh-1.0)*cff4
     &                         *rho0                                 !BF
#   endif
                Ksed = Ksed+1
             enddo
             frt = MAX(0.0,(cff1-tcb_top) ) /
     &            MAX( eps,
     &            ( tcb_bot-tcb_top ))
          endif
          bzactv = bzactv+frt*bed_thick(i,j,Ksed)
          bzactv = MAX( bzactv, 6.0*Ssize(i,j) ) !CRS
          bot_thick(i,j)=min(bot_thick(i,j),bzactv) !?CRS
!          bot_thick(i,j)=max(bot_thick(i,j),bzactv) !?CRS
#  endif /* defined COHESIVE_BED || defined MIXED_BED */
        enddo I_LOOP_CB
      enddo J_LOOP_CB

        J_LOOP2 : do j=Jstr,Jend
        do i=Istr,Iend

!
!          Calculate net deposition and erosion
           cff=0.0
           cff2=0.0
           do ised=1,NST
              cff=cff+settling_flux(i,j,ised)
              cff2=cff2+ero_flux(i,j,ised)
              dep_mass(i,ised)=0.0
              if ((ero_flux(i,j,ised)-settling_flux(i,j,ised)).lt.
     &             0.0) then
                 dep_mass(i,ised)=settling_flux(i,j,ised)-
     &                ero_flux(i,j,ised)
              endif
           enddo
           bot_dnet(i,j)=cff-cff2

           if ( cff-cff2.GT.0.0) then ! NET deposition
              !  Deposition. Determine if we need to create a new bed layer

#  if defined COHESIVE_BED || defined MIXED_BED
!
!  Calculate tau_crit of deposited bed
!
!  Calculate new mass in top layer
             bmz(1) = 0.0
             do ised=1,NST
                bmz(1) = bmz(1)+bed_mass(i,j,1,nnew,ised)
             enddo
             if (NLAY.GT.1) then
!  Average of cff1 and cff2, where
!    cff1 = linear extension of previous tcr slope to new surface
!    cff2 = minimum deposition
                cff1 = bed_bctr(i,j, 1) -
     &             bot_dnet(i,j) *
     &             (bed_bctr(i,j,2)-bed_bctr(i,j,1)) /
     &             (bmz(1)-bot_dnet(i,j))
                cff2 = MAX( rho0*tau_w(i,j) , tcr_min )
                bed_bctr(i,j,1) = MIN(bed_bctr(i,j,1),
     &                                 MAX( 0.5*(cff1+cff2), cff2 ))
!                cff1 = bot_dnet(i,j)/MAX(bmz(1),eps)
!                cff2 = bed_bctr(i,j,1)+bed_bctr(i,j,2)-2*tcr_min
!                bed_bctr(i,j,1) =bed_bctr(i,j,1) - cff1*cff2
             else
!  TODO: Not sure what mud tau_crit should be for single-layer bed.
!  Try weighted average of dep and old value
                cff1 = (bed_bctr(i,j,1)*(bmz(1)-bot_dnet(i,j))
     &                  + cff2*bot_dnet(i,j)) / bmz(1)
                bed_bctr(i,j,1) = MAX( cff1, cff2 )
             endif
#  endif
              ! (no test for age here)
              bed_age(i,j,1)=time
              if(bed_thick(i,j,1).gt.
     &             MAX(bot_thick(i,j),newlayer_thick)) then
                 ! Top layer is too thick
                 if (NLAY.gt.2) then
                    if(bed_thick(i,j,2).lt.minlayer_thick) then
                    ! Layer 2 is smaller than minimum size
                    ! Instead of pushing down all layers, just combine top 2 layers
                       cff=0.0
                       cff1=0.0
                       cff2=0.0
                       do ised=1,NST
                          cff =cff +dep_mass(i,ised)
                          cff1=cff1+bed_mass(i,j,1,nnew,ised)
                          cff2=cff2+bed_mass(i,j,2,nnew,ised)
                       enddo
#  if defined COHESIVE_BED || defined MIXED_BED
!
!   Assign new tau_crit at 2nd layer
!
                    bed_bctr(i,j,2)= bed_bctr(i,j,2) -
     &       (cff1-cff)/(cff1+cff2)*(bed_bctr(i,j,3)-bed_bctr(i,j,1))
!                  taucb(i,j,2) = 0.5*( taucb(i,j,1)+
!     &                                        taucb(i,j,2) )
#  endif
!

                       !  Update bed mass
                       do ised=1,NST
                          bed_mass(i,j,2,nnew,ised)=
     &                         MAX(bed_mass(i,j,2,nnew,ised)+
     &                         bed_mass(i,j,1,nnew,ised)-
     &                         dep_mass(i,ised),0.0)
                          bed_mass(i,j,1,nnew,ised)=dep_mass(i,ised)
                       enddo
                       ! ALA - average time and porosity
                       ! ALA CHECK WITH CRS cff1 or cff1-cff for first layer
                       bed_age(i,j,2)=(bed_age(i,j,1)*cff1+
     &                         bed_age(i,j,2)*cff2)/(cff1+cff2)
                       bed_age(i,j,1)=time
                       bed_poros(i,j,2)=(bed_poros(i,j,1)*cff1+
     &                         bed_poros(i,j,2)*cff2)/(cff1+cff2)
                       ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER
                       bed_poros(i,j,1)=bed_poros(i,j,1)
                    else
                    ! Layer 2 is > minlayer thick, need another layer
                    !  Combine bottom layers.
                       cff1=0.0
                       cff2=0.0
                       do ised=1,NST
                          cff1=cff1+bed_mass(i,j,NLAY-1,nnew,ised)
                          cff2=cff2+bed_mass(i,j,NLAY,nnew,ised)
                       enddo
                       bed_poros(i,j,NLAY)=
     &                      (bed_poros(i,j,NLAY-1)*cff1+
     &                      bed_poros(i,j,NLAY)*cff2)/(cff1+cff2)
                       bed_age(i,j,NLAY)=
     &                      (bed_age(i,j,NLAY-1)*cff1+
     &                      bed_age(i,j,NLAY)*cff2)/(cff1+cff2)
#  if defined COHESIVE_BED || defined MIXED_BED
!
!   Assign tcrit at top of new bottom bed tcrit for NLAY-1
              bed_bctr(i,j,NLAY)= bed_bctr(i,j,NLAY-1)
!              taucb(i,j,NLAY)=taucb(i,j,NLAY-1) +                &
!     &            cff2*(taucb(i,j,NLAY)-taucb(i,j,NLAY-1))/cff1
#  endif

                       do ised=1,NST
                          bed_mass(i,j,NLAY,nnew,ised)=
     &                         bed_mass(i,j,NLAY-1,nnew,ised)+
     &                         bed_mass(i,j,NLAY  ,nnew,ised)
                       enddo
                       !
                       !  Push layers down.
                       do k=NLAY-1,2,-1
                          bed_poros(i,j,k)=bed_poros(i,j,k-1)
                          bed_age(i,j,k)=bed_age(i,j,k-1)
                          do ised =1,NST
                             bed_mass(i,j,k,nnew,ised)=
     &                            bed_mass(i,j,k-1,nnew,ised)
                          enddo
#  if defined COHESIVE_BED || defined MIXED_BED
                bed_bctr(i,j,k)=bed_bctr(i,j,k-1)
#  endif
                       enddo
                       !  Set new top parameters for top 2 layers
#  if defined COHESIVE_BED || defined MIXED_BED
                       !  Tau_crit at top already determined
                       !  Set tau_crit at 2nd layer
                       cff=0.0
                       cff1=0.0
                       do ised=1,NST
                          cff =cff +dep_mass(i,ised)
                          cff1=cff1+bed_mass(i,j,1,nnew,ised)
                       enddo

                       cff2 = (bed_bctr(i,j,2)-bed_bctr(i,j,1))/cff1
                       bed_bctr(i,j,2) = bed_bctr(i,j,1)+cff*cff2
#  endif

                       do ised=1,NST
                          bed_mass(i,j,2,nnew,ised)=
     &                         MAX(bed_mass(i,j,2,nnew,ised)-
     &                         dep_mass(i,ised),0.0)
                          bed_mass(i,j,1,nnew,ised)=dep_mass(i,ised)
                       enddo
                    endif
                 elseIF (NLAY.eq.2) then
                  ! NBED=2
                  !
#  if defined COHESIVE_BED || defined MIXED_BED
                    !  Tau_crit at top already determined
                    !  Set tau_crit at 2nd layer
                    cff=0.0
                    cff1=0.0
                    do ised=1,NST
                       cff =cff +dep_mass(i,ised)
                       cff1=cff1+bed_mass(i,j,1,nnew,ised)
                    enddo
                    cff2 = (bed_bctr(i,j,2)-bed_bctr(i,j,1))/cff1
                    bed_bctr(i,j,2) = bed_bctr(i,j,1)+cff*cff2
#  endif
                    cff1=0.0
                    cff2=0.0
                    do ised=1,NST
                       cff1=cff1+bed_mass(i,j,1,nnew,ised)
                       cff2=cff2+bed_mass(i,j,2,nnew,ised)
                    enddo
                    do ised=1,NST
                       bed_mass(i,j,2,nnew,ised)=
     &                      MAX(bed_mass(i,j,2,nnew,ised)+
     &                      bed_mass(i,j,1,nnew,ised)-
     &                      dep_mass(i,ised),0.0)
                       bed_mass(i,j,1,nnew,ised)=dep_mass(i,ised)
                    enddo
                    ! ALA - average time and porosity
                    bed_age(i,j,2)=(bed_age(i,j,1)*cff1+
     &                      bed_age(i,j,2)*cff2)/(cff1+cff2)
                    bed_age(i,j,1)=time
                    bed_poros(i,j,2)=(bed_poros(i,j,1)*cff1+
     &                      bed_poros(i,j,2)*cff2)/(cff1+cff2)
                    ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER
                    bed_poros(i,j,1)=bed_poros(i,j,1)
                 else
                 ! NBED=1
                 endif
              else
                ! Net deposition has occurred, but no new bed layer was created
             endif
           else
                ! Net erosion occurred
#  if defined COHESIVE_BED || defined MIXED_BED
              bmz(1) = 0.0
              do ised=1,NST
                 bmz(1) = bmz(1)+bed_mass(i,j,1,nnew,ised)
              enddo
!  recalc tc for top of bed, based on linear
!  interpolation and mass removed / orig. mass in top layer
              bed_bctr(i,j,1)=bed_bctr(i,j,1)+
     &             MIN(1.0,-bot_dnet(i,j)/MAX(eps,bmz(1)))*
     &             MAX(0.0,(bed_bctr(i,j,2)-bed_bctr(i,j,1)))
#  endif
                bed_age(i,j,1)=time
                if (NLAY.eq.2) then
                 ! NBED=2
                    do ised=1,NST
                       bed_mass(i,j,2,nnew,ised)=
     &                      MAX(bed_mass(i,j,2,nnew,ised)+
     &                      bed_mass(i,j,1,nnew,ised)-
     &                      dep_mass(i,ised),0.0)
                       bed_mass(i,j,1,nnew,ised)=dep_mass(i,ised)
                    enddo
                elseif (NLAY.eq.1) then
                ! ALF NO NEED TO do ANYTHING
                else
                endif
           endif

          ! Recalculate thickness and fractions for all layers.
           do k=1,NLAY
              cff3=0.0
              do ised=1,NST
                 cff3=cff3+bed_mass(i,j,k,nnew,ised)
              enddo
              if (cff3.eq.0.0) then
                 cff3=eps
              endif
              bed_thick(i,j,k)=0.0
              do ised=1,NST
                 bed_frac(i,j,k,ised)=bed_mass(i,j,k,nnew,ised)/cff3
                 bed_thick(i,j,k)=MAX(bed_thick(i,j,k)+
     &                         bed_mass(i,j,k,nnew,ised)/
     &                         (Srho(ised)*
     &                (1.0-bed_poros(i,j,k))),0.0)
              enddo
           enddo
        enddo
      enddo J_LOOP2

      J_LOOP3 : do j=Jstr,Jend
        I_LOOP3 : do i=Istr,Iend
          if (bot_thick(i,j).gt.bed_thick(i,j,1)) then
            if (NLAY.eq.1) then
              bot_thick(i,j)=bed_thick(i,j,1)
            else
              thck_to_add=bot_thick(i,j)-bed_thick(i,j,1)
              thck_avail=0.0
              Ksed=1                                        ! initialize
              do k=2,NLAY
                if (thck_avail.lt.thck_to_add) then
                  thck_avail=thck_avail+bed_thick(i,j,k)
                  Ksed=k
                endif
              enddo
!
!  Catch here if there was not enough bed material.
!
              if (thck_avail.lt.thck_to_add) then
                bot_thick(i,j)=bed_thick(i,j,1)+thck_avail
                thck_to_add=thck_avail
              endif
!
!  Update bed mass of top layer and fractional layer.
!
              cff2=MAX(thck_avail-thck_to_add,0.0)/
     &             MAX(bed_thick(i,j,Ksed),eps)
              do ised=1,NST
                cff1=0.0
                do k=1,Ksed
                  cff1=cff1+bed_mass(i,j,k,nnew,ised)
                enddo
                cff3=cff2*bed_mass(i,j,Ksed,nnew,ised)
                bed_mass(i,j,1   ,nnew,ised)=cff1-cff3
                bed_mass(i,j,Ksed,nnew,ised)=cff3
              enddo
!
!  Update thickness of fractional layer ksource_sed.
!
              bed_thick(i,j,Ksed)=MAX(thck_avail-thck_to_add,0.0)

#  if defined COHESIVE_BED || defined MIXED_BED
!
! Update tau_cr of fractional layer   !RB debug
!           IF (Ksed < NLAY ) then
            bed_bctr(i,j,Ksed) = bed_bctr(i,j,Ksed+1)-
     &                 cff2*(bed_bctr(i,j,Ksed+1)-bed_bctr(i,j,Ksed))
!           endif
#  endif

!
!  Update bed fraction of top layer.
!
              cff3=0.0
              do ised=1,NST
                cff3=cff3+bed_mass(i,j,1,nnew,ised)
              enddo
              if (cff3.eq.0.0) then
                 cff3=eps
              endif
              do ised=1,NST
                bed_frac(i,j,1,ised)=bed_mass(i,j,1,nnew,ised)/cff3
              enddo
!
!  Upate bed thickness of top layer.
!
              bed_thick(i,j,1)=bot_thick(i,j)
!
!  Pull all layers closer to the surface.
!
              do k=Ksed,NLAY
                ks=Ksed-2
                bed_thick(i,j,k-ks)=bed_thick(i,j,k)
                bed_poros(i,j,k-ks)=bed_poros(i,j,k)
                bed_age(i,j,k-ks)=bed_age(i,j,k)

#  if defined COHESIVE_BED || defined MIXED_BED
                bed_bctr(i,j,k-ks)=bed_bctr(i,j,k)
#  endif

                do ised=1,NST
                  bed_frac(i,j,k-ks,ised)=bed_frac(i,j,k,ised)
                  bed_mass(i,j,k-ks,nnew,ised)=bed_mass(i,j,k,nnew,ised)
                enddo
              enddo
!
!  Add new layers onto the bottom. Split what was in the bottom layer to
!  fill these new empty cells. ("ks" is the number of new layers).
!
              ks=Ksed-2
              cff=1.0/REAL(ks+1)
#  if defined COHESIVE_BED || defined MIXED_BED
#   undef BF_TCR
#   if defined BF_TCR
          bmz(1) = 0.0
          do ised=1,NST
             bmz(1) = bmz(1)+bed_mass(i,j,1,nnew,ised)
          enddo
          do k=2,NLAY
             bmz(k) = bmz(k-1)
             do ised=1,NST
                bmz(k)=bmz(k)+bed_mass(i,j,k,nnew,ised)
             enddo
          enddo
          tcr(1) = tcr_min
          do k=2,NLAY
             tcr(k) = tcr_min
             if (bmz(k-1).GT.eps) then
!                tcr(k) = exp((log(bmz(k-1))-
!     &               bottom(i,j,idoff))/
!     &               bottom(i,j,idslp))
                tcr(k) = exp((log(bmz(k-1))-
     &                    tcr_off)/
     &                    tcr_slp)
             endif
             tcr(k) = MIN( MAX( tcr(k), tcr_min), tcr_max )
          enddo
#   else
!alf              cff1 = cff*(tcr_max-taucb(i,j,NLAY-ks))
!alf          Use the value from the bottom layer as max.
              cff1 = cff*(bed_bctr(i,j,NLAY)-bed_bctr(i,j,NLAY-ks))
#   endif

! Interpolate bottom layer to tau_crit_max
! TODO: should be the reference profile and not tau_crit_max
                  do k=NLAY-1,NLAY-ks,-1
#   if defined BF_TCR
                     bed_bctr(i,j,k)=tcr(k)
#   else
                     bed_bctr(i,j,k)=bed_bctr(i,j,NLAY-ks)+
     &                       REAL(k-NLAY+ks) * cff1
#   endif
                  enddo
#  endif
              ! ALA CHECK WITH CRS about bed_frac
              nnn=0
              do ised=1,NST
                 nlysm(ised)=newlayer_thick*REAL(ks+1)*
     &                  (Srho(ised)*
     &                  (1.0-bed_poros(i,j,NLAY-ks)))*
     &                  bed_frac(i,j,NLAY-ks,ised)
                  IF (ks.gt.0) then
                        if (bed_mass(i,j,NLAY-ks,nnew,ised).gt.
     &                       nlysm(ised)) then
                           nnn=nnn+1
                           nlysm(ised)=
     &                           newlayer_thick*REAL(ks)*
     &                          (Srho(ised)*
     &                          (1.0-bed_poros(i,j,NLAY-ks)))*
     &                          bed_frac(i,j,NLAY-ks,ised)
                        endif
                    endif
              enddo
              if (nnn.eq.NST) then
                  bed_thick(i,j,NLAY)=bed_thick(i,j,NLAY-ks)-
     &                                newlayer_thick*REAL(ks)
                  do ised=1,NST
                     bed_mass(i,j,NLAY,nnew,ised)=
     &                  bed_mass(i,j,NLAY-ks,nnew,ised)-nlysm(ised)
                  enddo
                  do k=NLAY-1,NLAY-ks,-1
                     bed_thick(i,j,k)=newlayer_thick
                     bed_age(i,j,k)=bed_age(i,j,NLAY-ks)
                     do ised=1,NST
                        bed_frac(i,j,k,ised)=bed_frac(i,j,NLAY-ks,ised)
                        bed_mass(i,j,k,nnew,ised)=
     &                           nlysm(ised)/REAL(ks)
                     enddo
                  enddo
              else
                 cff=1.0/REAL(ks+1)
                  do k=NLAY,NLAY-ks,-1
                     bed_thick(i,j,k)=bed_thick(i,j,NLAY-ks)*cff
                     bed_age(i,j,k)=bed_age(i,j,NLAY-ks)
                     do ised=1,NST
                        bed_frac(i,j,k,ised)=bed_frac(i,j,NLAY-ks,ised)
                        bed_mass(i,j,k,nnew,ised)=
     &                             bed_mass(i,j,NLAY-ks,nnew,ised)*cff
                     enddo
                  enddo
              endif
            endif  ! NLAY > 1
          endif  ! increase top bed layer
#  if defined MIXED_BED
!
!         Update cohesive property of seds in top layer
          cff1 = 0.0
          do ised=imud1,imud2
             cff1=cff1+bed_frac(i,j,1,ised)
          enddo
          bot_dprp(i,j)=min(max((cff1-transN)/
     &     (transC-transN),0.0),1.0)
#  endif
          enddo I_LOOP3

#  if defined SED_FLOCS && defined SED_DEFLOC
          I_LOOP_DFL: do i=Istr,Iend
!        Activate defloc in all layers but the top one
!            do k=2,NLAY
!        Activate defloc in all layers
            do k=1,NLAY
              alpha = MIN((time-bed_age(i,j,k))/t_dfloc,
     &                  1.0)
              cff3=0.0
              do ised=1,NST
                cff3=cff3+bed_mass(i,j,k,nnew,ised)
              enddo
              cff1 = 0.0
              do ised=imud1,imud2
                cff1 = cff1+bed_mass(i,j,k,nnew,ised)
              enddo
              do ised=imud1,imud2
                indx= imud1-NSAND
                masseq(indx)=mud_frac_eq(indx)*cff1
                bed_mass(i,j,k,nnew,ised)=alpha*masseq(indx)+
     &                  bed_mass(i,j,k,nnew,ised)*(1.0-alpha)
                bed_frac(i,j,k,ised)=bed_mass(i,j,k,nnew,ised)/cff3
              enddo
            enddo
          enddo I_LOOP_DFL
#  endif

#  if defined COHESIVE_BED || defined MIXED_BED
         I_LOOP_CB2: do i=Istr,Iend
!
!  Key to this algorithm: "mass available depth" (mad) [kg/m2]
!  present tau crit profile (tcp)
!  representative tau crit profile (tcr)
!  Compute depth and mass depth of sediment column
!  Compute cumulative depth
!  Compute mass depth
          bmz(1) = 0.0
          do ised=1,NST
             bmz(1) = bmz(1)+bed_mass(i,j,1,nnew,ised)
          enddo
          do k=2,NLAY
             bmz(k) = bmz(k-1)
             do ised=1,NST
                bmz(k)=bmz(k)+bed_mass(i,j,k,nnew,ised)
             enddo
          enddo

!  Calculate representative critical shear stress profile
!  Note that the values are for the TOP of the layer...we
!  assume the bottom of the bottom layer has tcr = tcr_max
          tcr(1) = tcr_min
          do k=2,NLAY
             tcr(k) = tcr_min
             if (bmz(k-1).GT.eps) then
!                tcr(k) = exp((log(bmz(k-1))-                            &
!     &               bottom(i,j,idoff))/                                &
!     &               bottom(i,j,idslp))
                tcr(k) = exp((log(bmz(k-1))-
     &                    tcr_off)/
     &                    tcr_slp)
             endif
             tcr(k) = MIN( MAX( tcr(k), tcr_min), tcr_max )
          enddo
!          ks=Ksed-2
!          do k=NLAY,NLAY-ks,-1
!             taucb(i,j,k)=tcr(k)
!          enddo
!  Relax tau crit profile taucb(i,j,k)
!  towards representative profile tcr...100 x slower for "swelling"
!  than consolidation
!  TODO - make the factor an input parameter
          if( bed_bctr(i,j,1).LE.tcr(1)) then
!             alpha = MIN(dt/bottom(i,j,idtim),1.0)
             alpha = MIN(dt/tcr_tim,1.0)
          else
!             alpha = MIN(dt/(100.0*bottom(i,j,idtim)),1.0)
             alpha = MIN(dt/(100.0*tcr_tim),1.0)
          endif
          do k=1,NLAY
             bed_bctr(i,j,k)=bed_bctr(i,j,k)+
     &            alpha*(tcr(k)-bed_bctr(i,j,k))
             bed_bctr(i,j,k)=
     &          MIN( MAX( bed_bctr(i,j,k), tcr_min), tcr_max )
          enddo
!ALA    Enforce last layer = tcr_max
!ALA          bed_bctr(i,j,NLAY)= tcr_max
        enddo I_LOOP_CB2
#  endif /* cohesive bed */

      enddo J_LOOP3
!
!-----------------------------------------------------------------------
!  Apply periodic or gradient boundary conditions to property arrays.
!-----------------------------------------------------------------------
!
! Set boundary conditions
!
#  ifndef EW_PERIODIC
      if (EASTERN_EDGE) then
        do j=Jstr,Jend
          do k=1,NLAY
            do ised=1,NST
            bed_mass(Iend+1,j,k,nnew,ised)=
     &                       bed_mass(Iend,j,k,nnew,ised)
            bed_frac(Iend+1,j,k,ised)=
     &                        bed_frac(Iend,j,k,ised)
             enddo
             bed_thick(Iend+1,j,k)=bed_thick(Iend,j,k)
             bed_poros(Iend+1,j,k)=bed_poros(Iend,j,k)
             bed_age(Iend+1,j,k)=bed_age(Iend,j,k)
#   if defined MIXED_BED || defined COHESIVE_BED
             bed_bctr(Iend+1,j,k)=bed_bctr(Iend,j,k)
#   endif
           enddo
           bot_thick(Iend+1,j)=bot_thick(Iend,j)
         enddo
        endif
      if (WESTERN_EDGE) then
        do j=Jstr,Jend
          do k=1,NLAY
            do ised=1,NST
            bed_mass(Istr-1,j,k,nnew,ised)=
     &                       bed_mass(Istr,j,k,nnew,ised)
            bed_frac(Istr-1,j,k,ised)=
     &                        bed_frac(Istr,j,k,ised)
             enddo
             bed_thick(Istr-1,j,k)=bed_thick(Istr,j,k)
             bed_poros(Istr-1,j,k)=bed_poros(Istr,j,k)
             bed_age(Istr-1,j,k)=bed_age(Istr,j,k)
#   if defined MIXED_BED || defined COHESIVE_BED
             bed_bctr(Istr-1,j,k)=bed_bctr(Istr,j,k)
#   endif
           enddo
           bot_thick(Istr-1,j)=bot_thick(Istr,j)
        enddo
      endif
#  endif
#  ifndef NS_PERIODIC
      if (NORTHERN_EDGE) then
        do i=Istr,Iend
          do k=1,NLAY
            do ised=1,NST
              bed_mass (i,Jend+1,k,nnew,ised)=
     &                            bed_mass(i,Jend,k,nnew,ised)
              bed_frac (i,Jend+1,k,ised)=
     &                            bed_frac(i,Jend,k,ised)
            enddo
            bed_thick(i,Jend+1,k)=bed_thick(i,Jend,k)
            bed_poros(i,Jend+1,k)=bed_poros(i,Jend,k)
            bed_age  (i,Jend+1,k)=bed_age(i,Jend,k)
#   if defined MIXED_BED || defined COHESIVE_BED
            bed_bctr (i,Jend+1,k)=bed_bctr(i,Jend,k)
#   endif
          enddo
          bot_thick  (i,Jend+1)=bot_thick(i,Jend)
        enddo
      endif
      if (SOUTHERN_EDGE) then
        do i=Istr,Iend
          do k=1,NLAY
            do ised=1,NST
              bed_mass (i,Jstr-1,k,nnew,ised)=
     &                            bed_mass(i,Jstr,k,nnew,ised)
              bed_frac (i,Jstr-1,k,ised)=
     &                            bed_frac(i,Jstr,k,ised)
            enddo
            bed_thick(i,Jstr-1,k)=bed_thick(i,Jstr,k)
            bed_poros(i,Jstr-1,k)=bed_poros(i,Jstr,k)
            bed_age  (i,Jstr-1,k)=bed_age(i,Jend,k)
#   if defined MIXED_BED || defined COHESIVE_BED
            bed_bctr (i,Jstr-1,k)=bed_bctr(i,Jend,k)
#   endif
          enddo
          bot_thick  (i,Jstr-1)=bot_thick(i,Jstr)
        enddo
      endif
#  endif
#  if !defined EW_PERIODIC && !defined NS_PERIODIC
      if (SOUTHERN_EDGE.and.WESTERN_EDGE) then
        bed_mass(Istr-1,Jstr-1,:,nnew,:)=
     &  0.5*(bed_mass(Istr,Jstr-1,:,nnew,:)
     &      +bed_mass(Istr-1,Jstr,:,nnew,:))
        bed_frac(Istr-1,Jstr-1,:,:)=
     &  0.5*(bed_frac(Istr,Jstr-1,:,:)
     &      +bed_frac(Istr-1,Jstr,:,:))
        bed_thick(Istr-1,Jstr-1,:)=
     &  0.5*(bed_thick(Istr,Jstr-1,:)
     &      +bed_thick(Istr-1,Jstr,:))
        bed_poros(Istr-1,Jstr-1,:)=
     &  0.5*(bed_poros(Istr,Jstr-1,:)
     &      +bed_poros(Istr-1,Jstr,:))
        bed_age(Istr-1,Jstr-1,:)=
     &  0.5*(bed_age(Istr,Jstr-1,:)
     &      +bed_age(Istr-1,Jstr,:))
#   if defined MIXED_BED || defined COHESIVE_BED
        bed_bctr(Istr-1,Jstr-1,:)=
     &  0.5*(bed_bctr(Istr,Jstr-1,:)
     &      +bed_bctr(Istr-1,Jstr,:))
#   endif
      endif
      if (SOUTHERN_EDGE.and.EASTERN_EDGE) then
        bed_mass(Iend+1,Jstr-1,:,nnew,:)=
     &  0.5*(bed_mass(Iend,Jstr-1,:,nnew,:)
     &      +bed_mass(Iend+1,Jstr,:,nnew,:))
        bed_frac(Iend+1,Jstr-1,:,:)=
     &  0.5*(bed_frac(Iend,Jstr-1,:,:)
     &      +bed_frac(Iend+1,Jstr,:,:))
        bed_thick(Iend+1,Jstr-1,:)=
     &  0.5*(bed_thick(Iend,Jstr-1,:)
     &      +bed_thick(Iend+1,Jstr,:))
        bed_poros(Iend+1,Jstr-1,:)=
     &  0.5*(bed_poros(Iend,Jstr-1,:)
     &      +bed_poros(Iend+1,Jstr,:))
        bed_age(Iend+1,Jstr-1,:)=
     &  0.5*(bed_age(Iend,Jstr-1,:)
     &      +bed_age(Iend+1,Jstr,:))
#   if defined MIXED_BED || defined COHESIVE_BED
        bed_bctr(Iend+1,Jstr-1,:)=
     &  0.5*(bed_bctr(Iend,Jstr-1,:)
     &      +bed_bctr(Iend+1,Jstr,:))
#   endif
      endif
      if (NORTHERN_EDGE.and.WESTERN_EDGE) then
        bed_mass(Istr-1,Jend+1,:,nnew,:)=
     &  0.5*(bed_mass(Istr,Jend+1,:,nnew,:)
     &      +bed_mass(Istr-1,Jend,:,nnew,:))
        bed_frac(Istr-1,Jend+1,:,:)=
     &  0.5*(bed_frac(Istr,Jend+1,:,:)
     &      +bed_frac(Istr-1,Jend,:,:))
        bed_thick(Istr-1,Jend+1,:)=
     &  0.5*(bed_thick(Istr,Jend+1,:)
     &      +bed_thick(Istr-1,Jend,:))
        bed_poros(Istr-1,Jend+1,:)=
     &  0.5*(bed_poros(Istr,Jend+1,:)
     &      +bed_poros(Istr-1,Jend,:))
        bed_age(Istr-1,Jend+1,:)=
     &  0.5*(bed_age(Istr,Jend+1,:)
     &      +bed_age(Istr-1,Jend,:))
#   if defined MIXED_BED || defined COHESIVE_BED
        bed_bctr(Istr-1,Jend+1,:)=
     &  0.5*(bed_bctr(Istr,Jend+1,:)
     &      +bed_bctr(Istr-1,Jend,:))
#   endif
      endif
      if (NORTHERN_EDGE.and.EASTERN_EDGE) then
        bed_mass(Iend+1,Jend+1,:,nnew,:)=
     &   0.5*(bed_mass(Iend,Jend+1,:,nnew,:)
     &       +bed_mass(Iend+1,Jend,:,nnew,:))
        bed_frac(Iend+1,Jend+1,:,:)=
     &   0.5*(bed_frac(Iend,Jend+1,:,:)
     &       +bed_frac(Iend+1,Jend,:,:))
        bed_thick(Iend+1,Jend+1,:)=
     &   0.5*(bed_thick(Iend,Jend+1,:)
     &       +bed_thick(Iend+1,Jend,:))
        bed_poros(Iend+1,Jend+1,:)=
     &   0.5*(bed_poros(Iend,Jend+1,:)
     &       +bed_poros(Iend+1,Jend,:))
        bed_age(Iend+1,Jend+1,:)=
     &   0.5*(bed_age(Iend,Jend+1,:)
     &       +bed_age(Iend+1,Jend,:))
#   if defined MIXED_BED || defined COHESIVE_BED
        bed_bctr(Iend+1,Jend+1,:)=
     &   0.5*(bed_bctr(Iend,Jend+1,:)
     &       +bed_bctr(Iend+1,Jend,:))
#   endif
      endif
#  endif

#  if defined BD_EXCHANGE
      do ised=1,NST
        do k=1,NLAY
          call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_frac(START_2D_ARRAY,k,ised))
          call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_mass(START_2D_ARRAY,k,nnew,ised))
         enddo
      enddo
      do k=1,NLAY
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_poros(START_2D_ARRAY,k))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_thick(START_2D_ARRAY,k))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_age(START_2D_ARRAY,k))
#   if defined MIXED_BED || defined COHESIVE_BED
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bed_bctr(START_2D_ARRAY,k))
#   endif
      enddo
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 bot_thick(START_2D_ARRAY))
#  endif
!
!-----------------------------------------------------------------------
!  Store total bed thickness
!-----------------------------------------------------------------------
!
#  ifdef MORPHODYN
      do j=Jstr,Jend
        do i=Istr,Iend
          bed_thick_tot(i,j,nnew)=0.
          do k=1,NLAY
            bed_thick_tot(i,j,nnew)=bed_thick_tot(i,j,nnew)+
     &                              bed_thick(i,j,k)
          enddo
        enddo
      enddo
!
! Set boundary conditions
!
#   ifndef EW_PERIODIC
      if (EASTERN_EDGE) then
        do j=Jstr,Jend
          bed_thick_tot(Iend+1,j,nnew)=bed_thick_tot(Iend,j,nnew)
         enddo
      endif
      if (WESTERN_EDGE) then
        do j=Jstr,Jend
          bed_thick_tot(Istr-1,j,nnew)=bed_thick_tot(Istr,j,nnew)
        enddo
      endif
#   endif
#   ifndef NS_PERIODIC
      if (NORTHERN_EDGE) then
        do i=Istr,Iend
          bed_thick_tot(i,Jend+1,nnew)=bed_thick_tot(i,Jend,nnew)
        enddo
      endif
      if (SOUTHERN_EDGE) then
        do i=Istr,Iend
          bed_thick_tot(i,Jstr-1,nnew)=bed_thick_tot(i,Jstr,nnew)
        enddo
      endif
#   endif
#   if !defined EW_PERIODIC && !defined NS_PERIODIC
      if (SOUTHERN_EDGE.and.WESTERN_EDGE) then
        bed_thick_tot(Istr-1,Jstr-1,nnew)=
     &  0.5*(bed_thick_tot(Istr,Jstr-1,nnew)
     &      +bed_thick_tot(Istr-1,Jstr,nnew))
      endif
      if (SOUTHERN_EDGE.and.EASTERN_EDGE) then
        bed_thick_tot(Iend+1,Jstr-1,nnew)=
     &  0.5*(bed_thick_tot(Iend,Jstr-1,nnew)
     &      +bed_thick_tot(Iend+1,Jstr,nnew))
      endif
      if (NORTHERN_EDGE.and.WESTERN_EDGE) then
        bed_thick_tot(Istr-1,Jend+1,nnew)=
     &  0.5*(bed_thick_tot(Istr,Jend+1,nnew)
     &      +bed_thick_tot(Istr-1,Jend,nnew))
      endif
      if (NORTHERN_EDGE.and.EASTERN_EDGE) then
        bed_thick_tot(Iend+1,Jend+1,nnew)=
     &   0.5*(bed_thick_tot(Iend,Jend+1,nnew)
     &       +bed_thick_tot(Iend+1,Jend,nnew))
      endif
#   endif
!
#   if defined BD_EXCHANGE
       call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                         bed_thick_tot(START_2D_ARRAY,nnew))
#   endif
!
!-----------------------------------------------------------------------
!  Store bed evolution
!-----------------------------------------------------------------------
!
      do j=Jstr-1,Jend+1
        do i=Istr-1,Iend+1
          dh(i,j)=bed_thick_tot(i,j,nstp)-bed_thick_tot(i,j,nnew)
        enddo
      enddo
#  endif
      return
      end

# endif
!=======================================================================
!                                                                      !
!  This routine computed sediment surface layer (sediment-water        !
!  interface) properties.                                              !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G.    !
!    Arango, 2008:  Development of a three-dimensional,  regional,     !
!    coupled wave, current, and sediment-transport model, Computers    !
!    & Geosciences, 34, 1284-1306.                                     !
!                                                                      !
!=======================================================================
!
      subroutine sed_surf_tile (Istr,Iend,Jstr,Jend)
!
      implicit none
# include "param.h"
# include "bbl.h"
# include "grid.h"
# include "scalars.h"
# include "sediment.h"
      integer Istr,Iend,Jstr,Jend
      integer i,j,ised
      real cff,cff1,cff2,cff3,cff4,eps
      parameter (eps=1.D-20)
# if defined BED_HIDEXP || defined BED_ARMOR
      integer ised0
      real ph,pe, phik,phim,phis, invlog2
# endif
!
# include "compute_auxiliary_bounds.h"

# define I_EXT_RANGE Istr-1,Iend+1
# define J_EXT_RANGE Jstr-1,Jend+1
!
!-----------------------------------------------------------------------
! Compute Hiding/Exposure and Armoring
!
! Hiding /Exposure:
! Resuspension of smaller size classes may be reduced due to hiding,
! and the coarser grains resuspend more easily due to exposure.
!
! Armoring:
! In addition, bed erosion may become limited when selective entrainment
! of the fine fraction causes "bed armoring": the remaining coarser
! grains inhibit resuspension of the fine grains.
!-----------------------------------------------------------------------
!
# ifdef BED_HIDEXP
#  ifdef BED_HIDEXP_WULIN
!
! Modulate critical stress tau_ce with hiding/exposure effects
! using hidden and exposed probabilities of each sediment particle
! in the bed material (Wu & Lin 2014)
!
      do j=J_EXT_RANGE
        do i=I_EXT_RANGE
          do ised=1,NST
            ph=0.
            pe=0.
            do ised0=1,NST
              cff =1./(Sd(ised)+Sd(ised0))
              ph=ph+bed_frac(i,j,1,ised0)*cff*Sd(ised0)  ! prob hiding
              pe=pe+bed_frac(i,j,1,ised0)*cff*Sd(ised )  ! prob exposure
            enddo
            tau_ce_2d(i,j,ised)=tau_ce(ised)*(pe/ph)**(-0.6)
          enddo
        enddo
      enddo
#  else
!
!  Hiding/exposure factor is function of sediment size:
!  (d(ised)/dm)^(-gamma), e.g., Wilcock & Crowe (2003)
!
      invlog2=1./log(2.)
      do j=J_EXT_RANGE
        do i=I_EXT_RANGE
          phim=0.
          do ised=1,NST
            phik=-log(Sd(ised))*invlog2
            phim=phim+phik*bed_frac(i,j,1,ised)
          enddo
          cff=2.**(phim)  ! 1/mean(Sd)
          do ised=1,NST
            cff1=1-0.67/(1+exp(1.5-Sd(ised)*cff)) ! or 1.
            tau_ce_2d(i,j,ised)=tau_ce(ised)*(Sd(ised)*cff)**(-cff1)
          enddo
        enddo
      enddo
#  endif
# else
      do ised=1,NST
        do j=J_EXT_RANGE
          do i=I_EXT_RANGE
            tau_ce_2d(i,j,ised)=tau_ce(ised)
          enddo
        enddo
      enddo
# endif
!
# ifdef BED_ARMOR
!
! Bed erosion may become limited when selective entrainment of the fine
! fraction causes "bed armoring": the remaining coarser grains inhibit
! resuspension of the fine grains.
! Formulation in Blaas et al. (2007) from Garcia and Parker (1991)
!
      invlog2=1./log(2.)
      do j=J_EXT_RANGE
        do i=I_EXT_RANGE
          phim=0.
          phis=0.
          do ised=1,NST
            phik=-log(Sd(ised))*invlog2
            phim=phim+phik*bed_frac(i,j,1,ised) ! mean
          enddo
          do ised=1,NST
            phik=-log(Sd(ised))*invlog2
            phis=phis+((phik-phim)**2)*bed_frac(i,j,1,ised) ! std
          enddo
          phis=SQRT(phis)
          cff=1./(1.-0.29*phis)**5. ! 1/lambda**5
          do ised=1,NST
            tau_ce_2d(i,j,ised)=tau_ce_2d(i,j,ised)*cff
          enddo
        enddo
      enddo
# endif

# undef DEBUG_ARMOR
# ifdef DEBUG_ARMOR
      do ised=1,NST
        MPI_master_only write(*,*) '================================'
        MPI_master_only write(*,*) ' ised      :',ised
        MPI_master_only write(*,*) ' Sd        :',1.e3*Sd(ised)
        MPI_master_only write(*,*) ' bed_frac  :',bed_frac(25,1,1,ised)
        MPI_master_only write(*,*) ' tau_ce_2d :',rho0
     &                                            *tau_ce_2d(25,1,ised)
        MPI_master_only write(*,*) '================================'
      enddo
# endif
!
!-----------------------------------------------------------------------
!  Update mean surface properties.
!  Sd50 must be positive definite, due to BBL routines.
!  Srho must be >1000, due to (s-1) in BBL routines
!-----------------------------------------------------------------------
!
      do j=Jstr-1,Jend+1
        do i=Istr-1,Iend+1
          cff1=1.
          cff2=1.
          cff3=1.
          cff4=1.
          do ised=1,NST
            cff1=cff1*tau_ce_2d(i,j,ised)**bed_frac(i,j,1,ised)
            cff2=cff2*Sd(ised)**bed_frac(i,j,1,ised)
            cff3=cff3*(wsed(ised)+eps)**bed_frac(i,j,1,ised)
            cff4=cff4*Srho(ised)**bed_frac(i,j,1,ised)
          enddo
          taucb(i,j)=cff1
          Ssize(i,j)=cff2
          w_set(i,j)=cff3
          Sdens(i,j)=MAX(cff4,1050.)
        enddo
      enddo
!
! Set boundary conditions
!
# ifndef EW_PERIODIC
      if (EASTERN_EDGE) then
        do j=Jstr,Jend
           taucb(Iend+1,j)=taucb(Iend,j)
           Ssize(Iend+1,j)=Ssize(Iend,j)
           w_set(Iend+1,j)=w_set(Iend,j)
           Sdens(Iend+1,j)=Sdens(Iend,j)
        enddo
      endif
      if (WESTERN_EDGE) then
        do j=Jstr,Jend
           taucb(Istr-1,j)=taucb(Istr,j)
           Ssize(Istr-1,j)=Ssize(Istr,j)
           w_set(Istr-1,j)=w_set(Istr,j)
           Sdens(Istr-1,j)=Sdens(Istr,j)
        enddo
      endif
# endif
# ifndef NS_PERIODIC
      if (NORTHERN_EDGE) then
        do i=Istr,Iend
          taucb  (i,Jend+1)=taucb(i,Jend)
          Ssize  (i,Jend+1)=Ssize(i,Jend)
          w_set  (i,Jend+1)=w_set(i,Jend)
          Sdens  (i,Jend+1)=Sdens(i,Jend)
        enddo
      endif
      if (SOUTHERN_EDGE) then
        do i=Istr,Iend
          taucb  (i,Jstr-1)=taucb(i,Jstr)
          Ssize  (i,Jstr-1)=Ssize(i,Jstr)
          w_set  (i,Jstr-1)=w_set(i,Jstr)
          Sdens  (i,Jstr-1)=Sdens(i,Jstr)
        enddo
      endif
# endif
# if !defined EW_PERIODIC && !defined NS_PERIODIC
      if (SOUTHERN_EDGE.and.WESTERN_EDGE) then
        taucb(Istr-1,Jstr-1)=0.5*(taucb(Istr,Jstr-1)+taucb(Istr-1,Jstr))
        Ssize(Istr-1,Jstr-1)=0.5*(Ssize(Istr,Jstr-1)+Ssize(Istr-1,Jstr))
        w_set(Istr-1,Jstr-1)=0.5*(w_set(Istr,Jstr-1)+w_set(Istr-1,Jstr))
        Sdens(Istr-1,Jstr-1)=0.5*(Sdens(Istr,Jstr-1)+Sdens(Istr-1,Jstr))
      endif
      if (SOUTHERN_EDGE.and.EASTERN_EDGE) then
        taucb(Iend+1,Jstr-1)=0.5*(taucb(Iend,Jstr-1)+taucb(Iend+1,Jstr))
        Ssize(Iend+1,Jstr-1)=0.5*(Ssize(Iend,Jstr-1)+Ssize(Iend+1,Jstr))
        w_set(Iend+1,Jstr-1)=0.5*(w_set(Iend,Jstr-1)+w_set(Iend+1,Jstr))
        Sdens(Iend+1,Jstr-1)=0.5*(Sdens(Iend,Jstr-1)+Sdens(Iend+1,Jstr))
      endif
      if (NORTHERN_EDGE.and.WESTERN_EDGE) then
        taucb(Istr-1,Jend+1)=0.5*(taucb(Istr,Jend+1)+taucb(Istr-1,Jend))
        Ssize(Istr-1,Jend+1)=0.5*(Ssize(Istr,Jend+1)+Ssize(Istr-1,Jend))
        w_set(Istr-1,Jend+1)=0.5*(w_set(Istr,Jend+1)+w_set(Istr-1,Jend))
        Sdens(Istr-1,Jend+1)=0.5*(Sdens(Istr,Jend+1)+Sdens(Istr-1,Jend))
      endif
      if (NORTHERN_EDGE.and.EASTERN_EDGE) then
        taucb(Iend+1,Jend+1)=0.5*(taucb(Iend,Jend+1)+taucb(Iend+1,Jend))
        Ssize(Iend+1,Jend+1)=0.5*(Ssize(Iend,Jend+1)+Ssize(Iend+1,Jend))
        w_set(Iend+1,Jend+1)=0.5*(w_set(Iend,Jend+1)+w_set(Iend+1,Jend))
        Sdens(Iend+1,Jend+1)=0.5*(Sdens(Iend,Jend+1)+Sdens(Iend+1,Jend))
      endif
# endif
# if defined BD_EXCHANGE
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 taucb(START_2D_ARRAY))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 Ssize(START_2D_ARRAY))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 w_set(START_2D_ARRAY))
         call exchange_r2d_tile (Istr,Iend,Jstr,Jend,
     &                 Sdens(START_2D_ARRAY))
# endif

# undef I_EXT_RANGE
# undef J_EXT_RANGE

      return
      end
!
!====================================================================
!
!  This routine reads in and processes sediment parameters and
!  some initial data from sediments.in file.
!
!=====================================================================
!
       subroutine init_sediment
# ifdef SED_FLOCS
        USE flocmod, ONLY : flocmod_alloc, flocmod_init
        USE flocmod, ONLY : f_ws
# endif
      implicit none
# include "param.h"
# include "grid.h"
# include "ocean3d.h"
# include "ncscrum.h"
# include "scalars.h"
# include "sediment.h"
!
      integer Nclass, i, icard, iunit, lstr, nl
      integer lenstr
      real rdum, cff1,cff2
      logical ldum
      parameter (iunit=50)
# ifdef SED_FLOCS
      real :: f_mneg_param, f_fp, f_fy, f_cfcst, f_dmin_frag
# endif
!
!---------------------------------------------------------------------
!  Read in initial float locations.
!---------------------------------------------------------------------
!
      lstr=lenstr(sedname)

      open(iunit,file=sedname(1:lstr),form='formatted',
     &  status='old', err=195)
!
!  Read input parameters according to their input card number.
!
      icard=0
      do while (icard.lt.99)
!
!  Read in sediment case title.
!
        if (icard.eq.1) then
          read(iunit,'(a)',err=60) Stitle
          lstr=lenstr(Stitle)
          write(stdout,10) Stitle(1:lstr)
  10      format(1x,'(',a,')',/)
!
!  Read in parameter and initial values per size class.
!
        elseif (icard.eq.2) then
          i=0
          Nclass=0
          do i=1,NST
           read(iunit,*,err=40) Sd(i), Csed(i), Srho(i), Wsed(i),
     &                          Erate(i), tau_ce(i), tau_cd(i),
     &                          (Bfr(nl,i),nl=1,NLAY)

!
! Re-compute critical shear stress tau_ce (Wu & Lin, 2014)
!
# ifdef TAU_CRIT_SHIELDS
           cff1=10*Sd(i)*(g*(Srho(i)/rho0 -1.))**0.33        ! Sd_star (m)
           cff2=0.3/(1.+1.2*cff1)+0.055*(1.-exp(-0.02*cff1)) ! Critical shield
           tau_ce(i)=cff2*Sd(i)/1000*g*(Srho(i)-rho0)        ! tau_ce (N/m2)
# elif defined TAU_CRIT_WULIN
           cff2=0.03                                         ! Critical shield
           tau_ce(i)=cff2*Sd(i)/1000*g*(Srho(i)-rho0)        ! tau_ce (N/m2)
# endif

           MPI_master_only write(stdout,'(/A,2x,i2)')
     &     'Sediment parameters for grain-size class:', i

           MPI_master_only write(stdout,'(7(/f10.5,2x,A))')
     &     Sd(i),    'Sd      Median sediment grain diameter (mm).',
     &     Csed(i),  'Csed    Analyt. init. of sed. conc. (kg/m3).',
     &     Srho(i),  'Srho    Sediment grain density (kg/m3).',
     &     Wsed(i),  'Wsed    Particle settling velocity (mm/s).',
     &     Erate(i), 'Erate   Surface erosion rate (kg/(m2s)).',
     &     tau_ce(i),'tau_ce  Critical shear for erosion (N/m2).',
     &     tau_cd(i),'tau_cd  Critical shear for deposition (N/m2).'

           do nl=1,NLAY
            MPI_master_only write(stdout,'(f10.5,2x,A,2x,i2)')
     &      Bfr(nl,i), 'bed_frac   Volume fraction, layer:',nl
           enddo

           Sd(i)=Sd(i)/1000.         ! mm --> m
           Wsed(i)=Wsed(i)/1000.     ! mm/s --> m/s
           tau_ce(i)=tau_ce(i)/rho0  ! N/m2 --> m2/s2
           tau_cd(i)=tau_cd(i)/rho0  ! N/m2 --> m2/s2
           Nclass=Nclass+1
          enddo

        elseif (icard.eq.3) then
           MPI_master_only write(stdout,'(/A)')
     &     'Bed parameters for all bed levels'

          read(iunit,*,err=60) (Bthk(nl),nl=1,NLAY)
           do nl=1,NLAY
            MPI_master_only write(stdout,'(f10.5,2x,A,2x,i2,2x,A)')
     &      Bthk(nl), 'Bthk   Init. bed thickness, layer:',nl,'(m).'
           enddo

        elseif (icard.eq.4) then

          read(iunit,*,err=60) (Bpor(nl),nl=1,NLAY)
           do nl=1,NLAY
            MPI_master_only write(stdout,'(f10.5,2x,A,2x,i2,2x,A)')
     &      Bpor(nl), 'Bpor   Init. bed porosity, layer:',nl,'(m).'
           enddo

        elseif (icard.eq.5) then
          read(iunit,*,err=60) Hrip
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      Hrip, 'Hrip   Init. ripple roughness (m).'

        elseif (icard.eq.6) then
          read(iunit,*,err=60) Lrip
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      Lrip, 'Lrip   Init. ripple height (m).'

        elseif (icard.eq.7) then
          read(iunit,*,err=60) bedload_coeff
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      bedload_coeff, 'bedload_coeff   Bedload rate.'

        elseif (icard.eq.8) then
          read(iunit,*,err=60) morph_fac
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      morph_fac, 'morph_fac   Morphological factor.'

# ifdef MIXED_BED
        elseif (icard.eq.9) then
          read(iunit,*,err=60) transC
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      transC, 'transC   Cohesive transition.'

        elseif (icard.eq.10) then
          read(iunit,*,err=60) transN
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      transN, 'transN   Non cohesive transition.'
# else
        elseif (icard.eq.9) then
          read(iunit,*,err=60) rdum
        elseif (icard.eq.10) then
          read(iunit,*,err=60) rdum
# endif
# if defined COHESIVE_BED || defined MIXED_BED
        elseif (icard.eq.11) then
          read(iunit,*,err=60) tcr_min
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      tcr_min, 'tcr_min   Mud min taucr.'

        elseif (icard.eq.12) then
          read(iunit,*,err=60) tcr_max
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      tcr_max, 'tcr_max   Mud max taucr.'

        elseif (icard.eq.13) then
          read(iunit,*,err=60) tcr_slp
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      tcr_slp, 'tcr_slp   Mud taucr slope.'

        elseif (icard.eq.14) then
          read(iunit,*,err=60) tcr_off
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      tcr_off, 'tcr_off   Mud profile offset.'

        elseif (icard.eq.15) then
          read(iunit,*,err=60) tcr_tim
            MPI_master_only write(stdout,'(f10.1,2x,A)')
     &      tcr_tim, 'tcr_tim   Mud consolidation rate.'

# endif
# ifdef SED_FLOCS
        elseif (icard.eq.16) then
          read(iunit,*,err=60) l_ADS,l_ASH,l_COLLFRAG,l_testcase
            MPI_master_only write(stdout,'(L10,2x,A)')
     &      l_ADS, 'l_ADS   Mud differential settling.'
            MPI_master_only write(stdout,'(L10,2x,A)')
     &      l_ASH, 'l_ASH   Mud shear aggregation.'
            MPI_master_only write(stdout,'(L10,2x,A,1x,A)')
     &      l_COLLFRAG, 'l_COLLFRAG   Mud collision-induced',
     &                   'fragmentation.'
            MPI_master_only write(stdout,'(L10,2x,A,1x,A)')
     &      l_testcase, 'l_testcase   Mud boolean for values',
     &                  'from Verney et al.'

        elseif (icard.eq.17) then
          read(iunit,*,err=60) f_dp0  ,f_nf   ,f_dmax,f_nb_frag,
     &                         f_alpha,f_beta ,f_ater,f_ero_frac,
     &                         f_ero_nbfrag,f_collfragparam,
     &                         f_clim,f_ero_iv
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      f_dp0, 'f_dp0   Mud Primary particle size(m).'
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      f_nf, 'f_nf   Mud floc fractal dimension.'
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      f_dmax, 'f_dmax   Mud maximum diameter.'
            MPI_master_only write(stdout,'(f10.5,2x,A,1x,A)')
     &      f_nb_frag, 'f_nb_frag   Mud Number of fragments by shear',
     &                             'erosion.'
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      f_alpha, 'f_alpha   Mud flocculation efficiency.'
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      f_beta, 'f_beta   Mud shear fragmentation rate.'
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      f_ater, 'f_ater   Mud ternary breakup.'
            MPI_master_only write(stdout,'(f10.5,2x,A,1x,A)')
     &      f_ero_frac, 'f_ero_frac   Mud fraction of shear',
     &             'fragmentation term transfered to shear erosion.'
            MPI_master_only write(stdout,'(f10.5,2x,A,1x,A)')
     &      f_ero_nbfrag, 'f_ero_nbfrag   Mud Number of fragments',
     &             'induced by shear erosion.'
            MPI_master_only write(stdout,'(f10.5,2x,A,1x,A)')
     &      f_collfragparam, 'f_collfragparam   Mud Fragmentation rate',
     &             'for collision-induced breakup.'
            MPI_master_only write(stdout,'(f10.5,2x,A,1x,A)')
     &      f_clim, 'f_clim   Mud Min concentration below which',
     &           'flocculation processes are not calc..'
            MPI_master_only write(stdout,'(I10,2x,A)')
     &      f_ero_iv, 'f_ero_iv   Mud Fragment size class.'

#  ifdef SED_DEFLOC
        elseif (icard.eq.18) then
          read(iunit,*,err=60) (mud_frac_eq(nl),nl=1,NMUD)
           do nl=1,NMUD
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      mud_frac_eq(nl),'mud_frac_eq   Equilibrium fractional class'
           enddo
        elseif (icard.eq.19) then
          read(iunit,*,err=60) t_dfloc
            MPI_master_only write(stdout,'(f10.5,2x,A)')
     &      t_dfloc, 't_dfloc   Floc decomposition rate.'
#  else
        elseif (icard.eq.18) then
          read(iunit,*,err=60) rdum
        elseif (icard.eq.19) then
          read(iunit,*,err=60) rdum
#  endif
# else
        elseif (icard.eq.16) then
          read(iunit,*,err=60) ldum
        elseif (icard.eq.17) then
          read(iunit,*,err=60) rdum
        elseif (icard.eq.18) then
          read(iunit,*,err=60) rdum
        elseif (icard.eq.19) then
          read(iunit,*,err=60) rdum
# endif
        endif ! icard
!
!  Read until last input card ID.
!
        read(iunit,*,err=60) icard
      enddo
      goto 90
!
!  Error while reading input parameters.
!
  40  write(stdout,50) icard, i, Nclass, sedname
  50  format(/,' INIT_SEDIMENTS - error reading variables in card: ',
     &             i2, ', entry: ',i3,/,15x,
     &             'nclass:',i3, 'input script: ',a)

  60  write(stdout,80) icard, sedname
  80  format(/,' INIT_SEDIMENTS - error while reading input card: ',
     &       i2,15x,'in input script: ',a)

 195  write(stdout,205) sedname
 205  format(/,'sediment inputfile : ',A,/,
     &   ' inputfile not found or bad sediment initialization')
      stop

  90  close(iunit)
      write(stdout,100) Nclass
 100  format(/,
     &       'Number of size classes in sediment computation:'i2)


       do i=1,NST
         tau_ce_2d(:,:,i)=tau_ce(i)
       enddo
# ifdef SED_FLOCS
        ! some parameters are not read in sediment.in
        ! f_mneg_param, f_fp, f_fy, f_cfcst, f_dmin_frag
        ! values from paraMUSTANG_default.txt
        ! && f_dmax not used but read in sediment.in
        f_mneg_param=0.001
        f_cfcst = 0.1875
        f_fp = 0.1
        f_fy = 1.0e-10
        f_dmin_frag = 0.00001
        CALL flocmod_alloc(NMUD)
        CALL flocmod_init(l_ADS, l_ASH, l_COLLFRAG,
     &      f_dp0, f_nf, f_nb_frag, f_alpha, f_beta, f_ater,
     &      f_ero_frac, f_ero_nbfrag, f_ero_iv, f_mneg_param,
     &      f_collfragparam, f_dmin_frag, f_cfcst, f_fp, f_fy,
     &      f_clim, Sd(imud1:imud1+NMUD-1) ,
     &      Srho(imud1:imud1+NMUD-1),
     &      rho0, l_testcase, stdout)
        Wsed(imud1:imud1+NMUD-1) = f_ws(1:NMUD)
# endif
      return
      end

# ifdef BEDLOAD_VANDERA
!======================================================================!
!======================================================================!
!                          VANDERA FUNCTIONS                           !
!                                                                      !
!   Subroutines and functions required for the sediment bedload        !
!   calculations using Van der A's formulations.                       !
!                                                                      !
! Tarandeep S. Kalra; Chris Sherwood; John C. Warner                   !
! ... adapted by P. Marchesiello (2021) for CROCO
!======================================================================!
!======================================================================!
!
      subroutine sandload_vandera( wavecycle,
     &                                   Hs, Td,  depth, RR,
     &                              d50, d50_mix, rhos, c_w,
     &                                             eta, dsf,
     &                       T_i, T_iu, uhat_i, mag_theta_i,
     &                                         om_ii, om_iy )
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real wavecycle
      real Hs, Td, depth, RR
      real d50, d50_mix, rhos, c_w
      real eta, dsf
      real T_i, T_iu
      real uhat_i, mag_theta_i
      real om_ii, om_iy
!
! local variables
!
      real m_fac,n_fac
      parameter(m_fac=11., n_fac=1.2)
      real alpha_fac
      parameter(alpha_fac=8.2)
      real xi
      parameter(xi=1.7) ! Based on Santoss_core.m
      real eps
      parameter(eps=1.0E-14)
      real eps_eff, P
      real om_i, theta_diff, theta_ieff, theta_cr
      real w_s, what, w_sc, delta
      real cff, cff1
      real w_s_calc, what_calc, theta_cr_calc
!
! Layer thickness
!
      if (eta.gt.0.) then
        delta=eta   ! For ripple regime
      else
        delta=dsf   ! For sheet flow regime
      endif
!
! Find settling velocity based on Soulsby (1997).
! VA2013 Use 0.8*d50 for settling velocity (text under equation 28).
!
      w_s=w_s_calc(0.8*d50, rhos)
!
! Correction by vertical advection (Stokes 2nd-order)
!    wavecycle=+1 => VA2013 Eq 29, for crest cycle
!    wavecycle=-1 => VA2013 Eq 30, for trough cycle
!
      what=what_calc(Hs, Td, depth, RR, delta)
      w_sc=MAX(w_s+wavecycle*what,eps)
!
! VA2013 Eq 27-28, Phase lag parameter
!
      cff=(1.-wavecycle*xi*uhat_i)/c_w
      cff1=1./MAX(2.*(T_i-T_iu)*w_sc, eps)
      P=alpha_fac*delta*cff*cff1
!
! Hiding/exposure correction factor
!
!      eps_eff=(d50/d50_mix)**0.25
      eps_eff=1.
!
! CRS for multiple sed types
!
      theta_ieff=eps_eff*mag_theta_i
!
! Find critical Shields parameters based on Soulsby (1997).
!
      theta_cr=theta_cr_calc(d50, rhos)
!
! VA2013 Eq 2: Sandload entrained during each half-cycle
!
      theta_diff=MAX((theta_ieff-theta_cr),0.)
      om_i=m_fac*(theta_diff)**n_fac
!
! VA2013 Eq 23-26: Sandload transported within the
!        half-cycle (om_ii) or in the next (om_iy)
!
      if (P.lt.1) then
        om_ii=om_i
        om_iy=0.
      else
        om_ii=om_i/P
        om_iy=om_i-om_ii
      endif
!
      return
      end
!
!======================================================================
!
      subroutine full_wave_cycle_stress_factors( d50, d90, osmgd,
     &                                            Td, c_w, depth,
     &                                               phi_curwave,
     &                                            RR, uhat, ahat,
     &                                                umax, umin,
     &                                                 mag_bstrc,
     &                                      T_c, T_t, T_cu, T_tu,
     &                                                ksd, ustrc,
     &                                          delw, udelta, fd,
     &                                   alpha, eta, ksw, tau_wRe )
!
! Iterative solution to obtain current and wave related bed roughness,
! computed from converged time-averaged Shields parameter (VA2013
! Apendix A). Also, compute full-cycle friction factor and wave Reynolds
! stress
!
!  This subroutine returns the following:
!  eta                 : ripple height
!  udelta              : current velocity at the wave boundary layer
!  fd                  : current friction factor
!  tau_wRe             : Wave averaged Reynolds stress
!  T_c, T_t, T_cu, T_tu: Updated time periods in half cycles
!                        based on current velocity
!
      implicit none
#  include "param.h"
#  include "scalars.h"
!
!  Imported variable declarations.
!
! Input the crest or trough half cycle velocity
! d50 -- grain size in meters
! Different for crest and trough half cycles
!
      real d50, d90, osmgd
      real Td, c_w, depth
      real phi_curwave
      real RR, uhat, ahat
      real umax, umin
      real mag_bstrc
      real T_c, T_t, T_cu, T_tu
      real udelta, delw, fd
      real alpha, eta, ksw, tau_wRe
!
!  Local variables
!
      integer iter
      integer total_iters
      parameter (total_iters=15)
      real tol,von_k
      parameter(tol=0.001, von_k=0.41)
      real eps
      parameter(eps=1.0E-14)
      real crs_fac
      parameter(crs_fac=1.)
      real theta_avg_old, theta_avg, theta_hat_i
      real psi  ! maximum mobility number
      real rlen ! ripple length
      real ksd, ustrc
      real fw
      real alpha_w, fwd
      real ustarw
      real ksw_calc, ksd_calc, fw_calc
      real fd_calc_madsen, fd_calc_santoss, mu_calc
!
! Maximum mobility number at crest and trough
! For irregular waves, use Rayleigh distributed maximum value
! - VA2013, text under equation Appendix B.4
!
      psi=(1.27*uhat)**2*osmgd
!
! Use Appendix B eqn B.1 and B.2 to get ripple height and length
!
      call ripple_dim(psi, d50, eta, rlen)
!
      eta=eta*ahat
      rlen=rlen*ahat
!
! Initiliaze with theta_avg=0 and theta_hat_i=theta_avg
!
      theta_avg=0.
      theta_avg_old=0.
!
! Calculate friction factor fd
!
#  ifdef BEDLOAD_VANDERA_FD_MADSEN
      fd=fd_calc_madsen(udelta, mag_bstrc)
#  elif defined BEDLOAD_VANDERA_ZEROCURR
      fd=0.
#  endif
!
      do iter=1,total_iters !---- ! Iterative solution
!
! Wave related bed roughness - VA2013 A.5
!
        ksw=ksw_calc(d50, mu_calc(d50), theta_avg, eta, rlen)
!
! Full-cycle wave friction factor - VA2013 Appendix Eq A.4
!
        fw=fw_calc(ahat, ksw)
!
#  if defined BEDLOAD_VANDERA_FD_SANTOSS &&\
     !defined BEDLOAD_VANDERA_ZEROCURR
!
! Current-related bed roughness - VA2013 Appendix A.1
!
        ksd=ksd_calc(d50, d90, mu_calc(d50), theta_avg, eta, rlen)
!
! Friction factor from ksd, udelta, delta
!
        fd=fd_calc_santoss(udelta, delw, ksd)
!
#  endif
!
! Time-averaged absolute Shields stress - VA2013 Appendix Eq. A.3
!
#  ifdef BEDLOAD_VANDERA_ZEROCURR
        theta_avg=osmgd*0.25*fw*uhat**2
#  else
        theta_avg=osmgd*(0.50*fd*udelta**2.+
     &                   0.25*fw*uhat**2.)
#  endif
!
        if (ABS(theta_avg-theta_avg_old).lt.tol) then
          EXIT
        endif
        theta_avg_old=theta_avg
!
      enddo  ! end iterations ----
!
! Full-cycle current friction factor - VA2013 Eq 20
!
#  ifdef BEDLOAD_VANDERA_ZEROCURR
      fwd=fw
      alpha=0.
#  else
      alpha=udelta/(udelta+uhat)
      fwd=alpha*fd+(1.-alpha)*fw
#  endif
!
! Wave Reynolds stress - VA2013 Eq 22
! (contribution associated with progressive surface waves)
!
#  ifdef BEDLOAD_VANDERA_WAVE_AVG_STRESS
      alpha_w=0.424
      tau_wRe=MAX((rho0*fwd*alpha_w*uhat**3./(2.*c_w)),eps)
#  else
      tau_wRe=0.
#  endif
!
! Compute the change in time period based on converged udelta
! (current velocity at wave boundary layer)
!
#  ifndef BEDLOAD_VANDERA_ZEROCURR
      call current_timeperiod(udelta, phi_curwave, umax, umin, RR,
     &                        T_c, T_t, Td)
#  endif
!
! Calculate the effect of surface waves
!
#  ifdef BEDLOAD_VANDERA_SURFACE_WAVE
      call surface_wave_mod(Td, depth, uhat, T_c, T_cu, T_t, T_tu)
#  endif
!
      return
      end
!
!======================================================================
!
      subroutine half_wave_cycle_stress_factors( T_iu, T_i, ahat, ksw,
     &                                                      fd, alpha,
     &                                                 alphac, alphaw,
     &                                                     d50, osmgd,
     &                              ui_r, uhat_i, udelta, phi_curwave,
     &                                                   tau_wRe, dsf,
     &                                theta_ix, theta_iy, mag_theta_i )
!
!  This subroutine returns the following:
!  dsf                 : sheetflow thickness
!  theta_ix, theta_iy  : Shields parameter in x and y dir.
!  mag_theta_i         : Magnitude of Shields parameter for half cycle
!
      implicit none
#  include "param.h"
#  include "scalars.h"
!
! Input the crest or trough half cycle velocity
! d50 -- grain size in meters
! Different for crest and trough half cycles
!
      real T_iu, T_i, ahat, ksw
      real fd, alpha
      real alphac, alphaw
      real d50, osmgd
      real ui_r, uhat_i, udelta, phi_curwave
      real tau_wRe
      real dsf, theta_ix, theta_iy, mag_theta_i
      real fwi_calc, dsf_calc
!
!  Local variables
!
      real eps
      parameter(eps = 1.0E-14)
      real fw_i, fwd_i
      real alpha_w, fwd, k, c_w
      real theta_hat_i
      real ui_rx, ui_ry, mag_ui
!
! Wave friction factor for wave and crest half cycle - VA2013 Eq 21
!
      fw_i=fwi_calc(T_iu, T_i, ahat, ksw)
!
! Wave current friction factor (Madsen and Grant) - VA2013 Eq 18
! Different for crest and trough
!
#  ifdef BEDLOAD_VANDERA_ZEROCURR
      fwd_i=fw_i
#  else
      fwd_i=alpha*fd+(1.-alpha)*fw_i
#  endif
!
! Magnitude of Shields parameter - VA2013 Eq 17
!
      theta_hat_i=0.5*fwd_i*uhat_i**2*osmgd
!
! Sheet flow thickness - VA2013 Appendix C C.1
!
      dsf=dsf_calc(d50, theta_hat_i) ! this dsf is in m
!
! Compute velocity magnitude based on representative velocities
!  - VA2013 Eq 12
!
! Get the representative trough half cycle water particle velocity
! as well as full cycle orbital velocity and excursion
!
#  ifdef BEDLOAD_VANDERA_ZEROCURR
      ui_rx=ui_r*alphaw
      ui_ry=0.
#  else
      ui_rx=udelta*COS(phi_curwave)*alphac+ui_r*alphaw
      ui_ry=udelta*SIN(phi_curwave)*alphac
#  endif
!
! mag_ui is set to a min value to avoid non-zero division
!
#  ifdef BEDLOAD_VANDERA_ZEROCURR
      mag_ui=MAX( SQRT(ui_rx*ui_rx), eps )
#  else
      mag_ui=MAX( SQRT(ui_rx*ui_rx+ui_ry*ui_ry), eps )
#  endif
!
! Magnitude of Shields parameter - VA2013 Eq 17
!
      mag_theta_i=MAX(0.5*fwd_i*osmgd*(mag_ui**2.), eps)
!
! Shields parameter in crest cycle
! rho0 is required for non-dimensionalizing
!
      theta_ix=ABS(mag_theta_i)*ui_rx/(mag_ui)+tau_wRe*osmgd/rho0
#  ifdef BEDLOAD_VANDERA_ZEROCURR
      theta_iy=0.
#  else
      theta_iy=ABS(mag_theta_i)*ui_ry/(mag_ui)
#  endif
!
! mag_theta_i is set to a min value to avoid non-zero division
!
#  ifdef BEDLOAD_VANDERA_ZEROCURR
      mag_theta_i=MAX( sqrt(theta_ix*theta_ix+theta_iy*theta_iy),eps )
#  else
      mag_theta_i=MAX( sqrt(theta_ix*theta_ix),eps )
#  endif
!
      return
      end
!
!======================================================================
!
      subroutine current_timeperiod(unet, phi_curwave, umax, umin,
     &                              RR, T_c, T_t, Td)
!
! This subroutine returns the following:
! T_c, T_t  : Time period in crest and trough cycle
!
! Modify the crest and trough time periods based on current velocity
! This function was developed by Chris Sherwood and Tarandeep Kalra
!
! The basis for these formulations are formed from Appendix A.3
! in SANTOSS report (Report nb: SANTOSS_UT_IR3 Date: January 2010)
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real unet, phi_curwave
      real umax, umin
      real RR, Td
      real T_c, T_t
!
!  Local variables
!
      real unet_xdir, udelta, delt
!
      unet_xdir=unet*cos(phi_curwave)
      if(RR.eq.0.5) then
        T_c=0.5*Td
        T_t=0.5*Td
        if(unet_xdir.ge.umax) then
          T_c=Td
          T_t=0.
        elseif(unet_xdir.le.umin) then
          T_c=0.
          T_t=Td
        elseif(unet_xdir.lt.0.0.and.unet_xdir.gt.umin) then
          delt=ASIN(-unet/umin)/pi
          T_t=T_t*(1.-2.*delt)
          T_c=Td-T_t
        elseif(unet_xdir.gt.0.0.and.unet_xdir.lt.umax) then
          delt=ASIN(unet_xdir/(-umax))/pi
          T_c=T_c*(1.-2.*delt)
          T_t=Td-T_c
        elseif(unet_xdir.eq.0.) then
          T_c=T_c
          T_t=T_t
        endif
      elseif(RR.gt.0.5) then
        T_c=T_c
        T_t=T_t
        if(unet_xdir.ge.umax) then
          T_c=Td
          T_t=0.
        elseif(unet_xdir<=umin) then
          T_c=0.
          T_t=Td
        elseif(unet_xdir.lt.0.0.and.unet_xdir.gt.umin) then
          delt=ASIN(-unet_xdir/umin)/pi
          T_t=T_t*(1.-2.*delt)
          T_c=Td-T_t
        elseif(unet_xdir.gt.0.0.and.unet_xdir.lt.umax) then
          delt=ASIN(unet_xdir/(-umax))/pi
          T_c=T_c*(1.-2.*delt)
          T_t=Td-T_c
        elseif(unet_xdir.eq.0.) then
          T_c=T_c
          T_t=T_t
        endif
      endif
      T_c=MAX(T_c,0.)
      T_t=MAX(T_t,0.)
!
      return
      end
!
!======================================================================
!
      subroutine surface_wave_mod(Td, depth, uhat, T_c, T_cu, T_t, T_tu)
!
! This subroutine returns the following:
! T_c, T_cu, T_t, T_tu  : Change in time period in crest and
!                         trough cycle due to particle displacement
!                         under surface waves.
!
! Crest period extension for horizontal particle displacement.
! Tuning factor eps_eff = 0.55 from measurements GWK Schretlen 2010
! Equations in Appendix B of SANTOSS Report (SANTOSS_UT_IR3, 2010)
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real Td, depth, uhat
      real T_c, T_cu, T_t, T_tu
!
!  Local variables
!
      real eps
      parameter(eps = 1.0E-14)
      real k_wn, eps_eff, c
      real delta_Tc, delta_Tt
      real T_c_new, T_cu_new
      real T_t_new, T_tu_new
      real kh_calc
!
      k_wn=kh_calc(Td,depth)/depth
      c=2.*pi/(k_wn*Td)
!
      eps_eff=0.55
!
      delta_Tc=eps_eff*uhat/(c*pi-eps_eff*2.*uhat)
      T_c_new=T_c+delta_Tc
!
! Avoid non zero values for time periods
!
      T_c_new=MAX( T_c_new, 0. )
      T_cu_new=MAX( T_cu*T_c_new/T_c, 0. )
!
      delta_Tt=eps_eff*uhat/(c*pi+eps_eff*2.*uhat)
      T_t_new=T_t-delta_Tt
      T_t_new=MAX( T_t_new, 0. )
      T_tu_new=MAX( T_tu*T_t_new/T_t, 0. )
!
      T_c=T_c_new
      T_cu=T_cu_new
      T_t=T_t_new
      T_tu=T_tu_new
!
      return
      end
!
!======================================================================
!
      subroutine ripple_dim(psi, d50, eta, rlen)
!
! This subroutine returns the following:
! eta, rlen : Ripple dimensions: (height and length)
!
! Calculate ripple dimensions of O'Donoghue et al. 2006
! based on VA2013 Appendix B
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real psi, d50
      real eta, rlen
      real d50_mm
      real m_eta, m_lambda, n_eta, n_lambda
      real eps
      parameter(eps=1.0E-14)
!
      d50_mm=1000.*d50
      if(d50_mm.lt.0.22) then
        m_eta=0.55
        m_lambda=0.73
      elseif(d50_mm.ge.0.22.and.d50_mm.lt.0.30) then
        m_eta=0.55+(0.45*(d50_mm-0.22)/(0.30-0.22))
        m_lambda=0.73+(0.27*(d50_mm-0.22)/(0.30-0.22))
      else
        m_eta=1.
        m_lambda=1.
      endif
!
! Smooth transition between ripple regime and bed sheet flow regime
!
      if(psi.le.190.) then
        n_eta=1.
      elseif(psi.gt.190..and.psi.lt.240.) then
        n_eta=0.5*(1.+cos(pi*(psi-190.)/(50.)))
      elseif(psi.ge.240.) then
        n_eta=0.
      endif
      n_lambda=n_eta
!
      eta=MAX(0.,m_eta*n_eta*(0.275-0.022*psi**0.42))
      rlen=MAX(eps,m_lambda*n_lambda*(1.97-0.44*psi**0.21))
!
      return
      end
!
!======================================================================
!
      subroutine skewness_params( H_s, T, depth, r, phi, Ur )
!
! Ruessink et al. (2012) provides equations for calculating skewness
! parameters. Uses Malarkey and Davies equations to get "bb" and "r",
! given input of H_s, T and depth
!
! r     - skewness/asymmetry parameter r=2b/(1+b^2)      [value]
! phi   - skewness/asymmetry parameter                   [value]
! Su    - umax/(umax-umin)                               [value]
! Au    - amax/(amax-amin)                               [value]
! alpha - tmax/pi                                        [value]
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real H_s, T, depth
      real Ur
      real r, phi
      real kh_calc
!
! Local variables
!
      real p1,p2,p3,p4,p5,p6
      parameter(p1=0.,p2=0.857,p3=-0.471,p4=0.297,p5=0.815,p6=0.672)
      real B, psi, bb
      real k_wn, cff
      real Su, Au
!
! Ruessink et al. (2012), Coastal Engineering 65:56-63.
!
! k is the local wave number computed with linear wave theory.
!
      k_wn=kh_calc(T,depth)/depth
!
! Ursell number defined as 3/(32*pi^2) * Hs*L^2/h^3
!               (Fourier coeff for NL terms)
!
      Ur=max(0.375*H_s/(k_wn*k_wn*depth**3.),1.e-12)
!
! Ruessink et al. (2012), Eq 9
!
      cff=EXP((p3-log10(Ur))/p4)
      B=p1+((p2-p1)/(1.+cff))
!
      psi=-90.*deg2rad*(1.-TANH(p5/Ur**p6))
!
#  ifdef BEDLOAD_VANDERA_ASYM_LIMITS
      B=MIN(B,0.8554) ! according to fig.2, max values w.r.t Ur=24
      psi=MAX(psi,-1.4233)
#  endif
!
! Markaley and Davies, equation provides bb which is "b" in paper
! Check from where CRS found these equations
!
      bb=sqrt(2.)*B/(sqrt(2.*B**2.+9.))
      r=2.*bb/(bb**2.+1.)
!
! Ruessink et al., 2012 under Equation 12.
!
      phi=-psi-0.5*pi
!
! Skewness Su and asymmetry Au
!
      Su=B*cos(psi)
      Au=B*sin(psi)
!
      return
      end
!
!======================================================================
!
      subroutine abreu_points( r, phi, Uw, T, T_c, T_t,
     &                            T_cu, T_tu, umax, umin, RR, beta )
!
!  Calculate umax, umin, and phases of asymmetrical wave orbital velocity
!
!  Use the asymmetry parameters from Ruessink et al, 2012
!  to get the umax, umin and phases of asymettrical wave
!  orbital velocity to be used by Van Der A.
!  T_c is duration of crest
!  T_cu Duration of accerating flow within crest half cycle
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real r, phi, Uw, T
      real T_c, T_t, T_cu, T_tu
      real umax, umin, RR, beta
!
! Local variables
!
      real b, c, ratio, tmt, tmc, tzd, tzu
      real omega, w, phi_new
      real P, F0, betar_0
!     real T_tu, T_cu, T_c, T_t
      real cff1, cff2, cff
      real Sk, Ak
!
      omega=2.*pi/T
!
      phi_new=-phi

! Malarkey and Davies (Under equation 16b)
      P=SQRT(1.-r*r)
!
! Malarkey and Davies (Under equation 16b)
!
      b=r/(1.+P)
!
! Appendix E of Malarkey and Davies
!
      c=b*SIN(phi_new)
!
      cff1=4.*c*(b*b-c*c)+(1.-b*b)*(1.+b*b-2.*c*c)
      cff2=(1.+b*b)**2.-4.*c*c
      ratio=cff1/cff2
!
! These "if" conditionals prevent ASIN to be between [-1,1] and prevent NaNs
! Not a problem in the MATLAB code
!
      if(ratio.gt.1.)then
        ratio=1.
      endif
      if(ratio.lt.-1.)then
        ratio=-1.
      endif
      tmc=ASIN(ratio)
!
      cff1=4.*c*(b*b-c*c)-(1.-b*b)*(1.+b*b-2.*c*c)
      cff2=(1.+b*b)**2.-4.*c*c
      ratio=cff1/cff2
      if(ratio.gt.1.)then
        ratio=1.
      endif
      if(ratio.lt.-1.)then
        ratio=-1.
      endif
      tmt=ASIN(ratio)
!
      if(tmc.lt.0.) then
        tmc=tmc+2.*pi
      endif
      if(tmt.lt.0.) then
        tmt=tmt+2.*pi
      endif
!
! Non dimensional umax and umin, under E5 in Malarkey and Davies
!
      umax=1.+c
      umin=umax-2.
!
!     Dimensionalize
!
      umax=umax*Uw
      umin=umin*Uw
!
! phase of zero upcrossing and downcrossing (radians)
!
      tzu=ASIN(b*SIN(phi_new))
      tzd=2.*ACOS(c)+tzu
!
! MD, equation 17
!
      RR=0.5*(1.+b*SIN(phi_new))
!
! MD, under equation 18
!
      if(r.le.0.5) then
        F0=1.0-0.27*(2.*r)**(2.1)
      else
        F0=0.59+0.14*(2.*r)**(-6.2)
      endif
!
! MD, Equation 15a,b
!
      if(r.ge.0..and.r.lt.0.5)then
        betar_0=0.5*(1.+r)
      elseif(r.gt.0.5.and.r.lt.1.)then
        cff1=4.*r*(1.+r)
        cff2=cff1+1.
        betar_0=cff1/cff2
      endif
!
! MD, Equation 18
!
      cff=SIN((0.5*pi-ABS(phi_new))*F0)/SIN(0.5*pi*F0)
      beta=0.5+(betar_0-0.5)*cff
!
! MD, Table 1, get asymmetry parameterization
! using GSSO (10a,b)
!
      cff=SQRT(2.*(1.+b*b)**3.)
      Sk=3.*SIN(phi_new)/cff
      Ak=-3.*COS(phi_new)/cff
!
! These are the dimensional fractions of wave periods needed by Van der A eqn.
!
      w=1./omega
      T_c=(tzd-tzu)*w
      T_t=T-T_c
      T_cu=(tmc-tzu)*w
      T_tu=(tmt-tzd)*w
!
      return
      end

!
!======================================================================
!
      FUNCTION kh_calc(Td,depth)
!
! Calculate wave number from Wave period and depth
!
!#  define KH_HUNT
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real kh_calc
      real Td, depth
      real omega, Khd
      real cff, x, y, t
#  ifdef KH_HUNT
      real K1, K2, K3, K4, K5, K6
      parameter (K1=0.6666666666, K2=0.3555555555,
     &           K3=0.1608465608, K4=0.0632098765,
     &           K5=0.0217540484, K6=0.0065407983)
#  endif
!
      omega=2.*pi/Td
!
#  ifdef KH_HUNT
!
! Hunt (1979) "Direct solutions of wave dispersion equation"
!
      Khd=depth*omega*omega/g
      kh_calc= SQRT( Khd*Khd+
     &         Khd/(1.+Khd*(K1+Khd*(K2+Khd*(K3+Khd*(K4+
     &         Khd*(K5+K6*Khd)))))) )
#  else
!
! Soulsby (2006) "Simplified calculation of wave orbital velocities"
!
      if(depth.lt.0.) then  ! avoid negative wavenumber
        x=0.
      else
        x=omega**2.*depth/g
      endif
!
      if(x.lt.1.) then
        y=SQRT(x)
      else
        y=x
      endif
!
      t=TANH(y)   ! Iteratively solving 3 times (Soulsby 1997)
      cff=(y*t-x)/(t+y*(1.-t*t))
      y=y-cff
!
      t=TANH(y)
      cff=(y*t-x)/(t+y*(1.-t*t))
      y=y-cff
!
      t=TANH(y)
      cff=(y*t-x)/(t+y*(1.-t*t))
      y=y-cff
!
      kh_calc=y
#  endif /* KH_HUNT */
!
      return
      end
!
!======================================================================
!
      FUNCTION w_s_calc(d50, rhos)
!
! Critical Shields parameter from Soulsby (1997).
! Dynamics of Marine Sands
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real w_s_calc
      real nu
      parameter(nu=1.36E-6)
      real d50, rhos
      real s, dstar
      real cff, cff1
!
      s=rhos/rho0
      dstar=(g*(s-1)/(nu*nu))**(1./3.)*d50
      cff=nu/d50
      cff1=10.36
      w_s_calc=cff*(sqrt(cff1*cff1+1.049*dstar**3.)-cff1)
!
      return
      end
!
!======================================================================
!
      FUNCTION what_calc(Hs, Td, depth, RR, zws)
!
! Second order Stokes theory to get vertical velocity of water particle
! at a given elevation based on santoss_core.m
!
      implicit none
      real what_calc
      real eps_inv
      parameter (eps_inv=1.0E14)
      real pi
      parameter (pi=3.14159265358979323846)
      real Hs, Td, depth, RR, zws
      real cff, worb1, worb2, worb
!
      worb1=pi*Hs*zws/(Td*depth)
      worb2=max(worb1*2.*(RR+RR-1.),1.e-14)
!
!  Using the SANTOSS model formulation
!
      cff=0.125
      worb=cff*worb1*SQRT(64.-(-worb1+
     &          SQRT(worb1**2.+32.*
     &               worb2**2.))**2./(worb2**2.))+
     &               worb2*SIN(2.*ACOS(cff*(-worb1+
     &          SQRT(worb1**2.+32.*worb2**2.))/worb2))
!
! Prevent worb from going to Infinity when worb2=0.0
!
      worb=MIN(worb, eps_inv)
      what_calc=worb
!
      return
      end
!
!======================================================================
!
      FUNCTION mu_calc(d50)
!
! Bed roughness factor based on grain size - VA2013 Appendix A
! required for current and wave related bed roughness.
!
      implicit none
      real mu_calc
      real d50, d50_mm
!
      d50_mm=d50*1000.
!
      if(d50_mm.le.0.15) then
        mu_calc=6.
      elseif(d50_mm.gt.0.15.and.d50_mm.lt.0.2) then
        mu_calc=6.-5.*((d50_mm-0.15)/(0.2-0.15))
      elseif(d50_mm.ge.0.2) then
        mu_calc=1.
      endif
!
      return
      end
!
!======================================================================
!
      FUNCTION ksd_calc(d50, d90, mu, theta_avg, eta, rlen)
!
! Current-related bed roughness - VA2013 Appendix A.1
!
      implicit none
      real ksd_calc
      real d50, d90, mu, theta_avg, eta, rlen
      real ripple_fac
!
      rlen=MAX(rlen,d50)
      ripple_fac=0.4*eta**2./rlen
      ksd_calc=MAX( 3.*d90, d50*(mu+6.*(theta_avg-1.)) )
     &         + ripple_fac
!
      return
      end
!
!======================================================================
!
      FUNCTION ksw_calc(d50, mu, theta_avg, eta, rlen)
!
! Wave related bed roughness - VA2013 Eq A.5
!
      implicit none
      real ksw_calc
      real d50, mu, theta_avg, eta, rlen
      real ripple_fac, ksw
!
      rlen=MAX(rlen,d50)
      ripple_fac=0.4*eta**2./rlen
      ksw_calc=MAX( d50, d50*(mu+6.*(theta_avg-1.)) )
     &         + ripple_fac
!
      return
      end
!
!======================================================================
!
      FUNCTION fw_calc(ahat, ksw)
!
! Full-cycle wave friction factor - VA2013 Eq A.4
!
      implicit none
      real fw_calc
      real ahat, ksw, ratio, fw
!
      ratio=ahat/ksw
      if(ratio.gt.1.587) then
        fw_calc=0.00251*EXP(5.21*(ratio)**(-0.19))
      else
        fw_calc=0.3
      endif
!
      return
      end
!
!======================================================================
!
      FUNCTION fd_calc_santoss(udelta, delta, ksd)
!
! Current friction factor Assuming logarithmic velocity profile.
! - VA2013 Eq 20
!
      implicit none
      real fd_calc_santoss
      real udelta, delta, ksd
      real min_udelta
      parameter(min_udelta=1.0E-4)
      real von_k
      parameter(von_k=0.41)

      if(udelta.lt.min_udelta) then
        fd_calc_santoss=0.
      else
        fd_calc_santoss=2.*(von_k/LOG(30.*delta/ksd))**2.
      endif

      return
      end
!
!======================================================================
!
      FUNCTION fd_calc_madsen(udelta, mag_bstrc)
!
! Current friction factor from current stresses.
!
      implicit none
#  include "param.h"
#  include "scalars.h"
!
      real fd_calc_madsen
      real eps
      parameter(eps=1.E-14)
      real min_udelta
      parameter(min_udelta=1.0E-4)
      real udelta, mag_bstrc
!
      if(udelta.lt.min_udelta) then
        fd_calc_madsen=0.
      else
        fd_calc_madsen=MAX((mag_bstrc/(0.5*udelta*udelta)),eps)
!        fd_calc_madsen=MIN(fd_calc_madsen,2.)
      endif
!
      return
      end
!
!======================================================================
!
      FUNCTION fwi_calc(T_iu, T_i, ahat, ksw)
!
! Wave friction factor for half cycles - VA2013 Eq 21
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real fwi_calc
      real eps
      parameter(eps = 1.0E-14)
      real T_iu, T_i, ahat, ksw
      real c1, ratio, fwi
      real cff1, cff2, cff3
!
      fwi_calc=0.3
!
      c1=2.6
      ratio=ahat/ksw
      if(ratio.gt.1.587) then
        cff1=MAX( (T_iu/T_i),0. )
        cff2=(2.*cff1)**c1
        cff3=cff2*ratio
!
! These "if" condition prevents arithematic overflow error
! when 0.0**-0.19
!
        if(cff3.le.0.) then
          fwi_calc=0.
        else
          fwi_calc=0.00251*EXP(5.21*(cff3)**(-0.19))
        end if
      end if
!
      return
      end
!
!======================================================================
!
      FUNCTION dsf_calc(d50, theta_i)
!
! Sheet flow thickness - VA2013 Appendix C.1.
!
      implicit none
      real dsf_calc
      real d50, theta_i
      real d50_mm
      real cff
!
      d50_mm=d50*1000.
      if(d50_mm.le.0.15)then
        cff=25.*theta_i
      elseif(d50_mm.gt.0.15.and.d50_mm.lt.0.2)then
        cff=25.-(12.*(d50_mm-0.15)/0.05)
      elseif(d50_mm.ge.0.2)then
        cff=13.*theta_i
      endif
      dsf_calc=MAX(d50*cff,d50)
!
      return
      end
!
!======================================================================
!
      FUNCTION theta_cr_calc(d50, rhos)
!
! Critical Shields parameter from Soulsby (1997).
!
      implicit none
#  include "param.h"
#  include "scalars.h"
      real theta_cr_calc
      real nu
      parameter(nu=1.36E-6)
      real d50, rhos
      real s, dstar
      real cff1, cff2
!
      s=rhos/rho0
      dstar=(g*(s-1)/(nu*nu))**(1./3.)*d50
      cff1=0.3/(1.+1.2*dstar)
      cff2=0.055*(1.-EXP(-0.02*dstar))
      theta_cr_calc=cff1+cff2
!
      return
      end
# endif /* BEDLOAD_VANDERA */
!======================================================================
!
# ifdef SED_FLOCS
      subroutine sed_flocmod_tile (Istr,Iend,Jstr,Jend)
        USE flocmod, ONLY : flocmod_main, flocmod_comp_g
      implicit none
#  include "param.h"
#  include "scalars.h"
#  include "ocean3d.h"
#  include "mixing.h"
#  include "bbl.h"
#  include "sediment.h"
#  include "forces.h"

      integer Istr,Iend,Jstr,Jend
      integer i, j, k
!
!  Variable declarations for floc model

      real Gval, diss, nueau
      real epsilon8, ustr2, effecz

      if (l_testcase) then
        CALL flocmod_comp_g(Gval, time)
      else
        epsilon8=epsilon(1.0)
      endif

      do j=Jstr,Jend
        do i=Istr,Iend
          do k=1,N
            if (.not. l_testcase) then

#  if defined FLOC_TURB_DISS && defined GLS_MIXING
!
! Dissipation from turbulence clossure
!
              if (k.eq.1) then
                diss = Eps_gls(i,j,k)
              elseif (k.eq.N) then
                diss = Eps_gls(i,j,k-1)
              else
                diss = 0.5*(Eps_gls(i,j,k-1)+Eps_gls(i,j,k))
              endif
#   ifdef FLOC_BBL_DISS
!
! Dissipation from wavecurrent bottom stress
! NOT READY FOR PRIME TIME: NEEDS VERTICAL DISTRIBUTION
! As first cut, use turbulence closure
!
              if (k.eq.1) then
#    ifdef BBL
                ustr2 =sqrt((bustrw(i,j)**2
     &                      +bvstrw(i,j)**2))
                effecz = (ustr2**0.5)*Pwave(i,j)*0.5/pi
                diss = MAX((ustr2**(1.5))/(vonKar*effecz),diss)
#    else
                ustr2 =sqrt((bustr(i,j)**2.0+bvstr(i,j)**2.0 ))
                diss = MAX((ustr2**(1.5))/(vonKar*
     &                       (z_r(i,j,1)-z_w(i,j,0))),diss)
#    endif
              endif
#   endif
#  else
!
! Constant tke dissipation
!
                diss = epsilon8
#  endif
                nueau = 1.5E-6
                Gval=sqrt(diss/nueau)
            endif ! not l_testcase

            CALL flocmod_main( dt,
     &            t(i,j,k,nnew,itrc_mud:itrc_mud+NMUD-1),
     &            Gval )

          enddo
        enddo
      enddo  !J_LOOP

      return
      end subroutine sed_flocmod_tile

# endif /* SED_FLOCS */
!
!====================================================================
!
!  UP5 interpolation for bedload fluxes
!
!=====================================================================
!
# ifdef BEDLOAD_UP5
      function flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua)
        implicit none
        REAL    :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2
        REAL    :: ua
        REAL    :: flux5, flux6

        flux6 = ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2)
     &       +(q_ip2+q_im3) )/60.0

        flux5= flux6-sign(1.,ua)*(
     &      (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0
        return
      end function flux5
# endif
#else
      subroutine sediment_empty
      end
#endif



