Manually merging lrdev_updates: WWMIII parts first.

schism-dev · Apr 27, 2022 · 40baa99 · 40baa99
1 parent 3fabce6
commit 40baa99
Show file tree

Hide file tree

Showing 12 changed files with 1,070 additions and 1,057 deletions.
diff --git a/src/WWMIII/wwm_babanin.F90 b/src/WWMIII/wwm_babanin.F90
@@ -26,7 +26,7 @@ subroutine calc_Sds(ip,nfreq,EDENS,f,Kds,ANAR_IN,LPOW,MPOW,a1,a2,ACLOC,IMATRA,IM
     real(rkind)             , INTENT(INOUT) :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)
     real(rkind)             , INTENT(OUT) ::  Kds(:)     ! Kds(f)=Sds(f)/E(f)
 
-    real(rkind)             ::  Sds(nfreq)     ! Sds(f), the source term
+    real(rkind)             ::  Sds_loc(nfreq)     ! Sds(f), the source term
     real(rkind)             ::  ndedens(nfreq) ! NDEDENS(f)=DEDENS(f)/EDENST(f)
     real(rkind)             ::  dedens(nfreq)  ! DEDENS(f)=EDENS(f)-EDENST(f)
     real(rkind)             ::  edenst(nfreq)  ! EDENST(f)=(2*g.^2)*((2*pi).^-4)*(f.^-5).*(A.^-1)*Bnt
@@ -146,7 +146,7 @@ subroutine calc_Sds(ip,nfreq,EDENS,f,Kds,ANAR_IN,LPOW,MPOW,a1,a2,ACLOC,IMATRA,IM
 ! 3/9/2009 Since we have already non-dimensionalized by using NDEDENS instead of DEDENS, we do it backwards: Sds=Kds*Edens
 
     do is=1,nfreq
-       Sds(is)=Kds(is)*Edens(is) ! note that in most cases, the calling routine will not use Sds
+       Sds_loc(is)=Kds(is)*Edens(is) ! note that in most cases, the calling routine will not use Sds
        ST1(is)=T1(is)*Edens(is)
        ST2(is)=T2(is)*Edens(is)
     end do
@@ -165,7 +165,7 @@ subroutine calc_Sds(ip,nfreq,EDENS,f,Kds,ANAR_IN,LPOW,MPOW,a1,a2,ACLOC,IMATRA,IM
 
     CALL integrate(ST1_INT,f,ST1,nfreq)
     CALL integrate(ST2_INT,f,ST2,nfreq)
-    CALL integrate(Sds_int,f,Sds,nfreq)
+    CALL integrate(Sds_int,f,Sds_loc,nfreq)
 
 ! (point output write location 3)
     if(.false. .AND. WRITEDBGFLAG == 1)then

diff --git a/src/WWMIII/wwm_breaking.F90 b/src/WWMIII/wwm_breaking.F90
@@ -25,7 +25,7 @@ SUBROUTINE SDS_SWB(IP, SME, KME, ETOT, HS, ACLOC, IMATRA, IMATDA, SSBR, DSSBR)
       REAL(rkind) :: BIPH,Ur,S0
       REAL(rkind) :: SBRD, WS, SURFA0, SURFA1
       REAL(rkind) :: ERF_BETA
-      REAL(rkind) :: WAVEDIR, DPDW
+      REAL(rkind) :: WAVEDIR, DEG, TANBETA
 
       INTEGER :: ie,inne,icount
       INTEGER :: IS, ID
@@ -71,9 +71,10 @@ SUBROUTINE SDS_SWB(IP, SME, KME, ETOT, HS, ACLOC, IMATRA, IMATDA, SSBR, DSSBR)
          ! (e.g. Sallenger and Holman, 1985)
          ! lower and upper limit values prescribed in wwminput.nml
          CALL PEAK_PARAMETER(IP,ACLOC,MSC,FPP,TPP,CPP,WNPP,CGPP,KPP,LPP,PEAKDSPR,PEAKDM,DPEAK,TPPD,KPPD,CGPD,CPPD)
-         WAVEDIR = PEAKDM * PI / 180.d0 + PI / 2.d0 !conversion naut. to math. and in rad.
-         DPDW = tanbeta_x(IP)*COS(WAVEDIR) + tanbeta_y(IP)*SIN(WAVEDIR)
-         BRCRIT(IP) = a_BRCR*dpdw + b_BRCR
+         CALL DEG2NAUT(PEAKDM,DEG,LNAUTOUT) !conversion naut. to math.
+         WAVEDIR = DEG * PI / 180.d0        !in rad
+         TANBETA = tanbeta_x(IP)*COS(WAVEDIR) + tanbeta_y(IP)*SIN(WAVEDIR)
+         BRCRIT(IP) = a_BRCR*TANBETA + b_BRCR
          IF (BRCRIT(IP) < min_BRCR) BRCRIT(IP) = min_BRCR
          IF (BRCRIT(IP) > max_BRCR) BRCRIT(IP) = max_BRCR
 
@@ -118,7 +119,7 @@ SUBROUTINE SDS_SWB(IP, SME, KME, ETOT, HS, ACLOC, IMATRA, IMATDA, SSBR, DSSBR)
            BR_COEF = B_ALP**3D0
         END SELECT
       ELSE IF (BR_COEF_METHOD == 2) THEN
-        ! Adaptive with bottom slope (Pezerat et al., 2020 - under review)
+        ! Adaptive with bottom slope (Pezerat et al., 2021 - OM)
         BR_COEF = A_BR_COEF(IP)
       ELSE
         CALL WWM_ABORT('BR_COEF_METHOD HAS A WRONG VALUE')
@@ -270,9 +271,21 @@ SUBROUTINE SDS_SWB(IP, SME, KME, ETOT, HS, ACLOC, IMATRA, IMATDA, SSBR, DSSBR)
            BETA = ZERO
          END IF
          ! Implicit solver
-         ! TO DO
+         ! Source terms are linearized using a Newton-Raphson approach
          IF (ICOMP .GE. 2) THEN
-           CALL WWM_ABORT('IBREAK=4 not yet implemented in implicit')
+           IF ( BRCRIT(IP) .GT. 0D0 ) THEN 
+             IF ( BRCRIT(IP) .LT. 1D0 ) THEN
+               WS   = 0.75D0 * BR_COEF * SME / SQRTPI * BRCRIT(IP)**2.5 * BETA
+             ELSE
+               WS   = 0.75D0 * BR_COEF * SME / SQRTPI * BETA
+             END IF
+             SbrD = 7.D0*WS/2.
+             SURFA0 = SbrD
+             SURFA1 = WS + SbrD
+           ELSE
+             SURFA0 = 0D0
+             SURFA1 = 0D0
+           END IF
          ! Explicit Solver
          ELSE
            IF ( BRCRIT(IP) .GT. 0D0 ) THEN