#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
! ----------------------------------------------------------------
! W3STRX Subr. Public User supplied triad interactions.
! INSTRX Subr. Public Corresponding initialization routine.
! ----------------------------------------------------------------
!
! 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, CG, WN, DEPTH, IX, S, D)
!/
!/ +-----------------------------------+
!/ | 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 )
!/ 05 Oct-2016 : Avoiding divide by zero for EMEAN ( version 5.15 )
!/
! 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
!/S USE W3SERVMD, ONLY: STRACE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(IN) :: CG(NK), WN(NK), DEPTH, A(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
!
!/S INTEGER, SAVE :: IENT = 0
INTEGER I1, I2, ID, IDDUM, II, IS, ISM, ISM1, ISMAX, &
ISP, ISP1, ITH, IK
REAL AUX1, AUX2, BIPH, C0, CM, DEP, DEP_2, DEP_3, E0, EM, HS, &
FT, RINT, SIGPICG, SINBPH, STRI, WISM, WISM1, WISP, &
WISP1, W0, WM, WN0, WNM, XIS, XISLN
REAL, ALLOCATABLE :: E(:), SA(:,:)
REAL :: EB(NK), EBAND, EMEAN, SIGM01
!----- Temp (to be moved) -----
REAL, ALLOCATABLE :: EF(:),SF(:)
REAL :: PTRIAD(5)
REAL :: URSELL
!------------------------------
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3STR1')
!
! 0. Initializations ------------------------------------------------ *
!
! **********************************************************
! *** The initialization routine should include all ***
! *** initialization, including reading data from files. ***
! **********************************************************
!
!> IF ( FIRST ) THEN
!> CALL INSTR1
!> FIRST = .FALSE.
!> END IF
!
! 1. .... ----------------------------------------------------------- *
!
!---- Compute SIGM01 (= 2pi/Tm01) for use in source term
!
! 1. Integral over directions
!
SIGM01 = 0.
EMEAN = 0.
! FMEAN = 0.
DO IK=1, NK
EB(IK) = 0.
DO ITH=1, NTH
EB(IK) = EB(IK) + A(ITH+(IK-1)*NTH)
END DO
END DO
!
! 2. Integrate over directions
!
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 = MAX ( 1.E-7 , SIGM01 ) / MAX(EMEAN,0.001)
!---- Temporary parameters (to be replaced by namelists) -----
PTRIAD(1) = 0.05
PTRIAD(2) = 2.5
PTRIAD(3) = 10.
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
!
! --- 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
ALLOCATE (E (1:NK))
ALLOCATE (SA(1:NTH,1:NK+ISP1))
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
!
! --- calculate biphase
!
BIPH = (0.5*PI)*(TANH(PTRIAD(4)/URSELL)-1.)
SINBPH = ABS( SIN(BIPH) )
!
ALLOCATE (EF (1:NK))
EF = 0.
DO ITH = 1, NTH
!
! --- initialize array with E(f) for the direction considered
! --- (convert from N(k) to E(f) using proper Jacobian)
!
DO IK = 1, NK
E(IK) = A(ITH+(IK-1)*NTH) * TPI * SIG(IK) / CG(IK)
!------------ Test ------------------------------------------
EF(IK) = EF(IK) + E(IK)
!------------------------------------------------------------
END DO
!
DO IK = 1, ISMAX
E0 = E(IK)
W0 = SIG(IK)
WN0 = WN(IK)
C0 = W0 / WN0
IF ( IK.GT.-ISM1 ) THEN
EM = WISM * E(IK+ISM1) + WISM1 * E(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.
WM = 0.
WNM = 0.
CM = 0.
END IF
AUX1 = WNM**2 * ( GRAV * DEP + 2.*CM**2 )
AUX2 = WN0 * DEP * ( GRAV * DEP + &
(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 ))
END DO
END DO
DEALLOCATE(EF)
!
! --- put source and diagonal terms together
! (using Jacobian for S(f) -> S(k))
!
ALLOCATE (SF (1:NK))
SF = 0.
DO IK = 1, NK
SIGPICG = SIG(IK) * 2. * PI / CG(IK)
DO ITH = 1, NTH
! --- Source term
S(ITH+(IK-1)*NTH) = 2.*( SA(ITH,IK) - &
2.*(WISP * SA(ITH,IK+ISP1) + &
WISP1 * SA(ITH,IK+ISP )) ) / &
SIGPICG
SF(IK) = 2.*( SA(ITH,IK) - &
2.*(WISP * SA(ITH,IK+ISP1) + &
WISP1 * SA(ITH,IK+ISP )) ) + SF(IK)
! --- Diagonal term
D = 0.
END DO
END DO
DEALLOCATE(SF)
ELSE
D = 0.
S = 0.
END IF
DEALLOCATE(E,SA)
RETURN
!/
!/ End of W3STRX ----------------------------------------------------- /
!/
END SUBROUTINE W3STR1
!/ ------------------------------------------------------------------- /
!/ SUBROUTINE INSTRX
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 23-Jun-2006 |
!/ +-----------------------------------+
!/
!/ 23-Jun-2006 : Origination. ( version 3.09 )
!/
! 1. Purpose :
!
! Initialization for source term routine.
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3STRX Subr. W3STRXMD Corresponding source term.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
! USE W3ODATMD, ONLY: NDSE
! USE W3SERVMD, ONLY: EXTCDE
!!/S USE W3SERVMD, ONLY: STRACE
!/
! IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!!/S CALL STRACE (IENT, 'INSTRX')
!
! 1. .... ----------------------------------------------------------- *
!
!/
!/ End of INSTRX ----------------------------------------------------- /
!/
! END SUBROUTINE INSTRX
!/
!/ End of module INSTRXMD -------------------------------------------- /
!/
END MODULE W3STR1MD