Modification of the bottom shear stress computation under waves+curre…

…nts in SED3D
schism-dev · Aug 5, 2020 · eaf9093 · eaf9093
1 parent 4e3c904
commit eaf9093
Show file tree

Hide file tree

Showing 4 changed files with 33 additions and 10 deletions.
diff --git a/src/Sediment/sed_bedload.F90 b/src/Sediment/sed_bedload.F90
@@ -46,6 +46,7 @@ SUBROUTINE sed_bedload_sd(ised,inea)
 !                         introduction of a limiter w_asym_max       !
 !                       > bed slope effects: now under option with   !
 !                         an external call to a new subroutine       !
+!          2020/07 - M.Pezerat,B.Mengual : tau_wc replaced by tau_m  !
 !                                                                    !
 !--------------------------------------------------------------------!
 
@@ -165,7 +166,7 @@ SUBROUTINE sed_bedload_sd(ised,inea)
       theta_w = 0.d0
 #endif
 
-      theta_m = tau_wc(inea)*osmgd ! component due to the waves+current 
+      theta_m = tau_m(inea)*osmgd ! component due to the waves+current 
       cff1 = theta_w*(1.0d0+w_asym)
       cff2 = theta_w*(1.0d0-w_asym)
 

diff --git a/src/Sediment/sed_friction.F90 b/src/Sediment/sed_friction.F90
@@ -108,13 +108,10 @@ SUBROUTINE sed_current_stress()
           z0s    = d50/12.0d0             ! Nikuradse roughness  (m)
 
           IF(hh>z0s) THEN
-          !IF(hh>Zob(i)) THEN
-            !IF(hh<=0.d0.OR.Zob(i)<=0.d0) THEN
             IF(hh<=0.OR.z0s<=0) THEN
               CALL parallel_abort('SED-current_stress: Cd failed')
             ENDIF
             cff1 = 1.0d0/LOG(hh/z0s)
-            !cff1 = 1.0d0/LOG(hh/Zob(i))
             cff2 = vonKar*vonKar*cff1*cff1
             wrk  = MIN(Cdb_max,MAX(Cdb_min,cff2))
           else
@@ -228,12 +225,17 @@ SUBROUTINE sed_wavecurrent_stress()
 !          2020/02 - X.Bertin, B.Mengual : fw based on Swart (1974)  !
 !                    using the skin friction                         !
 !                  - B.Mengual : tau_w limited by tau_max            !
+!          2020/07 - M.Pezerat, B.Mengual : tau_wc computed from     !
+!                    tau_m and tau_w, by including the angle between !
+!                    current and wave directions                     !
 !                                                                    !
 !--------------------------------------------------------------------!
 
       USE sed_mod,   ONLY: uorb,tp,Zob,bustr,bvstr,tau_c,      &
-     &                      tau_w,tau_wc,isd50,bottom,tau_max 
-      USE schism_glbl, ONLY: rkind,pi,nea,idry_e,errmsg,rho0
+     &                      tau_w,tau_wc,isd50,bottom,tau_max, &
+     &                      tau_m,ustress,vstress,wdir 
+      USE schism_glbl, ONLY: rkind,pi,nea,idry_e,errmsg,rho0, &
+     &                       elnode,i34,out_wwm
       USE schism_msgp, ONLY: myrank,parallel_abort,exchange_e2d
 
       IMPLICIT NONE
@@ -242,6 +244,8 @@ SUBROUTINE sed_wavecurrent_stress()
 !- Local variables --------------------------------------------------!
       INTEGER     :: i
       REAL(rkind) :: abskb,fw,ks,d50
+      ! Variables for tau_wc
+      REAL(rkind) :: phi,udir
 !- Start Statement --------------------------------------------------!
 
 !---------------------------------------------------------------------
@@ -308,20 +312,29 @@ SUBROUTINE sed_wavecurrent_stress()
         ENDIF
 
 !---------------------------------------------------------------------
-! - Computes mean wave-current bottom stress (Soulsby, 1997, eq. 69)
+! - Computes mean wave-current bottom stress, tau_m (Soulsby, 1997, eq. 69)
+! - MP,BM update: tau_wc is computed from tau_m and tau_w, by including 
+!   the angle between current and wave directions (Soulsby, 1997, eq. 70)
 !---------------------------------------------------------------------
 
+        ! Angle between current and wave directions (phi)                        
+        udir = ATAN2(vstress(i),ustress(i))
+        phi = wdir(i) - udir
+
         if(tau_c(i)+tau_w(i)==0) then
-          tau_wc(i)=0
+          tau_m(i) = 0
+          tau_wc(i) = 0
         else
-          tau_wc(i) = tau_c(i)*(1.0d0+1.2d0*(tau_w(i)/(tau_c(i)+tau_w(i)))**3.2d0)
+          tau_m(i)  = tau_c(i)*(1.0d0+1.2d0*(tau_w(i)/(tau_c(i)+tau_w(i)))**3.2d0)
+          tau_wc(i) = SQRT((tau_m(i) + tau_w(i)*ABS(COS(phi)))**2.0d0 + (tau_w(i)*SIN(phi))**2.0d0)
         endif
 
       ENDDO !i = 1,nea
 
 #else 
       ! No wave model
       DO i = 1,nea
+        tau_m(i) = tau_c(i)
         tau_wc(i)  = tau_c(i)
       ENDDO ! End loop nea
 #endif

diff --git a/src/Sediment/sed_init.F90 b/src/Sediment/sed_init.F90
@@ -65,6 +65,8 @@ SUBROUTINE sed_alloc()
 !                        porosity                                    !
 !                        > Different porosity initialization         !
 !                        depending on poro_option                    !
+!  2020/07 - M.Pezerat,B.Mengual : introduction of tau_m for         !
+!                        bottom shear stress computations            !
 !                                                                    !
 !--------------------------------------------------------------------!
 
@@ -132,6 +134,9 @@ SUBROUTINE sed_alloc()
         IF(i/=0) CALL parallel_abort('Sed: uorbp allocation failure')
       ALLOCATE(dirpeak(nea),stat=i) !Anouk
         IF(i/=0) CALL parallel_abort('Sed: dirpeak allocation failure')
+
+      ALLOCATE(tau_m(nea),stat=i)
+        IF(i/=0) CALL parallel_abort('Sed: tau_m allocation failure')
       ALLOCATE(tau_c(nea),stat=i)
         IF(i/=0) CALL parallel_abort('Sed: tau_c allocation failure')
       ALLOCATE(tau_w(nea),stat=i)
@@ -266,6 +271,7 @@ SUBROUTINE sed_alloc()
       wlpeak(:)    = IniVal
       uorb(:)      = IniVal
       uorbp(:)     = IniVal
+      tau_m(:)     = IniVal
       tau_c(:)     = IniVal
       tau_w(:)     = IniVal
       tau_wc(:)    = IniVal

diff --git a/src/Sediment/sed_mod.F90 b/src/Sediment/sed_mod.F90
@@ -57,6 +57,8 @@ MODULE sed_mod
 !                               > Bed compartment: top layer init.,  !
 !                                 porosity computations              !
 !                               > New outputs                        !
+!          2020/07 - M.Pezerat,B.Mengual : introduction of tau_m for !
+!                               bottom shear stress computations     !
 !                                                                    !
 !--------------------------------------------------------------------!
 !                                                                    !
@@ -253,9 +255,10 @@ MODULE sed_mod
        ! Variable shared by sed_friction subroutines
        REAL(rkind), ALLOCATABLE  :: bustr(:)       !Current-induced bottom stress in direction x (m^2/s/s)
        REAL(rkind), ALLOCATABLE  :: bvstr(:)       !Current-induced bottom stress in direction y (m^2/s/s)
+       REAL(rkind), ALLOCATABLE  :: tau_m(:)       !Wave-current mean bottom stress (m^2/s/s)
        REAL(rkind), ALLOCATABLE  :: tau_c(:)       !Current-induced bottom stress (m^2/s/s) =|(bustr,bvstr)|
        REAL(rkind), ALLOCATABLE  :: tau_w(:)       !Wave-induced bottom stres (m^2/s/s)
-       REAL(rkind), ALLOCATABLE  :: tau_wc(:)      !Wave-current mean bottom stress (m^2/s/s)
+       REAL(rkind), ALLOCATABLE  :: tau_wc(:)      !Wave-current max bottom stress (m^2/s/s)
 
        ! Bed variables defined at nodes
        REAL(rkind), ALLOCATABLE :: vc_area(:)      !Volume control area (vc_area(npa))