w3snlsmd.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3snlsmd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH-III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         13-Jul-2012 |
  !/                  +-----------------------------------+
  !/
  !/    04-Aug-2008 : Origination in research model.      ( version 3.13 )
  !/    27-Sep-2010 : Added to svn repository.            ( version 3.15 )
  !/    13-Jul-2012 : Moved from version 3.15 to 4.08.    ( version 4.08 )
  !/
  !/    Copyright 2009-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 :
  !
  !     Nonlinear interaction based `smoother' for high frequencies.
  !
  !  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.
  !      SITMIN    Real  Private  Minimum nondimensional radian
  !                               frequency in table.
  !      XSIT      Real  Private  Corresponding incremet factor.
  !      ABMAX     R.P.  Public   Maximum value of a34, b3 and b4.
  !     ----------------------------------------------------------------
  !
  !     Variables in W3GDATMD :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      CNLSA     Real  Public   a34 in quadruplet definition.
  !      CNLSC     Real  Public   C in Snl definition.
  !      CNLSFM    Real  Public   Maximum relative spectral change.
  !      CNLSC1/3  Real  Public   Constant in frequency filter.
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SNLS    Subr. Public   Nonlinear 'smoother' algorithm.
  !      EXPAND    Subr. W3SNLS   Expand spectrum for indirect address.
  !      INSNLS    Subr. Public   Initialization routine.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      WAVNU1    Subr. W3DISPMD Solve dispersion relation.
  !      WAVNU2    Subr. W3DISPMD Solve dispersion relation.
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !      EXTCDE    Subr. W3SERVMD Program abort.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !     !/S    Enable subroutine tracing.
  !     !/T    Enable test output.
  !
  !  7. Source code :
  !/
  !/ ------------------------------------------------------------------- /
  !/
  INTEGER, PRIVATE, PARAMETER :: NKD = 100
  REAL, PRIVATE, PARAMETER    :: KDMIN = 0.25 ,  kdmax = 10.
  REAL, PRIVATE               :: SITMIN, XSIT
  !
  REAL, PARAMETER             :: abmax = 0.25
  !
  PUBLIC
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3snls (  A, CG, WN, DEPTH, UABS, DT, SNL, AA )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH-III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         04-Aug-2008 |
    !/                  +-----------------------------------+
    !/
    !/    04-Aug-2008 : Origination.                        ( version 3.13 )
    !/
    !  1. Purpose :
    !
    !     High-frequeny filter based on the nonlinear interactions for
    !     an uresolved quadruplet.
    !
    !  2. Method :
    !
    !     Compute interactions for a quadruplet that is not resolved by
    !     the discrete spectral rsolution, and then reduces to a simple
    !     five-point stencil. Furthermore interactions are filtered by
    !     frequency to allow for high-frequency impact only, and the
    !     integration schem is embedded, and reduces to a filter technique
    !     for large time steps or strong 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.
    !       UABS    Real  I   Wind speed (m/s).
    !       DT      Real  I   Numerical time step (s).
    !       SNL     R.A.  O   Nonlinear source term.                (Opt)
    !       AA      R.A.  O   Averaged spectrum.                    (Opt)
    !     ----------------------------------------------------------------
    !       Note: A and AA may safely be same array/address.
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      WAVNU1    Subr. W3DISPMD Solve dispersion relation.
    !      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 :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S    Enable subroutine tracing.
    !     !/T    Enable test output.
    !     !/T1   Test output frequency filter.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nfr => nk, nth, sig, xfr, fachfa, dth,      &
         nthx, nfrx, nspl, nsph, snsst, cnlsc,       &
         cnlsfm, cnlsc1, cnlsc2, cnlsc3
    USE w3odatmd, ONLY: ndst, ndse
    !
    USE w3dispmd, ONLY: wavnu1
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_T2
    USE w3arrymd, ONLY: prt2ds
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)            :: A(NTH,NFR), CG(NFR), WN(NFR),    &
         DEPTH, UABS, DT
    REAL, INTENT(OUT), OPTIONAL :: SNL(NTH,NFR), AA(NTH,NFR)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IFR, IFRMIN, ITH, IFRMN2,            &
         IKD, JKD(0:NFR+2), ISPX0, ISPX
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: SIGP, CP, CM, XL, XH, EL, EH, DENOM, &
         SIT, XSITLN, MC, F3A,  F3B, F3C,     &
         F4A, F4B, F4C, F00, F31, F32, F41,   &
         F42, AUXB, AUX11, AUX21, AUX12,      &
         AUX22, FC1, FC2, FC3, FC4
    REAL                    :: XSI(NFR+2), XWN(NFR+2), XCG(NFR+2),  &
         UP(NSPL:NSPH), UN(NSPL:NSPH),        &
         E1(0:NFR+2), FILTFP(NFR+2),          &
         FPROP(NFR+2), DS1(NSPL:NSPH),        &
         DS2(NSPL:NSPH), DS3(NSPL:NSPH),      &
         DA1(NSPL:NSPH), DA2(NSPL:NSPH),      &
         DA3(NSPL:NSPH)
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SNLS')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) depth, uabs, dt
#endif
    !
    ! 1.  Initializations ------------------------------------------------ *
    ! 1.a Expanded frequency range
    !
    xsi(1:nfr) = sig(1:nfr)
    xwn(1:nfr) = wn
    xcg(1:nfr) = cg
    !
    xsi(nfr+1) = xsi(nfr) * xfr
    CALL wavnu1 ( xsi(nfr+1), depth, xwn(nfr+1), xcg(nfr+1) )
    xsi(nfr+2) = xsi(nfr+1) * xfr
    CALL wavnu1 ( xsi(nfr+2), depth, xwn(nfr+2), xcg(nfr+2) )
    !
    ! 1.b Expanded psuedo spectrum
    !
    CALL expand ( up, un )
    !
    ! 1.c Get relevant spectral peak frequency
    !
#ifdef W3_T1
    e1     = -1.
#endif
    sigp   = - tpi
    xl     = 1./xfr - 1.
    xh     =  xfr - 1.
    !
    ! 1.c.1 Wind too weak
    !
    IF ( uabs .LT. xsi(nfr)/xwn(nfr) ) THEN
      sigp   = grav / max( 0.01 , uabs )
    ELSE
      !
      ! 1.c.2 Compute 1D spectrum
      !
      e1(nfr+2) = sum(a(:,nfr)) * fachfa**2 * xsi(nfr+2)          &
           / xcg(nfr+2) * tpi * dth
      e1(nfr+1) = sum(a(:,nfr)) * fachfa    * xsi(nfr+1)          &
           / xcg(nfr+1) * tpi * dth
      !
      DO ifr=nfr, 1, -1
        e1(ifr) = sum(a(:,ifr)) * xsi(ifr) / xcg(ifr) * tpi * dth
        !
        ! 1.c.3 Reached PM frequency
        !
        IF ( uabs .LT. xsi(ifr)/xwn(ifr) ) THEN
          cp     = xsi(ifr)/xwn(ifr)
          cm     = xsi(ifr+1)/xwn(ifr+1)
          sigp   = xsi( ifr ) * (uabs-cm)/(cp-cm) +             &
               xsi(ifr+1) * (cp-uabs)/(cp-cm)
          EXIT
          !
        ELSE IF ( e1(ifr) .LT. e1(ifr+1) ) THEN
          !
          ! 1.c.4 Reached first peak
          !
          el     = e1(ifr  ) - e1(ifr+1)
          eh     = e1(ifr+2) - e1(ifr+1)
          denom  = xl*eh - xh*el
          sigp   = xsi(ifr+1) * (1.+0.5*(xl**2*eh-xh**2*el)     &
               / sign( max(abs(denom),1.e-15) , denom ) )
          EXIT
        ENDIF
        !
        ! ... End loop 1.c.2
        !
      END DO
      !
      ! 1.c.5 Nothing found
      !
      IF ( sigp .LT. 0. ) THEN
        !
        ! 1.c.5.a No energy there
        !
        IF ( e1(1) .EQ. 0. ) THEN
          sigp   = 2. * sig(nfr)
          !
          ! 1.c.5.b Peak at low boundary
          !
        ELSE
          sigp   = xsi(1)
        END IF
      END IF
      !
    END IF
    !
    ! 1.d Set up filter function etc.
    !
    xsitln = log(xsit)
    ifrmin =  1
    jkd    =  1
#ifdef W3_T1
    filtfp = -1.
#endif
    !
    DO ifr=nfr+2, 1, -1
      !
      filtfp(ifr) = exp(-cnlsc1/(xsi(ifr)/(cnlsc2*sigp))**cnlsc3)
      fprop(ifr) = filtfp(ifr) * cnlsc * xwn(ifr)**8 *             &
           xsi(ifr)**4 / tpi**9 / xcg(ifr)
      sit      = xsi(ifr) * sqrt(depth/grav)
      ikd      = 1 + nint( ( log(sit) - log(sitmin) ) / xsitln )
      jkd(ifr) = max( 1 , min(ikd,nkd) )
      !
      IF ( filtfp(ifr) .LT. 1.e-10 ) THEN
        ifrmin = ifr
        EXIT
      END IF
      !
    END DO
    !
    ifrmn2 = max( 1 , ifrmin - 1 )
    sit    = xsi(ifrmn2) * sqrt(depth/grav)
    ikd    = 1 + nint( ( log(sit) - log(sitmin) ) / xsitln )
    jkd(ifrmn2) = max( 1 , min(ikd,nkd) )
    !
#ifdef W3_T
    WRITE (ndst,9010) ifrmin, sigp * tpiinv
#endif
#ifdef W3_T1
    WRITE (ndst,9011)
    DO ifr=1, nfr
      WRITE (ndst,9012) ifr, xsi(ifr)/tpi, xsi(ifr)/xwn(ifr),   &
           e1(ifr), filtfp(ifr)
    END DO
#endif
    !
    ! 1.e Initialize arrays
    !
    !
    ! 2.  Compute base interactions -------------------------------------- *
    ! 2.a Loop over frequencies
    !
    DO ifr=ifrmin, nfr+1
      !
      ispx0  = (ifr-1)*nthx
      ikd    = jkd(ifr)
      !
      mc     = snsst( 1,ikd)
      f3a    = snsst( 2,ikd)
      f3b    = snsst( 3,ikd)
      f3c    = snsst( 4,ikd)
      f4a    = snsst( 5,ikd)
      f4b    = snsst( 6,ikd)
      f4c    = f3c
      !
      ! 2.b Loop over directions
      !
      DO ith=1, nth
        !
        ispx   = ispx0 + ith
        !
        f00    = up(ispx)
        f31    = up(ispx)*f3a + up(ispx+1)*f3b + up(ispx+nthx)*f3c
        f41    = up(ispx)*f4a + up(ispx-1)*f4b + up(ispx-nthx)*f4c
        f32    = up(ispx)*f3a + up(ispx-1)*f3b + up(ispx+nthx)*f3c
        f42    = up(ispx)*f4a + up(ispx+1)*f4b + up(ispx-nthx)*f4c
        !
        ds1(ispx) = fprop(ifr) * (f00**2*(f31+f41)-2.*f00*f31*f41)
        ds2(ispx) = fprop(ifr) * (f00**2*(f32+f42)-2.*f00*f32*f42)
        !
        aux11  = dt * ds1(ispx)
        aux21  = dt * ds2(ispx)
        auxb   = cnlsfm * filtfp(ifr) * max(1.e-10,un(ispx)) /   &
             max( 1.e-10 , abs(aux11)+abs(aux21) ) / mc
        aux12  = auxb * abs(aux11)
        aux22  = auxb * abs(aux21)
        !
        ! Expensive but more smooth limiter
        !
        !         DA1(ISPX) = AUX12 * TANH(AUX11/MAX(1.E-10,AUX12))
        !         DA2(ISPX) = AUX22 * TANH(AUX21/MAX(1.E-10,AUX22))
        !
        ! Crude but cheaper limiter
        !
        da1(ispx) = max( -aux12 , min( aux11 , aux12 ) )
        da2(ispx) = max( -aux22 , min( aux21 , aux22 ) )
        !
      END DO
      !
      ! ... End loop 2.b
      !
    END DO
    !
    ! 2.c Complete expanded arrays
    !
    ! ... End loop 2.a
    !
    ! 3.  Compute source term if requested ------------------------------- *
    ! 3.a Check for request
    !
    IF ( PRESENT(snl) ) THEN
#ifdef W3_T
      WRITE (ndst,9030) 'YES/--'
#endif
      !
      ! 3.b Initializations
      !
      snl(:,1:ifrmn2-1) = 0.
      !
      ds1(nspl:ifrmn2*nthx-1) = 0.
      ds2(nspl:ifrmn2*nthx-1) = 0.
      ds3(nspl:ifrmn2*nthx-1) = 0.
      !
      ispx  = ifrmn2*nthx
      ds1(ispx+nth+1:nsph:nthx) = ds1(ispx+ 1 :nsph:nthx)
      ds1(ispx      :nsph:nthx) = ds1(ispx+nth:nsph:nthx)
      ds2(ispx+nth+1:nsph:nthx) = ds2(ispx+ 1 :nsph:nthx)
      ds2(ispx      :nsph:nthx) = ds2(ispx+nth:nsph:nthx)
      ds3(ifrmn2*nthx:nsph)     = ds1(ifrmn2*nthx:nsph)  +        &
           ds2(ifrmn2*nthx:nsph)
      !
      ! 3.c Loop over frequencies
      !
      DO ifr=ifrmn2, nfr
        !
        ispx0  = (ifr-1)*nthx
        ikd    = jkd(ifr)
        !
        fc1    = - snsst(1,ikd)
        fc2    =   snsst(4,ikd)
        fc3    =   snsst(3,ikd)
        fc4    =   snsst(6,ikd)
        !
        ! 3.d Loop over directions
        !
        DO ith=1, nth
          ispx         = ispx0 + ith
          snl(ith,ifr) = fc1 *             ds3(   ispx  )         &
               +    fc2 * ( ds3(ispx-nthx) + ds3(ispx+nthx) )       &
               +    fc3 * ( ds1(ispx-  1 ) + ds2(ispx+  1 ) )       &
               +    fc4 * ( ds1(ispx+  1 ) + ds2(ispx-  1 ) )
          !
        END DO
        !
        ! ... End loop 3.d
        !
      END DO
      !
      ! ... End loop 3.c
      !
#ifdef W3_T
    ELSE
      WRITE (ndst,9030) '---/NO'
#endif
    END IF
    !
    ! 4.  Compute filtered spectrum if requested ------------------------- *
    ! 4.a Check for request
    !
    IF ( PRESENT(aa) ) THEN
#ifdef W3_T
      WRITE (ndst,9040) 'YES/--'
#endif
      !
      ! 4.b Initializations
      !
      aa(:,1:ifrmn2-1) = a(:,1:ifrmn2-1)
      !
      da1(nspl:ifrmn2*nthx-1) = 0.
      da2(nspl:ifrmn2*nthx-1) = 0.
      da3(nspl:ifrmn2*nthx-1) = 0.
      !
      ispx  = ifrmn2*nthx
      da1(ispx+nth+1:nsph:nthx) = da1(ispx+ 1 :nsph:nthx)
      da1(ispx      :nsph:nthx) = da1(ispx+nth:nsph:nthx)
      da2(ispx+nth+1:nsph:nthx) = da2(ispx+ 1 :nsph:nthx)
      da2(ispx      :nsph:nthx) = da2(ispx+nth:nsph:nthx)
      da3(ifrmn2*nthx:nsph)     = da1(ifrmn2*nthx:nsph)  +        &
           da2(ifrmn2*nthx:nsph)
      !
      ! 4.c Loop over frequencies
      !
      DO ifr=ifrmn2, nfr
        !
        ispx0  = (ifr-1)*nthx
        ikd    = jkd(ifr)
        !
        fc1    = - snsst(1,ikd)
        fc2    =   snsst(4,ikd)
        fc3    =   snsst(3,ikd)
        fc4    =   snsst(6,ikd)
        !
        ! 4.d Loop over directions
        !
        DO ith=1, nth
          ispx         = ispx0 + ith
          aa(ith,ifr) = max( 0. , a(ith,ifr) +                   &
               fc1 *   da3(ispx)                                  &
               + fc2 * ( da3(ispx-nthx) + da3(ispx+nthx) )          &
               + fc3 * ( da1(ispx-  1 ) + da2(ispx+  1 ) )          &
               + fc4 * ( da1(ispx+  1 ) + da2(ispx-  1 ) ) )
        END DO
        !
        ! ... End loop 4.d
        !
      END DO
      !
      ! ... End loop 4.c
      !
#ifdef W3_T
    ELSE
      WRITE (ndst,9040) '---/NO'
#endif
    END IF
    !
#ifdef W3_T
    stop
#endif
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (/' TEST W3SNLS: DEPTH, UABS, DT :',f9.2,f7.2,f7.2)
9010 FORMAT ( '              IFRMIN, FP  :',i4,f8.4)
9030 FORMAT ( ' TEST W3SNLS: SOURCE TERM REQUESTED : ',a)
9040 FORMAT ( ' TEST W3SNLS: AVERAGING REQUESTED   : ',a)
#endif
#ifdef W3_T1
9011 FORMAT ( ' TEST W3SNLS: IFR, FR, C, E1, FILT :')
9012 FORMAT (13x,i4,f10.4,2f10.2,f10.4)
#endif
    !/
    !/ Embedded subroutines
    !/
  CONTAINS
    !/ ------------------------------------------------------------------- /
    SUBROUTINE expand ( PSPC, SPEC )
      !/
      !/                  +-----------------------------------+
      !/                  | WAVEWATCH-III           NOAA/NCEP |
      !/                  |           H. L. Tolman            |
      !/                  |                        FORTRAN 90 |
      !/                  | Last update :         23-Jul-2008 |
      !/                  +-----------------------------------+
      !/
      !  1. Purpose :
      !
      !     Expand spectrum to simplify indirect addressing.
      !
      !  3. Parameters :
      !
      !     Parameter list
      !     ----------------------------------------------------------------
      !       PSPC    R.A.  O   Expanded spectrum.
      !       SPEC    R.A.  O   Expanded spectrum.
      !     ----------------------------------------------------------------
      !
      ! 10. Source code :
      !
      !/ ------------------------------------------------------------------- /
      IMPLICIT NONE
      !/
      !/ Parameter list
      !/
      REAL, INTENT(OUT)       :: PSPC(0:NTH+1,0:NFR+2),               &
           SPEC(0:NTH+1,0:NFR+2)
      !/
      !/ Local parameters
      !/
      INTEGER                 :: IFR, ITH
      !/
      !/ ------------------------------------------------------------------- /
      !
      spec(:,0) = 0.
      !
      spec(1:nth,1:nfr) = a
      spec(1:nth,nfr+1) = spec(1:nth,nfr) * fachfa
      spec(1:nth,nfr+2) = spec(1:nth,nfr+1) * fachfa
      !
      spec(nth+1,1:nfr+2) = spec( 1 ,1:nfr+2)
      spec(  0  ,1:nfr+2) = spec(nth,1:nfr+2)
      !
      DO ifr=1, nfr+2
        pspc(:,ifr) = spec(:,ifr) / xwn(ifr)
      END DO
      !
      RETURN
      !/
      !/ End of EXPAND ----------------------------------------------------- /
      !/
    END SUBROUTINE expand
    !/
    !/ End of W3SNLS ----------------------------------------------------- /
    !/
  END SUBROUTINE w3snls
  !/ ------------------------------------------------------------------- /
  SUBROUTINE insnls
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH-III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         04-Aug-2008 |
    !/                  +-----------------------------------+
    !/
    !/    04-Aug-2008 : Origination.                        ( version 3.13 )
    !/
    !  1. Purpose :
    !
    !     Initializations for the Snl / filter source term for high
    !     frequencies.
    !
    !  2. Method :
    !
    !     Precompute weight functions and store in array.
    !
    !  3. Parameters :
    !
    !     No parameter list.
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      WAVNU2    Subr. W3DISPMD Solve dispersion relation.
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !      EXTCDE    Subr. W3SERVMD Program abort.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3IOGR    Subr. W3IOGRMD Process model definition file.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     - Check a34, b4 and b5 against MAXAB to assure that the values
    !       are consistent with a reduced 5-point stencil for unresolved
    !       quadruplets. a34 is checked in ww3_grid, b3 and b4 are not.
    !
    !  7. Remarks :
    !
    !     - Small quadruplet compared to grid size reduces interactions
    !       so that distribution of results is purely local. This results
    !       in a much simpler model initialization than for the general
    !       MDIA.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S    Enable subroutine tracing.
    !     !/T    Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3odatmd, ONLY: ndst, ndse
    USE w3gdatmd,       nfr => nk, a34 => cnlsa
    !
    USE w3dispmd, ONLY: wavnu2
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IKD, IERR
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: DEPTH, SITMAX, OFF, S0, WN0, CG0,    &
         S3, WN3, CG3, S4, WN4, CG4, WN12,    &
         DT3, DT4, B3, B4
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'INSNLS')
#endif
    !
    ! 1.  Initializations ------------------------------------------------ *
    ! 1.a Set up relative depths
    !
    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
    !
    ! 1.b Set up quadruplet
    !
    off    = (xfr-1.) * a34
    !
    ! 1.c Set up storage
    !
    nthx   = nth + 2
    nfrx   = nfr + 2
    nspl   = - nthx
    nsph   = nfrx*nthx - 1
    !
    ALLOCATE ( mpars(igrid)%SNLPS%SNSST(6,nkd) )
    snsst  => mpars(igrid)%SNLPS%SNSST
    !
    ! 2.  Building quadruplet data base ---------------------------------- *
    !     For quadruplet and interaction strength evaluation
    !
    s0     = sitmin * sqrt( grav / depth )  / xsit
    !
    ! 2.a Loop over relative depths
    !
    DO ikd=1, nkd
      !
      ! 2.b Base quadruplet set up
      !
      s0     = s0 * xsit
      s3     = ( 1. + off ) * s0
      s4     = ( 1. - off ) * s0
      !
      CALL wavnu2 ( s0, depth, wn0, cg0, 1.e-6, 25, ierr)
      CALL wavnu2 ( s3, depth, wn3, cg3, 1.e-6, 25, ierr)
      CALL wavnu2 ( s4, depth, wn4, cg4, 1.e-6, 25, ierr)
      !
#ifdef W3_T
      WRITE (ndst,9020) ikd, wn0*depth, s0*tpiinv, depth
#endif
      !
      ! 2.c Offset angles
      !
      wn12   = 2. * wn0
      dt3    = acos( (wn3**2+wn12**2-wn4**2) / (2.*wn12*wn3) )
      dt4    = acos( (wn4**2+wn12**2-wn3**2) / (2.*wn12*wn4) )
      !
      b3     = dt3 / dth
      b4     = dt4 / dth
      !
#ifdef W3_T
      WRITE (ndst,9021) a34, b3, b4, dt3*rade, dt4*rade
#endif
      !
      IF ( a34.GT.abmax .OR. b3.GT.abmax .OR. b4.GT.abmax .OR.     &
           a34.LT.0. .OR. b3.LT.0. .OR. b4.LT.0. ) GOTO 801
      !
      ! 2.d Store weights
      !
      snsst( 1,ikd) = 2.*a34 + b3 + b4
      snsst( 2,ikd) = 1. - a34 - b3
      snsst( 3,ikd) = b3
      snsst( 4,ikd) = a34
      snsst( 5,ikd) = 1. - a34 - b4
      snsst( 6,ikd) = b4
      !
      ! ... End loop 2.a
      !
    END DO
    !
    RETURN
    !
    ! Error escape locations
    !
801 CONTINUE
    WRITE (ndse,1001) a34, b3, b4
    CALL extcde (1001)
    !
    ! Formats
    !
1001 FORMAT (/' *** WAVEWATCH-III ERROR IN INSNLS :'/                &
         '     PARAMETER FORCED OUT OF RANGE '/                 &
         '     A34, B3, B4 :', 3f10.4/)
    !
#ifdef W3_T
9010 FORMAT (/' TEST INSNLS: NKD, KDMIN/MAX/X :',i5,3f10.4)
9020 FORMAT ( '              IKD, KD, F, D    :',i5,3f10.4)
9021 FORMAT ( '              A34, B3,B4, TH3/4:',3f7.3,2f6.2)
#endif
    !/
    ! /End of INSNLS ------------------------------------------------------/
    !/
  END SUBROUTINE insnls
  !/
  !/ End of module W3SNLSMD -------------------------------------------- /
  !/
END MODULE w3snlsmd