swancom2.f90

Go to the documentation of this file.

   SUBROUTINE sbot (ABRBOT  ,DEP2    ,ECOS    ,ESIN    ,KWAVE   , &
                    SPCSIG  ,UBOT    ,UX2     ,UY2     ,IDCMIN  , &
            IDCMAX  ,ISSTOP  ,VARFR   ,FRCOEF  ,IG      , &
            IMATRA                             )
!  (This subroutine has been tested only for case IBOT = 1)
!
!****************************************************************
!
!     Computation of the source terms due to bottom friction
!
!  Method
!
!     In SWAN several bottom dissipation models are computed, i.e.:
!
!     IBOT = 1   Jonswap bottom friction model
!     IBOT = 2   Collins bottom friction model
!     IBOT = 3   Madsen bottom friction model (see Tolman)
!
!     Both methods are implemented in SWAN and the user has to make
!     a choice in the input file.
!
!     1. Jonswap model:
!     -----------------
!
!     The bottom interaction term SEbf(s,d) is supposed to take the
!     Jonswap form (Hasselman et al. 1973):
!                         2
!                    sigma  E(s,d)
!     SEbf = -GAMMA ----------------
!                     2     2
!                    g  sinh  (kD)
!                                                                 2 -3
!     where GAMMA is the decay parameter ,(default GAMMA = 0.067 m s  ).
!     In the Jonswap form the current velocities are not taken into
!     account.
!
!     2. COLLINS model:
!     -----------------
!
!     The energy dissipation due to bottom friction is modelled
!     according the quadratic friction law:
!                    2
!      SE = Tau * |U|
!
!      which for a spectrum can be written as:
!                          2
!                     sigma
!      SE(s,d)= - ---------------- * (Cfw.Ub + Cfc.Uc) * E(s,d)
!                       2
!                 g sinh (K(s) * D)
!
!     Ub is the velocity due to the wave at the bottom
!
!     The current velocity is Uc
!
!     2. MADSEN formulation:
!     ----------------------
!
!     The bottom dissipation applying Madsen formulation is as
!     follows:
!
!                          fw [n - 1/2] UBR E(s,d)
!     [1]    Sdsb(s,d) = -  ------------------------
!                                      D
!
!     in which :
!                            2
!                           s * D
!     [1a]   (n - 1/2) = -------------
!                                2
!                        2 g sinh (kD)
!
!     UBOT(IX,IY) is computed in the subroutine SINTGRL. The friction
!     factor fw is estimated using the formulation of Jonsson (1963,
!     1966a):
!
!                1                1                        Ab,r
!     [2]     -------- + log  { ---------- } = mf + log  { ----- }
!            4 sqrt(fw)     10  4 sqrt(fw)             10   Kn
!
!     with:
!
!               2        //      1
!     [3]   Ab,r  = 2 * // -------------- E(s,d) ds dd
!                      //      2
!                          sinh (kD)
!
!     with: Ab,r is the representative near bottom excursion
!                amplitude
!           Kn   equivalent roughness
!           mf   constant ( mf = -0.08) (determined by Jonsson
!                                        and Carlssen 1976 )
!
!     [2] is only valid for Ab,r/Kn larger than approximately 1.
!     For smaller values a constant value of fw is used (fw = 0.3
!     for Ab,r/Kn < 1.57 )
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4  
!   USE ALL_VARS, ONLY : MT,AC2                                                      
   USE vars_wave, ONLY : mt,ac2                                                      

   IMPLICIT NONE
   
   REAL    SPCSIG(MSC)                                                 
   INTEGER  ID,IS,ISSTOP,J,IDDUM,IG
   REAL :: XDUM,KD,SBOTEO,CFBOT,FACB,CFW,FW,CURR,UC,ABRBOT,ADUM,CDUM,DDUM
   LOGICAL  VARFR
   REAL :: DEP2(MT),ECOS(MDC),ESIN(MDC),IMATDA(MDC,MSC),IMATRA(MDC,MSC), &
           KWAVE(MSC,ICMAX),PLBTFR(MDC,MSC,NPTST),UBOT(MT),  &
       UX2(MT),UY2(MT),DISSC1(MDC,MSC),FRCOEF(MT)  
   INTEGER :: IDCMIN(MSC),IDCMAX(MSC)
   REAL    :: AKN

   IF(ibot >= 1 .AND. dep2(ig) > 0.)THEN
     IF(ibot == 1)THEN
!
!      *** Jonswap model ***
!
!      PBOT(3) = GAMMA (a) in the Jonswap equation,
!
       cfbot = pbot(3) / grav_w**2
       
!       CFBOT = CFBOT*10.0               !JQI added
     ELSE IF (ibot == 2)THEN
!
!      *** Collins model ***
!
!      PBOT(2) = [cfw]
!
       IF(varfr)THEN                                                
         cfw = frcoef(ig)
       ELSE
         cfw = pbot(2)
       ENDIF
       cfbot = cfw * ubot(ig) / grav_w
     ELSE IF(ibot == 3)THEN
!
!      *** Madsen model ***
!
       IF(varfr)THEN                                                 
         akn = frcoef(ig)
       ELSE
         akn = pbot(5)
       ENDIF
!
!        *** PBOT(4) = Mf                      ***
!        *** AKN = PBOT(5) = [kn]  (roughness) ***
!
       IF((abrbot/akn) > 1.57)THEN
         xdum = pbot(4) + log10( abrbot / akn )
!
!        *** solving the implicit equation using a Newton ***
!        *** Rapshon iteration proces : a + log a = b     ***
!        *** the start value for ADUM = 0.3 because 0.3626 ***
!        *** is the minimum value of ADUM with b=-0.08.    ***
!
         adum = 0.3
         DO j = 1, 50
           cdum  = adum
           ddum  = ( adum + log10(adum) - xdum ) / ( 1.+ ( 1. / adum) )
           adum  = adum - ddum
           IF(abs(cdum - adum) < 1.e-4) GOTO 29
         END DO    
         WRITE(*,*) ' error in iteration fw: Madsen formulation'
29       CONTINUE
!                                                 1               1
!        *** computation of FW -->  A = ----- --> FW = -----
!                                              4 uFW          16 A**2
         fw = 1. / (16. * adum**2)
       ELSE
         fw = 0.3
       ENDIF
       cfbot =  ubot(ig) * fw / (sqrt(2.) * grav_w)
     ENDIF

     DO is = 1, isstop
       kd      = kwave(is,1) * dep2(ig)
       IF(kd < 10.)THEN
         facb = cfbot * (spcsig(is) / sinh(kd)) **2                   
!
         DO iddum = idcmin(is) , idcmax(is)
           id = mod( iddum - 1 + mdc , mdc ) + 1
!
           sboteo = facb                                               
           IF(ibot == 2 .AND. icur == 1 .AND. pbot(1) > 0.)THEN
!            additional dissipation due to current, seldom used
             curr = ux2(ig)*ecos(id) + uy2(ig)*esin(id)
             uc   = abs(curr)
!             PBOT(1) = [cfc]
             sboteo = facb + pbot(1) * uc * (spcsig(is) / sinh(kd)) **2 
           END IF
!
!          *** store the results in the array IMATDA             ***
!
           imatra(id,is) = imatra(id,is) - sboteo*ac2(id,is,ig)
           imatda(id,is) = imatda(id,is) + sboteo
!           IF (TESTFL) PLBTFR(ID,IS,IPTST) = -1.* SBOTEO               
!           DISSC1(ID,IS) = DISSC1(ID,IS) + SBOTEO
     END DO   !IDDUM  
       END IF
     END DO   !IS  
!
   ENDIF
!
   RETURN
   END SUBROUTINE sbot
 
!
!****************************************************************
!
   SUBROUTINE frabre ( HM, ETOT, QBLOC )                               
!
!****************************************************************
!
!     to compute the fraction of breaking waves in point ix,iy
!     of the computational grid
!
!  Method (updated...)
!
!      The fraction of breaking waves in a point ix,iy is given by
!      the implicit relation:
!
!        1 - Qb        ETOT
!        ------ = -8 * -----
!        ln Qb         HM**2
!
!        from which Qb can be found by solving the equation:
!
!                         ETOT
!        F = 1 - Qb + 8 * ----  * ln(Qb) = 0.
!                           2
!                         HM
!
!        The following appproximation is applied:
!
!                            2
!  (1)|   B = sqrt( 8 ETOT/HM ), i.e. B = Hrms/HM
!
!
!     |   Qo = 0.                                      B <= 0.5
!  (2)|                 2
!     |   Qo = ( 2B -1 )                         0.5 < B <= 1
!
!
!     applying the Newton-Raphson procedure (for 0.2<B<1.0):
!
!     |   Qb = 0.                                      B <= 0.2
!     |
!     |                                   2
!     |               2  Qo - exp((Qo-1)/B )
!  (3)|   Qb = Qo  - B   ------------------      0.2 < B <  1.0
!     |                   2               2
!     |                  B  - exp((Qo-1)/B )
!     |
!     |
!     |   Qb = 1.                                      B >= 1.0
!     |
!
   USE swcomm4                                                         
   USE ocpcomm4                                                        
!
   IMPLICIT NONE

   REAL    :: ETOT,  HM,  QBLOC
   INTEGER :: IENT
   REAL    :: B,  B2,  QO,  Z

   IF((hm > 0.) .AND. (etot >= 0.))THEN
     b = sqrt(8. * etot / (hm*hm) )
   ELSE
     b = 0.0
   END IF
!
   IF(b <= 0.5)THEN
     qo = 0.
   ELSE IF(b <= 1.0)THEN
     qo = (2.*b - 1.)**2
   END IF
!
   IF(b <= 0.2)THEN
     qbloc = 0.0
   ELSE IF(b < 1.0)THEN
!
!    *** second iteration to find Qb ***
!
     b2 = b*b
     z  = exp((qo-1.)/b2)
     qbloc = qo - b2 * (qo-z)/(b2-z)
   ELSE
     qbloc = 1.0
   END IF
!
   RETURN
   END SUBROUTINE frabre
!
!****************************************************************
!
   SUBROUTINE ssurf (ETOT    ,HM      ,QB      ,SMEBRK  ,       &
                     IMATRA  ,IMATDA  ,IDCMIN  ,IDCMAX  ,       &
             PLWBRK  ,ISSTOP  ,IG               )
!
!****************************************************************
!     Computation of the source term due to wave breaking.
!     White capping is not taken into account
!
!  Method
!
!     The source term for surf breaking is implemented following
!     the approach of Battjes/Janssen (1978) for the energy dissipation:
!
!             Alpha      -     2                  -   SMEBRK
!     Dtot =  ----  Qb * f * Hm              with f = ------
!              4                                      2 * Pi
!
!     Now the source term is:
!
!                      SIGMA * AC2(ID,IS,IX,IY)
!     Sbr =  - Dtot *  ------------------------  =
!                              Etot
!
!
!              Alpha * SMEBRK * Qb * Hm * Hm    SIGMA * AC2(ID,IS,IX,IY)
!         =  - ------------------------------ * -------------------------
!                       8 * Pi                            Etot
!
!
!         =  WS * SIGMA * AC2(ID,IS,IX,IY)   =  WS * E
!
!
!     with
!
!     Alpha = PSURF(1)                            ;
!
!                   SMEBRK Qb
!     WS    = Alpha ------ --                     ;
!                     Pi   BB
!                        2
!     BB    = 8 Etot / Hm  = - (1 - Qb) / ln (Qb) ;
!
!
!     The local maximum wave height Hm and mean frequency SMEBRK are computed
!     in subroutine SINTGRL.
!     The fraction of breaking waves Qb is calculated in the subroutine FRABRE
!
!     The new value for the dissipation is computed implicitly using
!     the last computed value for the action density Nold (at the spatial
!     gridpoint under consideration).
!
!     Sbr = WS * N
!
!         = Sbr_new + (d Sbr/d N) (Nnew - Nold)
!
!         = WS * Nnew + SbrD * (Nnew - Nold)
!
!         = (WS + SbrD)* Nnew - SbrD * Nold
!
!         = SURFA1 * Nnew - SURFA0 * Nold
!
!     In order to do this we need the derivative
!     of the source term Sbr to the action density N
!
!             d Sbr     d WS
!     SbrD =  -----  =  ---- * N + WS
!             d N       d N
!
!     Since BB and SMEBRK * N are proportional, we have
!
!     d Sbr     d WS                   SMEBRK  (d Qb/ d BB) *BB - Qb
!     -----  =  ---- * BB + WS = Alpha ------  --------------------- * BB + WS =
!     d N       d BB                     Pi           sqr(BB)
!
!
!           SMEBRK d Qb
!     Alpha ------ ----
!            Pi    d BB
!
!     With:
!
!     d Qb         1
!     ---- = -------------                 ;
!     d BB   (d BB / d Qb)
!
!                      2
!     d Qb           ln (Qb)
!     ---- = ---------------------------
!     d BB   ln (Qb) + (1 - Qb) (1 / Qb)
!
!            Qb (1 - Qb)
!          = ------------                  ;
!            BB (BB - Qb)
!
!     ------------------------------------------------------------
!     Get HM, QB and ETOT from the subroutine SINTGRL
!     For spectral direction IS and ID do,
!       get the mean energy frequency average over the full spectrum
!       If ETOT > 0 then
!         compute source term for energy dissipation SURFA0 and SURFA1
!       Else
!         source term for wave breaking is 0.
!       End if
!       ----------------------------------------------------------
!       Compute source terms for energy averaged frequency
!       Store results in the arrays IMATDA and IMATRA
!     ------------------------------------------------------------
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4 
!   USE ALL_VARS, ONLY : MT,AC2                                                       
   USE vars_wave, ONLY : mt,ac2                                                       
!
   IMPLICIT NONE

   INTEGER :: ISSTOP,IDCMIN(MSC),IDCMAX(MSC)
   REAL    :: DISSC0(MDC,MSC),DISSC1(MDC,MSC),     &
              IMATDA(MDC,MSC),IMATRA(MDC,MSC),PLWBRK(MDC,MSC,NPTST)
   REAL    :: ETOT,HM,QB,SMEBRK                         
   INTEGER :: ID,IDDUM,IENT,IS,IG
   DOUBLE PRECISION BB,DIS0,SbrD,SURFA0,SURFA1,WS
!
!     ALFA = PSURF(1)   <default = 1.0>
!
   bb = 8. * dble(etot) / ( dble(hm)**2 )                              
!   SURFA0 = 0.
!   SURFA1 = 0.
!   IF(REAL(BB) > 0. .AND. REAL(ABS(BB - DBLE(QB))) > 0.)THEN  
!     IF(BB < 1.)THEN
!       WS  = ( DBLE(PSURF(1)) / DBLE(PI_W)) * DBLE(QB) * DBLE(SMEBRK) / BB 
!       SbrD = WS * (1. - DBLE(QB)) / (BB - DBLE(QB))                   
!     ELSE
!       WS  = ( DBLE(PSURF(1)) / DBLE(PI_W)) * DBLE(SMEBRK)               
!       SbrD = 0.
!     END IF
!     SURFA0 = SbrD
!     SURFA1 = WS + SbrD
!   ELSE
!     SURFA0 = 0.
!     SURFA1 = 0.
!   END IF
   IF(real(bb) > 0. .AND. real(abs(bb - dble(qb))) > 0.)THEN
       ws  = ( dble(psurf(1)) / dble(pi_w)) * dble(qb) * dble(smebrk) / bb
   ELSE
       ws = 0.
   ENDIF

!
!  *** store the results for surf wave breaking  ***
!  *** in the matrices IMATDA and IMATRA         ***
!
   DO is = 1, isstop
     DO iddum = idcmin(is), idcmax(is)
       id = mod( iddum - 1 + mdc , mdc ) + 1
       imatda(id,is) = imatda(id,is) + real(ws)
       dis0 = ws * dble(ac2(id,is,ig))
       imatra(id,is) = imatra(id,is) - real(dis0)
     END DO
   END DO


   RETURN
   END SUBROUTINE ssurf

!
!****************************************************************
!
   SUBROUTINE swcap  (SPCDIR  ,SPCSIG  ,KWAVE   ,IDCMIN  ,      &
                      IDCMAX  ,ISSTOP  ,ETOT    ,IMATDA  ,      &
              IMATRA  ,PLWCAP  ,CGO     ,UFRIC   ,      &
              DEP2    ,DISSIP  ,DISIMP  ,IG      )              
!  (This subroutine has not been fully tested. Only tested for 
!   case IWCAP = 2)
!
!****************************************************************
!
!     Calculates the dissipation due to whitecapping
!
!  Method
!
!     Whitecapping dissipation is formulated as follows:
!
!     S_wc(sig,th) = - C_wc E(sig,th)
!
!     where the coefficient C_wc has four basic forms:
!
!     C_wc1 = C_K  sig~ (k/k~): According to Komen (generalised)
!     C_wc2 = C_BJ sig~ (k/k~): According to Battjes-Janssen (modified)
!     C_wc3 = C_LH            : According to Longuet Higgins (1969)
!     C_wc4 = C_CSM           : According to cumulative steepness
!                               method (Ris et al. (1999) and Donelan (1999))
!
!     In these formulations C_K is defined as (Komen; 1994 p. 145):
!
!                                       n1           n2
!     C_K = C1 [(1-delta) + delta (k/k~)  ] (S~/S_PM)
!
!     where C1, delta, n1 and n2 can be varied
!
!
!     C_BJ is defined as:
!
!                    2
!            alpha Hm Qb
!     C_BJ = -----------
!              8 Pi m0
!
!     where alpha can be varied.
!
!     for Hrms > Hm the formulation changes in a limit to (Hrms->Hm; Qb->1):
!
!            alpha
!     C_BJ = -----
!             Pi
!
!     and C_LH is defined as (Komen; 1994 p. 151-152):
!
!                            4   2                        2
!     C_LH = C3 Sqrt[(m0 sig0 )/g ] exp(A) sig0 (sig/sig0)
!
!     where
!                    2    2         4                2       2
!     A = -1/8 [1-eps ] [g /(m0 sig0 )]  with  [1-eps ] = [m2 ] / [m0 m4]
!
!     and C3 can be varied
!
!
!     Cumulative steepness method. Basic idea is that dissipation depends
!     on the steepness of the wave spectrum at and below a particular
!     frequency. Details can be found in "SWAN fysica plus" (2002),
!     report H3937/A832, implemented by Delf Hydraulics and Alkyon by
!     order of RIKZ/RWS as a part of the project HR-Ontwikkeling.
!                                                                         40.41
!     Recent modifications are described in:                              40.41
!     Hurdle, D.P., and G.Ph. van Vledder, 2004: Improved spectral wave   40.41
!     modelling of white-capping dissipation in swell sea systems. Proc.  40.41
!     Offshore Mechanics, Arctic Engineering Conference, June 2004,       40.41
!     Vancouver, Canada,  paper OMAE2004-51562.                           40.41
!                                                                         40.41
!     Default values: C_st = 4                                            40.41
!                     m    = 2                                            40.41
!                                                                         40.41
!     C_CSM = A*C_st S(sig,th)                                            40.41
!
!     where
!
!     C_st is a tuneable coefficient and
!
!     S(sig,th) = integrals from 0 to sig and from 0 to 2pi of
!
!                 k^2 |cos(th-th')|^m E(sig,th) dsig dth'
!
!     where coefficient m controls the directional dependence in case
!     of straining mechanism (if dominant then m = 2 else m > 10) and
!     th is the actual direction where straining mechanism takes
!     place.
!                                                                         40.41
!     A is a normalisation coefficient computed                           40.41
!                                                                         40.41
!              1     GAMMA(0.5*m+1.)                                      40.41
!       A = -------- ----------------                                     40.41
!           SQRT(PI) GAMMA(0.5*m+0.5)                                     40.41
!
!     In these equations the variables have the following meaning:
!
!     Hm   : Maximum wave height
!     Hrms : Root mean square of the wave heights
!     eps^2: Measure for the spectral bandwidth
!     m0   : Total wave energy density (=ETOT)
!     m2   : Second moment of the variance spectrum (=ETOT2)
!     m4   : Fourth moment of the variance spectrum (=ETOT4)
!     k    : Wave number (=KWAVE(IS,1))
!     k~   : Mean wave number
!     Qb   : Fraction of breaking waves
!     sig  : Frequency (=SPCSIG(IS))
!     sig0 : Average zero crossing frequency
!     sig~ : Mean frequency
!     S~   : Overall steepness (STP_OV)
!     S_PM : Overall steepness for a Pierson-Moskowitz spectrum
!     th   : direction theta (=SPCDIR(ID))
!
!     ------------------------------------------------------------
!     Calculate needed parameters
!     ------------------------------------------------------------
!     If IWCAP = 1, 2, or 5; Calculate C_K
!     If IWCAP = 4 or 5; Calcualte C_BJ
!     If IWCAP = 3; Calculate C_LH
!     If IWCAP = 6; Calculate cumulative steepness spectrum
!     IF IWCAP = 7; Calculate with Alves and Banner (2003)
!     ------------------------------------------------------------
!     For frequency dependent part of the spectrum
!       Calculate dissipation term due to whitecapping
!     ------------------------------------------------------------
!     For the whole frequency domain
!       Fill the matrices and PLWCAP-array
!     ------------------------------------------------------------
!     End of SWCAP
!     ------------------------------------------------------------
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4                                                        
   USE m_wcap
!   USE ALL_VARS, ONLY : MT,AC2
   USE vars_wave, ONLY : mt,ac2

   IMPLICIT NONE

   INTEGER, INTENT(IN)  :: ISSTOP, IDCMIN(MSC), IDCMAX(MSC)
   REAL, INTENT(IN)     :: DEP2(MT)
   REAL, INTENT(IN)     :: ETOT
   REAL, INTENT(IN)     :: KWAVE(MSC,MICMAX)                           
   REAL, INTENT(IN)     :: SPCDIR(MDC,6), SPCSIG(MSC)
   REAL, INTENT(OUT)    :: PLWCAP(MDC,MSC,NPTST)
   REAL, INTENT(IN OUT) :: IMATDA(MDC,MSC), IMATRA(MDC,MSC)
   REAL, INTENT(IN OUT) :: DISSIP(MDC,MSC), DISIMP(MDC,MSC)            
   REAL, INTENT(IN)     :: UFRIC                                       
   REAL, INTENT(IN)     :: CGO(MSC,MICMAX)                             

   INTEGER, SAVE     :: IENT = 0
   INTEGER           :: ID, IDDUM, IS, ID1, ID2, IF, IL, MXWCP,IG
   REAL              :: A, C_BJ, HM, HRMS, N1, N2
   REAL              :: QB_WC, SIG0, STP_OV, STP_PM
   REAL              :: DDIF, DELTA, DSTEEP, EBIN, XFAC          
   REAL              :: CPOW, CTOT, GAMMA                              
   REAL              :: BINSIZE                                      
   REAL, ALLOCATABLE :: C_K(:), C_LH(:), WCAP(:), WCIMPL(:)
   REAL, ALLOCATABLE :: CUMSTP(:,:), DSIGMA(:)                         
   REAL, ALLOCATABLE :: FCOS(:)                                        
   REAL              :: B, P                                           
   REAL              :: EF(MSC)                                        
   CHARACTER*20 NUMSTR, CHARS                                          
   CHARACTER*80 MSGSTR                                                 

   mxwcp = 7
   IF(iwcap > mxwcp)THEN                                            
!
! Error message
!
     chars = numstr(mxwcp+1,rnan,'(I1)')
     CALL txpbla(chars,IF,il)
     msgstr = 'Value for IWCAP should be less than '// chars(if:il)
     CALL msgerr ( 4, msgstr )
     RETURN
   END IF
!
! Initialisation
!
   IF(etot    <= 0.) RETURN
   IF(etot2   <= 0.) RETURN
   IF(etot4   <= 0.) RETURN
   IF(actot   <= 0.) RETURN
   IF(edrktot <= 0.) RETURN
!
   ALLOCATE (c_k(msc), c_lh(msc), wcap(msc), wcimpl(msc))
   ALLOCATE (cumstp(0:msc,mdc), dsigma(1:msc))                         
   ALLOCATE (fcos(1:mdc))                                              
   wcimpl(1:msc) = 0.
!
! Calculate coefficients
!
   IF((iwcap == 1) .OR. (iwcap == 2) .OR. (iwcap == 5))THEN                 
!
! Calculate C_K
!
     stp_ov = km_wam * sqrt(etot)
     stp_pm = sqrt(pwcap(2))
     n1     = pwcap(11)
     n2     = 2. * pwcap(9)
     c_k(:) = pwcap(1) * (1. - pwcap(10) +               &
              pwcap(10) * (kwave(:,1) / km_wam)**n1) *   &
          (stp_ov / stp_pm)**n2
!
   ENDIF
!
   IF((iwcap == 4) .OR. (iwcap == 5))THEN
!
! Calculate values for Hm and Qb
!
     hrms   = sqrt(8. * etot)
     IF (iwcap.EQ.4) hm = pwcap(6) / km01
     IF (iwcap.EQ.5) hm = pwcap(6) / (pwcap(8) * km_wam)
     CALL frabre(hm, etot, qb_wc)
!
! Calculate C_BJ
!
     IF(hrms >= hm)THEN
       c_bj = pwcap(7)  /  pi_w
     ELSE IF(hrms > 0.)THEN
       c_bj = (pwcap(7) *  hm**2 * qb_wc) / (pi_w * hrms**2)
     ELSE
       c_bj = 0.
     END IF
   ENDIF
!
   IF(iwcap == 3)THEN
!
! Calculate C_LH
!
     sig0 = sqrt(etot2 / etot)
!
!       A = -(1./8.)*(ETOT2**2/(ETOT*ETOT4))*(GRAV_W**2/(ETOT*SIG0**4))
!       rewrite to prevent underflow
!
     a = -(1./8.) * grav_w**2 / etot4
     DO is=1, isstop
!          C_LH(IS) = PWCAP(5) * SQRT((ETOT * SIG0**4) / GRAV_W**2) *
!     &               EXP(A) * SIG0 * (SPCSIG(IS) / SIG0)**2
!          rewrite to prevent underflow:
!
       c_lh(is) = pwcap(5) * exp(a) * sqrt(etot2) * spcsig(is)**2 / grav_w
     END DO
   END IF
!
! In case of cumulative steepness method, calculate its spectrum          
!
   IF(iwcap == 6)THEN                                              
     xfac      = (spcsig(3)-spcsig(1))/(2.*spcsig(2))                 
     dsigma(:) = xfac*spcsig(:)                                       
     delta     = spcdir(2,1)-spcdir(1,1)                              
!
!    --- precompute cos dependence                                    
!                                                                     
     DO id1 = 1, mdc                                                  
       ddif = real(id1-1)*delta                                       
       fcos(id1) = (abs(cos(ddif)))**pwcap(13)                        
     END DO                                                           
!
!    --- compute normalisation coefficient                            
!
     cpow = pwcap(13)                                                 
     IF(cpow > 10)THEN                                              
       ctot = sqrt(cpow/(2.*pi_w))*(1.+0.25/cpow)                       
     ELSE                                                             
       ctot = 1./sqrt(pi_w)*gamma(0.5*cpow+1.)/gamma(0.5*cpow+0.5)      
     END IF                                                           
!
     cumstp = 0.                                                      
     DO is = 1,isstop                                                 
       binsize = spcsig(is)*delta*dsigma(is)                         
       DO id1 = 1,mdc                                                
         cumstp(is,id1) = cumstp(is-1,id1)                          
         DO id2 = 1,mdc                                             
           ebin= ac2(id2,is,ig)*binsize                      
!           EBIN= AC2(ID2,IS,kcgrd(1))*BINSIZE                      
           dsteep = kwave(is,1)**2*ebin*fcos(abs(id1-id2)+1)       
           cumstp(is,id1) = cumstp(is,id1) + dsteep                
         END DO                                                     
       END DO                                                        
     END DO                                                           
     cumstp = pwcap(12)*cumstp                                        
     cumstp = cumstp*ctot                                             
   END IF                                                              
!
! Calculate dissipation according to Alves & Banner (2003)                
!
   IF(iwcap == 7)THEN                                              
!                                                                         
!  Loop to calculate B(k)                                                 
!                                                                         
     DO is = 1, isstop                                                 
!                                                                         
!  Calculate E(f)                                                         
!                                                                         
       ef(is) = 0.                                                     
       DO id = 1,mdc                                                   
         ef(is) = ef(is) + ac2(id,is,ig)*spcsig(is)*pi2_w*ddir     
       END DO                                                          
!                                                                      
!  Calculate saturation spectrum B(k) from E(f)                          
!                                                                        
       b = (1./pi2_w) * cgo(is,1) * kwave(is,1)**3 * ef(is)             
!                                                                        
!  Calculate exponent P of the relative saturation (B/Br)                
!                                                                        
       pwcap(10)= 3. + tanh(25.76*(ufric*kwave(is,1)/spcsig(is)-0.1)) 
       p = 0.5*pwcap(10)*(1. + tanh( 10.*( (b/pwcap(12))**0.5 - 1.))) 
!                                                                        
!  Calculate WCAP(IS) from B(k) and P                                    
!                                                                        
       stp_ov = km_wam * sqrt(etot)                                   
       wcap(is) = pwcap(1)*(b/pwcap(12))**(p/2.) *                     &
                  stp_ov**pwcap(9) * (kwave(is,1)/km_wam)**pwcap(11) * &
          (grav_w**(0.5)*kwave(is,1)**(0.5)/spcsig(is))**(pwcap(10)/2-1) * &
          grav_w**(0.5)*kwave(is,1)**(0.5)                                 
     END DO                                                           
   END IF                                                             
!
   IF(iwcap < 6)THEN                                             
!
! Calculate the whitecapping source term WCAP(IS)
!
     DO is=1, isstop
       IF((iwcap == 1) .OR. (iwcap == 2) .OR.                         &
          ((iwcap == 5) .AND. (c_bj <= c_k(is))))THEN
         wcap(is) = c_k(is) * sigm_10 * (kwave(is,1) / km_wam)
       ELSE IF(iwcap == 3)THEN
         wcap(is) = c_lh(is)
       ELSE IF((iwcap == 4) .OR. ((iwcap == 5) .AND. (c_bj >= c_k(is))))THEN
         IF(iwcap == 4) wcap(is) = c_bj*sigm01 *(kwave(is,1)/km01  )
         IF(iwcap == 5) wcap(is) = c_bj*sigm_10*(kwave(is,1)/km_wam)
!
! Calculate a term that is added to both sides of the equation to compensate
! for the strong non-linearity in the fraction of breaking waves Qb
!
         IF(hrms < hm)THEN
           wcimpl(is)=wcap(is) * ((1.-qb_wc)/((hrms**2/hm**2)-qb_wc))
           wcap(is)  =wcap(is) + wcimpl(is)
         END IF
       ELSE
         CALL msgerr(2,'Whitecapping is inactive')
         WRITE (printf,*) 'Occurs in gridpoint: ', ig
       END IF
     END DO

   END IF
!
! Fill the diagonal of the matrix and the PLWCAP-array
!
   IF(iwcap /= 6)THEN                                              

     DO is=1, isstop
!
!      Only fill the values for the current sweep
!
       DO iddum = idcmin(is), idcmax(is)
         id = mod(iddum - 1 + mdc, mdc) + 1
         imatda(id,is) = imatda(id,is) + wcap(is)
         imatra(id,is) = imatra(id,is) - wcap(is) * ac2(id,is,ig)
     
         dissip(id,is) = dissip(id,is) + wcap(is)                      
         IF (testfl) plwcap(id,is,iptst) = -1.*(wcap(is)-wcimpl(is))
       END DO
     END DO

   ELSE

     DO is=1, isstop                                                   
!                                                                         
!      Only fill the values for the current sweep                      
!                                                                         
       DO iddum = idcmin(is), idcmax(is)                               
         id = mod(iddum - 1 + mdc, mdc) + 1                            
         imatda(id,is) = imatda(id,is) + cumstp(is,id)                 
         dissip(id,is) = dissip(id,is) + cumstp(is,id)                 
         IF (testfl) plwcap(id,is,iptst) = -1.*cumstp(is,id)           
       END DO                                                          
     END DO                                                            

   END IF
!
! Add the implicit part to the right-hand side
!
   IF((iwcap == 4) .OR. (iwcap == 5))THEN
     DO is=1, isstop
!
!    Only fill the values for the current sweep
!
       DO iddum = idcmin(is), idcmax(is)
         id = mod(iddum - 1 + mdc, mdc) + 1
         imatra(id,is) = imatra(id,is) + wcimpl(is) * ac2(id,is,ig)
         disimp(id,is) = disimp(id,is) + wcimpl(is) * ac2(id,is,ig)      
       END DO
     END DO
   END IF
!
   DEALLOCATE (c_k, c_lh, wcap, wcimpl)
   DEALLOCATE (cumstp, dsigma, fcos)                                   
!
   RETURN
   END SUBROUTINE swcap

!****************************************************************
!
!!$  SUBROUTINE BRKPAR(BRCOEF,ECOS,ESIN,AC2,SPCSIG,DEP2,RDX,RDY,IG)     
  SUBROUTINE brkpar(BRCOEF,ECOS,ESIN,AC2,SPCSIG,DEP2,IG)     
  
! (This subroutine has not been fully tested yet)    
!
!****************************************************************
!
      USE swcomm3                                                         
      USE swcomm4                                                         
      USE ocpcomm4 
      USE mod_prec
      USE all_vars, ONLY : mt,ntsn,nbsn,vx,vy,art2 
!
      IMPLICIT NONE                                                       
!
!  2. Purpose
!
!     Determine the bottom slope in upwave direction and calculate
!     the slope dependent breaking parameter according to Nelson (1987)
!     Note that Nelson (1987) is used here since in Nelson (1994a,1994b)
!     an error is present in the equation.
!
!  3. Method
!
!     The breaker parameter is given by:
!
!     Hm / d =  0.55 + 0.88 exp ( -0.012 * cot beta)
!
!     with beta the angle the bed makes with the horizontal. This
!     above equation is only valid for positive slopes (Negative
!     slopes were not considered by Nelson. For very steep slopes
!     (>0.05 say) a very large breaker parameter is obtained (>>1).
!
!     To ensure wave breaking in laboratory cases (with very steep
!     slopes an upper limit of 0.81 (which corresponds to a bottom
!     slope of 0.01) is imposed on the model of Nelson.
!
!     For negative bottom slopes (not considered by Nelson) a value
!     op 0.73 is imposed (which is the average value in Table 2 of
!     Battjes and Janssen (1978).
!
      REAL, INTENT(OUT) :: BRCOEF    ! variable breaker coefficient       

      REAL, INTENT(IN)  :: SPCSIG(MSC)                                    
      REAL, INTENT(IN)  :: AC2(MDC,MSC,0:MT)  ! action densities         
      REAL, INTENT(IN)  :: ECOS(MDC), ESIN(MDC)  ! Cos and Sin of Theta   
      REAL, INTENT(IN)  :: DEP2(MT)         ! depths at grid points    

!!$      REAL, INTENT(IN)  :: RDX(10), RDY(10)
      INTEGER, INTENT(IN) :: IG                               

!     9. STRUCTURE
!
!     ------------------------------------------------------------
!     Calculate total energy per direction for all frequencies
!       determine action density per direction weighted with cos/sin
!       determine mean propagation direction of energy
!     ------------------------------------------------------------
!     calculate the depth derivative in the mean wave direction
!      according to dd/ds (see also subroutine SPROSD)
!     calculate the slope dependend breaking coefficient
!     ------------------------------------------------------------
!     End of NELSON
!     -------------------------------------------------------------
!
!************************************************************************
!
  INTEGER :: ID    ,IS      ! counters 
  INTEGER :: J,I1,I2
  REAL(SP) :: F1    
!
!
  REAL  :: ETOTS,EEX,EEY,EAD,SIGMA1,COSDIR,SINDIR,DDDX,DDDY,DDDS,DETOT                                
!
!     *** determine the average wave direction ***
!
  eex   = 0.
  eey   = 0.
  etots = 0.
  DO id = 1, mdc
    ead = 0.
    DO is = 1, msc
      sigma1 = spcsig(is)                                            
      detot  = sigma1**2 * ac2(id,is,kcgrd(1))
      ead    = ead + detot
    END DO
    etots = etots + ead
    eex   = eex + ead * ecos(id)
    eey   = eey + ead * esin(id)
  END DO
!
  IF(etots > 0.)THEN
    cosdir = eex / etots
    sindir = eey / etots
  ELSE
    cosdir = 1.
    sindir = 0.
  END IF
!
! *** Determine bottom slope in average wave propagation direction ***
!
!  DDDX =  RDX(1) * (DEP2(KCGRD(1)) - DEP2(KCGRD(2)))              &
!        + RDX(2) * (DEP2(KCGRD(1)) - DEP2(KCGRD(3)))
!  DDDY =  RDY(1) * (DEP2(KCGRD(1)) - DEP2(KCGRD(2)))              &
!        + RDY(2) * (DEP2(KCGRD(1)) - DEP2(KCGRD(3)))
  dddx = 0.0_sp
  dddy = 0.0_sp
  
  DO j=1,ntsn(ig)-1
    i1 = nbsn(ig,j)
    i2 = nbsn(ig,j+1)
    f1 = 0.50_sp*(dep2(i1)+dep2(i2)) 
    dddx = dddx + f1*(vy(i1)-vy(i2))
    dddy = dddy + f1*(vx(i2)-vx(i1))
  END DO
  dddx = dddx/art2(ig)
  dddy = dddy/art2(ig)

!
  ddds = -1. * ( dddx * cosdir + dddy * sindir )
!
! *** calculate breaking coefficient according to Nelson (1987) ***
!
  IF(ddds >= 0.)THEN
    ddds   = max( 1.e-6 , ddds)
    brcoef = psurf(4) + psurf(7) * exp( -psurf(8) / ddds )           
    brcoef = min( psurf(5) , brcoef )                                
  ELSE
    brcoef = psurf(6)                                                 
  END IF
!
!  PSURF(2) = BRKVAR                              deleted             
!  PSURF(2) is no longer used to transmit br. coefficient             

  RETURN
  END SUBROUTINE brkpar