w3snl3md.F90

Go to the documentation of this file.

 
 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3snl3md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH-III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         13-Jul-2012 |
  !/                  +-----------------------------------+
  !/
  !/    21-Jul-2008 : Origination as NLX option.          ( version 3.13 )
  !/    03-Jan-2009 : Bug fixes INSNLX.                   ( version 3.13 )
  !/                  See remarks section for module.
  !/    25-Aug-2009 : Conversion to F(f,theta) form.      ( version 3.13 )
  !/    13-Nov-2009 : Bug fix DELTH in initialization.    ( version 3.13 )
  !/    01-Dec-2009 : Bug fix frequency filtering.        ( version 3.13 )
  !/    13-Aug-2010 : Move to NL3.                        ( version 3.15 )
  !/    13-Jul-2012 : Moved from version 3.15 to 4.08.    ( version 4.08 )
  !/
  !/    Copyright 2008-2012 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 :
  !
  !     Generalized and optimized multiple DIA implementation.
  !     Expressions in terms of original F(f,theta) spectrum.
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      NKD       I.P.  Private  Number of nondimensional depths in
  !                               storage array.
  !      KDMIN     R.P.  Private  Minimum relative depth in table.
  !      KDMAX     R.P.  Private  Maximum relative depth in table.
  !      LAMMAX    R.P.  Public   Maximum value for lambda or mu.
  !      DELTHM    R.P.  Public   Maximum angle gap (degree).
  !      SITMIN    Real  Private  Minimum nondimensional radian
  !                               frequency in table.
  !      XSIT      Real  Private  Corresponding increment factor.
  !     ----------------------------------------------------------------
  !
  !     See W3SNL3 and INSNL3 for documentation of variables in W3GDATMD
  !     as used here.
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SNL3    Subr. Public   Multiple DIA for arbitrary depth.
  !      EXPAND    Subr. W3SNL3   Expand spectrum for indirect address.
  !      EXPND2    Subr. W3SNL3   Expand Snl and D contributions.
  !      INSNL3    Subr. Public   Corresponding initialization routine.
  !      MINLAM    R.F.  INSNL3   Minimum lambda for quadruplet.
  !      MAXLAM    R.F.  INSNL3   Maximum lambda for quadruplet.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !      EXTCDE    Subr. W2SERVMD Program abort.
  !      WAVNU1    Subr. W3DISPMD Solve dispersion relation.
  !      WAVNU2    Subr. W3DISPMD Solve dispersion relation.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !     - Filtering techniques for computation of quadruplet spectral
  !       values and distribution in spectral space have been tested
  !       but were not found worth the large coding effort involved.
  !     - WAVNU1 is used in W3SNL3  for consistency with spectral grid
  !       description.
  !     - WAVNU2 is used in INSNL3  for accuracy in the computation of
  !       the layout of the quadruplets (higher computational cost is
  !       not an issue with initialization routine).
  !     - For large lambda or mu the original maximum kd = 10. still
  !       leads to significantly different quadruplet layout in
  !       secion 3. To remedy this, the orriginal settings of the
  !       lookup tables
  !
  !     INTEGER, PRIVATE, PARAMETER :: NKD = 250
  !     REAL, PRIVATE, PARAMETER    :: KDMIN = 0.025 ,  KDMAX = 10.
  !
  !       was reset to
  !
  !     INTEGER, PRIVATE, PARAMETER :: NKD = 275
  !     REAL, PRIVATE, PARAMETER    :: KDMIN = 0.025 ,  KDMAX = 20.
  !
  !       for the bug fix of 03-Jan-2009. Note that with this, the
  !       estimate of NTHMAX in INSNL3 also is needed to guarantee
  !       consistent NTHMAX and NTHM2 for any lambda and mu.
  !
  !  6. Switches :
  !
  !     !/S    Enable subroutine tracing.
  !     !/Tn   Test output (see main subroutines).
  !
  !  7. Source code :
  !/
  !/ ------------------------------------------------------------------- /
  !/
  INTEGER, PRIVATE, PARAMETER :: NKD = 275
  REAL, PRIVATE, PARAMETER    :: KDMIN = 0.025 ,  kdmax = 20.
  REAL, PUBLIC, PARAMETER     :: lammax = 0.49999
  REAL, PUBLIC, PARAMETER     :: delthm = 90.
  !
  REAL, PRIVATE               :: sitmin, xsit
  !
  PUBLIC
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE w3snl3 (  A, CG, WN, DEPTH, S, D )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH-III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         01-Dec-2009 |
    !/                  +-----------------------------------+
    !/
    !/    21-Jul-2008 : Origination as NLX option.          ( version 3.13 )
    !/    25-Aug-2009 : Conversion to F(f,theta) form.      ( version 3.13 )
    !/    01-Dec-2009 : Bug fix frequency filtering.        ( version 3.13 )
    !/
    !  1. Purpose :
    !
    !     Multiple Discrete Interaction Parameterization for arbitrary
    !     depths with generalized quadruplet layout.
    !
    !  2. Method :
    !
    !     This is a direct implementation of the Discrete Interaction
    !     Paramterization (DIA) with multiple representative quadruplets
    !     (MDIA) for arbitrary water depths.
    !
    !     The outer loop of the code is over quadruplet realizations,
    !     which implies two realizations for a conventional quadruplet
    !     definitions and four for extended definitions (with rescaling
    !     of the contants for consistency). Within this loop the compu-
    !     tations are performed in two stages. First, interactions
    !     contributions are computed for the entire spectral space,
    !     second all contributions are combined into the actual inter-
    !     actions and diagonal contributions.
    !
    !     Arbitrary depths are addressed by generating a lookup table
    !     for the relative depth. These tables are used for each discrete
    !     frequency separately. Efficient memory usages requires relative
    !     addressing to reduce the size of the lookup tables. To use this
    !     the spectral space is expanded to higher and lower frequencies
    !     as well as directional space is expanded/volded. This is done
    !     for the input (pseudo-) spectrum (action spectrum devided by the
    !     wavenumber) to determine spectral densities at the quadruplet
    !     components, and the spectral space describing individual contri-
    !     butions before they are combined into the actual interactions.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A       R.A.  I   Action spectrum A(ITH,IK) as a function of
    !                         direction (rad)  and wavenumber.
    !       CG      R.A.  I   Group velocities (dimension NK).
    !       WN      R.A.  I   Wavenumbers (dimension NK).
    !       DEPTH   Real  I   Water depth in meters.
    !       S       R.A.  O   Source term.
    !       D       R.A.  O   Diagonal term of derivative.
    !     ----------------------------------------------------------------
    !
    !     Variables describing the expanded frequency space from the
    !     dynamic storage in w3gdatmd.
    !
    !      Name      Type  Scope    Description
    !     ----------------------------------------------------------------
    !      NFR       Int.  Public   Number of frequencies or wavenumbers
    !                               in discrete spectral space (NFR=>NK).
    !      NFRMIN    Int.  Public   Minimum discrete frequency in the
    !                               expanded frequency space.
    !      NFRMAX    Int.  Public   Idem maximum for first part.
    !      NFRCUT    Int.  Public   Idem maximum for second part.
    !      NTHMAX    Int.  Public   Extension of directional space.
    !      NTHEXP    Int   Public   Number of bins in extended dir. space.
    !      NSPMIN, NSPMAX, NSPMX2
    !                Int.  Public   1D spectral space range.
    !      FRQ       R.A.  Public   Expanded frequency range (Hz).
    !      XSI       R.A.  Public   Expanded frequency range (rad/s).
    !     ----------------------------------------------------------------
    !
    !     Variables describing lookup tables.
    !
    !      Name      Type  Scope    Description
    !     ----------------------------------------------------------------
    !      NQA       Int.  Public   Number of actual quadruplets.
    !      QST1      I.A.  Public   Spectral offsets for compuation of
    !                               quadruplet spectral desnities.
    !      QST2      R.A.  Public   Idem weights.
    !      QST3      R.A.  Public   Norm. factors in product term and
    !                               in diagonal strength.
    !      QST4      I.A.  Public   Spectral offsets for combining of
    !                               interactions and diagonal.
    !      QST5      R.A.  Public   Idem weights for interactions.
    !      QST6      R.A.  Public   Idem weights for diagonal.
    !     ----------------------------------------------------------------
    !
    !     Variables describing model setup.
    !
    !      Name      Type  Scope    Description
    !     ----------------------------------------------------------------
    !      SNLMSC    Real  Public   Tuning power 'deep' scaling.
    !      SNLNSC    Real  Public   Tuning power 'shallow' scaling.
    !      SNLSFD    Real  Public   'Deep' nondimensional filer freq.
    !      SNLSFS    Real  Public   'Shallow' nondimensional filer freq.
    !     ----------------------------------------------------------------
    !
    !  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 :
    !
    !     - Note that this code uses explicit unroling of potential loop
    !       structures for optimization purposes.
    !     - Normalization with respect to the number of quadruplets is
    !       included in the proportionality constant.
    !     - Note that the outer loop in the routine considers one actual
    !       quadruplet realization per loop cycle. For the traditional
    !       quadruplet layout two realizations occure, for the expanded
    !       four realizations occur. For consistency, strength of a
    !       traditional layout is therefore doubled.
    !     - 1D representation is used of 2D spectral space for optimization
    !       purposes.
    !     - Contributions are first computed in the convetional spectral
    !       space and are then expancded "in place" into the expanded
    !       spectral space in EXPND2.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S    Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nfr => nk, nth, sig, fachfe, facti1, facti2,&
         nfrmin, nfrmax, nfrcut, nthmax, nthexp,     &
         nspmin, nspmax, nspmx2, frq, xsi, nqa,      &
         qst1, qst2, qst3, qst4, qst5, qst6, snlmsc, &
         snlnsc, snlsfd, snlsfs
    USE w3odatmd, ONLY: ndse, ndst
    !
    USE w3servmd, ONLY: extcde
    USE w3dispmd, ONLY: wavnu1, wavnu3
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)   :: A(NTH,NFR), CG(NFR), WN(NFR), DEPTH
    REAL, INTENT(OUT)  :: S(NTH,NFR), D(NTH,NFR)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER            :: IFR, IERR, IKD, JKD(NFRCUT), IQA, IF1MIN, &
         IF1MAX, IF2MIN, IF2MAX, ISP0, ISPX0, ITH, &
         ISP, ISPX
#ifdef W3_S
    INTEGER, SAVE      :: IENT = 0
#endif
    INTEGER            :: LQST1(16), LQST4(16)
    REAL               :: XSITLN, SIT, FPROP, FQ1, FQ2, FQ3, FQ4,   &
         AUX1, AUX2
    REAL               :: XWN(NFRMAX), XCG(NFRMAX), SCALE1(NFRCUT), &
         SCALE2(NFRCUT), LQST2(16), FACT(6),       &
         LQST5(16), LQST6(16)
    REAL               :: UE(NSPMIN:NSPMAX), DSB(NSPMIN:NSPMX2),    &
         DD1(NSPMIN:NSPMX2), DD2(NSPMIN:NSPMX2),   &
         DD3(NSPMIN:NSPMX2), DD4(NSPMIN:NSPMX2)
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SNL3')
#endif
    !
    ! 1.  Initialization ------------------------------------------------- *
    ! 1.a Constants and arrays
    !
    xsitln = log(xsit)
    !
    s      = 0.
    d      = 0.
    !     DSB    = 0.
    !     DD1    = 0.
    !     DD2    = 0.
    !     DD3    = 0.
    !     DD4    = 0.
    !
    ! 1.a Extended frequency range
    !
    xwn(1:nfr) = wn
    xcg(1:nfr) = cg
    !
    DO ifr = nfr+1, nfrmax
#ifdef W3_PDLIB
      CALL wavnu3(xsi(ifr), depth, xwn(ifr), xcg(ifr))
#else
      CALL wavnu1(xsi(ifr), depth, xwn(ifr), xcg(ifr))
#endif
    END DO
    !
    ! 1.b Expanded pseudo spetrum
    !
    CALL expand ( ue )
    !
    ! 1.c Set up scaling functions
    !
    aux1   = 1. / ( tpi**11 * grav**(4.-snlmsc) )
    aux2   = grav**2 / tpi**11
    !
    DO ifr=1, nfrcut
      scale1(ifr) = aux1 * xwn(ifr)**(4.+snlmsc) *                 &
           xsi(ifr)**(13.-2.*snlmsc) / xcg(ifr)**2
      scale2(ifr) = aux2 * xwn(ifr)**11 *                          &
           (xwn(ifr)*depth)**snlnsc / xcg(ifr)
    END DO
    !
    ! 1.d Set up depth scaling counters
    !
    DO ifr=1, nfrcut
      sit      = xsi(ifr) * sqrt(depth/grav)
      ikd      = 1 + nint( ( log(sit) - log(sitmin) ) / xsitln )
      jkd(ifr) = max( 1 , min(ikd,nkd) )
    END DO
    !
    ! 2.  Base loop over quadruplet realizations ------------------------- *
    !
    DO iqa=1 , nqa
      !
      ! 3.  Obtain quadruplet energies for all spectral bins --------------- *
      ! 3.a Set frequency ranges
      !
      aux1   = qst3(5,iqa,1)
      aux2   = qst3(6,iqa,1)
      !
      if1min = 1
      if1max = nfrcut
      if2min = 1
      if2max = nfr
      !
      IF ( aux1 .LE. 0. .AND. aux2 .LE. 0. ) THEN
        !
        cycle
        !
      ELSE IF ( aux2 .LE. 0. ) THEN
        !
        sit    = snlsfd * sqrt(grav/depth)
        ifr    = nint( facti2 + facti1*log(sit) )
        IF ( ifr .GT. nfr ) cycle
        !
        IF ( ifr .GT. 1 ) THEN
          if1min = max( 1 , ifr )
          if2min = max( 1 , if1min + nfrmin )
          dsb(1:(if1min-1)*nth) = 0.
          dd1(1:(if1min-1)*nth) = 0.
          dd2(1:(if1min-1)*nth) = 0.
          dd3(1:(if1min-1)*nth) = 0.
          dd4(1:(if1min-1)*nth) = 0.
        END IF
        !
      ELSE IF ( aux1 .LE. 0. ) THEN
        !
        sit    = snlsfs * sqrt(grav/depth)
        ifr    = nint( facti2 + facti1*log(sit) )
        IF ( ifr .LT. 1 ) cycle
        !
        IF ( ifr .LT. nfrcut ) THEN
          if1max = min( nfrcut, ifr )
          !               IF2MAX = NFR
          dsb(if1max*nth+1:nfrcut*nth) = 0.
          dd1(if1max*nth+1:nfrcut*nth) = 0.
          dd2(if1max*nth+1:nfrcut*nth) = 0.
          dd3(if1max*nth+1:nfrcut*nth) = 0.
          dd4(if1max*nth+1:nfrcut*nth) = 0.
        END IF
        !
      END IF
      !
      ! 3.b Loop over frequencies
      !
      DO ifr=if1min, if1max
        !
        ! 3.c Find discrete depths
        !
        ikd    = jkd(ifr)
        !
        ! 3.d Get offsets and weights
        !
        lqst1  = qst1(:,iqa,ikd)
        lqst2  = qst2(:,iqa,ikd)
        fact   = qst3(:,iqa,ikd)
        fact(1:4) = fact(1:4) * xcg(ifr) / ( xwn(ifr) *xsi(ifr) )
        fprop  = scale1(ifr)*fact(5) + scale2(ifr)*fact(6)
        !
        ! 3.e Loop over directions
        !
        isp0   = (ifr-1)*nth
        ispx0  = (ifr-1)*nthexp
        !
        DO ith=1, nth
          !
          isp    = isp0 + ith
          ispx   = ispx0 + ith
          !
          fq1    = ( ue(ispx+lqst1( 1)) * lqst2( 1) +               &
               ue(ispx+lqst1( 2)) * lqst2( 2) +               &
               ue(ispx+lqst1( 3)) * lqst2( 3) +               &
               ue(ispx+lqst1( 4)) * lqst2( 4) ) * fact(1)
          fq2    = ( ue(ispx+lqst1( 5)) * lqst2( 5) +               &
               ue(ispx+lqst1( 6)) * lqst2( 6) +               &
               ue(ispx+lqst1( 7)) * lqst2( 7) +               &
               ue(ispx+lqst1( 8)) * lqst2( 8) ) * fact(2)
          fq3    = ( ue(ispx+lqst1( 9)) * lqst2( 9) +               &
               ue(ispx+lqst1(10)) * lqst2(10) +               &
               ue(ispx+lqst1(11)) * lqst2(11) +               &
               ue(ispx+lqst1(12)) * lqst2(12) ) * fact(3)
          fq4    = ( ue(ispx+lqst1(13)) * lqst2(13) +               &
               ue(ispx+lqst1(14)) * lqst2(14) +               &
               ue(ispx+lqst1(15)) * lqst2(15) +               &
               ue(ispx+lqst1(16)) * lqst2(16) ) * fact(4)
          !
          aux1   = fq1 * fq2 * ( fq3 + fq4 )
          aux2   = fq3 * fq4 * ( fq1 + fq2 )
          dsb(isp) = fprop * ( aux1 - aux2 )
          !
          aux1   = fq3 + fq4
          aux2   = fq3 * fq4
          dd1(isp) = fprop * fact(1) * ( fq2 * aux1 - aux2 )
          dd2(isp) = fprop * fact(2) * ( fq1 * aux1 - aux2 )
          !
          aux1   = fq1 + fq2
          aux2   = fq1 * fq2
          dd3(isp) = fprop * fact(3) * ( aux2 - fq4*aux1 )
          dd4(isp) = fprop * fact(4) * ( aux2 - fq3*aux1 )
          !
          ! ... End loop 3.e
          !
        END DO
        !
        ! ... End loop 3.b
        !
      END DO
      !
      ! 3.e Expand arrays
      !
      CALL expnd2 ( dsb(1:nth*nfrcut), dsb )
      CALL expnd2 ( dd1(1:nth*nfrcut), dd1 )
      CALL expnd2 ( dd2(1:nth*nfrcut), dd2 )
      CALL expnd2 ( dd3(1:nth*nfrcut), dd3 )
      CALL expnd2 ( dd4(1:nth*nfrcut), dd4 )
      !
      ! 4.  Put it all together -------------------------------------------- *
      ! 4.a Loop over frequencies
      !
      DO ifr=if2min, if2max
        !
        ! 4.b Find discrete depths and storage
        !
        ikd    = jkd(ifr)
        !
        ! 4.c Get offsets and weights
        !
        lqst4  = qst4(:,iqa,ikd)
        lqst5  = qst5(:,iqa,ikd)
        lqst6  = qst6(:,iqa,ikd)
        !
        ! 4.d Loop over directions
        !
        ispx0  = (ifr-1)*nthexp
        !
        DO ith=1, nth
          !
          ispx   = ispx0 + ith
          !
          s(ith,ifr) = s(ith,ifr) + dsb(ispx+lqst4( 1)) * lqst5( 1) &
               + dsb(ispx+lqst4( 2)) * lqst5( 2) &
               + dsb(ispx+lqst4( 3)) * lqst5( 3) &
               + dsb(ispx+lqst4( 4)) * lqst5( 4) &
               + dsb(ispx+lqst4( 5)) * lqst5( 5) &
               + dsb(ispx+lqst4( 6)) * lqst5( 6) &
               + dsb(ispx+lqst4( 7)) * lqst5( 7) &
               + dsb(ispx+lqst4( 8)) * lqst5( 8) &
               + dsb(ispx+lqst4( 9)) * lqst5( 9) &
               + dsb(ispx+lqst4(10)) * lqst5(10) &
               + dsb(ispx+lqst4(11)) * lqst5(11) &
               + dsb(ispx+lqst4(12)) * lqst5(12) &
               + dsb(ispx+lqst4(13)) * lqst5(13) &
               + dsb(ispx+lqst4(14)) * lqst5(14) &
               + dsb(ispx+lqst4(15)) * lqst5(15) &
               + dsb(ispx+lqst4(16)) * lqst5(16)
          !
          d(ith,ifr) = d(ith,ifr) + dd1(ispx+lqst4( 1)) * lqst6( 1) &
               + dd1(ispx+lqst4( 2)) * lqst6( 2) &
               + dd1(ispx+lqst4( 3)) * lqst6( 3) &
               + dd1(ispx+lqst4( 4)) * lqst6( 4) &
               + dd2(ispx+lqst4( 5)) * lqst6( 5) &
               + dd2(ispx+lqst4( 6)) * lqst6( 6) &
               + dd2(ispx+lqst4( 7)) * lqst6( 7) &
               + dd2(ispx+lqst4( 8)) * lqst6( 8) &
               + dd3(ispx+lqst4( 9)) * lqst6( 9) &
               + dd3(ispx+lqst4(10)) * lqst6(10) &
               + dd3(ispx+lqst4(11)) * lqst6(11) &
               + dd3(ispx+lqst4(12)) * lqst6(12) &
               + dd4(ispx+lqst4(13)) * lqst6(13) &
               + dd4(ispx+lqst4(14)) * lqst6(14) &
               + dd4(ispx+lqst4(15)) * lqst6(15) &
               + dd4(ispx+lqst4(16)) * lqst6(16)
          !
          ! ... End loop 4.d
          !
        END DO
        !
        ! ... End loop 4.a
        !
      END DO
      !
      ! ... End of loop 2.
      !
    END DO
    !
    ! 5.  Convert back to wave action ------------------------------------ *
    !
    DO ifr=if2min, if2max
      s(:,ifr) = s(:,ifr) / xsi(ifr) * xcg(ifr) * tpiinv
    END DO
    !
    RETURN
    !/
    !/ Embedded subroutines
    !/
  CONTAINS
    !/ ------------------------------------------------------------------- /
 
    SUBROUTINE expand ( SPEC )
      !/
      !/                  +-----------------------------------+
      !/                  | WAVEWATCH-III           NOAA/NCEP |
      !/                  |           H. L. Tolman            |
      !/                  |                        FORTRAN 90 |
      !/                  | Last update :         21-Aug-2009 |
      !/                  +-----------------------------------+
      !/
      !/    03-Jul-2008 : Origination.                        ( version 3.13 )
      !/    21-Aug-2009 : Conversion to F(f,theta) form.      ( version 3.13 )
      !/
      !  1. Purpose :
      !
      !     Expand spectrum, subroutine used to simplify addressing.
      !
      !  3. Parameters :
      !
      !     Parameter list
      !     ----------------------------------------------------------------
      !       SPEC    R.A.  O   Expanded spectrum.
      !     ----------------------------------------------------------------
      !
      ! 10. Source code :
      !
      !/ ------------------------------------------------------------------- /
      IMPLICIT NONE
      !/
      !/ Parameter list
      !/
      REAL, INTENT(OUT)       :: SPEC(1-NTHMAX:NTH+NTHMAX,NFRMIN:NFRMAX)
      !/
      !/ Local parameters
      !/
      INTEGER                 :: IFR, ITH
      !/
      !/ ------------------------------------------------------------------- /
      !
      spec(:,nfrmin:0) = 0.
      !
      spec(1:nth,1:nfr) = a * tpi
      !
      DO ifr=1, nfr
        spec(1:nth,ifr) = spec(1:nth,ifr) * xsi(ifr) / xcg(ifr)
      END DO
      !
      DO ifr=nfr+1, nfrmax
        spec(1:nth,ifr) = spec(1:nth,ifr-1) * fachfe
      END DO
      !
      DO ith=1, nthmax
        spec(nth+ith,1:nfrmax) = spec(   ith   ,1:nfrmax)
        spec( 1 -ith,1:nfrmax) = spec(nth+1-ith,1:nfrmax)
      END DO
      !
      RETURN
      !/
      !/ End of EXPAND ----------------------------------------------------- /
      !/
    END SUBROUTINE expand
    !/ ------------------------------------------------------------------- /
 
    SUBROUTINE expnd2 ( ARIN, AROUT )
      !/
      !/                  +-----------------------------------+
      !/                  | WAVEWATCH-III           NOAA/NCEP |
      !/                  |           H. L. Tolman            |
      !/                  |                        FORTRAN 90 |
      !/                  | Last update :         16-Jul-2008 |
      !/                  +-----------------------------------+
      !/
      !/    16-Jul-2008 : Origination.                        ( version 3.13 )
      !/
      !  1. Purpose :
      !
      !     Expand spectrum to simplify indirect addressing.
      !     Done 'in place' with temporary array ( ARIN = AROUT )
      !
      !  3. Parameters :
      !
      !     Parameter list
      !     ----------------------------------------------------------------
      !       SPIN    R.A.  I   Input array.
      !       SPOUT   R.A.  I   Output array.
      !     ----------------------------------------------------------------
      !
      ! 10. Source code :
      !
      !/ ------------------------------------------------------------------- /
      IMPLICIT NONE
      !/
      !/ Parameter list
      !/
      REAL, INTENT(IN)      :: ARIN(NTH,NFRCUT)
      REAL, INTENT(OUT)     :: AROUT(1-NTHMAX:NTH+NTHMAX,NFRMIN:NFRCUT)
      !/
      !/ Local parameters
      !/
      INTEGER               :: IFR, ITH
      REAL                  :: TEMP(NTH,NFRCUT)
      !/
      !/ ------------------------------------------------------------------- /
      !
      temp   = arin
      !
      arout(:,nfrmin:0) = 0.
      !
      arout(1:nth,1:nfrcut) = temp
      !
      DO ith=1, nthmax
        arout(nth+ith,1:nfrcut) = arout(   ith   ,1:nfrcut)
        arout( 1 -ith,1:nfrcut) = arout(nth+1-ith,1:nfrcut)
      END DO
      !
      RETURN
      !/
      !/ End of EXPND2 ----------------------------------------------------- /
      !/
    END SUBROUTINE expnd2
    !/
    !/ End of W3SNL3 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3snl3
  !/ ------------------------------------------------------------------- /
  SUBROUTINE insnl3
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH-III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Nov-2009 |
    !/                  +-----------------------------------+
    !/
    !/    21-Jul-2008 : Origination as NLX option.          ( version 3.13 )
    !/    03-Jan-2009 : Bug fixes NTHMAX and NTHMX2.        ( version 3.13 )
    !/    21-Aug-2009 : Conversion to F(f,theta) form.      ( version 3.13 )
    !/    13-Nov-2009 : Harden DELTH computation.           ( version 3.13 )
    !/
    !  1. Purpose :
    !
    !     Initialization for generalized multiple DIA routine.
    !
    !  2. Method :
    !
    !     Fill storage aryays as described in the main subroutine with
    !     interpolation, model and distribution data.
    !
    !  3. Parameters :
    !
    !     Variables in W3GDATMD describing model setup.
    !
    !      Name      Type  Scope    Description
    !     ----------------------------------------------------------------
    !      SNLNQ     Int.  Public   Number of quadruplet definitions.
    !      SNLL      R.A.  Public   Array with lambda for quadruplet.
    !      SNLM      R.A.  Public   Array with mu for quadruplet.
    !      SNLT      R.A.  Public   Array with Dtheta for quadruplet.
    !      SNLCD     R.A.  Public   Array with Cd for quadruplet.
    !      SNLCS     R.A.  Public   Array with Cs for quadruplet.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !      EXTCDE    Subr. W3SERVMD Program abort.
    !      WAVNU2    Subr. W3DISPMD Solve dispersion relation.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3IOGR    Subr. W3IOGRMD Process model definiton file.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     See error escape location.
    !
    !  8. Remarks :
    !
    !     - Allocation of arrays directly done in data structure, using
    !       IGRID and resetting pointer of aliaases.
    !     - In the 03-Jan-2009 bug fix !/T3 error output was fixed, and
    !       NTHMAX is increased by 1 to assure that NTHMX2 .LE. NTHMAX
    !       for any lambda and mu. With this, the label 810 test is
    !       changed from equality testing to .LE. testing.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S    Enable subroutine tracing.
    !     !/T    General test output.
    !     !/T1   Filling of lookup table for quadruplet and interaction
    !            strength.
    !     !/T2   Filling of lookup table for combining interactions.
    !     !/T3   Display raw lookup table of second type.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3odatmd, ONLY: ndse, ndst
    USE w3gdatmd,       nfr => nk
    !
    USE w3dispmd, ONLY: wavnu2
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IFRMIN, IFRMAX, IKD, IERR, IQ, NQD,  &
         NQS, J, IFR, IQA, JJ, JF, NTHMX2,    &
         JIQ, JOF, JQR, IST
    INTEGER                 :: JFR(4), JFR1(4), JTH(4), JTH1(4)
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER, ALLOCATABLE    :: AST1(:,:,:), AST2(:,:,:)
    REAL                    :: SITMAX, XFRLN
    REAL                    :: OFF12, OFF34, TH12, DEPTH,           &
         S0, S1, S2, S3, S4, AUXFR(4),        &
         WN0, WN1, WN2, WN3, WN4,             &
         CG0, CG1, CG2, CG3, CG4, AUXF,       &
         AA, BB, CC, DELTH(4), AUX1, AUX2,    &
         WFR(4), WFR1(4), WTH(4), WTH1(4),    &
         WFROFF, SIOFF, WF
    !
    TYPE qst
      INTEGER               :: OFR(4), OFR1(4), OTH(4), OTH1(4)
      REAL                  :: HFR(4), HFR1(4), HTH(4), HTH1(4)
      REAL                  :: F1, F2, F3, F4, CQD, CQS
    END TYPE qst
    !
    TYPE(qst), ALLOCATABLE  :: TSTORE(:,:)
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'INSNL3')
#endif
    !
    ! 1.  Initialization ------------------------------------------------- *
    ! 1.a Checks
    !
    xfrln  = log(xfr)
    !
    IF ( lammax.LE.0. .OR. lammax.GT.0.5 .OR. delthm.LT.0. ) GOTO 800
    !
    ! 1.b Set up relative depths
    !
    ALLOCATE ( tstore(snlnq*4,1:nkd) )
    !
    depth  = 1.
    sitmin = sqrt( kdmin * tanh(kdmin) )
    sitmax = sqrt( kdmax * tanh(kdmax) )
    xsit   = (sitmax/sitmin)**(1./real(nkd-1))
    !
#ifdef W3_T
    WRITE (ndst,9010) nkd, kdmin, kdmax, xsit
#endif
    !
    ! 2.  Building quadruplet data base ---------------------------------- *
    !     For quadruplet and interaction strength evaluation
    !
    ifrmin = 0
    ifrmax = 0
    nthmax = 0
    !
    ! 2.a Loop over relative depths
    !
    s0     = sitmin * sqrt( grav / depth )  / xsit
    !
    DO ikd=1, nkd
      !
      s0     = s0 * xsit
      CALL wavnu2 ( s0, depth, wn0, cg0, 1.e-6, 25, ierr)
      !
      ! 2.b Loop over representative quadruplets
      !
      nqa    = 0
      nqd    = 0
      nqs    = 0
      !
      DO iq=1, snlnq
        !
#ifdef W3_T1
        WRITE (ndst,9020) ikd, iq, wn0*depth, s0*tpiinv, depth
#endif
        !
        off12  = snlm(iq)
        off34  = snll(iq)
        th12   = snlt(iq) * dera
        IF ( snlcd(iq) .GT. 0. ) nqd = nqd + 1
        IF ( snlcs(iq) .GT. 0. ) nqs = nqs + 1
        !
        IF ( th12 .LT. 0. ) THEN
          IF ( off12.LT.0. .OR. off12.GT.0.5 .OR.                 &
               off34.LT.0. .OR. off34.GT.0.5 ) GOTO 801
        ELSE
          IF ( snlt(iq).GT.delthm .OR. off12.LT.0. .OR.           &
               off12.GE.1.                                        &
               .OR.  off34.LT.minlam(off12,snlt(iq)) .OR.         &
               off34.GT.maxlam(off12,snlt(iq)) ) GOTO 802
        END IF
        !
#ifdef W3_T1
        WRITE (ndst,9021) snlt(iq), off12, off34,               &
             snlcd(iq), snlcs(iq)
#endif
        !
        ! 2.c Offset angles
        !
        s1     = s0 * ( 1. + off12 )
        CALL wavnu2 ( s1, depth, wn1, cg1, 1.e-6, 25, ierr)
        s2     = s0 * ( 1. - off12 )
        CALL wavnu2 ( s2, depth, wn2, cg2, 1.e-6, 25, ierr)
        s3     = s0 * ( 1. + off34 )
        CALL wavnu2 ( s3, depth, wn3, cg3, 1.e-6, 25, ierr)
        s4     = s0 * ( 1. - off34 )
        CALL wavnu2 ( s4, depth, wn4, cg4, 1.e-6, 25, ierr)
        !
        auxfr(1) = s1 / s0
        auxfr(2) = s2 / s0
        auxfr(3) = s3 / s0
        auxfr(4) = s4 / s0
        !
        IF ( th12 .LT. 0. ) THEN
          bb = 2. * wn0
        ELSE
          bb = wn1**2 + wn2**2 + 2.*wn1*wn2*cos(th12)
          bb = sqrt( max( bb , 0. ) )
        END IF
        !
        IF ( th12.LT.0. .AND. abs(off12).LE.1.e-4 ) THEN
          delth(1) = 0.
          delth(2) = 0.
        ELSE
          cc       = wn1
          aa       = wn2
          aux1     = (cc**2+bb**2-aa**2) / (2.*bb*cc)
          aux2     = (aa**2+bb**2-cc**2) / (2.*bb*aa)
          delth(1) = - acos( max( 0. , min( 1. , aux1 ) ) )
          delth(2) =   acos( max( 0. , min( 1. , aux2 ) ) )
        END IF
        cc       = wn3
        aa       = wn4
        aux1     = (cc**2+bb**2-aa**2) / (2.*bb*cc)
        aux2     = (aa**2+bb**2-cc**2) / (2.*bb*aa)
        delth(3) = - acos( max( 0. , min( 1. , aux1 ) ) )
        delth(4) =   acos( max( 0. , min( 1. , aux2 ) ) )
        !
#ifdef W3_T1
        WRITE (ndst,9022) delth(:) * rade
#endif
        !
        ! 2.d Frequency indices
        !
        DO j=1, 4
          jfr(j) = int( log(auxfr(j)) / xfrln )
          jfr1(j) = jfr(j) + 1 * sign(1.,auxfr(j)-1.)
          wfr(j) = (xfr**jfr1(j)-auxfr(j))/(xfr**jfr1(j)-xfr**jfr(j))
          wfr1(j) = 1. - wfr(j)
        END DO
        !
        ifrmin = min( ifrmin , minval(jfr1) )
        ifrmax = max( ifrmax , maxval(jfr1) )
        !
#ifdef W3_T1
        WRITE (ndst,9023) 1, jfr(1), jfr1(1), wfr(1), wfr1(1)
        DO, j=2, 4
          WRITE (ndst,9024) j, jfr(j), jfr1(j), wfr(j), wfr1(j)
        END DO
#endif
        !
        ! 2.e Directional indices
        !
        DO j=1, 4
          aux1    = delth(j) / dth
          jth(j) = int(aux1)
          jth1(j) = jth(j) + 1 * sign(1.,delth(j))
          wth1(j) = abs(aux1) - real(abs(jth(j)))
          wth(j) = 1. - wth1(j)
        END DO
        !
        nthmax = max( nthmax , maxval(abs(jth1)) )
        !
#ifdef W3_T1
        WRITE (ndst,9025) 1, jth(1), jth1(1), wth(1), wth1(1)
        DO, j=2, 4
          WRITE (ndst,9024) j, jth(j), jth1(j), wth(j), wth1(j)
        END DO
#endif
        !
        ! 2.f Temp storage of data
        !
        IF ( snlm(iq).EQ.0. .AND. snlt(iq).LT.0. ) THEN
          jj     = 2
        ELSE
          jj     = 4
        END IF
        !
        DO j=1, jj
          SELECT CASE (j)
          CASE (2)
            jth(3) = -jth(3)
            jth(4) = -jth(4)
            jth1(3) = -jth1(3)
            jth1(4) = -jth1(4)
          CASE (3)
            jth     = -jth
            jth1    = -jth1
          CASE (4)
            jth(3) = -jth(3)
            jth(4) = -jth(4)
            jth1(3) = -jth1(3)
            jth1(4) = -jth1(4)
          CASE DEFAULT
          END SELECT
          !
          nqa    = nqa + 1
          tstore(nqa,ikd)%OFR  = jfr
          tstore(nqa,ikd)%OFR1 = jfr1
          tstore(nqa,ikd)%HFR  = wfr
          tstore(nqa,ikd)%HFR1 = wfr1
          tstore(nqa,ikd)%OTH  = jth
          tstore(nqa,ikd)%OTH1 = jth1
          tstore(nqa,ikd)%HTH  = wth
          tstore(nqa,ikd)%HTH1 = wth1
          IF ( jj .EQ. 2 ) THEN
            tstore(nqa,ikd)%CQD  = snlcd(iq) * 2.
            tstore(nqa,ikd)%CQS  = snlcs(iq) * 2.
          ELSE
            tstore(nqa,ikd)%CQD  = snlcd(iq)
            tstore(nqa,ikd)%CQS  = snlcs(iq)
          END IF
          auxf                 = ( wn0 * s0 ) / cg0
          tstore(nqa,ikd)%F1   = auxf * cg1 / ( wn1 * s1 )
          tstore(nqa,ikd)%F2   = auxf * cg2 / ( wn2 * s2 )
          tstore(nqa,ikd)%F3   = auxf * cg3 / ( wn3 * s3 )
          tstore(nqa,ikd)%F4   = auxf * cg4 / ( wn4 * s4 )
          !
        END DO
        !
        ! ... End loop 2.b
        !
      END DO
      !
      ! ... End loop 2.a
      !
    END DO
    !
#ifdef W3_T1
    WRITE (ndst,*)
#endif
#ifdef W3_T
    WRITE (ndst,9026) nqa, snlnq*4, nqd, nqs
#endif
    !
    ! 2.g Expanded spectral range
    !
    nthmax = nthmax + 1
    !
    nfrmin =  1  + ifrmin
    nfrmax = nfr + ifrmax - ifrmin
    nfrcut = nfr          - ifrmin
    nthexp = nth + 2*nthmax
    !
    nspmin = 1 + (nfrmin-1)*nthexp - nthmax
    nspmax = nfrmax * nthexp - nthmax
    nspmx2 = nfrcut * nthexp - nthmax
    !
#ifdef W3_T
    WRITE (ndst,9027) nfr, nfrmin, nfrmax, nfrcut, nth,          &
         1-nthmax, nth+nthmax, nthexp
#endif
    !
    ALLOCATE ( mpars(igrid)%SNLPS%FRQ(nfrmax),                      &
         mpars(igrid)%SNLPS%XSI(nfrmax) )
    frq     => mpars(igrid)%SNLPS%FRQ
    xsi     => mpars(igrid)%SNLPS%XSI
    !
    xsi(1:nfr) = sig(1:nfr)
    DO ifr=nfr+1, nfrmax
      xsi(ifr) = xsi(ifr-1) * xfr
    END DO
    frq    = xsi * tpiinv
    !
    ! 2.h Final storage
    !
    ALLOCATE ( mpars(igrid)%SNLPS%QST1(16,nqa,nkd),                 &
         mpars(igrid)%SNLPS%QST3(6,nqa,nkd),                  &
         mpars(igrid)%SNLPS%QST2(16,nqa,nkd) )
    qst1   => mpars(igrid)%SNLPS%QST1
    qst2   => mpars(igrid)%SNLPS%QST2
    qst3   => mpars(igrid)%SNLPS%QST3
    !
    ! 2.h.1 Basic data
    !
    DO ikd=1, nkd
      DO iqa=1, nqa
        !
        DO j=1, 4
          !
          qst1((j-1)*4+1,iqa,ikd) = tstore(iqa,ikd)%OTH (j) +       &
               tstore(iqa,ikd)%OFR (j) * nthexp
          qst1((j-1)*4+2,iqa,ikd) = tstore(iqa,ikd)%OTH1(j) +       &
               tstore(iqa,ikd)%OFR (j) * nthexp
          qst1((j-1)*4+3,iqa,ikd) = tstore(iqa,ikd)%OTH (j) +       &
               tstore(iqa,ikd)%OFR1(j) * nthexp
          qst1((j-1)*4+4,iqa,ikd) = tstore(iqa,ikd)%OTH1(j) +       &
               tstore(iqa,ikd)%OFR1(j) * nthexp
          !
          qst2((j-1)*4+1,iqa,ikd) = tstore(iqa,ikd)%HFR (j) *       &
               tstore(iqa,ikd)%HTH (j)
          qst2((j-1)*4+2,iqa,ikd) = tstore(iqa,ikd)%HFR (j) *       &
               tstore(iqa,ikd)%HTH1(j)
          qst2((j-1)*4+3,iqa,ikd) = tstore(iqa,ikd)%HFR1(j) *       &
               tstore(iqa,ikd)%HTH (j)
          qst2((j-1)*4+4,iqa,ikd) = tstore(iqa,ikd)%HFR1(j) *       &
               tstore(iqa,ikd)%HTH1(j)
          !
        END DO
        !
        qst3(1,iqa,ikd) = tstore(iqa,ikd)%F1
        qst3(2,iqa,ikd) = tstore(iqa,ikd)%F2
        qst3(3,iqa,ikd) = tstore(iqa,ikd)%F3
        qst3(4,iqa,ikd) = tstore(iqa,ikd)%F4
        qst3(5,iqa,ikd) = tstore(iqa,ikd)%CQD
        qst3(6,iqa,ikd) = tstore(iqa,ikd)%CQS
        !
      END DO
    END DO
    !
    IF ( nqd .GT. 0 ) qst3(5,:,:) = qst3(5,:,:) / real(nqd)
    IF ( nqs .GT. 0 ) qst3(6,:,:) = qst3(6,:,:) / real(nqs)
    !
    DEALLOCATE ( tstore )
    !
    ! 3.  Building quadruplet data base ---------------------------------- *
    !     For constructing interactions and diagonal from contributions
    !
    nthmx2 = 0
    ALLOCATE ( mpars(igrid)%SNLPS%QST4(16,nqa,nkd),                 &
         mpars(igrid)%SNLPS%QST5(16,nqa,nkd),                 &
         mpars(igrid)%SNLPS%QST6(16,nqa,nkd) )
    qst4   => mpars(igrid)%SNLPS%QST4
    qst5   => mpars(igrid)%SNLPS%QST5
    qst6   => mpars(igrid)%SNLPS%QST6
    ALLOCATE ( ast1(16,nqa,nkd), ast2(16,nqa,nkd) )
    !
    ! 3.a Loop over relative depths
    !
    s0     = sitmin * sqrt( grav / depth )  / xsit
    !
    DO ikd=1, nkd
      !
      s0     = s0 * xsit
      CALL wavnu2 ( s0, depth, wn0, cg0, 1.e-6, 25, ierr)
      !
      ! 3.b Loop over representative quadruplets
      !
      nqa    = 0
      !
      DO iq=1, snlnq
        !
#ifdef W3_T2
        WRITE (ndst,9030) ikd, iq, wn0*depth, s0*tpiinv, depth
#endif
        !
        off12  = snlm(iq)
        off34  = snll(iq)
        th12   = snlt(iq) * dera
        !
#ifdef W3_T2
        WRITE (ndst,9031) snlt(iq), off12, off34
#endif
        !
        ! 3.c Frequency indices
        !
        auxfr(1) = ( 1. + off12 )
        auxfr(2) = ( 1. - off12 )
        auxfr(3) = ( 1. + off34 )
        auxfr(4) = ( 1. - off34 )
        !
        DO j=1, 4
          jfr(j) = int( log(auxfr(j)) / xfrln )
          jfr1(j) = jfr(j) + 1 * sign(1.,auxfr(j)-1.)
          wfr(j) = (xfr**jfr1(j)-auxfr(j))/(xfr**jfr1(j)-xfr**jfr(j))
          wfr1(j) = 1. - wfr(j)
        END DO
        !
#ifdef W3_T2
        WRITE (ndst,9032) 1, jfr(1), jfr1(1), wfr(1), wfr1(1)
        DO, j=2, 4
          WRITE (ndst,9033) j, jfr(j), jfr1(j), wfr(j), wfr1(j)
        END DO
#endif
        !
        ! 3.d Loop over quadruplet components
        !
        DO jiq=1, 4
          !
          IF ( jiq .LE. 2 ) THEN
            wf     = -1.
          ELSE
            wf     =  1.
          END IF
          !
          ! 3.e Loop over frequency offsets, get directional offsets
          !
          DO jof=1, 2
            !
            IF ( jof .EQ. 1 ) THEN
              ifr    = -jfr(jiq)
              wfroff =  wfr(jiq)
            ELSE
              ifr    = -jfr1(jiq)
              wfroff =  wfr1(jiq)
            END IF
            !
            sioff  = s0 * xfr**ifr
            CALL wavnu2 ( sioff, depth, wn0, cg0, 1.e-6, 25, ierr)
            s1     = sioff * ( 1. + off12 )
            CALL wavnu2 ( s1, depth, wn1, cg1, 1.e-6, 25, ierr)
            s2     = sioff * ( 1. - off12 )
            CALL wavnu2 ( s2, depth, wn2, cg2, 1.e-6, 25, ierr)
            s3     = sioff * ( 1. + off34 )
            CALL wavnu2 ( s3, depth, wn3, cg3, 1.e-6, 25, ierr)
            s4     = sioff * ( 1. - off34 )
            CALL wavnu2 ( s4, depth, wn4, cg4, 1.e-6, 25, ierr)
            !
#ifdef W3_T2
            WRITE (ndst,9034) jiq, jof, ifr, wfroff, sioff/s0
#endif
            !
            IF ( th12 .LT. 0. ) THEN
              bb = 2. * wn0
            ELSE
              bb = wn1**2 + wn2**2 + 2.*wn1*wn2*cos(th12)
              bb = sqrt( max( bb , 0. ) )
            END IF
            !
            IF ( th12.LT.0. .AND. abs(off12).LE.1.e-4 ) THEN
              delth(1) = 0.
              delth(2) = 0.
            ELSE
              cc       = wn1
              aa       = wn2
              aux1     = (cc**2+bb**2-aa**2) / (2.*bb*cc)
              aux2     = (aa**2+bb**2-cc**2) / (2.*bb*aa)
              delth(1) = - acos( max( 0. , min( 1. , aux1 ) ) )
              delth(2) =   acos( max( 0. , min( 1. , aux2 ) ) )
            END IF
            cc       = wn3
            aa       = wn4
            aux1     = (cc**2+bb**2-aa**2) / (2.*bb*cc)
            aux2     = (aa**2+bb**2-cc**2) / (2.*bb*aa)
            delth(3) = - acos( max( 0. , min( 1. , aux1 ) ) )
            delth(4) =   acos( max( 0. , min( 1. , aux2 ) ) )
            !
#ifdef W3_T2
            WRITE (ndst,9035) delth(:) * rade
#endif
            !
            aux1    = delth(jiq) / dth
            jth(jiq) = int(aux1)
            jth1(jiq) = jth(jiq) + 1 * sign(1.,delth(jiq))
            wth1(jiq) = abs(aux1) - real(abs(jth(jiq)))
            wth(jiq) = 1. - wth1(jiq)
            !
            nthmx2 = max( nthmx2 , abs(jth1(jiq)) )
            !
#ifdef W3_T2
            WRITE (ndst,9036) jiq, jth(jiq), jth1(jiq),         &
                 wth(jiq), wth1(jiq)
#endif
            !
            ! 3.f Loop over quadruplet realizations
            !
            IF ( snlm(iq).EQ.0. .AND. snlt(iq).LT.0. ) THEN
              jj     = 2
            ELSE
              jj     = 4
            END IF
            !
            DO jqr=1, jj
              !
              SELECT CASE (jqr)
              CASE (2)
                jth(3) = -jth(3)
                jth(4) = -jth(4)
                jth1(3) = -jth1(3)
                jth1(4) = -jth1(4)
              CASE (3)
                jth     = -jth
                jth1    = -jth1
              CASE (4)
                jth(3) = -jth(3)
                jth(4) = -jth(4)
                jth1(3) = -jth1(3)
                jth1(4) = -jth1(4)
              CASE DEFAULT
                jth     = -jth
                jth1    = -jth1
              END SELECT
              !
              ist    = (jiq-1)*4 + (jof-1)*2 + 1
              ast1(ist,nqa+jqr,ikd) = ifr
              ast2(ist,nqa+jqr,ikd) = jth(jiq)
              qst5(ist,nqa+jqr,ikd) = wf * ( wfroff * wth(jiq) )
              qst6(ist,nqa+jqr,ikd) = wf * ( wfroff * wth(jiq) )**2
              ist    = ist + 1
              ast1(ist,nqa+jqr,ikd) = ifr
              ast2(ist,nqa+jqr,ikd) = jth1(jiq)
              qst5(ist,nqa+jqr,ikd) = wf * ( wfroff * wth1(jiq) )
              qst6(ist,nqa+jqr,ikd) = wf * ( wfroff * wth1(jiq) )**2
              !
              ! ... End loop 3.f
              !
            END DO
            !
            ! ... End loop 3.e
            !
          END DO
          !
          ! ... End loop 3.d
          !
        END DO
        !
#ifdef W3_T3
        DO jqr=1, jj
          WRITE (ndst,9037) ikd, nqa+jqr
          DO ist=1, 16
            WRITE (ndst,9038) ist, ast1(ist,nqa+jqr,ikd),       &
                 ast2(ist,nqa+jqr,ikd),       &
                 qst5(ist,nqa+jqr,ikd),       &
                 qst6(ist,nqa+jqr,ikd)
          END DO
        END DO
#endif
        !
        ! ... End loop 3.b
        !
        nqa    = nqa + jj
        !
      END DO
      !
      ! ... End loop 3.a
      !
    END DO
    !
    ! 3.g Finalize storage
    !
    qst4 = ast1*nthexp + ast2
    !
    IF ( nthmax .LT. nthmx2 ) GOTO 810
    IF ( nqa .NE. SIZE(ast1(1,:,1)) ) GOTO 811
    !
    DEALLOCATE ( ast1, ast2 )
    !
    RETURN
    !
    ! Error escape locations
    !
800 CONTINUE
    WRITE (ndse,1000) lammax, delthm
    CALL extcde ( 1000 )
    !
801 CONTINUE
    WRITE (ndse,1001) off12, off34
    CALL extcde ( 1001 )
    !
802 CONTINUE
    WRITE (ndse,1002) off12, off34, snlt(iq),                       &
         minlam(off12,snlt(iq)), maxlam(off12,snlt(iq))
    CALL extcde ( 1002 )
    !
810 CONTINUE
    WRITE (ndse,1010) nthmax, nthmx2
    CALL extcde ( 1010 )
    !
811 CONTINUE
    WRITE (ndse,1011) nqa, SIZE(ast1(1,:,1))
    CALL extcde ( 1011 )
    !
    RETURN
    !
    ! Formats
    !
1000 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/                &
         '     PARAMETER OUT OF RANGE '/                        &
         '     LAMMAX, DELTHM :', 2e12.4/)
1001 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/                &
         '     PARAMETER OUT OF RANGE '/                        &
         '     MU, LAMBDA :', 2e12.4/)
1002 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/                &
         '     PARAMETER OUT OF RANGE '/                        &
         '     MU, LAMBDA, TH12 :',3e12.4/                      &
         '     LAMBDA RANGE     :',2e12.4)
1010 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/                &
         '     NTHMAX LESS THAN NTHMX2 :', 2i8/)
1011 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNL3 :'/                &
         '     NQA INCONSISTENT :', 2i8/)
    !
#ifdef W3_T
9010 FORMAT (/' TEST INSNL3: NKD, KDMIN/MAX/X : ',i8,3f10.4)
#endif
    !
#ifdef W3_T1
9020 FORMAT (/' TEST INSNL3: IKD, IQ, KD, F, D: ',2i8,2f10.4,f10.2)
9021 FORMAT (/' TEST INSNL3: TH12             : ',3x,f8.2/       &
         '              OFF12, OFF34     : ',3x,2f8.2/      &
         '              CD, CS           : ',3x,2e10.2)
9022 FORMAT ( '              ANGLES (DEGR)    : ',1x,4f8.2)
9023 FORMAT ( '              FREQUENCY IND.   : ',1x,3i4,2f6.2)
9024 FORMAT ( '                               : ',1x,3i4,2f6.2)
9025 FORMAT ( '              DIRECTIONAL IND. : ',1x,3i4,2f6.2)
#endif
#ifdef W3_T
9026 FORMAT ( ' TEST INSNL3: FILLING FIRST DATA TABLES :'/        &
         '              NQA AND MAXIMUM  : ',2i8/            &
         '              NQD AND NQS      : ',2i8)
9027 FORMAT ( '              NFR, MIN/MAX/CUT : ',4i8/            &
         '              NTH, MIN/MAX/EXP : ',4i8)
#endif
    !
#ifdef W3_T2
9030 FORMAT (/' TEST INSNL3: IKD, IQ, KD, F, D: ',2i8,2f10.4,f10.2)
9031 FORMAT (/' TEST INSNL3: TH12             : ',3x,f8.2/       &
         '              OFF12, OFF34     : ',3x,2f8.2)
9032 FORMAT ( '              FREQUENCY IND.   : ',1x,3i4,2f6.2)
9033 FORMAT ( '                               : ',1x,3i4,2f6.2)
9034 FORMAT ( '              J,J,J, W, SIn    : ',1x,3i4,2f6.2)
9035 FORMAT ( '              ANGLES (DEGR)    : ',3x,4f8.2)
9036 FORMAT ( '              DIRECTIONAL IND. : ',1x,3i4,2f6.2)
#endif
#ifdef W3_T3
9037 FORMAT (/' TEST INSNL3: STORAGE ARRAYS FOR IKD, IQA =',2i6)
9038 FORMAT (23x,3i4,3f8.3)
#endif
    !/
    !/ Embedded subroutines
    !/
  CONTAINS
    !/ ------------------------------------------------------------------- /
 
    REAL FUNCTION MINLAM ( MU, THETA )
      !/
      !/                  +-----------------------------------+
      !/                  | WAVEWATCH-III           NOAA/NCEP |
      !/                  |           H. L. Tolman            |
      !/                  |                        FORTRAN 90 |
      !/                  | Last update :         28-Jan-2004 |
      !/                  +-----------------------------------+
      !/
      !/    28-Jan-2009 : Origination.
      !/
      !  1. Purpose :
      !
      !     Calculate minimum allowed lambda for quadruplet configuration.
      !
      !  3. Parameters :
      !
      !     Parameter list
      !     ----------------------------------------------------------------
      !       MU, THETA  Real   Quadruplet parameters, theta in degree.
      !     ----------------------------------------------------------------
      !
      ! 10. Source code :
      !
      !/ ------------------------------------------------------------------- /
      IMPLICIT NONE
      !/
      !/ Parameter list
      !/
      REAL, INTENT(IN)        :: mu, theta
      !/
      !/ Local parameters
      !/
      REAL                    :: muloc, thetar, bb, aux
      !/
      !/ ------------------------------------------------------------------- /
      !/
      IF ( theta .LT. 0. ) THEN
        minlam = 0.
      ELSE
        muloc  = max( 0. , min( 1., mu ) )
        thetar = theta * atan(1.) / 45.
        bb     = (1.+muloc)**4 + (1.-muloc)**4 +                    &
             2. * (1.+muloc)**2 * (1.-muloc)**2 * cos(thetar)
        bb     = sqrt( max( bb , 0. ) )
        aux    = max( 0. , 0.5*bb-1. )
        minlam = sqrt( aux )
      END IF
      !
      RETURN
      !/
      !/ End of MINLAM ----------------------------------------------------- /
      !/
    END FUNCTION minlam
    !/ ------------------------------------------------------------------- /
 
    REAL function maxlam ( mu, theta )
      !/
      !/                  +-----------------------------------+
      !/                  | WAVEWATCH-III           NOAA/NCEP |
      !/                  |           H. L. Tolman            |
      !/                  |                        FORTRAN 90 |
      !/                  | Last update :         28-Jan-2004 |
      !/                  +-----------------------------------+
      !/
      !/    28-Jan-2009 : Origination.
      !/
      !  1. Purpose :
      !
      !     Calculate minimum allowed lambda for quadruplet configuration.
      !
      !  3. Parameters :
      !
      !     Parameter list
      !     ----------------------------------------------------------------
      !       MU, THETA  Real   Quadruplet parameters, theta in degree.
      !     ----------------------------------------------------------------
      !
      ! 10. Source code :
      !
      !/ ------------------------------------------------------------------- /
      IMPLICIT NONE
      !/
      !/ Parameter list
      !/
      REAL, INTENT(IN)        :: mu, theta
      !/
      !/ Local parameters
      !/
      REAL                    :: muloc, thetar, bb, aux
      !/
      !/ ------------------------------------------------------------------- /
      !/
      IF ( theta .LT. 0. ) THEN
        maxlam = 0.5
      ELSE
        muloc  = max( 0. , min( 1., mu ) )
        thetar = theta * atan(1.) / 45.
        bb     = (1.+muloc)**4 + (1.-muloc)**4 +                    &
             2. * (1.+muloc)**2 * (1.-muloc)**2 * cos(thetar)
        bb     = sqrt( max( bb , 0. ) )
        maxlam = 0.25 * bb
      END IF
      !
      RETURN
      !/
      !/ End of MAXLAM ----------------------------------------------------- /
      !/
    END FUNCTION maxlam
    !/
    !/ End of INSNL3 ----------------------------------------------------- /
    !/
  END SUBROUTINE insnl3
  !/
  !/ End of module W3SNL3MD -------------------------------------------- /
  !/
END MODULE w3snl3md