Module for inclusion of triad nonlinear interaction according to Eldeberky's (1996) Lumped Triad Interaction (LTA) source term. More...

Functions/Subroutines
subroutine	w3str1 (A, AOLD, CG, WN, DEPTH, IX, S, D)
	Triad interaction source term computed using the Lumped Triad Appproximation (LTA) of Eldeberky (1996). More...

Detailed Description

Module for inclusion of triad nonlinear interaction according to Eldeberky's (1996) Lumped Triad Interaction (LTA) source term.

Author: A. J. van der Westhuysen

Date: 13-Jan-2013

Function/Subroutine Documentation

◆ w3str1()

subroutine w3str1md::w3str1	(	real, dimension(nspec), intent(in)	A,
		real, dimension(nspec), intent(in)	AOLD,
		real, dimension(nk), intent(in)	CG,
		real, dimension(nk), intent(in)	WN,
		real, intent(in)	DEPTH,
		integer, intent(in)	IX,
		real, dimension(nspec), intent(out)	S,
		real, dimension(nspec), intent(out)	D
	)

Triad interaction source term computed using the Lumped Triad Appproximation (LTA) of Eldeberky (1996).

     The parametrized biphase is given by:

                                  0.2
     beta = - pi/2 + pi/2 tanh ( ----- )
                                   Ur

     where Ur is the Ursell number.

     The source term as function of frequency p is:

             +      -
     S(p) = S(p) + S(p)

     in which

      +
     S(p) = alpha Cp Cg,p (R(p/2,p/2))**2 sin (|beta|) ( E(p/2)**2 -2 E(p) E(p/2) )

      -          +
     S(p) = - 2 S(2p)

     with alpha a tunable coefficient and R(p/2,p/2) is the interaction
     coefficient of which the expression can be found in Eldeberky (1996).

     Note that a slightly adapted formulation of the LTA is used in
     in the SWAN model:

     - Only positive contributions to higher harmonics are considered
       here (no energy is transferred to lower harmonics).

     - The mean frequency in the expression of the Ursell number
       is calculated according to the first order moment over the
       zeroth order moment (personal communication, Y.Eldeberky, 1997).

     - The interactions are calculated up to 2.5 times the mean
       frequency only.

     - Since the spectral grid is logarithmically distributed in frequency
       space, the interactions between central bin and interacting bin
       are interpolated such that the distance between these bins is
       factor 2 (nearly).

     - The interactions are calculated in terms of energy density
       instead of action density. So the action density spectrum
       is firstly converted to the energy density grid, then the
       interactions are calculated and then the spectrum is converted
       to the action density spectrum back.

Parameters

[in]	A	Action density spectrum (1-D)
[in]	CG	Group velocities.
[in]	WN	Wavenumbers.
[in]	DEPTH	Mean water depth.
[in]	IX
[out]	S	Source term (1-D version).
[out]	D	Diagonal term of derivative (1-D version).

Author: A. J. van der Westhuysen

Date: 13-Jan-2013

Definition at line 184 of file w3str1md.F90.

     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III           NOAA/NCEP |
     !/                  |     A. J. van der Westhuysen      |
     !/                  |     A. Roland                     |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         13-Jan-2013 |
     !/                  +-----------------------------------+
     !/
     !/    13 Jan-2013 : Origination, based on SWAN v40.91 code ( version 4.08 )
     !/    05 Oct-2016 : Avoiding divide by zero for EMEAN      ( version 5.15 )
     !/    28 Feb-2023 : Improvement of efficiency and stability ( version 7.xx)
     !/
     !  1. Purpose :
     !
     !     Triad interaction source term computed using the Lumped Triad
     !     Appproximation (LTA) of Eldeberky (1996).
     !
     !  2. Method :
     !
     !     (Taken from SWAN v40.91, based on code by Marcel Zijlema, TU Delft)
     !
     !     The parametrized biphase is given by:
     !
     !                                  0.2
     !     beta = - pi/2 + pi/2 tanh ( ----- )
     !                                   Ur
     !
     !     where Ur is the Ursell number.
     !
     !     The source term as function of frequency p is:
     !
     !             +      -
     !     S(p) = S(p) + S(p)
     !
     !     in which
     !
     !      +
     !     S(p) = alpha Cp Cg,p (R(p/2,p/2))**2 sin (|beta|) ( E(p/2)**2 -2 E(p) E(p/2) )
     !
     !      -          +
     !     S(p) = - 2 S(2p)
     !
     !     with alpha a tunable coefficient and R(p/2,p/2) is the interaction
     !     coefficient of which the expression can be found in Eldeberky (1996).
     !
     !     Note that a slightly adapted formulation of the LTA is used in
     !     in the SWAN model:
     !
     !     - Only positive contributions to higher harmonics are considered
     !       here (no energy is transferred to lower harmonics).
     !
     !     - The mean frequency in the expression of the Ursell number
     !       is calculated according to the first order moment over the
     !       zeroth order moment (personal communication, Y.Eldeberky, 1997).
     !
     !     - The interactions are calculated up to 2.5 times the mean
     !       frequency only.
     !
     !     - Since the spectral grid is logarithmically distributed in frequency
     !       space, the interactions between central bin and interacting bin
     !       are interpolated such that the distance between these bins is
     !       factor 2 (nearly).
     !
     !     - The interactions are calculated in terms of energy density
     !       instead of action density. So the action density spectrum
     !       is firstly converted to the energy density grid, then the
     !       interactions are calculated and then the spectrum is converted
     !       to the action density spectrum back.
     !
     !  3. Parameters :
     !
     !     Parameter list
     !     ----------------------------------------------------------------
     !       A       R.A.  I   Action density spectrum (1-D)
     !       CG      R.A.  I   Group velocities.
     !       WN      R.A.  I   Wavenumbers.
     !       DEPTH   Real  I   Mean water depth.
     !       EMEAN   Real  I   Mean wave energy.
     !       FMEAN   Real  I   Mean wave frequency.
     !       S       R.A.  O   Source term (1-D version).
     !       D       R.A.  O   Diagonal term of derivative (1-D version).
     !     ----------------------------------------------------------------
     !
     !  4. Subroutines used :
     !
     !      Name      Type  Module   Description
     !     ----------------------------------------------------------------
     !      STRACE    Subr. W3SERVMD Subroutine tracing.
     !     ----------------------------------------------------------------
     !
     !  5. Called by :
     !
     !      Name      Type  Module   Description
     !     ----------------------------------------------------------------
     !      W3SRCE    Subr. W3SRCEMD Source term integration.
     !      W3EXPO    Subr.   N/A    Point output post-processor.
     !      GXEXPO    Subr.   N/A    GrADS point output post-processor.
     !     ----------------------------------------------------------------
     !
     !  6. Error messages :
     !
     !       None.
     !
     !  7. Remarks :
     !
     !  8. Structure :
     !
     !     Determine resonance condition and the maximum discrete freq.
     !     for which the interactions are calculated.
     !
     !     If Ursell number larger than prescribed value compute interactions
     !        Calculate biphase
     !        Do for each direction
     !           Convert action density to energy density
     !           Do for all frequencies
     !             Calculate interaction coefficient and interaction factor
     !             Compute interactions and store results in matrix
     !
     !  9. Switches :
     !
     !     !/S  Enable subroutine tracing.
     !
     ! 10. Source code :
     !
     !/ ------------------------------------------------------------------- /
     USE constants, ONLY: grav, pi, tpi
     USE w3gdatmd, ONLY: nk, nth, nspec, dth, sig, dden, fte, ftf
     USE w3odatmd, ONLY: ndse
     USE w3servmd, ONLY: extcde
 #ifdef W3_S
     USE w3servmd, ONLY: strace
 #endif
     !/
     IMPLICIT NONE
     !/
     !/ ------------------------------------------------------------------- /
     !/ Parameter list
     !/
     REAL, INTENT(IN)        :: CG(NK), WN(NK), DEPTH, A(NSPEC), AOLD(NSPEC) 
     INTEGER, INTENT(IN)     :: IX
     REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC)
     !/
     !/ ------------------------------------------------------------------- /
     !/ Local parameters
     !/
     !     AUX1  :     auxiliary real
     !     AUX2  :     auxiliary real
     !     BIPH  :     parameterized biphase of the spectrum
     !     C0    :     phase velocity at central bin
     !     CM    :     phase velocity at interacting bin
     !     DEP   :     water depth
     !     DEP_2 :     water depth to power 2
     !     DEP_3 :     water depth to power 3
     !     E     :     energy density as function of frequency
     !     E0    :     energy density at central bin
     !     EM    :     energy density at interacting bin
     !     HS    :     significant wave height
     !     FT    :     auxiliary real indicating multiplication factor
     !                 for triad contribution
     !     I1    :     auxiliary integer
     !     I2    :     auxiliary integer
     !     ID    :     counter
     !     IDDUM :     loop counter in direction space
     !     IENT  :     number of entries
     !     II    :     loop counter
     !     IS    :     loop counter in frequency space
     !     ISM   :     negative range for IS
     !     ISM1  :     negative range for IS
     !     ISMAX :     maximum of the counter in frequency space for
     !                 which the triad interactions are calculated (cut-off)
     !     ISP   :     positive range for IS
     !     ISP1  :     positive range for IS
     !     RINT  :     interaction coefficient
     !     SA    :     interaction contribution of triad
     !     SIGPICG :   sigma times 2pi/Cg (a Jacobian for E(f) -> E(k))
     !     SINBPH:     absolute sine of biphase
     !     STRI  :     total triad contribution
     !     WISM  :     interpolation weight factor corresponding to lower harmonic
     !     WISM1 :     interpolation weight factor corresponding to lower harmonic
     !     WISP  :     interpolation weight factor corresponding to higher harmonic
     !     WISP1 :     interpolation weight factor corresponding to higher harmonic
     !     W0    :     radian frequency of central bin
     !     WM    :     radian frequency of interacting bin
     !     WN0   :     wave number at central bin
     !     WNM   :     wave number at interacting bin
     !     XIS   :     rate between two succeeding frequency counters
     !     XISLN :     log of XIS
     !
 #ifdef W3_S
     INTEGER, SAVE :: IENT = 0
 #endif
     INTEGER :: I1, I2, ID, IDDUM, II, IS, ISM, ISM1, ISMAX
     INTEGER :: ISP, ISP1, ITH, IK
     REAL    :: AUX1, AUX2, BIPH, C0, CM, DEP, DEP_2, DEP_3, E0, EM, HS
     REAL    :: FT, RINT, SIGPICG, SINBPH, STRI, WISM, WISM1, WISP
     REAL    :: WISP1, W0, WM, WN0, WNM, XIS, XISLN, EDM, ED0, G9DEP, STRI2
     REAL    :: E(NK), SA(NTH,100), SA2(NTH,100), A2(NSPEC), A3(NSPEC), HMAX
     REAL    :: EB(NK), EBAND, EMEAN, SIGM01, ED(NK)
 !----- Temp (to be moved) -----
     REAL    :: EF(NK), JACEPS, DIFFSTR
     REAL    :: PTRIAD(5)
     REAL    :: URSELL, ALPHAR
 !------------------------------
 !/
 !/ ------------------------------------------------------------------- /
 !/
 #ifdef W3_S
     CALL strace (ient, 'W3STR1')
 #endif
  
 !AR: todo: check all PRX routines for differences, check original thesis of elderberky. 
 !
 ! 1.  Integral over directions
 !
     sigm01 = 0.
     emean  = 0.
     jaceps = 1e-12
  
     hmax   = depth * 0.73
  
     DO ik=1, nk
       eb(ik) = 0.
       ed(ik) = 0.
       DO ith=1, nth
         eb(ik) = eb(ik) + a(ith+(ik-1)*nth)
         ed(ik) = ed(ik) + a(ith+(ik-1)*nth) * dden(ik) / cg(ik)
       END DO
     END DO
 !
 ! 2.  Integrate over frequencies. 
 !
     DO ik=1, nk
       eb(ik) = eb(ik) * dden(ik) / cg(ik)
       emean  = emean  + eb(ik)
       sigm01  = sigm01  + eb(ik)*sig(ik)
     END DO
 !
 ! 3.  Add tail beyond discrete spectrum
 !     ( DTH * SIG(NK) absorbed in FTxx )
 !
     eband  = eb(nk) / dden(nk)
     emean  = emean  + eband * fte
     sigm01 = sigm01 + eband * ftf
 !
 ! 4.  Final processing
 !
     sigm01 = sigm01 / emean
  
 !---- Temporary parameters (to be replaced by namelists) -----
  
     ptriad(1) = 1. 
     ptriad(2) = 10.
     ptriad(3) = 10. ! not used 
     ptriad(4) = 0.2
     ptriad(5) = 0.01
  
     hs = 4.*sqrt( max(0.,emean) )
     ursell = (grav*hs)/(2.*sqrt(2.)*sigm01**2*depth**2)
 !---------------------------------------------
  
     dep   = depth
     dep_2 = dep**2
     dep_3 = dep**3
     g9dep = grav * dep
 !
 !     --- compute some indices in sigma space
 !
     i2     = int(float(nk) / 2.)
     i1     = i2 - 1
     xis    = sig(i2) / sig(i1)
     xisln  = log( xis )
  
     isp    = int( log(2.) / xisln )
     isp1   = isp + 1
     wisp   = (2. - xis**isp) / (xis**isp1 - xis**isp)
     wisp1  = 1. - wisp
  
     ism    = int( log(0.5) / xisln )
     ism1   = ism - 1
     wism   = (xis**ism -0.5) / (xis**ism - xis**ism1)
     wism1  = 1. - wism
  
     e  = 0.
     sa = 0.
 !
 !     --- compute maximum frequency for which interactions are calculated
 !
     ismax = 1
     DO ik = 1, nk
      IF ( sig(ik) .LT. ( ptriad(2) * sigm01) ) THEN
         ismax = ik
      ENDIF
     ENDDO
     ismax = max( ismax , isp1 )
 !
 !     --- compute 3-wave interactions
 !
       IF (ursell.GE.ptriad(5) ) THEN ! AR: No need for switching it off from my point of view!
 !
 !       --- calculate biphase
 !
         biph   = (0.5*pi)*(tanh(ptriad(4)/ursell)-1.)
         sinbph = abs(sin(biph) )
         ef     = 0.
  
         DO ith = 1, nth
           DO ik = 1, nk
             e(ik)  = a(ith+(ik-1)*nth) * tpi * sig(ik) / cg(ik)
             ef(ik) = ef(ik) + e(ik)        
           END DO
           DO ik = 1, ismax
             e0  = e(ik)
             ed0 = eb(ik) 
             w0  = sig(ik)
             wn0 = wn(ik)
             c0  = w0 / wn0
             IF ( ik.GT.-ism1 ) THEN
                em  = wism * e(ik+ism1)   + wism1 * e(ik+ism)
                edm = wism * eb(ik+ism1)  + wism1 * eb(ik+ism)
                wm  = wism * sig(ik+ism1) + wism1 * sig(ik+ism)
                wnm = wism * wn(ik+ism1)  + wism1 * wn(ik+ism)
                cm  = wm / wnm
             ELSE
                em  = 0.
                edm = 0.
                wm  = 0.
                wnm = 0.
                cm  = 0.
             END IF
             aux1 = wnm**2 * ( g9dep + 2*cm**2 ) 
             aux2 = wn0*dep* (g9dep+(2./15.)*grav*dep_3*wn0**2-(2./5.)*w0**2*dep_2)
             rint = aux1 / aux2
             ft   = ptriad(1) * c0 * cg(ik) * rint**2 * sinbph 
             sa(ith,ik) = max(0.,ft * ( em * em - 2. * em * e0)) ! 1/(m²*s²) * m4 = m²/s² !!! [m²/s]
           END DO
         END DO
  
         DO ik = 1, nk - 1 
           sigpicg = sig(ik)*tpi/cg(ik) ! 1/s * s/m = 1/m
           DO ith = 1, nth
             stri = sa(ith,ik) - 2 * (wisp *  sa(ith,ik+isp1) + wisp1 *  sa(ith,ik+isp))
             IF (a(ith+(ik-1)*nth) .gt. jaceps) THEN
               d(ith+(ik-1)*nth) = stri / ((a(ith+(ik-1)*nth)) * sigpicg) 
               s(ith+(ik-1)*nth) = stri / sigpicg 
             ELSE
               d(ith+(ik-1)*nth) = 0.
               s(ith+(ik-1)*nth) = 0.
             ENDIF
           END DO
         END DO
       ELSE
         d = 0.
         s = 0.
       END IF
  
     !/
     !/ End of W3STR1 ----------------------------------------------------- /
     !/

References w3gdatmd::dden, w3gdatmd::dth, w3servmd::extcde(), w3gdatmd::fte, w3gdatmd::ftf, constants::grav, w3odatmd::ndse, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, constants::pi, w3gdatmd::sig, w3servmd::strace(), and constants::tpi.

Referenced by w3srcemd::w3srce().

Functions/Subroutines

Detailed Description

Function/Subroutine Documentation

◆ w3str1()