w3snl1md.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
 
MODULE w3snl1md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         28-Feb-2023 |
  !/                  +-----------------------------------+
  !/
  !/    04-Feb-2000 : Origination.                        ( version 2.00 )
  !/    09-May-2002 : Switch clean up.                    ( version 2.21 )
  !/    24-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
  !/    03-Sep-2012 : Clean up of test output T0, T1      ( version 4.07 )
  !/    28-Feb-2023 : Adds GQM separate routines          ( version 7.07 )
  !/
  !/    Copyright 2009 National Weather Service (NWS),
  !/       National Oceanic and Atmospheric Administration.  All rights
  !/       reserved.  WAVEWATCH III is a trademark of the NWS.
  !/       No unauthorized use without permission.
  !/
  !  1. Purpose :
  !
  !     Bundles routines calculate nonlinear wave-wave interactions
  !     according to the Discrete Interaction Approximation (DIA) of
  !     Hasselmann et al. (JPO, 1985).
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SNL1    Subr. Public   Calculate interactions.
  !      INSNL1    Subr. Public   Initialization routine.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !     See subroutine documentation.
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !       !/S      Enable subroutine tracing.
  !       !/T(n)   Test output, see subroutines.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  !/
  !/
  PUBLIC
  !/
  !/ These are the arrays and variables used for GQM method
  !/
  INTEGER              :: nconf
  INTEGER, ALLOCATABLE :: k_if2 (:,:,:) , k_if3 (:,:,:) , k_1p2p(:,:,:) , &
       k_1p3m(:,:,:) , k_1p2m(:,:,:) , k_1p3p(:,:,:) , &
       k_1m2p(:,:,:) , k_1m3m(:,:,:) , k_1m2m(:,:,:) , &
       k_1m3p(:,:,:)
  INTEGER, ALLOCATABLE :: f_poin(:) , t_poin(:) , k_if1(:) , k_1p(:,:) ,  &
       k_1m(:,:) , idconf(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: f_coef(:) , f_proj(:) , tb_sca(:) , tb_v14(:)
  DOUBLE PRECISION, ALLOCATABLE :: tb_v24(:,:,:) , tb_v34(:,:,:) ,        &
       tb_tpm(:,:,:) , tb_tmp(:,:,:) , tb_fac(:,:,:)
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE w3snl1 (A, CG, KDMEAN, S, D)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         06-Jun-2018 |
    !/                  +-----------------------------------+
    !/
    !/    12-Jun-1996 : Final FORTRAN 77                    ( version 1.18 )
    !/    04-Feb-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    09-May-2002 : Switch clean up.                    ( version 2.21 )
    !/    24-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    03-Sep-2012 : Clean up of test output T0, T1      ( version 4.07 )
    !/    06-Jun-2018 : Add optional DEBUGSRC               ( version 6.04 )
    !/
    !  1. Purpose :
    !
    !     Calculate nonlinear interactions and the diagonal term of
    !     its derivative.
    !
    !  2. Method :
    !
    !     Discrete interaction approximation. (Hasselmann and Hasselmann
    !     1985; WAMDI group 1988)
    !
    !     The DIA is applied to the energy spectrum (instead of the action
    !     spectrum), for which is was originally developped. Because the
    !     frequency grid is invariant, the nonlinear interactions are
    !     calculated for the frequency spectrum, as in WAM. This requires
    !     only a single set of interpolation data which can be applied
    !     throughout the spatial domain. For deep water this is idenitical
    !     to a direct application to the wavenumber spectrum, for shallow
    !     water it is not. As the shallow water correction is nothing but
    !     a crude approximation, the choice between spectra is expected to
    !     be irrelevant.
    !
    !     The nonlinear interactions are calculated for two "mirror image"
    !     quadruplets as described in the manual. The central bin of these
    !     quadruples is placed on the discrete complonents of the spectrum,
    !     which requires interpolation to obtain other eneregy densities.
    !     The figure below defines the diferent basic counters and weights
    !     necessary for this interpolation.
    !
    !
    !               IFRM1  IFRM
    !                5        7    T |
    !          ITHM1  +------+     H +
    !                 |      |     E |      IFRP      IFRP1
    !                 |   \  |     T |       3           1
    !           ITHM  +------+     A +        +---------+  ITHP1
    !                6       \8      |        |         |
    !                                |        |  /      |
    !                           \    +        +---------+  ITHP
    !                                |      /4           2
    !                              \ |  /
    !          -+-----+------+-------#--------+---------+----------+
    !                              / |  \        FREQ.
    !                                |      \4           2
    !                           /    +        +---------+  ITHP
    !                                |        |  \      |
    !                6       /8      |        |         |
    !           ITHM  +------+       +        +---------+  ITHP1
    !                 |   \  |       |       3           1
    !                 |      |       |      IFRP      IFRP1
    !          ITHM1  +------+       +
    !                5        7      |
    !
    !     To create long vector loops and to efficiently deal with the
    !     closed nature of the directional space, the relative counters
    !     above are replaced by complete addresses stored in 32 arrays
    !     (see section 3 and INSNL1). The interaction are furthermore
    !     calucated for an extended spectrum, making it unnecessary to
    !     introduce extra weight factors for low and high frequencies.
    !     Therefore low and high frequencies are added to the local
    !     (auxiliary) spectrum as illustraed below.
    !
    !
    !              ^  +---+---------------------+---------+- NTH
    !              |  |   :                     :         |
    !                 |   :                     :         |
    !              d  | 2 :  original spectrum  :    1    |
    !              i  |   :                     :         |
    !              r  |   :                     :         |
    !                 +---+---------------------+---------+-  1
    !                            Frequencies -->     ^
    !         IFR =   0   1                    NFR   |  NFRHGH
    !                                                |
    !                                             NFRCHG
    !
    !     where : 1 : Extra tail added beyond NFR
    !             2 : Empty bins at low frequencies
    !
    !             NFRHGH = NFR + IFRP1 - IFRM1
    !             NFRCHG = NFR - IFRM1
    !
    !     All counters and arrays are set in INSNL1. See also section 3
    !     and section 8.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A       R.A.  I   Action spectrum A(ISP) as a function of
    !                         direction (rad)  and wavenumber.
    !       CG      R.A.  I   Group velocities (dimension NK).
    !       KDMEAN  Real  I   Mean relative depth.
    !       S       R.A.  O   Source term.                           *)
    !       D       R.A.  O   Diagonal term of derivative.           *)
    !     ----------------------------------------------------------------
    !                             *) 1-D array with dimension NTH*NK
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !      PRT2DS    Subr. W3ARRYMD Print plot of spectra.
    !      OUTMAT    Subr. W3WRRYMD Print out 2D matrix.
    !     ----------------------------------------------------------------
    !
    !  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 :
    !
    !       None.
    !
    !  8. Structure :
    !
    !     -------------------------------------------
    !      1.  Calculate proportionality constant.
    !      2.  Prepare auxiliary spectrum
    !      3.  Calculate (unfolded) interactions
    !        a Energy at interacting bins
    !        b Contribution to interactions
    !        c Fold interactions to side angles
    !      4.  Put source and diagonal term together
    !     -------------------------------------------
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable general test output.
    !     !/T0  2-D print plot of source term.
    !     !/T1  Print arrays.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE constants
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, fachfe,                &
         kdcon, kdmn, snlc1, snls1, snls2, snls3
    USE w3adatmd, ONLY: nfr, nfrhgh, nfrchg, nspecx, nspecy,        &
         ip11, ip12, ip13, ip14, im11, im12, im13, im14,   &
         ip21, ip22, ip23, ip24, im21, im22, im23, im24,   &
         ic11, ic12, ic21, ic22, ic31, ic32, ic41, ic42,   &
         ic51, ic52, ic61, ic62, ic71, ic72, ic81, ic82,   &
         dal1, dal2, dal3, af11,                           &
         awg1, awg2, awg3, awg4, awg5, awg6, awg7, awg8,   &
         swg1, swg2, swg3, swg4, swg5, swg6, swg7, swg8
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
#ifdef W3_T1
    USE w3odatmd, ONLY: ndst
#endif
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_T0
    USE w3arrymd, ONLY: prt2ds
#endif
#ifdef W3_T1
    USE w3arrymd, ONLY: outmat
#endif
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: A(NSPEC), CG(NK), KDMEAN
    REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IFR, ISP
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: X, X2, CONS, CONX, FACTOR,           &
         E00, EP1, EM1, EP2, EM2,             &
         SA1A, SA1B, SA2A, SA2B
#ifdef W3_T0
    REAL                    :: SOUT(NK,NFR), DOUT(NK,NFR)
#endif
    REAL               ::  UE  (1-NTH:NSPECY), SA1 (1-NTH:NSPECX),  &
         SA2 (1-NTH:NSPECX), DA1C(1-NTH:NSPECX),  &
         DA1P(1-NTH:NSPECX), DA1M(1-NTH:NSPECX),  &
         DA2C(1-NTH:NSPECX), DA2P(1-NTH:NSPECX),  &
         DA2M(1-NTH:NSPECX), CON (      NSPEC )
    !/
    !/ ------------------------------------------------------------------- /
    !/
    ! initialisations
    !
#ifdef W3_S
    CALL strace (ient, 'W3SNL1')
#endif
    !
    ! 1.  Calculate prop. constant --------------------------------------- *
    !
    x      = max( kdcon*kdmean , kdmn )
    x2     = max( -1.e15, snls3*x)
    cons   = snlc1 * ( 1. + snls1/x * (1.-snls2*x) * exp(x2) )
    !
#ifdef W3_T
    WRITE (ndst,9000) kdmean, cons
#endif
    !
    ! 2.  Prepare auxiliary spectrum and arrays -------------------------- *
    !
    DO ifr=1, nfr
      conx = tpiinv / sig(ifr) * cg(ifr)
      DO ith=1, nth
        isp       = ith + (ifr-1)*nth
        ue(isp) = a(isp) / conx
        con(isp) = conx
      END DO
    END DO
    !
    DO ifr=nfr+1, nfrhgh
      DO ith=1, nth
        isp      = ith + (ifr-1)*nth
        ue(isp) = ue(isp-nth) * fachfe
      END DO
    END DO
    !
    DO isp=1-nth, 0
      ue(isp) = 0.
      sa1(isp) = 0.
      sa2(isp) = 0.
      da1c(isp) = 0.
      da1p(isp) = 0.
      da1m(isp) = 0.
      da2c(isp) = 0.
      da2p(isp) = 0.
      da2m(isp) = 0.
    END DO
    !
    ! 3.  Calculate interactions for extended spectrum ------------------- *
    !
    DO isp=1, nspecx
      !
      ! 3.a Energy at interacting bins
      !
      e00    =        ue(isp)
      ep1    = awg1 * ue(ip11(isp)) + awg2 * ue(ip12(isp))        &
           + awg3 * ue(ip13(isp)) + awg4 * ue(ip14(isp))
      em1    = awg5 * ue(im11(isp)) + awg6 * ue(im12(isp))        &
           + awg7 * ue(im13(isp)) + awg8 * ue(im14(isp))
      ep2    = awg1 * ue(ip21(isp)) + awg2 * ue(ip22(isp))        &
           + awg3 * ue(ip23(isp)) + awg4 * ue(ip24(isp))
      em2    = awg5 * ue(im21(isp)) + awg6 * ue(im22(isp))        &
           + awg7 * ue(im23(isp)) + awg8 * ue(im24(isp))
      !
      ! 3.b Contribution to interactions
      !
      factor = cons * af11(isp) * e00
      !
      sa1a   = e00 * ( ep1*dal1 + em1*dal2 )
      sa1b   = sa1a - ep1*em1*dal3
      sa2a   = e00 * ( ep2*dal1 + em2*dal2 )
      sa2b   = sa2a - ep2*em2*dal3
      !
      sa1(isp) = factor * sa1b
      sa2(isp) = factor * sa2b
      !
      da1c(isp) = cons * af11(isp) * ( sa1a + sa1b )
      da1p(isp) = factor * ( dal1*e00 - dal3*em1 )
      da1m(isp) = factor * ( dal2*e00 - dal3*ep1 )
      !
      da2c(isp) = cons * af11(isp) * ( sa2a + sa2b )
      da2p(isp) = factor * ( dal1*e00 - dal3*em2 )
      da2m(isp) = factor * ( dal2*e00 - dal3*ep2 )
      !
    END DO
    !
    ! 4.  Put source and diagonal term together -------------------------- *
    !
    DO isp=1, nspec
      !
      s(isp) = con(isp) * ( - 2. * ( sa1(isp) + sa2(isp) )       &
           + awg1 * ( sa1(ic11(isp)) + sa2(ic12(isp)) )    &
           + awg2 * ( sa1(ic21(isp)) + sa2(ic22(isp)) )    &
           + awg3 * ( sa1(ic31(isp)) + sa2(ic32(isp)) )    &
           + awg4 * ( sa1(ic41(isp)) + sa2(ic42(isp)) )    &
           + awg5 * ( sa1(ic51(isp)) + sa2(ic52(isp)) )    &
           + awg6 * ( sa1(ic61(isp)) + sa2(ic62(isp)) )    &
           + awg7 * ( sa1(ic71(isp)) + sa2(ic72(isp)) )    &
           + awg8 * ( sa1(ic81(isp)) + sa2(ic82(isp)) ) )
      !
      d(isp) =  - 2. * ( da1c(isp) + da2c(isp) )                 &
           + swg1 * ( da1p(ic11(isp)) + da2p(ic12(isp)) )     &
           + swg2 * ( da1p(ic21(isp)) + da2p(ic22(isp)) )     &
           + swg3 * ( da1p(ic31(isp)) + da2p(ic32(isp)) )     &
           + swg4 * ( da1p(ic41(isp)) + da2p(ic42(isp)) )     &
           + swg5 * ( da1m(ic51(isp)) + da2m(ic52(isp)) )     &
           + swg6 * ( da1m(ic61(isp)) + da2m(ic62(isp)) )     &
           + swg7 * ( da1m(ic71(isp)) + da2m(ic72(isp)) )     &
           + swg8 * ( da1m(ic81(isp)) + da2m(ic82(isp)) )
      !
    END DO
    !
    ! ... Test output :
    !
#ifdef W3_T0
    DO ifr=1, nfr
      DO ith=1, nth
        isp          = ith + (ifr-1)*nth
        sout(ifr,ith) = s(isp) * tpi * sig(ifr) / cg(ifr)
        dout(ifr,ith) = d(isp)
      END DO
    END DO
    CALL prt2ds (ndst, nk, nk, nth, sout, sig(1:), '  ', 1.,  &
         0.0, 0.001, 'Snl(f,t)', ' ', 'NONAME')
    CALL prt2ds (ndst, nk, nk, nth, dout, sig(1:), '  ', 1.,  &
         0.0, 0.001, 'Diag Snl', ' ', 'NONAME')
#endif
    !
#ifdef W3_T1
    CALL outmat (ndst, s, nth, nth, nk, 'Snl')
    CALL outmat (ndst, d, nth, nth, nk, 'Diag Snl')
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3SNL1 : KDMEAN, CONS :',f8.2,f8.1)
#endif
    !/
    !/ End of W3SNL1 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3snl1
!/ ------------------------------------------------------------------- /
  SUBROUTINE insnl1 ( IMOD )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         24-Dec-2004 |
    !/                  +-----------------------------------+
    !/
    !/    19-Oct-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    04-Feb-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    09-May-2002 : Switch clean up.                    ( version 2.21 )
    !/    24-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/
    !  1. Purpose :
    !
    !     Preprocessing for nonlinear interactions (weights).
    !
    !  2. Method :
    !
    !     See W3SNL1.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       IMOD    Int.  I   Model number.
    !     ----------------------------------------------------------------
    !
    !     Local variables
    !     ----------------------------------------------------------------
    !       ITHxn   Real  Directional indices.                 (relative)
    !       IFRxn   Real  Frequency indices.                   (relative)
    !       IT1     R.A.  Directional indices.                      (1-D)
    !       IFn     R.A.  Frequency indices.                        (1-D)
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3IOGR    Subr. W3IOGRMD Model definition file processing.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     - Check on array dimensions for local arrays in W3SNL.
    !
    !  7. Remarks :
    !
    !     - Test output is generated through W3IOGR.
    !     - No testing of IMOD ir resetting of pointers.
    !
    !  8. Structure :
    !
    !     - See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nth, nspec, dth, xfr, sig, lam
    USE w3adatmd, ONLY: w3dmnl
    USE w3adatmd, ONLY: nfr, nfrhgh, nfrchg, nspecx, nspecy,        &
         ip11, ip12, ip13, ip14, im11, im12, im13, im14,   &
         ip21, ip22, ip23, ip24, im21, im22, im23, im24,   &
         ic11, ic12, ic21, ic22, ic31, ic32, ic41, ic42,   &
         ic51, ic52, ic61, ic62, ic71, ic72, ic81, ic82,   &
         dal1, dal2, dal3, af11,                           &
         awg1, awg2, awg3, awg4, awg5, awg6, awg7, awg8,   &
         swg1, swg2, swg3, swg4, swg5, swg6, swg7, swg8
    USE w3odatmd, ONLY: ndst, ndse
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: IMOD
    !/
    !/ Local parameters
    !/
    INTEGER                 :: IFR, ITH, ISP, ITHP, ITHP1, ITHM,    &
         ITHM1,IFRP, IFRP1, IFRM, IFRM1
    INTEGER, ALLOCATABLE    :: IF1(:), IF2(:), IF3(:), IF4(:),      &
         IF5(:), IF6(:), IF7(:), IF8(:),      &
         IT1(:), IT2(:), IT3(:), IT4(:),      &
         IT5(:), IT6(:), IT7(:), IT8(:)
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: DELTH3, DELTH4, LAMM2, LAMP2, CTHP,  &
         WTHP, WTHP1, CTHM, WTHM, WTHM1,      &
         XFRLN, WFRP, WFRP1, WFRM, WFRM1, FR, &
         AF11A
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'INSNL1')
#endif
#ifdef W3_T
    WRITE (ndst,9000) imod
#endif
    !
    nfr     = nk
    !
    ! 1.  Internal angles of quadruplet.
    !
    lamm2  = (1.-lam)**2
    lamp2  = (1.+lam)**2
    delth3 = acos( (lamm2**2+4.-lamp2**2) / (4.*lamm2) )
    delth4 = asin(-sin(delth3)*lamm2/lamp2)
    !
    ! 2.  Lambda dependend weight factors.
    !
    dal1   = 1. / (1.+lam)**4
    dal2   = 1. / (1.-lam)**4
    dal3   = 2. * dal1 * dal2
    !
    ! 3.  Directional indices.
    !
    cthp   = abs(delth4/dth)
    ithp   = int(cthp)
    ithp1  = ithp + 1
    wthp   = cthp - real(ithp)
    wthp1  = 1.- wthp
    !
    cthm   = abs(delth3/dth)
    ithm   = int(cthm)
    ithm1  = ithm + 1
    wthm   = cthm - real(ithm)
    wthm1  = 1.- wthm
    !
    ! 4.  Frequency indices.
    !
    xfrln  = log(xfr)
    !
    ifrp   = int( log(1.+lam) / xfrln )
    ifrp1  = ifrp + 1
    wfrp   = (1.+lam - xfr**ifrp) / (xfr**ifrp1 - xfr**ifrp)
    wfrp1  = 1. - wfrp
    !
    ifrm   = int( log(1.-lam) / xfrln )
    ifrm1  = ifrm - 1
    wfrm   = (xfr**ifrm -(1.-lam)) / (xfr**ifrm - xfr**ifrm1)
    wfrm1  = 1. - wfrm
    !
    ! 5.  Range of calculations
    !
    nfrhgh = nfr + ifrp1 - ifrm1
    nfrchg = nfr - ifrm1
    nspecy = nfrhgh * nth
    nspecx = nfrchg * nth
    !
    ! 6.  Allocate arrays or check array sizes
    !
    CALL w3dmnl ( imod, ndse, ndst, nspec, nspecx )
    !
    ALLOCATE ( if1(nfrchg), if2(nfrchg), if3(nfrchg), if4(nfrchg),  &
         if5(nfrchg), if6(nfrchg), if7(nfrchg), if8(nfrchg),  &
         it1(nth), it2(nth), it3(nth), it4(nth),              &
         it5(nth), it6(nth), it7(nth), it8(nth) )
    !
    ! 7.  Spectral addresses
    !
    DO ifr=1, nfrchg
      if1(ifr) =           ifr+ifrp
      if2(ifr) =           ifr+ifrp1
      if3(ifr) = max( 0 , ifr+ifrm  )
      if4(ifr) = max( 0 , ifr+ifrm1 )
      if5(ifr) = max( 0 , ifr-ifrp  )
      if6(ifr) = max( 0 , ifr-ifrp1 )
      if7(ifr) =           ifr-ifrm
      if8(ifr) =           ifr-ifrm1
    END DO
    !
    DO ith=1, nth
      it1(ith) = ith + ithp
      it2(ith) = ith + ithp1
      it3(ith) = ith + ithm
      it4(ith) = ith + ithm1
      it5(ith) = ith - ithp
      it6(ith) = ith - ithp1
      it7(ith) = ith - ithm
      it8(ith) = ith - ithm1
      IF ( it1(ith).GT.nth) it1(ith) = it1(ith) - nth
      IF ( it2(ith).GT.nth) it2(ith) = it2(ith) - nth
      IF ( it3(ith).GT.nth) it3(ith) = it3(ith) - nth
      IF ( it4(ith).GT.nth) it4(ith) = it4(ith) - nth
      IF ( it5(ith).LT. 1 ) it5(ith) = it5(ith) + nth
      IF ( it6(ith).LT. 1 ) it6(ith) = it6(ith) + nth
      IF ( it7(ith).LT. 1 ) it7(ith) = it7(ith) + nth
      IF ( it8(ith).LT. 1 ) it8(ith) = it8(ith) + nth
    END DO
    !
    DO isp=1, nspecx
      ifr       = 1 + (isp-1)/nth
      ith       = 1 + mod(isp-1,nth)
      ip11(isp) = it2(ith) + (if2(ifr)-1)*nth
      ip12(isp) = it1(ith) + (if2(ifr)-1)*nth
      ip13(isp) = it2(ith) + (if1(ifr)-1)*nth
      ip14(isp) = it1(ith) + (if1(ifr)-1)*nth
      im11(isp) = it8(ith) + (if4(ifr)-1)*nth
      im12(isp) = it7(ith) + (if4(ifr)-1)*nth
      im13(isp) = it8(ith) + (if3(ifr)-1)*nth
      im14(isp) = it7(ith) + (if3(ifr)-1)*nth
      ip21(isp) = it6(ith) + (if2(ifr)-1)*nth
      ip22(isp) = it5(ith) + (if2(ifr)-1)*nth
      ip23(isp) = it6(ith) + (if1(ifr)-1)*nth
      ip24(isp) = it5(ith) + (if1(ifr)-1)*nth
      im21(isp) = it4(ith) + (if4(ifr)-1)*nth
      im22(isp) = it3(ith) + (if4(ifr)-1)*nth
      im23(isp) = it4(ith) + (if3(ifr)-1)*nth
      im24(isp) = it3(ith) + (if3(ifr)-1)*nth
    END DO
    !
    DO isp=1, nspec
      ifr       = 1 + (isp-1)/nth
      ith       = 1 + mod(isp-1,nth)
      ic11(isp) = it6(ith) + (if6(ifr)-1)*nth
      ic21(isp) = it5(ith) + (if6(ifr)-1)*nth
      ic31(isp) = it6(ith) + (if5(ifr)-1)*nth
      ic41(isp) = it5(ith) + (if5(ifr)-1)*nth
      ic51(isp) = it4(ith) + (if8(ifr)-1)*nth
      ic61(isp) = it3(ith) + (if8(ifr)-1)*nth
      ic71(isp) = it4(ith) + (if7(ifr)-1)*nth
      ic81(isp) = it3(ith) + (if7(ifr)-1)*nth
      ic12(isp) = it2(ith) + (if6(ifr)-1)*nth
      ic22(isp) = it1(ith) + (if6(ifr)-1)*nth
      ic32(isp) = it2(ith) + (if5(ifr)-1)*nth
      ic42(isp) = it1(ith) + (if5(ifr)-1)*nth
      ic52(isp) = it8(ith) + (if8(ifr)-1)*nth
      ic62(isp) = it7(ith) + (if8(ifr)-1)*nth
      ic72(isp) = it8(ith) + (if7(ifr)-1)*nth
      ic82(isp) = it7(ith) + (if7(ifr)-1)*nth
    END DO
    !
    DEALLOCATE ( if1, if2, if3, if4, if5, if6, if7, if8,  &
         it1, it2, it3, it4, it5, it6, it7, it8 )
    !
    ! 8.  Fill scaling array (f**11)
    !
    DO ifr=1, nfr
      af11a  = (sig(ifr)*tpiinv)**11
      DO ith=1, nth
        af11(ith+(ifr-1)*nth) = af11a
      END DO
    END DO
    !
    fr     = sig(nfr)*tpiinv
    DO ifr=nfr+1, nfrchg
      fr     = fr * xfr
      af11a  = fr**11
      DO ith=1, nth
        af11(ith+(ifr-1)*nth) = af11a
      END DO
    END DO
    !
    ! 9.  Interpolation weights
    !
    awg1   = wthp  * wfrp
    awg2   = wthp1 * wfrp
    awg3   = wthp  * wfrp1
    awg4   = wthp1 * wfrp1
    awg5   = wthm  * wfrm
    awg6   = wthm1 * wfrm
    awg7   = wthm  * wfrm1
    awg8   = wthm1 * wfrm1
    !
    swg1   = awg1**2
    swg2   = awg2**2
    swg3   = awg3**2
    swg4   = awg4**2
    swg5   = awg5**2
    swg6   = awg6**2
    swg7   = awg7**2
    swg8   = awg8**2
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST INSNL1 : IMOD :',i4)
#endif
    !/
    !/ End of INSNL1 ----------------------------------------------------- /
    !/
  END SUBROUTINE insnl1
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3snlgqm(A,CG,WN,DEPTH,TSTOTn,TSDERn)
    ! This and the following routines are adapted to WW3 from TOMAWAC qnlin3.f
    !***********************************************************************
    ! TOMAWAC   V6P1                                   24/06/2011
    !***********************************************************************
    !
    !brief    COMPUTES THE CONTRIBUTION OF THE NON-LINEAR INTERACTIONS
    !+                SOURCE TERM BETWEEN QUADRUPLETS USING THE GQM METHOD
    !+                ("GAUSSIAN QUADRATURE METHOD") PROPOSED BY LAVRENOV
    !+                (2001)
    !+
    !+            PROCEDURE SPECIFIC TO THE CASE WHERE THE FREQUENCIES
    !+                FOLLOW A GEOMETRICAL PROGRESSION AND THE DIRECTIONS
    !+                ARE EVENLY DISTRIBUTED OVER [0;2.PI].
    !
    !note     THIS SUBROUTINE USES THE OUTPUT FROM 'PRENL3' TO OPTIMISE
    !+          THE COMPUTATIONS FOR DIA.
    !
    !reference  LAVRENOV, I.V. (2001):
    !+           "EFFECT OF WIND WAVE PARAMETER FLUCTUATION ON THE NONLINEAR
    !+           SPECTRUM EVOLUTION". J. PHYS. OCEANOGR. 31, 861-873.
    !
    !history  E. GAGNAIRE-RENOU
    !+        04/2011
    !+        V6P1
    !+   CREATED
    !
    !history  G.MATTAROLO (EDF - LNHE)
    !+        24/06/2011
    !+        V6P1
    !+   Translation of French names of the variables in argument
 
    !
    !/ Warning, contrary to the DIA routine, there is no extension to frequencies below IK=1
    !/ as a result the first two frequencies are not fully treated.
    !==================================================================================
    !     This subroutine is same as qnlin3 in TOMWAC
    USE constants, ONLY: tpi
    USE w3gdatmd,  ONLY: sig, nk ,  nth , dth, xfr, fr1, gqthrsat, gqamp
 
    IMPLICIT NONE
 
    REAL, intent(in) :: A(NTH,NK), CG(NK), WN(NK)
    REAL, intent(in) :: DEPTH
    REAL, intent(out) :: TSTOTn(NTH,NK), TSDERn(NTH,NK)
 
    INTEGER          :: ITH,IK,NT,NF
    REAL             :: q_dfac, SATVAL(NK), SUME, ACCVAL, ACCMAX, AMPFAC
    DOUBLE PRECISION :: RAISF, FREQ(NK)
    DOUBLE PRECISION :: TSTOT(NTH,NK) , TSDER(NTH,NK), F(NTH,NK)
    DOUBLE PRECISION :: TEMP
 
    !.....LOCAL VARIABLES
    INTEGER             JF    , JT    , JF1   , JT1  , IQ_OM2 &
         , JFM0  , JFM1  , JFM2  , JFM3  , IXF1 , IXF2   &
         , IXF3  , JFMIN , JFMAX , ICONF , LBUF
    INTEGER            KT1P  , KT1M  , JT1P  , JT1M  , KT1P2P, KT1P2M &
         , KT1P3P, KT1P3M, KT1M2P, KT1M2M, KT1M3P, KT1M3M &
         , JT1P2P, JT1P2M, JT1P3P, JT1P3M, JT1M2P, JT1M2M &
         , JT1M3P, JT1M3M
    DOUBLE PRECISION  V1_4  , V2_4  , V3_4  , Q_2P3M, Q_2M3P, FACTOR &
         , T_2P3M, T_2M3P, S_2P3M, S_2M3P, SCAL_T, T2P3M &
         , T2M3P , SP0   , SP1P  , SP1M  , SP1P2P, SP1P2M &
         , SP1P3P, SP1P3M, SP1M2P, SP1M2M, SP1M3P, SP1M3M &
         , CF0   , CP0   , CF1   , CP1   , CF2   , CP2   &
         , CF3   , CP3   , Q2PD0 , Q2PD1 , Q2PD2P, Q2PD3M &
         , Q2MD0 , Q2MD1 , Q2MD2M, Q2MD3P ,AUX00 , AUX01  &
         , AUX02 , AUX03 , AUX04 , AUX05 , SEUIL  &
         , AUX06 , AUX07 , AUX08 , AUX09 , AUX10 , FSEUIL
 
    nt = nth
    nf = nk
    lbuf = 500
    seuil = 0.
    raisf = xfr
 
    DO ik = 1,nk
      freq(ik) = fr1*raisf**(ik-1)
    ENDDO
 
    DO ith = 1,nth
      DO ik = 1,nk
        ! F is the E(f,theta) spectrum ...
        f(ith,ik) = dble(a(ith,ik)*sig(ik))*dble(tpi)/dble(cg(ik))
      ENDDO
    ENDDO
    !   CALL INSNLGQM
    ! it returns: F_POIN , T_POIN , F_COEF , F_PROJ, TB_SCA , K_IF1, K_1P, k_1M , K_IF2
    !             K_IF3, K_1P2P , K_1P3M , K_1P2M , K_1P3P , K_1M2P , K_1M3M ,  K_1M2M
    !             K_1M3P , TB_V14 , TB_FAC , TB_V24 , TB_V34 , TB_TMP , TB_TPM , IDCONF, NCONF
    !=======================================================================
    !     COMPUTES THE GENERALIZED MIN AND MAX FREQUENCIES : INSTEAD OF GOING
    !     FROM 1 TO NF IN FREQ(JF) FOR THE MAIN FREQUENCY, IT GOES FROM JFMIN
    !     TO JFMAX
    !     JFMIN IS GIVEN BY Fmin=FREQ(1) /Gamma_min
    !     JFMAX IS GIVEN BY Fmax=FREQ(NF)*Gamma_max
    !     TESTS HAVE SHOWN THAT IT CAN BE ASSUMED Gamma_min=1. (JFMIN=1) AND
    !     Gamma_max=1.3 (JFMAX>NF) TO OBTAIN IMPROVED RESULTS
    !     Note by Fabrice Ardhuin: this appears to give the difference in tail benaviour with Gerbrant's WRT
    !=======================================================================
    jfmin=max(1-int(log(1.0d0)/log(raisf)),1)
    jfmax=min(nf+int(log(1.3d0)/log(raisf)),nk)
    !
    !=======================================================================
    !     COMPUTES THE SPECTRUM THRESHOLD VALUES (BELOW WHICH QNL4 IS NOT
    !     CALCULATED). THE THRESHOLD IS SET WITHIN 0 AND 1.
    ! This was commented by FA
    !=======================================================================
    !        AUX00=0.0D0
    !        DO JF=1,NF
    !          DO JT=1,NT
    !            IF (F(JT,JF).GT.AUX00) AUX00=F(JT,JF)
    !          ENDDO
    !        ENDDO
    !        FSEUIL=AUX00*SEUIL
 
    tstot = 0.
    tsder = 0.
    !=======================================================================
    accmax=0.
    DO jf=jfmin,jfmax
      sume=sum(f(:,jf))*dth
      satval(jf) = sume*freq(jf)**5
      accval = sume*freq(jf)**4
      IF (accval.GT.accmax) accmax=accval
    END DO
 
 
    !     ==================================================
    !     STARTS LOOP 1 OVER THE SELECTED CONFIGURATIONS
    !     ==================================================
    DO iconf=1,nconf
      !       ---------selected configuration characteristics
      jf1   =idconf(iconf,1)
      jt1   =idconf(iconf,2)
      iq_om2=idconf(iconf,3)
      !
      !       ---------Recovers V1**4=(f1/f0)**4
      v1_4  =tb_v14(jf1)
      !       ---------Recovers the shift of the frequency index on f1
      ixf1  =k_if1(jf1)
      !       ---------Recovers the direction indexes for Delat1
      kt1p  =k_1p(jt1,jf1)
      kt1m  =k_1m(jt1,jf1)
      !       ---------Recovers V2**4=(f2/f0)**4 and V3**4=(f3/f0)**4
      v2_4  =tb_v24(iq_om2,jt1,jf1)
      v3_4  =tb_v34(iq_om2,jt1,jf1)
      !       ---------Recovers the frequency indexes shift on f2 and f3
      ixf2  =k_if2(iq_om2,jt1,jf1)
      ixf3  =k_if3(iq_om2,jt1,jf1)
      !       ---------Recovers the direction indexes shift
      kt1p2p=k_1p2p(iq_om2,jt1,jf1)
      kt1p2m=k_1p2m(iq_om2,jt1,jf1)
      kt1p3p=k_1p3p(iq_om2,jt1,jf1)
      kt1p3m=k_1p3m(iq_om2,jt1,jf1)
      kt1m2p=k_1m2p(iq_om2,jt1,jf1)
      kt1m2m=k_1m2m(iq_om2,jt1,jf1)
      kt1m3p=k_1m3p(iq_om2,jt1,jf1)
      kt1m3m=k_1m3m(iq_om2,jt1,jf1)
      !       ---------Recovers the coupling coefficients
      t2p3m =tb_tpm(iq_om2,jt1,jf1)
      t2m3p =tb_tmp(iq_om2,jt1,jf1)
      !       ---------Recovers the multiplicative factor of QNL4
      factor=tb_fac(iq_om2,jt1,jf1)
 
      !       = = = = = = = = = = = = = = = = = = = = = = = = =
      !       STARTS LOOP 2 OVER THE SPECTRUM FREQUENCIES
      !       = = = = = = = = = = = = = = = = = = = = = = = = =
      DO jf=jfmin,jfmax
        IF (satval(jf).GT.gqthrsat) THEN
          !
          !.........Recovers the coefficient for the coupling factor
          !.........Computes the coupling coefficients for the case +Delta1 (SIG=1)
          scal_t=tb_sca(lbuf+jf)*factor
          t_2p3m=t2p3m*scal_t
          t_2m3p=t2m3p*scal_t
          !
          !.........Frequency indexes and coefficients
          jfm0=f_poin(jf+lbuf)
          cf0 =f_coef(jf+lbuf)
          cp0 =f_proj(jf+lbuf)
          jfm1=f_poin(jf+ixf1)
          cf1 =f_coef(jf+ixf1)
          cp1 =f_proj(jf+ixf1)
          jfm2=f_poin(jf+ixf2)
          cf2 =f_coef(jf+ixf2)
          cp2 =f_proj(jf+ixf2)
          jfm3=f_poin(jf+ixf3)
          cf3 =f_coef(jf+ixf3)
          cp3 =f_proj(jf+ixf3)
          !
          !         -------------------------------------------------
          !         STARTS LOOP 3 OVER THE SPECTRUM DIRECTIONS
          !         -------------------------------------------------
          DO jt=1,nt
            !
            !...........Direction indexes
            !           direct config (+delta1) (sig =1)
            jt1p  =t_poin(jt+kt1p)
            jt1p2p=t_poin(jt+kt1p2p)
            jt1p2m=t_poin(jt+kt1p2m)
            jt1p3p=t_poin(jt+kt1p3p)
            jt1p3m=t_poin(jt+kt1p3m)
            !           image config (-delta1)
            jt1m  =t_poin(jt+kt1m)
            jt1m2p=t_poin(jt+kt1m2p)
            jt1m2m=t_poin(jt+kt1m2m)
            jt1m3p=t_poin(jt+kt1m3p)
            jt1m3m=t_poin(jt+kt1m3m)
            !
            !           - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            !           STARTS LOOP 4 OVER THE MESH NODES
            !           - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
            !
            sp0=f(jt,jfm0)*cf0
            !
            !              IF (SP0.GT.FSEUIL) THEN
            !
            !               Config. +Delta1 (SIG=1)
            !               =======================
            !...............Computes the spectrum values in 1, 2, 3
            sp1p  =f(jt1p  ,jfm1)*cf1
            sp1p2p=f(jt1p2p,jfm2)*cf2
            sp1p3m=f(jt1p3m,jfm3)*cf3
            sp1p2m=f(jt1p2m,jfm2)*cf2
            sp1p3p=f(jt1p3p,jfm3)*cf3
            !
            !...............Computes auxiliary products and variables
            aux01=sp0*v1_4+sp1p
            aux02=sp0*sp1p
            aux03=sp1p2p*sp1p3m
            aux04=sp1p2p*v3_4+sp1p3m*v2_4
            aux05=sp1p2m*sp1p3p
            aux06=sp1p2m*v3_4+sp1p3p*v2_4
            aux07=aux02*v3_4
            aux08=aux02*v2_4
            !
            !...............Computes the components of the transfer term
            s_2p3m=aux03*aux01-aux02*aux04
            s_2m3p=aux05*aux01-aux02*aux06
            q_2p3m=t_2p3m*s_2p3m
            q_2m3p=t_2m3p*s_2m3p
            aux00 =q_2p3m+q_2m3p
            !
            !...............Computes the components of the derived terms (dQ/dF)
            q2pd0 =t_2p3m*(aux03*v1_4   - sp1p*aux04)*cf0
            q2pd1 =t_2p3m*(aux03        - sp0 *aux04)*cf1
            q2pd2p=t_2p3m*(aux01*sp1p3m - aux07     )*cf2
            q2pd3m=t_2p3m*(aux01*sp1p2p - aux08     )*cf3
            q2md0 =t_2m3p*(aux05*v1_4   - sp1p*aux06)*cf0
            q2md1 =t_2m3p*(aux03        - sp0 *aux06)*cf1
            q2md2m=t_2m3p*(aux01*sp1p3p - aux07     )*cf2
            q2md3p=t_2m3p*(aux01*sp1p2m - aux08     )*cf3
            aux09=q2pd0+q2md0
            aux10=q2pd1+q2md1
            !
            !...............Sum of Qnl4 term in the table TSTOT
            tstot(jt,jfm0    )=tstot(jt,jfm0    )+aux00 *cp0
            tstot(jt1p,jfm1  )=tstot(jt1p,jfm1  )+aux00 *cp1
            tstot(jt1p2p,jfm2)=tstot(jt1p2p,jfm2)-q_2p3m*cp2
            tstot(jt1p2m,jfm2)=tstot(jt1p2m,jfm2)-q_2m3p*cp2
            tstot(jt1p3m,jfm3)=tstot(jt1p3m,jfm3)-q_2p3m*cp3
            tstot(jt1p3p,jfm3)=tstot(jt1p3p,jfm3)-q_2m3p*cp3
            !
            !...............Sum of the term dQnl4/dF in the table TSDER
            tsder(jt,jfm0)=tsder(jt,jfm0)+aux09 *cp0
            tsder(jt1p,jfm1)=tsder(jt1p,jfm1)+aux10 *cp1
            tsder(jt1p2p,jfm2)=tsder(jt1p2p,jfm2)-q2pd2p*cp2
            tsder(jt1p2m,jfm2)=tsder(jt1p2m,jfm2)-q2md2m*cp2
            tsder(jt1p3m,jfm3)=tsder(jt1p3m,jfm3)-q2pd3m*cp3
            tsder(jt1p3p,jfm3)=tsder(jt1p3p,jfm3)-q2md3p*cp3
#ifdef W3_TGQM
            ! Test output to set up triplet method ...
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt,    jfm0,aux00 *cp0, f(jt,jfm0),tstot(jt    ,jfm0)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1p,  jfm1,aux00 *cp1, f(jt1p,jfm1),tstot(jt1p,jfm1)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1p2p,jfm2,-q_2p3m*cp2,f(jt1p2p,jfm2),tstot(jt1p2p,jfm2)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1p2m,jfm2,-q_2m3p*cp2,f(jt1p2m,jfm2),tstot(jt1p2m,jfm2)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1p3m,jfm2,-q_2p3m*cp3,f(jt1p3m,jfm3),tstot(jt1p3m,jfm3)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1p3p,jfm2,-q_2m3p*cp3,f(jt1p3p,jfm3),tstot(jt1p3p,jfm3)
            temp=(tb_tpm(iq_om2,jt1,jf1)*(( f(jt1p2p,jfm2)*cf2 *f(jt1p3m,jfm3)*cf3)* &
                 (f(jt,jfm0    )*cf0*tb_v14(jf1)+f(jt1p  ,jfm1)*cf1) &
                 -sp0*sp1p*(sp1p2p*v3_4+sp1p3m*v2_4))+t_2m3p*(aux05*aux01-aux02*aux06)) *cp0
            WRITE(995,'(3I3,3E12.3)') iconf,jf,jt, f(jt,jfm0)
            temp=(q_2p3m+q_2m3p) *cp1
            WRITE(995,'(5I3,3E12.3)') iconf,jf,jt,jt1p,  jfm1,aux00 *cp1, f(jt1p,jfm1),tstot(jt1p,jfm1)
            WRITE(995,'(5I3,3E12.3)') iconf,jf,jt,jt1p2p,jfm2,-q_2p3m*cp2,f(jt1p2p,jfm2),tstot(jt1p2p,jfm2)
            WRITE(995,'(5I3,3E12.3)') iconf,jf,jt,jt1p2m,jfm2,-q_2m3p*cp2,f(jt1p2m,jfm2),tstot(jt1p2m,jfm2)
            WRITE(995,'(5I3,3E12.3)') iconf,jf,jt,jt1p3m,jfm2,-q_2p3m*cp3,f(jt1p3m,jfm3),tstot(jt1p3m,jfm3)
            WRITE(995,'(5I3,3E12.3)') iconf,jf,jt,jt1p3p,jfm2,-q_2m3p*cp3,f(jt1p3p,jfm3),tstot(jt1p3p,jfm3)
#endif
            !
            !               Config. -Delta1 (SIG=-1)
            !               ========================
            !...............Computes the spectrum values in 1, 2, 3
            sp1m  =f(jt1m  ,jfm1)*cf1
            sp1m2p=f(jt1m2p,jfm2)*cf2
            sp1m3m=f(jt1m3m,jfm3)*cf3
            sp1m2m=f(jt1m2m,jfm2)*cf2
            sp1m3p=f(jt1m3p,jfm3)*cf3
            !
            !...............Computes auxiliary products and variables
            aux01=sp0*v1_4+sp1m
            aux02=sp0*sp1m
            aux03=sp1m2p*sp1m3m
            aux04=sp1m2p*v3_4+sp1m3m*v2_4
            aux05=sp1m2m*sp1m3p
            aux06=sp1m2m*v3_4+sp1m3p*v2_4
            aux07=aux02*v3_4
            aux08=aux02*v2_4
            !
            !...............Computes the transfer term components
            s_2p3m=aux03*aux01-aux02*aux04
            s_2m3p=aux05*aux01-aux02*aux06
            q_2p3m=t_2m3p*s_2p3m
            q_2m3p=t_2p3m*s_2m3p
            aux00 =q_2p3m+q_2m3p   ! Same as in +Delta1, can be commented out
            !
            !...............Computes the derived terms components (dQ/dF)
            q2pd0 =t_2p3m*(aux03*v1_4   - sp1m*aux04)*cf0
            q2pd1 =t_2p3m*(aux03        - sp0 *aux04)*cf1
            q2pd2p=t_2p3m*(aux01*sp1m3m - aux07     )*cf2
            q2pd3m=t_2p3m*(aux01*sp1m2p - aux08     )*cf3
            q2md0 =t_2m3p*(aux05*v1_4   - sp1m*aux06)*cf0
            q2md1 =t_2m3p*(aux03        - sp0 *aux06)*cf1
            q2md2m=t_2m3p*(aux01*sp1m3p - aux07     )*cf2
            q2md3p=t_2m3p*(aux01*sp1m2m - aux08     )*cf3
            aux09=q2pd0+q2md0
            aux10=q2pd1+q2md1
            !
            !...............Sum of Qnl4 term in the table TSTOT
            tstot(jt    ,jfm0)=tstot(jt    ,jfm0)+aux00 *cp0
            tstot(jt1m  ,jfm1)=tstot(jt1m  ,jfm1)+aux00 *cp1
            tstot(jt1m2p,jfm2)=tstot(jt1m2p,jfm2)-q_2p3m*cp2
            tstot(jt1m2m,jfm2)=tstot(jt1m2m,jfm2)-q_2m3p*cp2
            tstot(jt1m3m,jfm3)=tstot(jt1m3m,jfm3)-q_2p3m*cp3
            tstot(jt1m3p,jfm3)=tstot(jt1m3p,jfm3)-q_2m3p*cp3
            !
            !...............Sum of the term dQnl4/dF in the table TSDER
            tsder(jt    ,jfm0)=tsder(jt    ,jfm0)+aux09 *cp0
            tsder(jt1m  ,jfm1)=tsder(jt1m  ,jfm1)+aux10 *cp1
            tsder(jt1m2p,jfm2)=tsder(jt1m2p,jfm2)-q2pd2p*cp2
            tsder(jt1m2m,jfm2)=tsder(jt1m2m,jfm2)-q2md2m*cp2
            tsder(jt1m3m,jfm3)=tsder(jt1m3m,jfm3)-q2pd3m*cp3
            tsder(jt1m3p,jfm3)=tsder(jt1m3p,jfm3)-q2md3p*cp3
            !
#ifdef W3_TGQM
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt,    jfm0,aux00 *cp0, f(jt,jfm0),tstot(jt    ,jfm0)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1m,  jfm1,aux00 *cp1, f(jt1m,jfm1),tstot(jt1m,jfm1)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1m2p,jfm2,-q_2p3m*cp2,f(jt1m2p,jfm2),tstot(jt1m2p,jfm2)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1m2m,jfm2,-q_2m3p*cp2,f(jt1m2m,jfm2),tstot(jt1m2m,jfm2)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1m3m,jfm2,-q_2p3m*cp3,f(jt1m3m,jfm3),tstot(jt1m3m,jfm3)
            WRITE(994,'(5I3,3E12.3)') iconf,jf,jt,jt1m3p,jfm2,-q_2m3p*cp3,f(jt1m3p,jfm3),tstot(jt1m3p,jfm3)
#endif
            !
            !              ENDIF ! this was the test on SEUIL
            !
          ENDDO
          !         -------------------------------------------------
          !         END OF LOOP 3 OVER THE SPECTRUM DIRECTIONS
          !         -------------------------------------------------
          !
        ENDIF ! End of test on saturation level
      ENDDO
      !       = = = = = = = = = = = = = = = = = = = = = = = = =
      !       END OF LOOP 2 OVER THE SPECTRUM FREQUENCIES
      !       = = = = = = = = = = = = = = = = = = = = = = = = =
      !
    ENDDO
    !     ==================================================
    !     END OF LOOP 1 OVER THE SELECTED CONFIGURATIONS
    !     ==================================================
    ! Applying WAM DEPTH SCALING ! to be added later ...
    !      CALL q_dscale(F,WN,SIG,DTH,NK,NTH,DEPTH,q_dfac)
    q_dfac=1
 
    ! Amplification inspired by Lavrenov 2001, eq 10.
    ampfac=gqamp(4)*min(max(accmax/gqamp(2),1.)**gqamp(1),gqamp(3))
    !WRITE(991,*) ACCMAX,q_dfac,AMPFAC,GQAMP(1:3),SATVAL(10),SATVAL(30)
 
    ! Replacing Double Precision with Simple Real and scaling
    tstotn = tstot*q_dfac*ampfac
    tsdern = tsder*q_dfac*ampfac
 
 
    ! Converting Snl(theta,f) to Snl(theta,k)/sigma
    DO ith = 1,nt
      DO ik = 1,nf
        tstotn(ith,ik) = tstotn(ith,ik)*cg(ik)/(tpi*sig(ik))
      ENDDO
    ENDDO
    !CLOSE(994)
    !STOP
  END SUBROUTINE w3snlgqm
 
  !/ ------------------------------------------------------------------- /
  FUNCTION couple(XK1 ,YK1 ,XK2 ,YK2 ,XK3 ,YK3 ,XK4 ,YK4)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  | M. Benoit & E. Gagnaire-Renou     |
    !/                  | Last update :         20-Nov-2022 |
    !/                  +-----------------------------------+
    !/
    !/    19-Nov-2022 : Transfer from TOMAWAC code          ( version 7.xx )
    !/
    !  1. Purpose :
    !
    !     Computes the 4-wave coupling coefficient used in Snl4
    !
    !  2. Method :
    !
    !     Uses theoretical expression by Webb (1978)
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       XK1     Real  I   x component of k1 wavenumber ...
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      INNSLGQM  Subr. W3SNL2   Prepares source term integration.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: grav
    !
    IMPLICIT NONE
    
    DOUBLE PRECISION, INTENT(IN)    :: xk1   , yk1   , xk2   , yk2
    DOUBLE PRECISION, INTENT(IN)    :: xk3   , yk3
    DOUBLE PRECISION, INTENT(IN)    :: xk4   , yk4
    DOUBLE PRECISION couple
    !
    !.....LOCAL VARIABLES
    !     """"""""""""""""""
    DOUBLE PRECISION rk1   , rk2   , rk3   , rk4   , wk1   , wk2
    DOUBLE PRECISION wk3   , wk4   , s12   , s13   , s14   , s23
    DOUBLE PRECISION s24   , s34   , w1p2  , q12   , w1m3  , q13
    DOUBLE PRECISION w1m4  , q14   , ddd   , coef  , deno13, nume13
    DOUBLE PRECISION deno14, nume14, zero, pi
 
    !
    pi = acos(-1.)
    coef=pi*grav*grav/4.d0
    zero=1.d-10
    !
    rk1=sqrt(xk1*xk1+yk1*yk1)
    rk2=sqrt(xk2*xk2+yk2*yk2)
    rk3=sqrt(xk3*xk3+yk3*yk3)
    rk4=sqrt(xk4*xk4+yk4*yk4)
    !
    wk1=sqrt(rk1)
    wk2=sqrt(rk2)
    wk3=sqrt(rk3)
    wk4=sqrt(rk4)
    !
    s12=xk1*xk2+yk1*yk2
    s13=xk1*xk3+yk1*yk3
    s14=xk1*xk4+yk1*yk4
    s23=xk2*xk3+yk2*yk3
    s24=xk2*xk4+yk2*yk4
    s34=xk3*xk4+yk3*yk4
    !
    w1p2=sqrt((xk1+xk2)*(xk1+xk2)+(yk1+yk2)*(yk1+yk2))
    w1m3=sqrt((xk1-xk3)*(xk1-xk3)+(yk1-yk3)*(yk1-yk3))
    w1m4=sqrt((xk1-xk4)*(xk1-xk4)+(yk1-yk4)*(yk1-yk4))
    q12=(wk1+wk2)*(wk1+wk2)
    q13=(wk1-wk3)*(wk1-wk3)
    q14=(wk1-wk4)*(wk1-wk4)
    !
    !.....COMPUTES THE D COEFFICIENT OF WEBB (1978)
    !     """"""""""""""""""""""""""""""""""""""
    ddd=2.00d0*q12*(rk1*rk2-s12)*(rk3*rk4-s34)/(w1p2-q12) &
         +0.50d0*(s12*s34+s13*s24+s14*s23) &
         +0.25d0*(s13+s24)*q13*q13 &
         -0.25d0*(s12+s34)*q12*q12 &
         +0.25d0*(s14+s23)*q14*q14 &
         +2.50d0*rk1*rk2*rk3*rk4 &
         +q12*q13*q14*(rk1+rk2+rk3+rk4)
 
    deno13=w1m3-q13
    nume13=2.00d0*q13*(rk1*rk3+s13)*(rk2*rk4+s24)
    IF (abs(deno13).LT.zero) THEN
      IF (abs(nume13).LT.zero) THEN
        WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-3)  0/0 !'
      ELSE
        WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-3) inifinte value'
      ENDIF
      WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-3) term not used'
    ELSE
      ddd=ddd+nume13/deno13
    ENDIF
    deno14=w1m4-q14
    nume14=2.00d0*q14*(rk1*rk4+s14)*(rk2*rk3+s23)
    IF (abs(deno14).LT.zero) THEN
      IF (abs(nume14).LT.zero) THEN
        WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-4)  0/0 !'
      ELSE
        WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-4) inifinte value'
      ENDIF
      WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-4) term not used'
    ELSE
      ddd=ddd+nume14/deno14
    ENDIF
 
    couple=coef*ddd*ddd/(wk1*wk2*wk3*wk4)
    !      RETURN
  END FUNCTION couple
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE gauleg (W_LEG ,X_LEG ,NPOIN)
    !/ ------------------------------------------------------------------- /
    !.....VARIABLES IN ARGUMENT
    !     """"""""""""""""""""
    IMPLICIT NONE
    INTEGER ,         INTENT(IN)    :: NPOIN
    DOUBLE PRECISION ,INTENT(INOUT) :: W_LEG(NPOIN) , X_LEG(NPOIN)
    !
    !.....LOCAL VARIABLES
    !     """""""""""""""""
    INTEGER           I, M, J
    DOUBLE PRECISION  EPS, Z, P1, P2, P3, PP, Z1, PI
    parameter(eps=3.d-14)
    !
    pi = acos(-1.)
    m=(npoin+1)/2
    DO i=1,m
      z=cos(pi*(dble(i)-0.25d0)/(dble(npoin)+0.5d0))
1     CONTINUE
      p1=1.0d0
      p2=0.0d0
      DO j=1,npoin
        p3=p2
        p2=p1
        p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)
      ENDDO
      pp=dble(npoin)*(z*p1-p2)/(z*z-1.d0)
      z1=z
      z=z-p1/pp
      IF (abs(z-z1).GT.eps) GOTO 1
      x_leg(i)=-z
      x_leg(npoin+1-i)=z
      w_leg(i)=2.d0/((1.d0-z**2)*pp**2)
      w_leg(npoin+1-i)=w_leg(i)
    ENDDO
  END SUBROUTINE gauleg
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE f1f1f1(F1SF,NF1,IQ_OM1)
    ! TOMAWAC   V6P3                                   15/06/2011
    !***********************************************************************
    !
    !brief   SUBROUTINE CALLED BY PRENL3
    !+         COMPUTES VALUES OF RATIO F1/F AS FUNCTION OF THE IQ_OM1
    !+         INDICATOR
    !
    !history  E. GAGNAIRE-RENOU
    !+        04/2011
    !+        V6P1
    !+   CREATED
    !
    !history  G.MATTAROLO (EDF - LNHE)
    !+        15/06/2011
    !+        V6P1
    !+   Translation of French names of the variables in argument
    !
    !history  E. GAGNAIRE-RENOU
    !+        12/03/2013
    !+        V6P3
    !+   Better formatted: WRITE(LU,*), etc.
    !/ ------------------------------------------------------------------- /
    IMPLICIT NONE
    INTEGER,          INTENT(IN)    :: IQ_OM1
    INTEGER,          INTENT(INOUT) :: NF1
    DOUBLE PRECISION, INTENT(INOUT) :: F1SF(*)
    !
    INTEGER I,M
    DOUBLE PRECISION RAISON
    !
    IF(iq_om1.EQ.1) THEN
      IF(nf1.NE.14) THEN
        WRITE(*,*) '#1 Incorrect value for NF1',nf1
      ENDIF
      f1sf( 1)=0.30d0
      f1sf( 2)=0.40d0
      f1sf( 3)=0.50d0
      f1sf( 4)=0.60d0
      f1sf( 5)=0.70d0
      f1sf( 6)=0.80d0
      f1sf( 7)=0.90d0
      f1sf( 8)=1.00d0
      f1sf( 9)=1.11d0
      f1sf(10)=1.25d0
      f1sf(11)=1.42d0
      f1sf(12)=1.67d0
      f1sf(13)=2.00d0
      f1sf(14)=2.50d0
      f1sf(15)=3.30d0
    ELSEIF(iq_om1.EQ.2) THEN
      IF (nf1.NE.26) THEN
        WRITE(*,*) '#2 Incorrect value for NF1', nf1
      ENDIF
      f1sf( 1)=0.32d0
      f1sf( 2)=0.35d0
      f1sf( 3)=0.39d0
      f1sf( 4)=0.44d0
      f1sf( 5)=0.50d0
      f1sf( 6)=0.56d0
      f1sf( 7)=0.63d0
      f1sf( 8)=0.70d0
      f1sf( 9)=0.78d0
      f1sf(10)=0.86d0
      f1sf(11)=0.92d0
      f1sf(12)=0.97d0
      f1sf(13)=1.00d0
      f1sf(14)=1.03d0
      f1sf(15)=1.08d0
      f1sf(16)=1.13d0
      f1sf(17)=1.20d0
      f1sf(18)=1.28d0
      f1sf(19)=1.37d0
      f1sf(20)=1.48d0
      f1sf(21)=1.50d0
      f1sf(22)=1.65d0
      f1sf(23)=1.85d0
      f1sf(24)=2.10d0
      f1sf(25)=2.40d0
      f1sf(26)=2.70d0
      f1sf(27)=3.20d0
    ELSEIF(iq_om1.EQ.3) THEN
      IF(nf1.NE.11) THEN
        WRITE(*,*) 'Incorrect value for NF1', nf1
      ENDIF
      f1sf( 1)=0.30d0
      f1sf( 2)=0.48d0
      f1sf( 3)=0.64d0
      f1sf( 4)=0.78d0
      f1sf( 5)=0.90d0
      f1sf( 6)=1.00d0
      f1sf( 7)=1.12d0
      f1sf( 8)=1.28d0
      f1sf( 9)=1.50d0
      f1sf(10)=1.80d0
      f1sf(11)=2.40d0
      f1sf(12)=3.40d0
    ELSEIF(iq_om1.EQ.4) THEN
      IF(nf1.NE.40) THEN
        WRITE(*,*) 'Incorrect value for NF1', nf1
      ENDIF
      nf1=20
      m=10
      raison=9.d0**(1.d0/dble(nf1))
      f1sf(m+1)=1.0d0/3.0d0
      nf1=2*m+nf1
      DO i=m+2,nf1+1
        f1sf(i)=f1sf(i-1)*raison
      ENDDO
      DO i=m,1,-1
        f1sf(i)=f1sf(i+1)/raison
      ENDDO
    ELSEIF(iq_om1.EQ.5) THEN
      raison=9.d0**(1.d0/dble(nf1))
      f1sf(1)=1.d0/3.d0
      DO i=2,nf1+1
        f1sf(i)=f1sf(i-1)*raison
      ENDDO
    ELSEIF(iq_om1.EQ.6) THEN
      raison=(3.d0-1.d0/3.d0)/dble(nf1)
      f1sf(1)=1.d0/3.d0
      DO i=2,nf1+1
        f1sf(i)=f1sf(i-1)+raison
      ENDDO
    ELSEIF(iq_om1.EQ.7) THEN
      IF(nf1.NE.20) THEN
        WRITE(*,*) 'Incorrect value for NF1', nf1
      ENDIF
      f1sf( 1)=1.d0/3.d0
      f1sf( 2)=0.40d0
      f1sf( 3)=0.46d0
      f1sf( 4)=0.52d0
      f1sf( 5)=0.60d0
      f1sf( 6)=0.70d0
      f1sf( 7)=0.79d0
      f1sf( 8)=0.86d0
      f1sf( 9)=0.92d0
      f1sf(10)=0.97d0
      f1sf(11)=1.00d0
      f1sf(12)=1.04d0
      f1sf(13)=1.10d0
      f1sf(14)=1.18d0
      f1sf(15)=1.28d0
      f1sf(16)=1.42d0
      f1sf(17)=1.60d0
      f1sf(18)=1.84d0
      f1sf(19)=2.14d0
      f1sf(20)=2.52d0
      f1sf(21)=3.00d0
    ENDIF
    !
  END SUBROUTINE f1f1f1
  !/ ------------------------------------------------------------------- /
  SUBROUTINE insnlgqm
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |       E. Gagnaire-Renou &         |
    !/                  |       M. Benoit                   |
    !/                  |       S. Mostafa Siadatamousavi   |
    !/          |       M. Beyramzadeh              |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Nov-2022 |
    !/                  +-----------------------------------+
    !/
    !/    20-Nov-2022 : Merging with NL2 in WW3.            ( version 7.00 )
    !/
    !  1. Purpose :
    !
    !     Preprocessing for nonlinear interactions (Xnl).
    !
    !  2. Method :
    !
    !     See Xnl documentation.
    !
    !  3. Parameters :
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module      Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD    Subroutine tracing.
    !                Subr. GAULEG      Gauss-Legendre weights
    !      xnl_init  Subr. m_constants Xnl initialization routine.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3IOGR    Subr. W3IOGRMD Model definition file management.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     - See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: grav
    USE w3gdatmd,  ONLY: nk , nth , xfr , fr1, gqnf1, gqnt1, gqnq_om2, nltail, gqthrcou
 
#ifdef W3_S
    CALL strace (ient, 'INSNLGQM')
#endif
    IMPLICIT NONE
    !.....LOCAL VARIABLES
    INTEGER           JF    , JT    , JF1   , JT1   , NF1P1 , IAUX , NT , NF , IK
    INTEGER           IQ_TE1 , IQ_OM2 , LBUF , DIMBUF , IQ_OM1 , NQ_TE1 , NCONFM
 
    DOUBLE PRECISION  EPSI_A, AUX   , CCC   , DENO  , AAA   , DP2SG , TAILF
    DOUBLE PRECISION  V1    , V1_4  , DV1   , DTETAR , ELIM , RAISF
    DOUBLE PRECISION  V2    , V2_4  , V3    , V3_4
    DOUBLE PRECISION  W2    , W2_M  , W2_1  , W_MIL , W_RAD
    DOUBLE PRECISION  RK0   , XK0   , YK0   , RK1   , XK1   , YK1
    DOUBLE PRECISION  RK2   , XK2P  , YK2P  , XK2M  , YK2M
    DOUBLE PRECISION  RK3   , XK3P  , YK3P  , XK3M  , YK3M
    DOUBLE PRECISION  D01P  , C_D01P, S_D01P, D0AP  , C_D0AP, S_D0AP
    DOUBLE PRECISION  GA2P  , C_GA2P, S_GA2P, GA3P  , C_GA3P, S_GA3P, TWOPI, PI, SEUIL1 , SEUIL2 , SEUIL
    !
    !.....Variables related to the Gaussian quadratures
    DOUBLE PRECISION  W_CHE_TE1, W_CHE_OM2, C_LEG_OM2
    !
    !.....Variables related to the configuration selection
    DOUBLE PRECISION  TEST1 , TEST2
    DOUBLE PRECISION :: FREQ(NK)
    DOUBLE PRECISION, ALLOCATABLE :: F1SF(:) , X_CHE_TE1(:) , X_CHE_OM2(:) , X_LEG_OM2(:) , W_LEG_OM2(:) &
         ,  MAXCLA(:)
 
    pi = acos(-1.)
    lbuf = 500
    dimbuf = 2*lbuf+200
    twopi  = 2.*pi
    !
    ! Defines some threshold values for filtering (See Gagnaire-Renou Thesis,  p 52)
    !
    seuil1 = 1e10
    seuil2 = gqthrcou
 
    IF(gqnf1.EQ.14) iq_om1=1
    IF(gqnf1.EQ.26) iq_om1=2
    IF(gqnf1.EQ.11) iq_om1=3
    IF(gqnf1.EQ.40) iq_om1=4
    IF(gqnf1.EQ.11) iq_om1=3
    IF(gqnf1.EQ.40) iq_om1=4
    IF(gqnf1.EQ.20) iq_om1=7
    !
    ! Note by FA: not sure what the 5 and 6 cases correspond to
    !
    nq_te1 = gqnt1/2
    nconfm = gqnf1*gqnt1*gqnq_om2
 
    raisf = xfr
    nt = nth
    nf = nk
    dtetar = twopi/dble(nt)
 
    DO ik = 1,nk
      freq(ik) = fr1*raisf**(ik-1)
    ENDDO
 
    tailf = -nltail
 
    !===============ALLOCATE MATRICES=============================================
    if (Allocated(k_if2) ) then
      deallocate(k_if2)
    endif
    ALLOCATE(k_if2(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_if3) ) then
      deallocate(k_if3)
    endif
    ALLOCATE(k_if3(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_1p2p) ) then
      deallocate(k_1p2p)
    endif
    ALLOCATE(k_1p2p(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_1p3m) ) then
      deallocate(k_1p3m)
    endif
    ALLOCATE(k_1p3m(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_1p2m) ) then
      deallocate(k_1p2m)
    endif
    ALLOCATE(k_1p2m(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_1p3p) ) then
      deallocate(k_1p3p)
    endif
    ALLOCATE(k_1p3p(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_1m2p) ) then
      deallocate(k_1m2p)
    endif
    ALLOCATE(k_1m2p(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_1m3m) ) then
      deallocate(k_1m3m)
    endif
    ALLOCATE(k_1m3m(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_1m2m) ) then
      deallocate(k_1m2m)
    endif
    ALLOCATE(k_1m2m(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_1m3p) ) then
      deallocate(k_1m3p)
    endif
    ALLOCATE(k_1m3p(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(tb_v24) ) then
      deallocate(tb_v24)
    endif
    ALLOCATE(tb_v24(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(tb_v34) ) then
      deallocate(tb_v34)
    endif
    ALLOCATE(tb_v34(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(tb_tpm) ) then
      deallocate(tb_tpm)
    endif
    ALLOCATE(tb_tpm(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(tb_tmp) ) then
      deallocate(tb_tmp)
    endif
    ALLOCATE(tb_tmp(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(tb_fac) ) then
      deallocate(tb_fac)
    endif
    ALLOCATE(tb_fac(gqnq_om2,gqnt1,gqnf1))
 
    if (Allocated(k_if1) ) then
      deallocate(k_if1)
    endif
    ALLOCATE(k_if1(gqnf1))
 
    if (Allocated(k_1p) ) then
      deallocate(k_1p)
    endif
    ALLOCATE(k_1p(gqnt1,gqnf1))
 
    if (Allocated(k_1m) ) then
      deallocate(k_1m)
    endif
    ALLOCATE(k_1m(gqnt1,gqnf1))
 
    if (Allocated(tb_v14) ) then
      deallocate(tb_v14)
    endif
    ALLOCATE(tb_v14(gqnf1))
 
    if (Allocated(idconf) ) then
      deallocate(idconf)
    endif
    ALLOCATE(idconf(nconfm,3))
 
    !=======================================================================
    !     INITIALISATION OF AUXILIAIRY TABLES FOR SPECTRUM INTERPOLATION
    !=======================================================================
    if (Allocated(f_poin) ) then
      deallocate(f_poin)
    endif
    ALLOCATE(f_poin(dimbuf))
 
    if (Allocated(t_poin) ) then
      deallocate(t_poin)
    endif
    ALLOCATE(t_poin(dimbuf))
 
    if (Allocated(f_coef) ) then
      deallocate(f_coef)
    endif
    ALLOCATE(f_coef(dimbuf))
 
    if (Allocated(f_proj) ) then
      deallocate(f_proj)
    endif
    ALLOCATE(f_proj(dimbuf))
 
    if (Allocated(tb_sca) ) then
      deallocate(tb_sca)
    endif
    ALLOCATE(tb_sca(dimbuf))
 
 
    f_poin(:)=0
    t_poin(:)=0
    f_coef(:)=0.d0
    f_proj(:)=0.d0
    tb_sca(:)=0.0d0
 
    DO jf=1,lbuf
      f_poin(jf)=1
      f_coef(jf)=0.0d0
      f_proj(jf)=0.0d0
    ENDDO
    DO jf=1,nf
      iaux=lbuf+jf
      f_poin(iaux)=jf
      f_coef(iaux)=1.0d0
      f_proj(iaux)=1.0d0
    ENDDO
    aux=1.d0/raisf**tailf
    DO jf=1,lbuf
      iaux=lbuf+nf+jf
      f_poin(iaux)=nf
      f_coef(iaux)=aux**jf
      f_proj(iaux)=0.0d0
    ENDDO
    !
    DO jt=lbuf,1,-1
      t_poin(jt)=nt-mod(lbuf-jt,nt)
    ENDDO
    DO jt=1,nt
      t_poin(lbuf+jt)=jt
    ENDDO
    DO jt=1,lbuf
      t_poin(lbuf+nt+jt)=mod(jt-1,nt)+1
    ENDDO
    !======================================================================
    !
    !=======================================================================
    !     COMPUTES SCALE COEFFICIENTS FOR THE COUPLING COEFFICIENT
    !     Would be easier to pass these on from W3SRCE ???
    !=======================================================================
    dp2sg=twopi*twopi/grav
    DO jf=1,lbuf
      aux=freq(1)/raisf**(lbuf-jf+1)
      tb_sca(jf)=(dp2sg*aux**2)**6/(twopi**3*aux)
    ENDDO
    DO jf=1,nf
      tb_sca(lbuf+jf)=(dp2sg*freq(jf)**2)**6/(twopi**3*freq(jf))
    ENDDO
    DO jf=1,lbuf
      iaux=lbuf+nf+jf
      aux=freq(nf)*raisf**jf
      tb_sca(iaux)=(dp2sg*aux**2)**6/(twopi**3*aux)
    ENDDO
    !=======================================================================
    !
    !=======================================================================
    !     COMPUTES VALUES FOR GAUSSIAN QUADRATURES
    !=======================================================================
    if (Allocated(x_che_te1) ) then
      deallocate(x_che_te1)
    endif
    ALLOCATE(x_che_te1(1:nq_te1),x_che_om2(1:gqnq_om2))
 
    if (Allocated(x_leg_om2) ) then
      deallocate(x_leg_om2)
    endif
    ALLOCATE(x_leg_om2(1:gqnq_om2),w_leg_om2(1:gqnq_om2))
    !
    !.....Abscissa and weight (constant) for Gauss-Chebyshev
    DO iq_te1=1,nq_te1
      x_che_te1(iq_te1)=cos(pi*(dble(iq_te1)-0.5d0)/dble(nq_te1))
    ENDDO
    w_che_te1=pi/dble(nq_te1)
    DO iq_om2=1,gqnq_om2
      x_che_om2(iq_om2)=cos(pi*(dble(iq_om2)-0.5d0)/dble(gqnq_om2))
    ENDDO
    w_che_om2=pi/dble(gqnq_om2)
    !
    !.....Abscissa et weight for Gauss-Legendre
    CALL gauleg( w_leg_om2 , x_leg_om2 , gqnq_om2 )
    DO iq_om2=1,gqnq_om2
      x_leg_om2(iq_om2)=0.25d0*(1.d0+x_leg_om2(iq_om2))**2
    ENDDO
    !=======================================================================
    !
    !
    !=======================================================================
    !     COMPUTES VALUES OF RATIO F1/F AS FUNCTION OF THE IQ_OM1 INDICATOR
    !=======================================================================
    nf1p1=gqnf1+1
    if (Allocated(f1sf) ) then
      deallocate(f1sf)
    endif
    ALLOCATE(f1sf(1:nf1p1))
 
    CALL f1f1f1 ( f1sf  , gqnf1   , iq_om1)
    !=======================================================================
    !
    !     ==================================================
    !     STARTS LOOP 1 OVER THE RATIOS F1/F0
    !     ==================================================
    DO jf1=1,gqnf1
      !       ---------Computes and stores v1=f1/f0 and v1**4
      v1=(f1sf(jf1+1)+f1sf(jf1))/2.d0
      k_if1(jf1)=nint(dble(lbuf)+log(v1)/log(raisf))
      v1_4=v1**4
      tb_v14(jf1)=v1_4
      !       ---------Computes and stores dv1=df1/f0
      dv1=f1sf(jf1+1)-f1sf(jf1)
      !       ---------Computes the A parameter
      aaa=((1.d0+v1)**4-4.d0*(1.d0+v1_4))/(8.d0*v1**2)
      !
      !       =================================================
      !       STARTS LOOP 2 OVER THE DELTA_1+ VALUES
      !       =================================================
      DO jt1=1,gqnt1
        !
        !......Computes the Delta1+ values (=Theta_1-Theta_0) between 0 and Pi.
        IF (jt1.LE.nq_te1) THEN
          !           ---------First interval : X from -1 to A
          iq_te1=jt1
          c_d01p=(-1.d0+aaa)/2.d0+(1.d0+aaa)/2.d0*x_che_te1(iq_te1)
          ccc=dv1*sqrt((aaa-c_d01p)/(1.d0-c_d01p))*w_che_te1
        ELSE
          !           ---------Second interval : X from A to 1
          iq_te1=jt1-nq_te1
          c_d01p=( 1.d0+aaa)/2.d0+(1.d0-aaa)/2.d0*x_che_te1(iq_te1)
          ccc=dv1*sqrt((c_d01p-aaa)/(1.d0+c_d01p))*w_che_te1
        ENDIF
        s_d01p=sqrt(1.d0-c_d01p*c_d01p)
        d01p  =acos(c_d01p)
        k_1p(jt1,jf1)=lbuf+nint(d01p/dtetar)
        k_1m(jt1,jf1)=lbuf-nint(d01p/dtetar)
        !
        !         ---------Computes Epsilon_a
        epsi_a=2.d0*sqrt(1.d0+v1_4+2.d0*v1*v1*c_d01p)/(1.d0+v1)**2
        !         ---------Computes Delta_A+ and its cosinus
        c_d0ap=(1.d0-v1_4+0.25d0*epsi_a**2*(1.d0+v1)**4) &
             /(epsi_a*(1.d0+v1)**2)
        s_d0ap=sqrt(1.0d0-c_d0ap*c_d0ap)
        d0ap  = acos(c_d0ap)
        !
        !.......Integration over OMEGA2 depending on EPS_A
        IF (epsi_a.LT.1.d0) THEN
          !        - - - - - - - - - - - - - - - - - - - - - - - - - - - -
          !........Case of a single singularity (in OMEGA2-)
          !        - - - - - - - - - - - - - - - - - - - - - - - - - - - -
          w2_m=0.5d0*(1.d0-epsi_a/2.d0)
          w2_1=0.5d0
          !
          w_rad=w2_1-w2_m
          c_leg_om2=sqrt(w_rad)
          !
          !        ----------------------------------------------------
          !........STARTS LOOP 3 OVER OMEGA_2 (CASE Epsilon_A < 1)
          !........Case of a single singularity (in OMEGA2-)
          !........Integration over OMEGA2 via GAUSS-LEGENDRE quadrature
          !        ----------------------------------------------------
          DO iq_om2=1,gqnq_om2
            !             ---------Computes W2, V2, and V3
            w2=w2_m+w_rad*x_leg_om2(iq_om2)
            v2=w2*(1.d0+v1)
            v2_4=v2**4
            tb_v24(iq_om2,jt1,jf1)=v2_4
            k_if2(iq_om2,jt1,jf1) = nint(dble(lbuf) &
                 + log(v2)/log(raisf))
            v3=1.d0+v1-v2
            v3_4=v3**4
            tb_v34(iq_om2,jt1,jf1)=v3_4
            k_if3(iq_om2,jt1,jf1) = nint(dble(lbuf) &
                 + log(v3)/log(raisf))
            !             ---------Computes Gamma_2+ et Gamma_3+ angles
            c_ga2p=(epsi_a**2/4.d0+w2**4-(1.d0-w2)**4)/(epsi_a*w2*w2)
            c_ga2p=max(min(c_ga2p,1.d0),-1.d0)
            s_ga2p=sqrt(1.d0-c_ga2p*c_ga2p)
            ga2p  =acos(c_ga2p)
            c_ga3p=(epsi_a**2/4.d0-w2**4+(1.d0-w2)**4)/epsi_a &
                 /(1.d0-w2)**2
            c_ga3p=max(min(c_ga3p,1.d0),-1.d0)
            s_ga3p=sqrt(1.d0-c_ga3p*c_ga3p)
            ga3p  =acos(c_ga3p)
            !             Shifting of the direction indexes - Config. +Delta1 (SIG=1)
            k_1p2p(iq_om2,jt1,jf1)=nint(( d0ap+ga2p)/dtetar &
                 +dble(lbuf))
            k_1p3m(iq_om2,jt1,jf1)=nint(( d0ap-ga3p)/dtetar &
                 +dble(lbuf))
            k_1p2m(iq_om2,jt1,jf1)=nint(( d0ap-ga2p)/dtetar &
                 +dble(lbuf))
            k_1p3p(iq_om2,jt1,jf1)=nint(( d0ap+ga3p)/dtetar &
                 +dble(lbuf))
            !             Shifting of the direction indexes - Config. -Delta1 (SIG=-1)
            k_1m2p(iq_om2,jt1,jf1)=nint((-d0ap+ga2p)/dtetar &
                 +dble(lbuf))
            k_1m3m(iq_om2,jt1,jf1)=nint((-d0ap-ga3p)/dtetar &
                 +dble(lbuf))
            k_1m2m(iq_om2,jt1,jf1)=nint((-d0ap-ga2p)/dtetar &
                 +dble(lbuf))
            k_1m3p(iq_om2,jt1,jf1)=nint((-d0ap+ga3p)/dtetar &
                 +dble(lbuf))
            !
            !.........Computes the coupling coefficients (only for Delta_1+ )
            rk0=1.d0
            rk1=v1*v1
            rk2=v2*v2
            rk3=(1.d0+v1-v2)**2
            xk0  = rk0
            yk0  = 0.0d0
            xk1  = rk1*c_d01p
            yk1  = rk1*s_d01p
            xk2p = rk2*(c_d0ap*c_ga2p-s_d0ap*s_ga2p)
            yk2p = rk2*(s_d0ap*c_ga2p+c_d0ap*s_ga2p)
            xk2m = rk2*(c_d0ap*c_ga2p+s_d0ap*s_ga2p)
            yk2m = rk2*(s_d0ap*c_ga2p-c_d0ap*s_ga2p)
            xk3p = rk3*(c_d0ap*c_ga3p-s_d0ap*s_ga3p)
            yk3p = rk3*(s_d0ap*c_ga3p+c_d0ap*s_ga3p)
            xk3m = rk3*(c_d0ap*c_ga3p+s_d0ap*s_ga3p)
            yk3m = rk3*(s_d0ap*c_ga3p-c_d0ap*s_ga3p)
            tb_tpm(iq_om2,jt1,jf1)=couple( xk0   , yk0   , xk1   , yk1   , xk2p  , yk2p  , xk3m  , yk3m)
            tb_tmp(iq_om2,jt1,jf1)=couple( xk0   , yk0   , xk1   , yk1   , xk2m  , yk2m  , xk3p  , yk3p)
            !
            !.........Computes the multiplicative coefficient for QNL4
            deno=2.d0*sqrt( (0.5d0*(1.d0+epsi_a/2.d0)-w2) &
                 *((w2-0.5d0)**2+0.25d0*(1.d0+epsi_a)) &
                 *((w2-0.5d0)**2+0.25d0*(1.d0-epsi_a)) )
            tb_fac(iq_om2,jt1,jf1)=1.d0/(deno*v1*w2*(1.d0-w2)) &
                 /(1.d0+v1)**5 * w_leg_om2(iq_om2)*c_leg_om2* ccc
          ENDDO
          !        -----------------------------------------------
          !........END OF THE LOOP 3 OVER OMEGA_2 (CASE Epsilon_A < 1)
          !        -----------------------------------------------
          !
        ELSE
          !        - - - - - - - - - - - - - - - - - - - - - - - - - - - -
          !........STARTS LOOP 3 OVER OMEGA_2 (CASE Epsilon_A > 1)
          !........Case of two singularities (in OMEGA2- and OMEGA2_1)
          !........Integration over OMEGA2 via GAUSS-CHEBYSCHEV quadrature
          !        - - - - - - - - - - - - - - - - - - - - - - - - - - - -
          w2_m=0.5d0*(1.d0-epsi_a/2.d0)
          w2_1=0.5d0*(1.d0-sqrt(epsi_a-1.d0))
          !
          w_mil=(w2_m+w2_1)/2.d0
          w_rad=(w2_1-w2_m)/2.d0
          !
          DO iq_om2=1,gqnq_om2
            !             ---------Computes W2, V2, and V3
            w2=w_mil+w_rad*x_che_om2(iq_om2)
            v2=w2*(1.d0+v1)
            v2_4=v2**4
            tb_v24(iq_om2,jt1,jf1)=v2_4
            k_if2(iq_om2,jt1,jf1)=nint(dble(lbuf) &
                 +log(v2)/log(raisf))
            v3=1.d0+v1-v2
            v3_4=v3**4
            tb_v34(iq_om2,jt1,jf1)=v3_4
            k_if3(iq_om2,jt1,jf1)=nint(dble(lbuf) &
                 +log(v3)/log(raisf))
            !             ---------Computes Gamma_2+ et Gamma_3+ angles
            c_ga2p=(epsi_a**2/4.d0+w2**4-(1.d0-w2)**4)/(epsi_a*w2*w2)
            c_ga2p=max(min(c_ga2p,1.d0),-1.d0)
            s_ga2p=sqrt(1.d0-c_ga2p*c_ga2p)
            ga2p  =acos(c_ga2p)
            c_ga3p=(epsi_a**2/4.d0-w2**4+(1.d0-w2)**4)/epsi_a &
                 /(1.d0-w2)**2
            c_ga3p=max(min(c_ga3p,1.d0),-1.d0)
            s_ga3p=sqrt(1.d0-c_ga3p*c_ga3p)
            ga3p  =acos(c_ga3p)
            !             Shifts the direction indexes - Config. +Delta1 (SIG=1)
            k_1p2p(iq_om2,jt1,jf1)=nint(( d0ap+ga2p)/dtetar &
                 +dble(lbuf))
            k_1p3m(iq_om2,jt1,jf1)=nint(( d0ap-ga3p)/dtetar &
                 +dble(lbuf))
            k_1p2m(iq_om2,jt1,jf1)=nint(( d0ap-ga2p)/dtetar &
                 +dble(lbuf))
            k_1p3p(iq_om2,jt1,jf1)=nint(( d0ap+ga3p)/dtetar &
                 +dble(lbuf))
            !             Shifts the direction indexes - Config. -Delta1 (SIG=-1)
            k_1m2p(iq_om2,jt1,jf1)=nint((-d0ap+ga2p)/dtetar &
                 +dble(lbuf))
            k_1m3m(iq_om2,jt1,jf1)=nint((-d0ap-ga3p)/dtetar &
                 +dble(lbuf))
            k_1m2m(iq_om2,jt1,jf1)=nint((-d0ap-ga2p)/dtetar &
                 +dble(lbuf))
            k_1m3p(iq_om2,jt1,jf1)=nint((-d0ap+ga3p)/dtetar &
                 +dble(lbuf))
            !
            !.........Computes the coupling coefficients (only for Delta_1+ )
            rk0=1.d0
            rk1=v1*v1
            rk2=v2*v2
            rk3=(1.d0+v1-v2)**2
            xk0  = rk0
            yk0  = 0.0d0
            xk1  = rk1*c_d01p
            yk1  = rk1*s_d01p
            xk2p = rk2*(c_d0ap*c_ga2p-s_d0ap*s_ga2p)
            yk2p = rk2*(s_d0ap*c_ga2p+c_d0ap*s_ga2p)
            xk2m = rk2*(c_d0ap*c_ga2p+s_d0ap*s_ga2p)
            yk2m = rk2*(s_d0ap*c_ga2p-c_d0ap*s_ga2p)
            xk3p = rk3*(c_d0ap*c_ga3p-s_d0ap*s_ga3p)
            yk3p = rk3*(s_d0ap*c_ga3p+c_d0ap*s_ga3p)
            xk3m = rk3*(c_d0ap*c_ga3p+s_d0ap*s_ga3p)
            yk3m = rk3*(s_d0ap*c_ga3p-c_d0ap*s_ga3p)
            tb_tpm(iq_om2,jt1,jf1)=couple( xk0   , yk0   , xk1   , yk1   , xk2p  , yk2p  , xk3m  , yk3m)
            tb_tmp(iq_om2,jt1,jf1)=couple( xk0   , yk0   , xk1   , yk1   , xk2m  , yk2m  , xk3p  , yk3p)
            !
            !.........Computes the multiplicative coefficient for QNL4
            deno=2.d0*sqrt( (0.5d0*(1.d0+epsi_a/2.d0)-w2) &
                 *((w2-0.5d0)**2+0.25d0*(1.d0+epsi_a)) &
                 *(0.5d0*(1.d0+sqrt(epsi_a-1.d0))-w2) )
            tb_fac(iq_om2,jt1,jf1)=1.d0/(deno*v1*w2*(1.d0-w2)) &
                 /(1.d0+v1)**5 * w_che_om2* ccc
            !
          ENDDO
          !        -----------------------------------------------
          !........END OF LOOP 3 OVER OMEGA_2 (CASE Epsilon_A > 1)
          !        -----------------------------------------------
          !
        ENDIF
      ENDDO
      !       =================================================
      !       END OF LOOP 2 OVER THE DELTA_1+ VALUES
      !       =================================================
      !
    ENDDO
    !     ==================================================
    !     END OF LOOP 1 OVER THE F1/F0 RATIOS
    !     ==================================================
    DEALLOCATE(f1sf)
    DEALLOCATE(x_che_te1)
    DEALLOCATE(x_che_om2)
    DEALLOCATE(x_leg_om2)
    DEALLOCATE(w_leg_om2)
 
    !     ===========================================================
    !     POST-PROCESSING TO ELIMINATE PART OF THE CONFIGURATIONS
    !     ===========================================================
    !
    !.....It looks, for every value of the ratio V1, for the maximum value
    !.....of FACTOR*COUPLING : it is stored in the local table NAXCLA(.)
    !     """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
    ALLOCATE(maxcla(1:gqnf1))
    DO jf1=1,gqnf1
      aux=0.0d0
      DO jt1=1,gqnt1
        DO iq_om2=1,gqnq_om2
          aux=max(aux,tb_fac(iq_om2,jt1,jf1)*tb_tpm(iq_om2,jt1,jf1),tb_fac(iq_om2,jt1,jf1)*tb_tmp(iq_om2,jt1,jf1))
        ENDDO
      ENDDO
      maxcla(jf1)=aux
    ENDDO
    !
    !.....It looks for the max V1 value
    !     """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""
    aux=0.0d0
    DO jf1=1,gqnf1
      IF (maxcla(jf1).GT.aux) aux=maxcla(jf1)
    ENDDO
 
    test1=seuil1*aux
    !
    !.....Set to zero the coupling coefficients not used
    !     """""""""""""""""""""""""""""""""""""""""""""""""""""
    nconf=0
    DO jf1=1,gqnf1
      test2 =seuil2*maxcla(jf1)
      DO jt1=1,gqnt1
        DO iq_om2=1,gqnq_om2
          aaa=tb_fac(iq_om2,jt1,jf1)*tb_tpm(iq_om2,jt1,jf1)
          ccc=tb_fac(iq_om2,jt1,jf1)*tb_tmp(iq_om2,jt1,jf1)
          IF ((aaa.GT.test1.OR.aaa.GT.test2).OR. &
               (ccc.GT.test1.OR.ccc.GT.test2)) THEN
            nconf=nconf+1
            idconf(nconf,1)=jf1
            idconf(nconf,2)=jt1
            idconf(nconf,3)=iq_om2
          ENDIF
#ifdef W3_TGQM
          WRITE(993,*) nconf,jf1,jt1,iq_om2,aaa,ccc,(aaa.GT.test1.OR.aaa.GT.test2), &
               (ccc.GT.test1.OR.ccc.GT.test2)
#endif
        ENDDO
      ENDDO
    ENDDO
    DEALLOCATE(maxcla)
    !
    !..... counts the fraction of the eliminated configurations
    elim=(1.d0-dble(nconf)/dble(nconfm))*100.d0
#ifdef W3_TGQM
    WRITE(994,*) 'NCONF, ELIM FRACTION:',nconf,elim
#endif
  END SUBROUTINE insnlgqm
  !/
  !/ End of module W3SNL1MD -------------------------------------------- /
  !/
END MODULE w3snl1md