WW3/w3str1md_8F90_source.html

#include "w3macros.h"

!/ ------------------------------------------------------------------- /

MODULE w3str1md

  !/

  !/                  +-----------------------------------+

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |     A. J. van der Westhuysen      |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         13-Jan-2013 |

  !/                  +-----------------------------------+

  !/

  !/    13 Jan-2013 : Origination, based on SWAN v40.91 code ( version 4.08 )

  !/

  !/    Copyright 2009 National Weather Service (NWS),

  !/       National Oceanic and Atmospheric Administration.  All rights

  !/       reserved.  WAVEWATCH III is a trademark of the NWS.

  !/       No unauthorized use without permission.

  !/

  !  1. Purpose :

  !

  !     Module for inclusion of triad nonlinear interaction according to

  !     Eldeberky's (1996) Lumped Triad Interaction (LTA) source term.

  !

  !  2. Variables and types :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !     ----------------------------------------------------------------

  !

  !  3. Subroutines and functions :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !      W3STR1    Subr. Public   User supplied triad interactions.

  !     ----------------------------------------------------------------

  !

  !  4. Subroutines and functions used :

  !

  !      Name      Type  Module   Description

  !     ----------------------------------------------------------------

  !      STRACE    Subr. W3SERVMD Subroutine tracing.

  !     ----------------------------------------------------------------

  !

  !  5. Remarks :

  !

  !     WAVEWATCH III is designed as a highly plug-compatible code.

  !     Source term modules can be included as self-contained modules,

  !     with limited changes needed to the interface of routine calls

  !     in W3SRCE, and in the point postprocessing programs only.

  !     Codes submitted for inclusion in WAVEWATCH III should be

  !     self-contained in the way described below, and might be

  !     provided with distributions fully integrated in the data

  !     structure, or as an optional version of this module to be

  !     included by the user.

  !

  !     Rules for preparing a module to be included in or distributed

  !     with WAVEWATCH III :

  !

  !      - Fully document the code following the outline given in this

  !        file, and according to all other WAVEWATCH III routines.

  !      - Provide a file with necessary modifications to W3SRCE and

  !        all other routines that require modification.

  !      - Provide a test case with expected results.

  !      - It is strongly recommended that the programming style used

  !        in WAVEWATCH III is followed, in particular

  !          a) for readability, write as if in fixed FORTRAN format

  !             regarding column use, even though all files are F90

  !             free format.

  !          b) I prefer upper case programming for permanent code,

  !             as I use lower case in debugging and temporary code.

  !

  !     This module needs to be self-contained in the following way.

  !

  !      a) All saved variables connected with this source term need

  !         to be declared in the module header. Upon acceptance as

  !         permanent code, they will be converted to the WAVEWATCH III

  !         dynamic data structure.

  !      b) Provide a separate computation and initialization routine.

  !         In the submission, the initialization should be called

  !         from the computation routine upon the first call to the

  !         routine. Upon acceptance as permanent code, the

  !         initialization routine will be moved to a more appropriate

  !         location in the code (i.e., being absorbed in ww3_grid or

  !         being moved to W3IOGR).

  !

  !     See notes in the file below where to add these elements.

  !

  !  6. Switches :

  !

  !     !/S  Enable subroutine tracing.

  !

  !  7. Source code :

  !/

  !/ ------------------------------------------------------------------- /

  !/

  !     *****************************************

  !     ***    Declare saved variables here   ***

  !     ***  public or private as appropriate ***

  !     *****************************************

  !

  PUBLIC

  !/

CONTAINS

  !/ ------------------------------------------------------------------- /

  SUBROUTINE w3str1 (A, AOLD, CG, WN, DEPTH, IX, S, D)

    !/

    !/                  +-----------------------------------+

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

    !/

  END SUBROUTINE w3str1

  !/ ------------------------------------------------------------------- /

END MODULE w3str1md