swancom2.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	sbot (ABRBOT, DEP2, ECOS, ESIN, KWAVE, SPCSIG, UBOT, UX2, UY2, IDCMIN, IDCMAX, ISSTOP, VARFR, FRCOEF, IG, IMATRA)
subroutine	frabre (HM, ETOT, QBLOC)
subroutine	ssurf (ETOT, HM, QB, SMEBRK, IMATRA, IMATDA, IDCMIN, IDCMAX, PLWBRK, ISSTOP, IG)
subroutine	swcap (SPCDIR, SPCSIG, KWAVE, IDCMIN, IDCMAX, ISSTOP, ETOT, IMATDA, IMATRA, PLWCAP, CGO, UFRIC, DEP2, DISSIP, DISIMP, IG)
subroutine	brkpar (BRCOEF, ECOS, ESIN, AC2, SPCSIG, DEP2, IG)

Function/Subroutine Documentation

brkpar()

subroutine brkpar	(	real, intent(out)	BRCOEF,
		real, dimension(mdc), intent(in)	ECOS,
		real, dimension(mdc), intent(in)	ESIN,
		real, dimension(mdc,msc,0:mt), intent(in)	AC2,
		real, dimension(msc), intent(in)	SPCSIG,
		real, dimension(mt), intent(in)	DEP2,
		integer, intent(in)	IG
	)

Definition at line 886 of file swancom2.f90.

  
! (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

frabre()

subroutine frabre	(	real	HM,
		real	ETOT,
		real	QBLOC
	)

Definition at line 225 of file swancom2.f90.

!
!****************************************************************
!
!     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

Here is the caller graph for this function:

sbot()

subroutine sbot	(	real	ABRBOT,
		real, dimension(mt)	DEP2,
		real, dimension(mdc)	ECOS,
		real, dimension(mdc)	ESIN,
		real, dimension(msc,icmax)	KWAVE,
		real, dimension(msc)	SPCSIG,
		real, dimension(mt)	UBOT,
		real, dimension(mt)	UX2,
		real, dimension(mt)	UY2,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		integer	ISSTOP,
		logical	VARFR,
		real, dimension(mt)	FRCOEF,
		integer	IG,
		real, dimension(mdc,msc)	IMATRA
	)

Definition at line 16 of file swancom2.f90.

!  (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

ssurf()

subroutine ssurf	(	real	ETOT,
		real	HM,
		real	QB,
		real	SMEBRK,
		real, dimension(mdc,msc)	IMATRA,
		real, dimension(mdc,msc)	IMATDA,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real, dimension(mdc,msc,nptst)	PLWBRK,
		integer	ISSTOP,
		integer	IG
	)

Definition at line 314 of file swancom2.f90.

!
!****************************************************************
!     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

swcap()

subroutine swcap	(	real, dimension(mdc,6), intent(in)	SPCDIR,
		real, dimension(msc), intent(in)	SPCSIG,
		real, dimension(msc,micmax), intent(in)	KWAVE,
		integer, dimension(msc), intent(in)	IDCMIN,
		integer, dimension(msc), intent(in)	IDCMAX,
		integer, intent(in)	ISSTOP,
		real, intent(in)	ETOT,
		real, dimension(mdc,msc), intent(inout)	IMATDA,
		real, dimension(mdc,msc), intent(inout)	IMATRA,
		real, dimension(mdc,msc,nptst), intent(out)	PLWCAP,
		real, dimension(msc,micmax), intent(in)	CGO,
		real, intent(in)	UFRIC,
		real, dimension(mt), intent(in)	DEP2,
		real, dimension(mdc,msc), intent(inout)	DISSIP,
		real, dimension(mdc,msc), intent(inout)	DISIMP,
		integer	IG
	)

Definition at line 482 of file swancom2.f90.

!  (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

Here is the call graph for this function: