!======================================================================
! 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
!======================================================================
!
#include "cppdefs.h"
#ifdef BBL
      subroutine bblm (tile)
!====================================================================!
!  Compute bottom stresses for combined waves & currents             !
!  using the parametric approximation by Soulsby 1997:               !
!  t_cw = t_c[1+1.2(t_w/(t_c+t_w))^3.2]                              !
!  in which                                                          !
!     t_cw = the combined wave-averaged stress (in current dir)      !
!     t_c  = stress due to currents if waves would be absent         !
!     t_w  = amplitude of stress due to waves without currents       !
!  and                                                               !
!     t_cw_max = SQRT([t_cw+t_w cos(phi_cw)]^2 + [t_w sin(phi_cw)]^2)!
!     in which                                                       !
!       t_cw_max ~ the maximum combined wave-averaged stress         !
!       phi_cw   = the angle between current and waves               !
!--------------------------------------------------------------------!
! References:
!
!  - Dyer 1986, Coastal & Estuarine Sediment Dynamics, Wiley, 342 pp.
!  - Harris & Wiberg 2001, Comp. & Geosci. 27, 675-690
!  - Li & Amos 2001, Comp. & Geosci. 27, 619-645
!  - Soulsby 1997, Dynamics of Marine Sands, Telford  Publ., 249 pp.
!  - Soulsby 1995, Bed shear-stresses due to combined waves and currents,
!    in: Stive et al: Advances in Coastal Morphodynamics, Wiley, 4.20-4.23
!  - Wiberg & Harris 1994, J. Geophys. Res. 99(C4), 775-789
!
! First Implementation: Meinte Blaas 2002
! Changes:
!   P. Marchesiello & R. Benshila 2013-2014: WKB wave model;
!        bottom stress stability limit; 2D application
!   P. Marchesiello 2020: updates for BEDLOAD_VANDERA
!        and new SANDBAR test case
!   P. Marchesiello 2021: clear separation between form and skin stress
! with sediment modeling only using skin stress (no bedload contribution);
! addition of wave-breaking turbulence effect in skin stress.
!
!====================================================================!
      implicit none
      integer  tile, trd, omp_get_thread_num
# include "param.h"
# include "private_scratch.h"
!
# include "compute_tile_bounds.h"
!
      trd=omp_get_thread_num()
      call bblm_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))
      return
      end

!********************************************************************!
      SUBROUTINE bblm_tile (Istr,Iend,Jstr,Jend,Ub,Zr,Ur,Vr,Umag,
     &                                            taucw,taucwmax)
!********************************************************************!
!  output: bu/vstr = tauc = t_cw
!     effective wave-averaged bottom stress
!     (after applying apparent roughness)
!     bu/vstrw = tauw = t_cw_max
!     maximum skin-frictional bottom stress
!     rheight = ripple height (m)
!     rlength = ripple length (m)
!
!  input: Awave, Pwave, Dwave, h, Zob, d50, u,v
!         NB! Dwave is assumed to be in Cartesian convention (CROCO):
!             the direction to where the vector points, measured
!             counterclockwise from the positive x-axis of the
!             model grid (must account for grid angle)
!********************************************************************!
!
      implicit none
# include "param.h"
# include "bbl.h"
# include "forces.h"
# include "grid.h"
# ifdef SOLVE3D
#  include "ocean3d.h"
# endif
# include "ocean2d.h"
# include "scalars.h"
# include "sediment.h"
# if defined M3FAST && defined BSTRESS_FAST
#  include "nbq.h"
# endif
# include "mixing.h"
# ifdef WKB_WWAVE
#  include "wkb_wwave.h"
# endif
# ifdef SANDBAR
#  undef  BBL_WAVE_SKEWNESS
#  define BBL_BREAKING_STIR
# endif
      logical SAND, SILT
      integer Iend, Istr, Jend, Jstr, i, ised, j, k
      integer ISTRm,ISTRp,IENDp,JSTRm,JSTRp,JENDp
      real    Ab, anglec, anglew, cff, cff1, cff2,
     &        Fwave, Kbh, Kbh2, Kdh, phic, phicw,
     &        d50, rhosed, rhow, smgd, osmgd, D, mu_w,
     &        tau_cb, tau_ex, tau_bf, tau_up, tau_en,
     &        tau_w, tau_ws, tau_c, tau_cs, tau_cw, tau_cws,
     &        thetw, twopi, Ucur, Vcur, Ubc, visk,
     &        znots, znot, znota, znot_bl, znot_rip,
     &        eps, zomin, zomax, depth, Ut
      parameter
     &       (Ubc   = 1.e-2,  ! minimum orbital velocity for waves
     &        eps   = 1.e-10, ! small number to avoid zero division
     &        zomin = 1.e-5,  ! minimum roughness
     &        zomax = 3.e-2)  ! maximum roughness
# ifdef Z0_RIP
      real    rhgt, rlen, rhmax, rhmin, rlmin, psi,
     &        rhbio, rlbio, rhbiomax, rhfac
      parameter (
     &        rhmin = 1.e-4,      ! minimum ripple height [m]
     &        rlmin = rhmin/0.12, ! minimum ripple length [m]
     &        rhbiomax = 6.e-3,   ! maximum biogenic ripple height [m]
     &        rhfac = 0.0164)     ! biogenic ripple factor
# endif
      real       taucw(PRIVATE_2D_SCRATCH_ARRAY),
     &        taucwmax(PRIVATE_2D_SCRATCH_ARRAY),
     &              Ub(PRIVATE_2D_SCRATCH_ARRAY),
     &            Umag(PRIVATE_2D_SCRATCH_ARRAY),
     &              Ur(PRIVATE_2D_SCRATCH_ARRAY),
     &              Vr(PRIVATE_2D_SCRATCH_ARRAY),
     &              Zr(PRIVATE_2D_SCRATCH_ARRAY)
      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)
      real scf1, scf2, scf3, scf4, scf5
      parameter (scf1 = 0.5 * 1.39, scf2 = 0.52,
     &           scf3 = 2.0 - scf2, scf4 = 1.2,
     &           scf5 = 3.2)
# ifdef BBL_WAVE_SKEWNESS
      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 Ursell, cff3,cff4,cff5,cff6,cff7
# endif
!
# include "compute_auxiliary_bounds.h"
!
# ifdef MASKING
#  define SWITCH *
# else
#  define SWITCH !
# endif
# ifdef WET_DRY
#  define SWITCH_WET *
# else
#  define SWITCH_WET !
# endif
!
      twopi=2.*pi
!
!-----------------------------------------------------------------------
!  Initalize stresses due to currents and waves.
!-----------------------------------------------------------------------
!
      do j=JstrV-1,Jend
        do i=IstrU-1,Iend
          taucw(i,j)=0.        ! [m^2/s^2]
          taucwmax(i,j)=0.
        enddo
      enddo
!
!-----------------------------------------------------------------------
!  Set currents above bed.
!-----------------------------------------------------------------------
!
      DO j=JstrV-1,Jend+1
        DO i=IstrU-1,Iend+1
# ifdef SOLVE3D
          Zr(i,j)=max(z_r(i,j,1)-z_w(i,j,0),Zob(i,j)+1.E-4)
#  if defined M3FAST && defined BSTRESS_FAST
          Ur(i,j)=2.*qdmu_nbq(i,j,1)/(Hz(i,j,1)+Hz(i-1,j,1))
          Vr(i,j)=2.*qdmv_nbq(i,j,1)/(Hz(i,j,1)+Hz(i,j-1,1))
#  else
          Ur(i,j)=u(i,j,1,nrhs)
          Vr(i,j)=v(i,j,1,nrhs)
#  endif
# else
          Zr(i,j)=0.5*(h(i,j)+zeta(i,j,kstp))
          Ur(i,j)=ubar(i,j,kstp)
          Vr(i,j)=vbar(i,j,kstp)
# endif
        ENDDO
      ENDDO
!
!=======================================================================
!  Compute bottom stresses: start main i,j loop
!=======================================================================
!
      DO j=JstrV-1,Jend
        DO i=IstrU-1,Iend
!
# ifdef Z0_RIP
          rhbio = 0.
          rlbio = 0.
          rlen = Lripple(i,j)
          rhgt = Hripple(i,j)
# endif
# ifdef SOLVE3D
          rhow = rho(i,j,1)+rho0
          visk = 1.3e-3/rhow         ! kinem. viscosity
# else
          visk = 1.3e-3/rho0
# endif
!
!---------------------------------------------------------------------
! Compute bed wave orbital velocity (m/s), excursion amplitude (m)
! from wind-induced waves, and angle between current and waves.
!---------------------------------------------------------------------
!
! Use Dean & Dalrymple 1991 6th-degree polynomial to approximate
! wavenumber on shoaling water.
!
          Fwave=twopi/Pwave(i,j)
          depth=h(i,j)+zeta(i,j,kstp)
          Kdh=depth*Fwave*Fwave/g
          Kbh2=Kdh*Kdh+Kdh/(1.+Kdh*(K1+Kdh*(K2+Kdh*(K3+Kdh*(K4+
     &                         Kdh*(K5+K6*Kdh))))))
          Kbh = SQRT(Kbh2)
!
!  Compute bed wave orbital velocity and excursion amplitude.
!
          Ab=Awave(i,j)/SINH(Kbh) + eps
# ifdef OW_COUPLING_FULL
          Ub(i,j) = ubr(i,j) SWITCH rmask(i,j)
     &                       SWITCH_WET rmask_wet(i,j)
# else
          Ub(i,j)=Fwave*Ab SWITCH rmask(i,j)
     &                     SWITCH_WET rmask_wet(i,j)
!
!  Correction by skewness factor
!  (Abreu et al., 2010; Malarkey and Davis, 2012)
!
#  ifdef BBL_WAVE_SKEWNESS
          Ursell=1.0607*Awave(i,j)/(depth*kbh2)
          cff1=exp((p3-log10(Ursell))/p4)
          cff2=p1+(p2-p1)/(1+cff1)
          cff3=1.4142*cff2/(sqrt(2.*cff2**2+9.))
          cff4=-0.5*pi*(1-tanh(p5/Ursell**p6))+0.5*pi
          cff5=2.*cff3/(cff3**2+1.)
          cff6=cff5/(1.+sqrt(1.-cff5**2))
          cff7=cff6*sin(cff4)
          Ub(i,j)=Ub(i,j)*0.9*(1+cff7)  SWITCH rmask(i,j)
     &                                  SWITCH_WET rmask_wet(i,j)
#  endif
# endif
!
! Stirring velocity from wave breaking turbulence
!  (Roelvink & Reniers, 1995)
!
# if defined MRL_WCI && defined BBL_BREAKING_STIR
#  if defined WKB_WWAVE && defined WAVE_ROLLER
          cff=wepr(i,j)+(1-wkb_roller)*wepb(i,j)
#  elif defined WAVE_ROLLER
          cff=wepr(i,j)
#  else
          cff=wepb(i,j)
#  endif
!          Ut=cff/(exp(depth/whrm(i,j))-1)
          Ut=cff/(exp(depth/max(whrm(i,j),1.e-3))-1)
     &                              SWITCH rmask(i,j)
     &                              SWITCH_WET rmask_wet(i,j)
          Ut=500.*Ut  ! * efficient factor
# else
          Ut=0.
# endif
!
!  Compute bottom current magnitude at RHO-points.
!
          Ucur=0.5*(Ur(i,j)+Ur(i+1,j))
          Vcur=0.5*(Vr(i,j)+Vr(i,j+1))
          Umag(i,j)=SQRT(Ucur*Ucur+Vcur*Vcur)
     &                             SWITCH rmask(i,j)
     &                             SWITCH_WET rmask_wet(i,j)
     &                             + eps
!
!  Compute angle between currents and waves (radians)
!
          if (Ucur .ne. 0.) then
            phic=ATAN2(Vcur,Ucur)
          else
            phic=0.5*pi*SIGN(1.,Vcur)
          endif
          phicw=Dwave(i,j)-phic  ! Dwave in CROCO grid convention
!
!---------------------------------------------------------------------
! Sediment parameters:
! Establish local median grain size for all calculations in this
! subroutine. Since most parameterizations have been derived ignoring
! multiple grain sizes, we apply this single d50 also in the case of
! mixed beds.
!---------------------------------------------------------------------
!
# if defined ANA_BSEDIM || defined SEDIMENT
          d50 = Ssize(i,j)             ! [m]
          tau_cb = taucb(i,j)          ! [m^2/s^2]
          rhosed = Sdens(i,j)/rhow     ! [ ] relative density
# else
          d50 = 0.16e-3       ! default values: fine sand
          tau_cb = 0.17/rhow
          rhosed = 2650./rhow
# endif /* ANA_BSEDIM || SEDIMENT */
!
          smgd=(rhosed-1.)*g*d50
          osmgd=1./smgd
!
!---------------------------------------------------------------------
! Determine grain roughness from sediment size
!---------------------------------------------------------------------
!
          znots = d50/12.
          znota = znots
!
!---------------------------------------------------------------------
! Determine critical stresses for ripple formation
!---------------------------------------------------------------------
!
# ifdef Z0_RIP
!
! Transition to sheet flow (Li & Amos, 2001)
          tau_up = 0.03*smgd*d50**(-0.376)                ! [m^2/s^2]
!
! Break off (Grant & Madsen,1982)
          tau_bf = 0.78*(d50/visk)**0.6*smgd**0.3*tau_cb  ! [m^2/s^2]
# endif
!
!---------------------------------------------------------------------
! If significant waves (Ub > Ubc= 0.01 m/s):
! Wave-current interaction case according to Soulsby 1995.
! Otherwise: tauw = tauc for sediment purposes
!---------------------------------------------------------------------
!
          if(Ub(i,j).gt.Ubc) then   !<======
!
!---------------------------------------------------------------------
! Determine skin shear stress for combined flow [m^2/s^2]
!
!  1- Soulsby wave skin stress tau_ws with friction factor fw
!    (tau_ws = 0.5*fw*Ub**2 with fw=1.39*(znots/Ab)**0.52)
!  2- Soulsby combined wave-current skin stress tau_cws
!  3- Soulsby Maximum of combined wave-current skin stress
!---------------------------------------------------------------------
!
! Pure wave skin stress
            tau_ws = scf1*((znots*Fwave)**scf2)*(Ub(i,j)+Ut)**scf3

! current-only skin stress
            cff1 = vonKar/LOG(Zr(i,j)/znots)
            tau_cs = cff1*cff1*Umag(i,j)*Umag(i,j)

!  combined wave-current skin stress
            tau_cws = tau_cs*(1.+scf4*((tau_ws/(tau_ws+tau_cs))**scf5))
         !   taucw(i,j) = tau_cws

!  Maximum of combined wave-current skin stress (for sediment model)
            taucwmax(i,j) = SQRT((tau_cws + tau_ws*COS(phicw))**2
     &                                    +(tau_ws*SIN(phicw))**2)
            tau_cw = taucwmax(i,j)   ! from RUTGERS
!
!---------------------------------------------------------------------
! Compute apparent bedload roughness
!---------------------------------------------------------------------
!
# ifdef Z0_BL
            tau_ex=max(tau_cws-tau_cb,0.)
            znot_bl=17.4*d50*(tau_ex*osmgd)**0.75
            znota = znots + znot_bl
!
!---------------------------------------------------------------------
! Determine bedform roughness ripple height (m) and ripple length (m)
! for sandy beds.
!---------------------------------------------------------------------
!
#  ifdef Z0_RIP
#   undef Z0_RIP_VANDERA
#   ifdef Z0_RIP_VANDERA
!  O'Donoghue et al. (2006) - Van der A et al. (2013)
!
            psi=osmgd*(1.27*Ub(i,j))**2          ! Max mobility number
            CALL rip_dim(psi, d50, rhgt, rlen)   ! Ripple dimensions
            rhgt=max(rhmin,rhgt*Ab)              ! height (m)
            rlen=max(rlmin,rlen*Ab)              ! length (m)
#   else
!  Li & Amos (2001) - Blaas et al. (2007)
!
            if (d50 .ge. 0.063e-3) then
              rhmax=0.25*rlen
              rhgt=max(min(rhmax,rhgt),rhmin)
              tau_en=tau_cws*max(1.,(rlen/(rlen-pi*rhgt))**2)
              if (tau_cws.lt.tau_cb .and. tau_en.ge.tau_cb) then
                rhgt = (19.6*SQRT(tau_cws/tau_cb)+20.9)*d50
                rlen = rhgt/0.12       ! local transport
              elseif (tau_cws.ge.tau_cb .and. tau_cw.lt.tau_bf) then
                if (tau_ws.ge.1.56*tau_cs) then
                  rhgt = (27.14*SQRT(tau_cw/tau_cb)+16.36)*d50
                  rlen = rhgt/0.15     ! wave bedload regime
                else
                  rhgt = (22.15*SQRT(tau_cw/tau_cb)+6.38)*d50
                  rlen = rhgt/0.12     ! curr bedload regime
                endif
              elseif (tau_cw.ge.tau_bf .and. tau_cw.lt.tau_up) then
                rlen = 535.*d50        ! break-off regime
                rhgt = 0.15*rlen*(SQRT(tau_up)-SQRT(tau_cw))/
     &                           (SQRT(tau_up)-SQRT(tau_bf ))
              elseif (tau_cw.ge.tau_up) then
                rlen = 0.              ! sheet flow (plane bed)
                rhgt = 0.
              else
                rhgt=Hripple(i,j)      ! tau_en < tau_cb: no transport,
                rlen=Lripple(i,j)      ! pre-existing conditions
              endif !tau_cws
            endif !d50
#   endif
#  endif /* Z0_RIP */
!
!---------------------------------------------------------------------
! Determine (biogenic) bedform roughness, ripple height (m)
! and ripple length (m) for silty beds, using Harris & Wiberg 2001.
!---------------------------------------------------------------------
!
#  ifdef Z0_BIO
            if (d50.lt.0.063e-3) then
              rlbio = 0.1 ! biogenic ripple length (Wheatcroft 1994)
              thetw = tau_cws*osmgd
              rhbio = (thetw**(-1.67))*rlbio*rhfac
              rhgt = min(rhbio,rhbiomax)
              rlen = rlbio
            endif
#  endif /* Z0_BIO */
!
!---------------------------------------------------------------------
! Total roughness znota consists of grain-scale roughness +
! bedload thickness + sandy ripple-enhanced roughness + silty
! biogenic ripple roughness.
!---------------------------------------------------------------------
!
#  if defined Z0_RIP || defined Z0_BIO
            znot_rip = 0.267*rhgt*rhgt/(max(rlen,rlmin)) ! Nielsen (1992)
            znota = znota + znot_rip
#  endif /* Z0_RIP || Z0_BIO */
!
!---------------------------------------------------------------------
! Limit total roughness between zomin or Zob and zomax
!---------------------------------------------------------------------
!
            znota = max(zomin,min(zomax,max(znota,Zob(i,j))))
!
!---------------------------------------------------------------------
! Compute bottom stress components based on total roughness [m/s]^2
! and store for use in computing flow drag and eddy diffusivity/viscosity
!---------------------------------------------------------------------
!
            cff1 = vonKar/LOG(Zr(i,j)/znota)
            cff2 = MIN(Cdb_max,MAX(Cdb_min,cff1*cff1))
            tau_c = cff2*Umag(i,j)*Umag(i,j)
            tau_w = scf1*((znota*Fwave)**scf2)*(Ub(i,j)**scf3)
            tau_cw = tau_c*(1.+scf4*((tau_w/(tau_w+tau_c))**scf5))
            taucw(i,j) = tau_cw
# endif /* Z0_BL */
!
          else ! (Ub < Ubc)  <================
!
!---------------------------------------------------------------------
! If current-only: tauw = tauc(skin) for use in sediment.F (ifdef BBL)
! tauc for current still depending on roughness due to current ripples
! (ifdef Z0_RIP)
!---------------------------------------------------------------------
!
# ifdef Z0_RIP
            if(tau_cs.gt.tau_up) then
              rhgt=0.
              rlen=0.
            else if(tau_cs.lt.tau_cb) then
              rhgt=Hripple(i,j)
              rlen=Lripple(i,j)
            else
              rlen=1000.*d50            ! Yalin (1964)
              rhgt=0.117*rlen**1.19     ! Allen (1970)
            endif
            znota = znots + 0.267*rhgt*rhgt/(max(rlen,rlmin))
# else
            znota = znots
# endif
            znota = max(zomin,min(zomax,max(znota,Zob(i,j))))
            cff1=vonKar/LOG(Zr(i,j)/znota)
            cff2=MIN(Cdb_max,MAX(Cdb_min,cff1*cff1))
            tau_c=cff2*Umag(i,j)*Umag(i,j)
            taucw(i,j) = tau_c

            cff1 = vonKar/LOG(Zr(i,j)/znots)
            tau_cs = cff1*cff1*Umag(i,j)*Umag(i,j)
            taucwmax(i,j) = tau_cs

          endif ! (Ub < Ubc)  <================
!
!---------------------------------------------------------------------
!  Load variables for output purposes.
!---------------------------------------------------------------------
!
          Abed(i,j)    = Ab            SWITCH rmask(i,j)
         Zbnot(i,j)    = znots         SWITCH rmask(i,j)
         Zbapp(i,j)    = znota         SWITCH rmask(i,j)
# ifdef Z0_RIP
         Hripple(i,j)  = rhgt          SWITCH rmask(i,j)
         Lripple(i,j)  = rlen          SWITCH rmask(i,j)
# endif
       enddo
      enddo  ! end main loop ==========
!
!=======================================================================
!  Finalize computation and store in global arrays
!=======================================================================
!
!  Compute kinematic bottom stress components (m2/s2) for flow
!  due to combined current and wind-induced waves.
!
      do j=Jstr,Jend
        do i=IstrU,Iend
          cff1=0.5*(taucw(i-1,j)+taucw(i,j))          ! at u point
          anglec=Ur(i,j)/(0.5*(Umag(i-1,j)+Umag(i,j)))
          bustr(i,j)=cff1*anglec SWITCH umask(i,j)
        enddo
      enddo

      do j=JstrV,Jend
        do i=Istr,Iend
          cff1=0.5*(taucw(i,j-1)+taucw(i,j))          ! at v point
          anglec=Vr(i,j)/(0.5*(Umag(i,j-1)+Umag(i,j)))
          bvstr(i,j)=cff1*anglec SWITCH vmask(i,j)
        enddo
      enddo
!
!  Compute form and skin stress components at rho points
!
      do j=Jstr,Jend
        do i=Istr,Iend
          Ucur=0.5*(Ur(i,j)+Ur(i+1,j))
          Vcur=0.5*(Vr(i,j)+Vr(i,j+1))
!
          anglew=cos(Dwave(i,j))
          bustrw(i,j)=taucwmax(i,j)*anglew ! total skin stress
!          Ubot(i,j)=Ub(i,j)*anglew
!          Ur(i,j)=Ucur
!
          anglew=sin(Dwave(i,j))
          bvstrw(i,j)=taucwmax(i,j)*anglew
!          Vbot(i,j)=Ub(i,j)*anglew
!          Vr(i,j)=Vcur
        enddo
      enddo

# ifdef LIMIT_BSTRESS
!---------------------------------------------------------------------
!  From J. Warner's code:
!  Set limiting factor for bottom stress. The bottom stress is adjusted
!  to not change the direction of momentum.  It only should slow down
!  to zero.  The value of 0.75 is arbitrary limitation assigment.
!---------------------------------------------------------------------
!
      cff=0.75/dt
      do j=Jstr,Jend
        do i=IstrU,Iend
          cff1=cff*0.5*(Zr(i-1,j)+Zr(i,j))
          bustr(i,j)=SIGN(1.0, bustr(i,j))*
     &               MIN(ABS(bustr(i,j)),
     &                   ABS(Ur(i,j))*cff1)
        enddo
      enddo
      do j=JstrV,Jend
        do i=Istr,Iend
          cff1=cff*0.5*(Zr(i,j-1)+Zr(i,j))
          bvstr(i,j)=SIGN(1.0, bvstr(i,j))*
     &               MIN(ABS(bvstr(i,j)),
     &                   ABS(Vr(i,j))*cff1)
        enddo
      enddo
# endif
!
!---------------------------------------------------------------------
! Set boundary conditions
!---------------------------------------------------------------------
!
# if defined M3FAST && defined BSTRESS_FAST
      if (LAST_FAST_STEP) then
# endif

# ifndef EW_PERIODIC
      IF (EASTERN_EDGE) THEN
        DO j=Jstr,Jend
          bustr(Iend+1,j)=bustr(Iend,j)
        END DO
        DO j=JstrV,Jend
          bvstr(Iend+1,j)=bvstr(Iend,j)
        END DO
        !
        DO j=Jstr,Jend
          bustrw(Iend+1,j)=bustrw(Iend,j)
          bvstrw(Iend+1,j)=bvstrw(Iend,j)
        END DO
      END IF
      IF (WESTERN_EDGE) THEN
        DO j=Jstr,Jend
          bustr(IstrU-1,j)=bustr(IstrU,j)
        END DO
        DO j=JstrV,Jend
          bvstr(Istr-1,j)=bvstr(Istr,j)
        END DO
        DO j=Jstr,Jend
          bustrw(Istr-1,j)=bustrw(Istr,j)
          bvstrw(Istr-1,j)=bvstrw(Istr,j)
        END DO
      END IF
# endif
# ifndef NS_PERIODIC
      IF (NORTHERN_EDGE) THEN
        DO i=IstrU,Iend
          bustr(i,Jend+1) =bustr(i,Jend)
        END DO
        DO i=Istr,Iend
          bvstr(i,Jend+1) =bvstr(i,Jend)
        END DO
        DO i=Istr,Iend
          bustrw(i,Jend+1) =bustrw(i,Jend)
          bvstrw(i,Jend+1) =bvstrw(i,Jend)
        END DO
      END IF
      IF (SOUTHERN_EDGE) THEN
        DO i=IstrU,Iend
          bustr(i,Jstr-1)=bustr(i,Jstr)
        END DO
        DO i=Istr,Iend
          bvstr(i,JstrV-1)=bvstr(i,JstrV)
        END DO
        DO i=Istr,Iend
          bustrw(i,Jstr-1)=bustrw(i,Jstr)
          bvstrw(i,Jstr-1)=bvstrw(i,Jstr)
        END DO
      END IF
# endif
# if !defined EW_PERIODIC && !defined NS_PERIODIC
      ISTRm=Istr-1
      ISTRp=Istr+1
      IENDp=Iend+1
      JSTRm=Jstr-1
      JSTRp=Jstr+1
      JENDp=Jend+1
      IF (SOUTHERN_EDGE.and.WESTERN_EDGE) THEN
        bustr(Istr,JSTRm)  =0.5*(bustr(ISTRp,JSTRm)+bustr(Istr,Jstr))
        bvstr(ISTRm,Jstr)  =0.5*(bvstr(Istr,Jstr)+bvstr(ISTRm,JSTRp))
        bustrw(ISTRm,JSTRm)=0.5*(bustrw(Istr,JSTRm)+bustrw(ISTRm,Jstr))
        bvstrw(ISTRm,JSTRm)=0.5*(bvstrw(Istr,JSTRm)+bvstrw(ISTRm,Jstr))
      ENDIF
      IF (SOUTHERN_EDGE.and.EASTERN_EDGE) THEN
        bustr(IENDp,JSTRm) =0.5*(bustr(IENDp,Jstr)+bustr(Iend,JSTRm))
        bvstr(IENDp,Jstr)  =0.5*(bvstr(IENDp,JSTRp)+bvstr(Iend,Jstr))
        bustrw(IENDp,JSTRm)=0.5*(bustrw(IENDp,Jstr)+bustrw(Iend,JSTRm))
        bvstrw(IENDp,JSTRm)=0.5*(bvstrw(IENDp,Jstr)+bvstrw(Iend,JSTRm))
      ENDIF
      IF (NORTHERN_EDGE.and.WESTERN_EDGE) THEN
        bustr(Istr,JENDp)  =0.5*(bustr(Istr,Jend)+bustr(ISTRp,JENDp))
        bvstr(ISTRm,JENDp) =0.5*(bvstr(ISTRm,Jend)+bvstr(Istr,JENDp))
        bustrw(ISTRm,JENDp)=0.5*(bustrw(ISTRm,Jend)+bustrw(Istr,JENDp))
        bvstrw(ISTRm,JENDp)=0.5*(bvstrw(ISTRm,Jend)+bvstrw(Istr,JENDp))
      ENDIF
      IF (NORTHERN_EDGE.and.EASTERN_EDGE) THEN
        bustr(IENDp,JENDp) =0.5*(bustr(IENDp,Jend)+bustr(Iend,JENDp))
        bvstr(IENDp,JENDp) =0.5*(bvstr(IENDp,Jend)+bvstr(Iend,JENDp))
        bustrw(IENDp,JENDp)=0.5*(bustrw(IENDp,Jend)+bustrw(Iend,JENDp))
        bvstrw(IENDp,JENDp)=0.5*(bvstrw(IENDp,Jend)+bvstrw(Iend,JENDp))
      ENDIF
# endif
# if defined BD_EXCHANGE
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend, Abed)
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend, Hripple)
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend, Lripple)
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend, Zbnot)
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend, Zbapp)
      call exchange_u2d_tile (Istr,Iend,Jstr,Jend, bustr)
      call exchange_v2d_tile (Istr,Iend,Jstr,Jend, bvstr)
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend, bustrw)
      call exchange_r2d_tile (Istr,Iend,Jstr,Jend, bvstrw)
# endif

# if defined M3FAST && defined BSTRESS_FAST
      endif
# endif

      return
      end
!
!======================================================================
!
      SUBROUTINE rip_dim(psi, d50, rhgt, rlen)
!
! This subroutine returns rhgt, rlen: ripple 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, rhgt, rlen
      real m_eta, m_lam, n_eta, n_lam
      real d50_mm, eps
      parameter(eps=1.e-14)
!
      d50_mm=0.001*d50
      IF (d50_mm .lt. 0.22) THEN
        m_eta=0.55
        m_lam=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_lam=0.73+(0.27*(d50_mm-0.22)/(0.30-0.22))
      ELSE
        m_eta=1.
        m_lam=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_lam=n_eta
!
      rhgt=MAX(0., m_eta*n_eta*(0.275-0.022*psi**0.42))
      rlen=MAX(eps,m_lam*n_lam*(1.970-0.440*psi**0.21))
!
      RETURN
      END
!
!======================================================================
!
#else
      subroutine bblm_empty
      end
#endif

