#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)
#endif
!
#ifdef W3_T
 9010 FORMAT ( '              IFRMIN, FP  :',I4,F8.4)
#endif
#ifdef W3_T1
 9011 FORMAT ( ' TEST W3SNLS: IFR, FR, C, E1, FILT :')
 9012 FORMAT (13X,I4,F10.4,2F10.2,F10.4)
#endif
!
#ifdef W3_T
 9030 FORMAT ( ' TEST W3SNLS: SOURCE TERM REQUESTED : ',A)
 9040 FORMAT ( ' TEST W3SNLS: AVERAGING REQUESTED   : ',A)
#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