w3canomd.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3canomd
  !/
  !/                  +-----------------------------------+
  !/                  |                                   |
  !/                  |   P.A.E.M. Janssen                |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         21-Aug-2014 |
  !/                  +-----------------------------------+
  !/
  !/    XX-Jul-2010 : Origination by  PAEM JANSSEN
  !/    18-Oct-2012 : Adapted to WAVEWATCH III: F. Ardhuin( version 4.07 )
  !/    21-Aug-2014 : Bug corrected: only first call wasOK( version 5.01 )
  !/
  !  0. Note by F. Ardhuin:
  !     In adapting the orginal program to be a WAVEWATCH module, I
  !     have so far strived to keep the original code. As a result
  !     some routines are unnecessarily duplicated (e.g. the calculation of
  !     the group velocity ...). But this can improve the traceability of
  !     the code.
  !     The first spectrum (JONSWAP) has been removed from the code
  !
  !  1. Purpose :
  !
  !
  !     CALCULATION OF THE SECOND ORDER CORRECTION TO THE SURFACE GRAVITY
  !     WAVE SPECTRUM
  !
  !     DOCUMENTATION.
  !     -------------
  !
  !     PRESENTLY, THE SOFTWARE IS SET UP TO DO FOR A GIVING FIRST-ORDER
  !     SPECTRUM AT A GIVEN DEPTH THE DETERMINATION OF THE SECOND ORDER
  !     CORRECTION (INCLUDING SECOND-HARMONICS, WAVE SET DOWN AND DOPPLER SHIFT
  !     OWING TO THE STOKES FREQUENCY CORRECTION.
  !
  !     EVALUATION OF THE INTERACTION COEFFICIENTS FOR ARBRITRARY DEPTH WOULD
  !     BE VERY TIME CONSUMING. THEREFORE, THE APPROACH IN THE WAM MODEL IS
  !     FOLLOWED, WHERE TABLES ARE GENERATED FOR A LOGARITHMIC DEPTH TABLE
  !
  !               D(JD) = DEPTHA*DEPTHD**(JD-1)
  !
  !     WITH JD AN INTEGER. IN THE PRESENT OPERATIONAL VERSION OF ECWAM JD
  !     RANGES FROM 1 TO NDEPTH = 74, WHILE DEPTHA = 1. AND DEPTHD = 1.1
  !
  !     FINALLY, THIS IS A VERY TIME-CONSUMING CALCULATION, AT LEAST FOR AN
  !     OPERATIONAL MODEL. i HAVE THEREFORE INTRODUCED THE OPTION THAT THE SECOND-ORDER
  !     SPECTRUM IS CALCULATED ON A LOWER RESOLUTION GRID (TYPICALLY HALF THE
  !     RESOLUTION) WHILE THE INFORMATION CONTAINED IN THE FIRST-ORDER SPECTRUM
  !     IS KEPT ON THE ORIGINAL SPECTRAL GRID.
  !
  ! ----------------------------------------------------------------------
  !
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SREF    Subr. Public   Reflection of waves (shorline, islands...)
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !
  !  6. Switches :
  !
  !     !/S  Enable subroutine tracing.
  !
  !  7. Source code :
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !
  !/
  !/ Public variables
  !/
 
  REAL                     :: g, pi, zpi, rad, deg
  INTEGER                  :: ndepth
  REAL                     :: deptha         ! first depth in table
  REAL, SAVE , PRIVATE, ALLOCATABLE   :: omega(:)
#ifdef W3_OMPG
  !$omp threadprivate( OMEGA )
#endif
  INTEGER, SAVE , PRIVATE             :: counter = 0
#ifdef W3_OMPG
  !$omp threadprivate( COUNTER )
#endif
  ! Tables for non-linear coefficients ...
  REAL, SAVE , PRIVATE, ALLOCATABLE   :: ta(:,:,:,:),tb(:,:,:,:),tc_ql(:,:,:,:),&
       tt_4m(:,:,:,:),tt_4p(:,:,:,:),tfakh(:,:),    &
       tfak(:,:)
#ifdef W3_OMPG
  !$omp threadprivate( TA, TB, TC_QL, TT_4M, TT_4P, TFAKH, TFAK )
#endif
  INTEGER, SAVE, PRIVATE, ALLOCATABLE :: im_p(:,:),im_m(:,:)
#ifdef W3_OMPG
  !$omp threadprivate( IM_P, IM_M )
#endif
 
 
  !
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3add2ndorder(E,DEPTH,WN,CG,IACTION)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |            F. Ardhuin             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         19-Oct-2012 |
    !/                  +-----------------------------------+
    !/
    !/    19-Oct-2012 : Origination                         ( version 4.08 )
    !/
    !  1. Purpose :
    !
    !     Adds second order spectrum on top of first order spectrum
    !
    !  2. Method :
    !
    !     Uses P. Janssen's code for the inverse canonical transform
    !
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A         R.A.  I   Action density spectrum (1-D)
    !       CG        R.A.  I   Group velocities.
    !       WN        R.A.  I   Wavenumbers.
    !       DEPTH     Real  I   Mean water depth.
    !       S         R.A.  O   Source term (1-D version).
    !       D         R.A.  O   Diagonal term of derivative (1-D version).
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !       E         R.A. I/O   Energy density spectrum (1-D), f-theta
    !       DEPTH     Real I     Water depth
    !       WN        R.A.       wavenumbers
    !       CG        R.A.       group velocities
    !       IACTION   Int  I     Switch to specify if the input spectrum
    !                            is E(f,theta) or A(k,theta)
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SREF    Subr. W3REF1MD Shoreline reflection source term
    !      W3EXPO    Subr.   N/A    Point output post-processor.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: grav
    USE w3dispmd
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, th, dth, igpars
 
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(INOUT)     :: E(NSPEC)
    REAL, INTENT(IN)        :: DEPTH
    REAL, INTENT(IN)        :: WN(NK)
    REAL, INTENT(IN)        :: CG(NK)
    INTEGER, INTENT(IN)     :: IACTION
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER           :: ISPEC, IK, ITH, M
    REAL              :: CO1, ATOE, DPTH
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    LOGICAL, SAVE     :: FIRST = .true.
#ifdef W3_OMPG
    !$omp threadprivate( FIRST )
#endif
    REAL, ALLOCATABLE, SAVE  :: FR(:), DFIM(:)
    REAL, ALLOCATABLE, SAVE  :: F1(:,:), F3(:,:)
#ifdef W3_OMPG
    !$omp threadprivate( FR, DFIM, F1, F3 )
#endif
    INTEGER, SAVE     ::  NFRE, NANG
    INTEGER, SAVE     ::  NFREH, NANGH
#ifdef W3_OMPG
    !$omp threadprivate( NFRE, NANG, NFREH, NANGH )
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3ADD2NDORDER')
#endif
    !
    ! 0.  Initializations ------------------------------------------------ *
    !
    IF (first) THEN
      first=.false.
      nfre=nk
      nang=nth
      nfreh=nk
      nangh=nth
      g=grav
      pi = 4.*atan(1.)
      zpi=2*pi
      rad = pi/180.
      deg = 180./pi
      ALLOCATE(fr(nfre), dfim(nfre))
      fr(1:nfre)=sig(1:nk)/zpi
      ! The following can be replaced using DSIP from WWATCH
      co1 = 0.5*dth
      dfim(1)= co1*(fr(2)-fr(1))
      DO m=2,nfre-1
        dfim(m)=co1*(fr(m+1)-fr(m-1))
      ENDDO
      dfim(nfre)=co1*(fr(nfre)-fr(nfre-1))
      !
      ALLOCATE(f1(nang,nfre), f3(nang,nfre))
      ndepth=igpars(6)
      deptha=igpars(7)
    END IF
    dpth = depth
 
    DO ik=1,nk
      IF (iaction.EQ.0) THEN
        atoe=1
      ELSE
        atoe=sig(ik)*zpi / cg(ik)
      END IF
      DO ith=1,nth
        ispec=ith+(ik-1)*nth
        f1(ith,ik)=e(ispec)*atoe
      END DO
      !WRITE(100,'(100G16.8)') SIG(IK)*ZPI,(F1(ITH,IK),ITH=1,NTH)
 
    END DO
    !
    ! 1. DETERMINE SECOND-ORDER SPECTRUM.
    !
 
    CALL cal_sec_order_spec(f1,f3,nfre,nang,fr,dfim,th,   &
         dth,dpth,+1., nfreh, nangh)
 
    !
    ! 2. Adds 2nd order spectrum to 1st order
    !
    DO ik=1,nk
      IF (iaction.EQ.0) THEN
        atoe=1
      ELSE
        atoe=sig(ik)*zpi / cg(ik)
      END IF
      DO ith=1,nth
        ispec=ith+(ik-1)*nth
        e(ispec)=f3(ith,ik)/atoe
      END DO
      !WRITE(101,'(I3,100G16.8)') SIG(IK)*ZPI,(F3(ITH,IK),ITH=1,NTH)
    END DO
 
#ifdef W3_T
    print*,' END CAL_SEC_ORDER_SPEC'
#endif
    RETURN
 
  END SUBROUTINE w3add2ndorder
  !/ ------------------------------------------------------------------- /
 
  !-----------------------------------------------------------------------
  !
  SUBROUTINE cal_sec_order_spec(F1,F3,NFRE,NANG,FR,DFIM,TH,DELTH, &
       DPTH,SIGM, NFREH, NANGH)
    !
    !***  *CAL_SEC_ORDER_SPEC*   DETERMINES SECOND_ORDER SPECTRUM
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              DETERMINATION OF SECOND-ORDER SPECTRUM
    !
    !     INTERFACE.
    !     ----------
    !              *CALL*  *CAL_SEC_ORDER_SPEC(F1,F3,NFRE,NANG,FR,
    !                                DFIM,TH,DELTH,DPTH,SIGM)*
    !
    !                       INPUT:
    !                            *F1*    - 2-D FREE WAVE SPECTRUM
    !                            *NFRE*  - NUMBER OF FREQUENCIES
    !                            *NANG*  - NUMBER OF DIRECTIONS
    !                            *FR*    - FREQUENCIES
    !                            *DFIM*  - FREQUENCY INCREMENT
    !                            *TH*    - DIRECTIONAL ARRAY
    !                            *DELTH* - DIRECTIONAL INCREMENT
    !                            *DPTH*  - DEPTH ARRAY
    !                            *SIGM*   - FOR SIGM = 1 FORWARD MAPPING
    !                                      WHILE FOR SIGM = -1 INVERSE
    !                                      MAPPING.
    !
    !                       OUTPUT:
    !                            *F3*    - 2-D SPECTRUM INCLUDING SECOND-ORDER
    !                                      CORRECTION
    !
    !     METHOD.
    !     -------
    !              IS DESCRIBED IN JANSSEN (2009), JFM, 637, 1-44.
    !
    !     EXTERNALS.
    !     ----------
    !              NONE
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
 
    REAL, INTENT(IN)        :: F1(NANG,NFRE)
    REAL, INTENT(OUT)       :: F3(NANG,NFRE)
 
    INTEGER, INTENT(IN)     :: NFRE,NANG,NFREH, NANGH
 
    REAL, INTENT(IN)        :: DFIM(NFRE),FR(NFRE), TH(NANG), DELTH
    REAL, INTENT(IN)        :: DPTH, SIGM
 
    LOGICAL FRSTIME,DOUBLEP
 
    INTEGER MDW,M,K, K0,M0,MP,KP,MM,KM,KL,KLL,ML,JD
    INTEGER, SAVE :: MR, MA,NMAX
#ifdef W3_OMPG
    !$omp threadprivate( MR, MA, NMAX )
#endif
 
    !      PARAMETER (NFREH=32,NANGH=36)
 
    INTEGER, SAVE            :: INDEP
#ifdef W3_OMPG
    !$omp threadprivate( INDEP )
#endif
    REAL,ALLOCATABLE         :: PF1(:,:),PF3(:,:)
 
 
    REAL DEPTH,ALPHA,GAM_J,DEPTHD
    REAL OM0,AA1,BB1,&
         F,EPSMIN,DELFF,SPEC1,SQRTK
    REAL FRAC,DEL,DELF,D1,D2,D3,D4,C1,&
         C2,XM,XK
    REAL, SAVE :: OMSTART
    REAL, SAVE :: XMR,XMA, DELTHH, CO1
#ifdef W3_OMPG
    !$omp threadprivate( OMSTART, XMR,XMA, DELTHH, CO1 )
#endif
    REAL :: F13(NFREH,NANGH)
    REAL :: SUM0,AKMEAN
    REAL :: DELOM(NFREH),THH(NANGH),DFDTH(NFREH)
 
    DATA frstime/.true./
 
    common/const/depth,alpha,mdw,gam_j,depthd
    common/precis/doublep
 
    !
    !***  2. DETERMINE SECOND ORDER CORRECTION TO THE SPECTRUM
    !     ----------------------------------------------------
 
    !
#ifdef W3_T
    print*,' START SECOND-ORDER CALC.'
#endif
 
    doublep = .true.
    !
    !***  2.1 SET UP OF LOW-RESOLUTION CALCULATION GRID.
    !     ---------------------------------------------
    !
    epsmin = 1.0e-4
    frac = 0.1
    omstart = zpi*fr(1)
    mr = max(1,nfre/nfreh)
    xmr = 1./float(mr)
    ma = nang/nangh
    xma = 1./float(ma)
    delthh = float(ma)*delth
 
    IF (frstime) THEN
      ! IF (COUNTER.GT.0) THEN
      !  DEALLOCATE(OMEGA,TFAK,TA,TB,TC_QL,TT_4M,TT_4P,IM_P,IM_M,TFAKH)
      ! ENDIF
      ALLOCATE(omega(nfreh))
      ALLOCATE(tfak(nfre,ndepth))
      ALLOCATE(ta(nangh,nfreh,nfreh,ndepth))
      ALLOCATE(tb(nangh,nfreh,nfreh,ndepth))
      ALLOCATE(tc_ql(nangh,nfreh,nfreh,ndepth))
      ALLOCATE(tt_4m(nangh,nfreh,nfreh,ndepth))
      ALLOCATE(tt_4p(nangh,nfreh,nfreh,ndepth))
      ALLOCATE(im_p(nfreh,nfreh))
      ALLOCATE(im_m(nfreh,nfreh))
      ALLOCATE(tfakh(nfreh,ndepth))
 
      DO m=1,nfreh
        omega(m) = zpi*fr(mr*m)
      ENDDO
 
      DO k=1,nangh
        k0 = ma*k+1
        IF (k0.GT.nang) k0 = k0-nang
        thh(k) = th(k0)
      ENDDO
 
      co1   = 1./2.*delthh
      delom(1) = co1*(omega(2)-omega(1))
      DO m=2,nfreh-1
        delom(m)=co1*(omega(m+1)-omega(m-1))
      ENDDO
      delom(nfreh)=co1*(omega(nfreh)-omega(nfreh-1))
      !
      dfdth = delom/zpi
      !
      !***  2.2 INITIALISE TABLES
      !     ---------------------
      !
      nmax = xmr*(1+nint(log(2.*omega(nfreh)/omstart)/log(1.+frac)))
      nmax = nmax+1
#ifdef W3_T
      print*,' NMAX = ',nmax
#endif
 
      depthd = 1.1
 
      DO jd=1,ndepth
        depth = deptha*depthd**(jd-1)
        DO m=1,nfre
          om0 = zpi*fr(m)
          tfak(m,jd) = aki(om0,depth)
        ENDDO
      ENDDO
 
      indep = 1+nint(log(dpth/deptha)/log(depthd))
      indep = min(ndepth,indep)
      indep = max(1,indep)
 
      CALL tables_2nd(nfreh,nangh,ndepth,deptha,omstart,frac,xmr,&
           dfdth,omega,thh)
      print*, '2ND ORDER TABLES GENERATED:',ndepth,deptha, delthh
 
      frstime = .false.
    ENDIF ! end of test on FRSTIME
    !
    counter=counter+1
    !
    !***  DETERMINE SOME MOMENTS.
    !     ----------------------
    !
    sum0 = 0.
    akmean = 0.
    DO m=1,nfre
      DO k=1,nang
        sqrtk=sqrt(tfak(m,indep))
        sum0 = sum0+f1(k,m)*dfim(m)
        akmean = akmean+f1(k,m)*dfim(m)/sqrtk
      ENDDO
    ENDDO
    !
    ! NB: AKMEAN is the mean wavenumber corresponding to Tm0,-1 in deep water
    !
    akmean = (sum0/akmean)**2
 
    !
    !***  2.2 INTERPOLATION OR NOT.
    !     ------------------------
    !
    IF (mr.EQ.1 .AND. ma.EQ.1) THEN
      !
      !***     2.21 NO INTERPOLATION.
      !        ----------------------
      !
#ifdef W3_T
      print*,' NO THINNING AND INTERPOLATION'
      print*,'nanG:',nang,nmax,nfre,ndepth,deptha,depthd,dpth,'##',delth,delthh
#endif
 
      CALL secspom(f1,f3,nfre,nang,nmax,ndepth,&
           deptha,depthd,omstart,frac,mr,dfdth,omega,&
           dpth,akmean,ta,tb,tc_ql,tt_4m,tt_4p,&
           im_p,im_m,counter)
      DO m=1,nfre
        DO k=1,nang
          delf = f3(k,m)
          f3(k,m)=max(0.00000001,f1(k,m)+sigm*delf)
        ENDDO
      ENDDO
 
    ELSE
 
      !
      !***     2.22 ENERGY CONSERVING INTERPOLATION SCHEME
      !        -------------------------------------------
      !
      print*,' !THINNING AND INTERPOLATION!'
      ALLOCATE(pf1(nangh,nfreh))
      ALLOCATE(pf3(nangh,nfreh))
 
      pf1 = 0.
      DO m=1,nfreh
        DO k=1,nangh
          m0 = mr*m
          mp = m0+1
          mp = min(nfre,mp)
          mm = m0-1
 
          k0 = ma*k+1
          kp = k0+1
          km = k0-1
          delff = 0.
          DO kl = km,kp
            kll = kl
            IF (kll.GT.nang) kll = kll-nang
            IF (kll.LT.1) kll = kll+nang
            DO ml = mm,mp
              del = dfim(ml)
              delff = delff+del
              spec1 = f1(kll,ml)
              pf1(k,m)=pf1(k,m)+spec1*del
            ENDDO
          ENDDO
          pf1(k,m) =pf1(k,m)/delff
        ENDDO
      ENDDO
      !
      !***     2.23 DETERMINE SECOND-ORDER SPEC
      !        --------------------------------
      !
      CALL secspom(pf1,pf3,nfreh,nangh,nmax,ndepth,&
           deptha,depthd,omstart,frac,mr,dfdth,omega,&
           dpth,akmean,ta,tb,tc_ql,tt_4m,tt_4p,&
           im_p,im_m,counter)
      !
      !***     2.24 INTERPOLATE TOWARDS HIGH-RES GRID
      !        --------------------------------------
      !
      DO m=1,nfre
        DO k=1,nang
          xm = real(m/mr)
          xk = real((k-1)/ma)
 
          m0 = max(1,int(xm))
          k0 = int(xk)
 
          d1 = real(m)/real(mr)-xm
          d2 = 1.-d1
          d3 = real(k-1)/real(ma)-xk
          d4 = 1.-d3
 
          IF (k0.LT.1) k0 = k0+nangh
          mp = min(nfreh,m0+1)
          kp = k0+1
          IF (kp.GT.nangh) kp = kp-nangh
 
          c1 = pf3(k0,m0)*d4+pf3(kp,m0)*d3
          c2 = pf3(kp,mp)*d3+pf3(k0,mp)*d4
 
          delf = c1*d2+c2*d1
          f3(k,m)=max(0.00000001,f1(k,m)+sigm*delf)
        ENDDO
      ENDDO
 
    ENDIF
 
 
    IF (mr.GT.1 .OR. ma.GT.1 ) THEN
      DO m=1,nfreh
        aa1 = 0.
        DO k=1,nangh
          aa1 = aa1+pf1(k,m)*delthh
        ENDDO
        aa1 = max(aa1,epsmin)
 
        bb1 = 0.
        DO k=1,nangh
          bb1 = bb1+(pf1(k,m)+pf3(k,m))*delthh
        ENDDO
        bb1 = max(bb1,epsmin)
        f   = omega(m)/zpi
 
#ifdef W3_T
        WRITE(6,62) m,f,aa1,bb1,delthh
        WRITE(80,62) m,f,aa1,bb1,delthh
#endif
      ENDDO
 
      DO m=1,nfreh
        DO k=1,nangh
          f13(m,k)=pf1(k,m)+pf3(k,m)
        ENDDO
      ENDDO
    ENDIF
 
    !
#ifdef W3_T
62  FORMAT(i4,9f16.9)
#endif
    !
    RETURN
  END SUBROUTINE cal_sec_order_spec
  !
  !--------------------------------------------------------------------
  !
  SUBROUTINE tables_2nd(NFRE,NANG,NDEPTH,DEPTHA,OMSTART,FRAC,XMR,&
       DFDTH,OMEGA,TH)
    !
    !--------------------------------------------------------------------
    !
    !*****TABLES** COMPUTES TABLES FOR SECOND ORDER SPECTRUM IN FREQUENCY SPACE.
    !
    !     P.JANSSEN DECEMBER 2008
    !
    !     PURPOSE
    !     -------
    !             DETERMINES TABLES, BASED ON JANSSEN (2008)
    !             THERE ARE THREE CORRECTIONS:
    !                   1) GENERATION OF SECOND-HARMONICS
    !                   2) QUASI-LINEAR EFFECT
    !                   3) SHIFT OF SPECTRUM BECAUSE OF STOKES FREQUENCY
    !                      CORRECTION.
    !
    !     INTERFACE
    !     ---------
    !             *CALL* *TABLES(NFRE,NANG,NDEPTH,OMSTART,FRAC,XMR,
    !                            OMEGA,TA,TB,TC_QL,TT_4M,TT_4P,IM_P,IM_M,
    !                            TFAK)*
    !
    !
    !     PARAMETER   TYPE      PURPOSE.
    !     ---------   ----      -------
    !
    !       NFRE      INTEGER   NUMBER OF FREQUENCIES
    !       NANG      INTEGER   NUMBER OF DIRECTIONS
    !       NDEPTH    INTEGER   NUMBER OF ENTRIES IN THE DEPTH TABLE
    !       OMSTART   REAL      START FREQUENCY
    !       FRAC      REAL      FRACTIONAL INCREASE IN FREQUENCY SPACE
    !       XMR       REAL      INVERSE OF THINNING FACTOR IN FREQUENCY SPACE
    !       OMEGA     REAL      ANGULAR FREQUENCY ARRAY
    !       DFDTH     REAL      PRODUCT OF INCREMENT IN FREQUENCY AND DIRECTION
    !       TH        REAL      DIRECTION ARRAY
    !       TA        REAL      TABLE FOR MINUS INTERACTIONS
    !       TB        REAL      TABLE FOR PLUS INTERACTIONS
    !       TC_QL     REAL      TABLE FOR QUASI-LINEAR INTERACTIONS
    !       TT_4M     REAL      TABLE FOR STOKES FREQUENCY CORRECTION
    !       TT_4P     REAL      TABLE FOR STOKES FREQUENCY CORRECTION
    !       IM_P      INTEGER   TABLE FOR WAVENUMBER M2 PLUS
    !       IM_M      INTEGER   TABLE FOR WAVENUMBER M2 MIN
    !       TFAK      REAL      WAVENUMBER TABLE
    !
    !
    !     METHOD
    !     ------
    !
    !     EXTERNALS
    !     ---------
    !             NONE
    !
    !     REFERENCES
    !     ----------
    !             V.E. ZAKHAROV, HAMILTONIAN APPROACH (1968)
    !             M.A. SROKOSZ, J.G.R.,91,995-1006 (1986)
    !             P.A.E.M. JANSSEN, ECMWF TECH MEMO (2008),JFM PAPER (2009)
    !
    !
    !--------------------------------------------------------------------
    !
    !
    !
    IMPLICIT NONE
 
    INTEGER NFRE,NANG,NDEPTH,MDW,JD,M,K,M1,K1,MP,MM,L
 
    REAL DEPTH,ALPHA,GAM_J,DEPTHA,DEPTHD
    REAL OM0,TH0,XK0,OM1,TH1,XK1,OM2,XK2,OM0P,XK0P,OM0M,XK0M,OMSTART,&
         FRAC,XMR,XM2,FAC
    REAL OMEGA(NFRE),TH(NANG),DFDTH(NFRE)
 
    common/const/depth,alpha,mdw,gam_j,depthd
    !
    !     1. COMPUTATION OF WAVENUMBER ARRAY TFAK
    !     ---------------------------------------
    !
    !
    DO jd=1,ndepth
      depth = deptha*depthd**(jd-1)
      DO m=1,nfre
        om0 = omega(m)
        tfak(m,jd) = aki(om0,depth)
      ENDDO
      WRITE(6,*) 'GENERATING TABLES FOR DEPTH:',jd,depth,deptha,ndepth
      !
      !     2. COMPUTATION OF THE 2nd ORDER COEFFICIENTS.
      !     ---------------------------------------------
      !
      !
      k1 = 0
      th1 = th(nang)
      DO m=1,nfre
        om0 = omega(m)
        xk0 = tfak(m,jd)
 
        mp   = min(m+1,nfre)
        om0p = omega(mp)
        xk0p = tfak(mp,jd)
 
        mm   = max(m-1,1)
        om0m = omega(mm)
        xk0m = tfak(mm,jd)
 
        DO m1=1,nfre
 
          om1 = omega(m1)
 
          DO l=1,nang
            !
            !              XK0-XK1 CASE
            !
            k = k1+l
            th0 = th(k)
            om2 = om0-om1
 
 
            IF (abs(om1).LT.om0/2.) THEN
              xm2  = log(om2/omstart)/log(1.+frac)
              im_m(m1,m) = nint(xmr*(xm2+1.))
              xk1 = tfak(m1,jd)
              xk2 = aki(om2,depth)
 
              ta(l,m1,m,jd) = dfdth(m1)*a(xk1,xk2,th1,th0)**2
            ELSE
              ta(l,m1,m,jd) = 0.
              im_m(m1,m) = 1
            ENDIF
            !
            !              XK1+XK0 CASE
            !
            om2 = om1+om0
            xm2  = log(om2/omstart)/log(1.+frac)
            im_p(m1,m) = nint(xmr*(xm2+1.))
            xk1 = tfak(m1,jd)
            xk2 = aki(om2,depth)
 
            tb(l,m1,m,jd) = dfdth(m1)*b(xk1,xk2,th1,th0)**2
            !
            !              QUASI-LINEAR EFFECT
            !
            !
            tc_ql(l,m1,m,jd) = dfdth(m1)*c_ql(xk0,xk1,th0,th1)
            !
            !              STOKES-FREQUENCY CORRECTION
            !
            !
            fac = 2.*g/om1*dfdth(m1)
            tt_4m(l,m1,m,jd) = &
                 fac*(w2(xk0m,xk1,xk1,xk0m,th0,th1,th1,th0)+&
                 v2(xk0m,xk1,xk1,xk0m,th0,th1,th1,th0))
            tt_4p(l,m1,m,jd) = &
                 fac*(w2(xk0p,xk1,xk1,xk0p,th0,th1,th1,th0)+&
                 v2(xk0p,xk1,xk1,xk0p,th0,th1,th1,th0))
            ! Table identical to Janssen: verified.
            !              IF (JD.EQ.1) WRITE(998,'(F4.1,3I3,5G11.3)') DEPTH,M,M1,L, TB(L,M1,M,JD),  &
            !                           TC_QL(L,M1,M,JD) , FAC, TT_4M(L,M1,M,JD), TT_4P(L,M1,M,JD)
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    !
    !
    !--------------------------------------------------------------------
    !
    RETURN
  END SUBROUTINE tables_2nd
  !
  !--------------------------------------------------------------------
  !
  SUBROUTINE secspom(F1,F3,NFRE,NANG,NMAX,NDEPTH,&
       DEPTHA,DEPTHD,OMSTART,FRAC,MR,DFDTH,OMEGA,&
       DEPTH,AKMEAN,TA,TB,TC_QL,TT_4M,TT_4P,&
       IM_P,IM_M,COUNTER)
    !
    !--------------------------------------------------------------------
    !
    !*****SECSPOM** COMPUTES SECOND ORDER SPECTRUM IN FREQUENCY SPACE.
    !
    !     P.JANSSEN JULY 2008
    !
    !     PURPOSE
    !     -------
    !             DETERMINES SECOND-ORDER SPECTRUM, BASED ON JANSSEN (2008)
    !             THERE ARE THREE CORRECTIONS:
    !                   1) GENERATION OF SECOND-HARMONICS
    !                   2) QUASI-LINEAR EFFECT
    !                   3) SHIFT OF SPECTRUM BECAUSE OF STOKES FREQUENCY
    !                      CORRECTION.
    !
    !     INTERFACE
    !     ---------
    !             *CALL* *SECSPOM(F1,F3,NFRE,NANG,NMAX,NDEPTH,
    !                             DEPTHA,DEPTHD,OMSTART,FRAC,MR,DFDTH,OMEGA,
    !                             DEPTH,AKMEAN,TA,TB,TC_QL,TT_4M,TT_4P,
    !                             IM_P,IM_M)*
    !
    !
    !     PARAMETER   TYPE      PURPOSE.
    !     ---------   ----      -------
    !
    !       F1        REAL      2D FREE WAVE SPECTRUM (INPUT)
    !       F3        REAL      BOUND WAVES SPECTRUM (OUTPUT)
    !       NFRE      INTEGER   NUMBER OF FREQUENCIES
    !       NANG      INTEGER   NUMBER OF DIRECTIONS
    !       NMAX      INTEGER   MAXIMUM INDEX CORRESPONDS TO TWICE THE CUT-OFF
    !                           FREQUENCY
    !       NDEPTH    INTEGER   NUMBER OF ENTRIES IN DEPTH TABLE
    !       DEPTHA    REAL      START VALUE DEPTH ARRAY
    !       DEPTHD    REAL      INCREMENT DEPTH ARRAY
    !       OMSTART   REAL      START VALUE ANG. FREQUENCY ARRAY
    !       FRAC      REAL      FRACTIONAL INCREASE IN FREQUENCY SPACE
    !       MR        INTEGER   THINNING FACTOR IN FREQUENCY SPACE
    !       OMEGA     REAL      ANGULAR FREQUENCY ARRAY
    !       DEPTH     REAL      DEPTH ARRAY
    !       AKMEAN    REAL      MEAN WAVENUMBER ARRAY
    !       TA        REAL      TABLE FOR MINUS INTERACTIONS
    !       TB        REAL      TABLE FOR PLUS INTERACTIONS
    !       TC_QL     REAL      TABLE FOR QUASI-LINEAR INTERACTIONS
    !       TT_4M     REAL      TABLE FOR STOKES FREQUENCY CORRECTION
    !       TT_4P     REAL      TABLE FOR STOKES FREQUENCY CORRECTION
    !       IM_P      INTEGER   TABLE FOR WAVENUMBER M2 PLUS
    !       IM_M      INTEGER   TABLE FOR WAVENUMBER M2 MIN
    !
    !
    !
    !     METHOD
    !     ------
    !             EVALUATE SECOND ORDER SPECTRUM IN FREQUENCY BASED ON
    !             KRASITSKII'S CANONICAL TRANSFORMATION.
    !
    !     EXTERNALS
    !     ---------
    !             NONE
    !
    !     REFERENCES
    !     ----------
    !             V.E. ZAKHAROV, HAMILTONIAN APPROACH (1968)
    !             M.A. SROKOSZ, J.G.R.,91,995-1006 (1986)
    !             P.A.E.M. JANSSEN, JFM (2009)
    !
    !
    !--------------------------------------------------------------------
    !
    !
    !
    USE w3gdatmd, ONLY: igpars
    IMPLICIT NONE
 
    INTEGER NFRE,NANG,NDEPTH,M,K,M1,K1,M2_M,M2_P,K2,MP,&
         MM,L,MR,NMAX,JD,COUNTER
    INTEGER IM_P(NFRE,NFRE),IM_M(NFRE,NFRE),IL(NANG,NANG)
 
    REAL OM0,OM0H,OM1,OM0P,OM0M,&
         OMSTART,FRAC,XINCR1,XINCR2,XINCR3,XINCR4,FAC1,FAC2,&
         FAC3,T_4M,T_4P,F2K,F2KP,F2KM,F2K1,F2K2,DELM1,DEPTHA,DEPTHD,&
         XD,X_MIN
    REAL OMEGA(NFRE), DFDTH(NFRE), OMEGAHF(NFRE+1:NMAX)
    REAL TA(NANG,NFRE,NFRE,NDEPTH),TB(NANG,NFRE,NFRE,NDEPTH),&
         TC_QL(NANG,NFRE,NFRE,NDEPTH),TT_4M(NANG,NFRE,NFRE,NDEPTH),&
         TT_4P(NANG,NFRE,NFRE,NDEPTH)
    REAL F1(NANG,NFRE),F3(NANG,NFRE),DEPTH
    REAL AKMEAN
    REAL G1(NANG,NMAX),G3(NANG,NFRE)
 
    LOGICAL :: LL2H
 
    !
    !***  1. COMPUTATION OF TAIL OF THE SPECTRUM AND INDEX JD
    !     ---------------------------------------------------
    !
    !
    x_min = igpars(9)   ! this was 1.1 in Janssen's original code
 
    DO m=nfre+1,nmax
      omegahf(m) = omstart*(1.+frac)**(mr*m-1)
    ENDDO
 
    DO k=1,nang
      DO k1=1,nang
        l = k-k1
        IF (l.GT.nang) l=l-nang
        IF (l.LT.1) l=l+nang
        il(k,k1) = l
      ENDDO
    ENDDO
 
 
    !   This was Janssen's version ... limited to kD > X_MIN ... (here set to 1.1)
    xd = max(x_min/akmean,depth)    ! note by FA: why do we have X_MIN/AKMEAN??!
    xd = depth
    xd = log(xd/deptha)/log(depthd)+1.
    jd = nint(xd)
    jd = max(jd,1)
    jd = min(jd,ndepth)
 
    DO m=1,nfre
      DO k=1,nang
        g1(k,m) = f1(k,m)
        g3(k,m) = 0.
      ENDDO
    ENDDO
 
    DO m=nfre+1,nmax
      DO k=1,nang
        g1(k,m) = omega(nfre)**5*g1(k,nfre)/omegahf(m)**5
      ENDDO
    ENDDO
    !
    !
    !
    !
    !***  2. COMPUTATION OF THE 2nd ORDER FREQUENCY SPECTRUM.
    !     ---------------------------------------------------
    !
    !
    DO m=1,nfre
      om0 = omega(m)
      om0h = om0/2.
      mp   = min(m+1,nfre)
      om0p = omega(mp)
      mm   = max(m-1,1)
      om0m = omega(mm)
      delm1 = 1./(om0p-om0m)
      DO k=1,nang
        k2 = k
        f2k = g1(k,m)
        f2kp = g1(k,mp)
        f2km = g1(k,mm)
        DO m1=1,nfre
          om1 = omega(m1)
          ll2h = (abs(om1).LT.om0h)
          m2_m = im_m(m1,m)
          m2_p = im_p(m1,m)
          DO k1=1,nang
            f2k1 = g1(k1,m1)
            l = il(k,k1)
            !
            !                   2.1 OM0-OM1 CASE: SECOND HARMONICS
            !                   OM2 = OM0-OM1
            !
            IF (ll2h) THEN
              f2k2 = g1(k2,m2_m)
              fac1 = ta(l,m1,m,jd)
              fac2 = f2k1*f2k2+g1(k2,m1)*g1(k1,m2_m)
 
              xincr1 = fac1*fac2
              g3(k,m) = g3(k,m)+xincr1
            ENDIF
            !
            !                   2.2 OM1+OM0 CASE: INFRA-GRAVITY WAVES
            !                    OM2 = OM1+OM0
            !
            f2k2 = g1(k2,m2_p)
            fac3 = 2.*tb(l,m1,m,jd)
            xincr2 = fac3*f2k2
            !
            !                   2.3 QUASI-LINEAR EFFECT
            !
            xincr3 = tc_ql(l,m1,m,jd)*f2k
            !
            !                   2.4 STOKES-FREQUENCY CORRECTION
            !
            t_4m = tt_4m(l,m1,m,jd)
            t_4p = tt_4p(l,m1,m,jd)
            xincr4 = -(f2kp*t_4p-f2km*t_4m)*delm1
 
            g3(k,m) = g3(k,m)+f2k1*(xincr2+xincr3+xincr4)
 
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    !
    DO m=1,nfre
      DO k=1,nang
        f3(k,m) = g3(k,m)
      ENDDO
    ENDDO
    !
    !--------------------------------------------------------------------
    !
    RETURN
  END SUBROUTINE secspom
 
  !
  REAL FUNCTION A(XI,XJ,THI,THJ)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *A(XI,XJ,THI,THJ)
    !
    !-----------------------------------------------------------------------
    !
    !***  *A*  DETERMINES THE MINUS INTERACTIONS.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE
    !              WAVE INTERACTIONS OF GRAVITY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV)
    !
    !     INTERFACE.
    !     ----------
    !              *A(XI,XJ)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL ri,rj,rk,xi,xj,thi,thj,thk,oi,oj,ok,fi,fj,fk
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
 
    ri = xi
    rj = xj
    rk  = vabs(ri,rj,thi,thj)
    thk = vdir(ri,rj,thi,thj)
 
    oi=omeg(ri)
    oj=omeg(rj)
    ok=omeg(rk)
 
    fi = sqrt(oi/(2.*g))
    fj = sqrt(oj/(2.*g))
    fk = sqrt(ok/(2.*g))
 
 
    a = fk/(fi*fj)*(a1(rk,ri,rj,thk,thi,thj)+&
         a3(rk,ri,rj,thk-pi,thi,thj))
 
    RETURN
  END FUNCTION a
  !
  REAL function b(xi,xj,thi,thj)
    !***  *REAL FUNCTION* *B(XI,XJ,THI,THJ)
    !
    !-----------------------------------------------------------------------
    !
    !***  *B*  DETERMINES THE PLUS INTERACTION COEFFICIENTS.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE
    !              WAVE INTERACTIONS OF GRAVITY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV)
    !
    !     INTERFACE.
    !     ----------
    !              *B(XI,XJ)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL del,ri,rj,rk,xi,xj,thi,thj,thk,oi,oj,ok,fi,fj,fk
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    del = 0.
    ri = xi
    rj = xj
    rk  = vabs(rj,ri,thj,thi-pi)
    thk = vdir(rj,ri,thj,thi-pi)
 
    oi=omeg(ri)+del
    oj=omeg(rj)+del
    ok=omeg(rk)+del
 
    fi = sqrt(oi/(2.*g))
    fj = sqrt(oj/(2.*g))
    fk = sqrt(ok/(2.*g))
 
    b = 0.5*fk/(fi*fj)*(a2(rk,ri,rj,thk,thi,thj)+&
         a2(rk,rj,ri,thk-pi,thj,thi))
 
    RETURN
  END FUNCTION b
  !
  REAL function c_ql(xk0,xk1,th0,th1)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *C_QL(XK0,XK1,TH0,TH1)
    !
    !-----------------------------------------------------------------------
    !
    !***  *A*  DETERMINES THE QUASI-LINEAR TERM.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              DETERMINE CONTRIBUTION BY QUASI-LINEAR TERMS
    !
    !     INTERFACE.
    !     ----------
    !              *C_QL(XK0,XK1)*
    !                        *XK0*  - WAVE NUMBER
    !                        *XK1*  - WAVE NUMBER
    !     METHOD.
    !     -------
 
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL xk0,xk1,th0,th1,om1,f1
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    om1 = omeg(xk1)
    f1  = sqrt(om1/(2.*g))
 
    c_ql = 2./f1**2*(b2(xk0,xk1,xk1,xk0,th0,th1,th1,th0)+&
         b3(xk0,xk0,xk1,xk1,th0-pi,th0,th1,th1))
 
    RETURN
  END FUNCTION c_ql
 
  !
  !
  REAL function vplus(xi,xj,xk,thi,thj,thk)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *VPLUS(XI,XJ,XK,THI,THJ,THK)
    !
    !-----------------------------------------------------------------------
    !
    !***  *VPLUS*  DETERMINES THE SECOND-ORDER TRANSFER COEFFICIENT
    !              FOR THREE WAVE INTERACTIONS OF GRAVITY WAVES.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE
    !              WAVE INTERACTIONS OF GRAVITY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV)
    !
    !     INTERFACE.
    !     ----------
    !              *VPLUS(XI,XJ,XK)*
    !                      *XI*   - WAVE NUMBER
    !                      *XJ*   - WAVE NUMBER
    !                      *XK*   - WAVE NUMBER
    !                      *THI*  - WAVE DIRECTION
    !                      *THJ*  - WAVE DIRECTION
    !                      *THK*  - WAVE DIRECTION
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL del1,ri,rj,rk,xi,xj,xk,thi,thj,thk,oi,oj,ok,qi,qj,qk,&
         rij,rik,rjk,sqijk,sqikj,sqjki,zconst
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    del1 = 10.**(-12)
    zconst=1./(4*sqrt(2.))
 
    ri = xi
    rj = xj
    rk = xk
 
    oi=omeg(ri)+del1
    oj=omeg(rj)+del1
    ok=omeg(rk)+del1
 
    qi=oi**2/g
    qj=oj**2/g
    qk=ok**2/g
 
    rij = ri*rj*cos(thj-thi)
    rik = ri*rk*cos(thk-thi)
    rjk = rj*rk*cos(thk-thj)
 
    sqijk=sqrt(g*ok/(oi*oj))
    sqikj=sqrt(g*oj/(oi*ok))
    sqjki=sqrt(g*oi/(oj*ok))
 
    vplus=zconst*( (rij+qi*qj)*sqijk + (rik+qi*qk)*sqikj&
         + (rjk+qj*qk)*sqjki )
    RETURN
  END FUNCTION vplus
  !
  REAL function vmin(xi,xj,xk,thi,thj,thk)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *VMIN(XI,XJ,XK,THI,THJ,THK)
    !
    !-----------------------------------------------------------------------
    !
    !***  *VMIN*  DETERMINES THE SECOND-ORDER TRANSFER COEFFICIENT FOR
    !             THREE WAVE INTERACTIONS OF GRAVITY WAVES.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE
    !              WAVE INTERACTIONS OF GRAVITY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV)
    !
    !     INTERFACE.
    !     ----------
    !              *VMIN(XI,XJ,XK)*
    !                      *XI*   - WAVE NUMBER
    !                      *XJ*   - WAVE NUMBER
    !                      *XK*   - WAVE NUMBER
    !                      *THI*  - WAVE DIRECTION
    !                      *THJ*  - WAVE DIRECTION
    !                      *THK*  - WAVE DIRECTION
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL del1,ri,rj,rk,xi,xj,xk,thi,thj,thk,oi,oj,ok,qi,qj,qk,&
         rij,rik,rjk,sqijk,sqikj,sqjki,zconst
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    del1 = 10.**(-12)
    zconst=1./(4*sqrt(2.))
 
    ri = xi
    rj = xj
    rk = xk
 
    oi=omeg(ri)+del1
    oj=omeg(rj)+del1
    ok=omeg(rk)+del1
 
    qi=oi**2/g
    qj=oj**2/g
    qk=ok**2/g
 
    rij = ri*rj*cos(thj-thi)
    rik = ri*rk*cos(thk-thi)
    rjk = rj*rk*cos(thk-thj)
 
    sqijk=sqrt(g*ok/(oi*oj))
    sqikj=sqrt(g*oj/(oi*ok))
    sqjki=sqrt(g*oi/(oj*ok))
 
    vmin=zconst*( (rij-qi*qj)*sqijk + (rik-qi*qk)*sqikj&
         + (rjk+qj*qk)*sqjki )
    RETURN
  END FUNCTION vmin
  !
  REAL function u(xi,xj,xk,xl,thi,thj,thk,thl)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *U(XI,XJ,XK,XL,THI,THJ,THK,THL)
    !
    !-----------------------------------------------------------------------
    !
    !***  *U*  DETERMINES THE THIRD-ORDER TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY WAVES.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV,AND CRAWFORD ET AL)
    !
    !     INTERFACE.
    !     ----------
    !              *U(XI,XJ,XK,XL)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !                      *XL*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL xi,xj,xk,xl,thi,thj,thk,thl,oi,oj,ok,ol,xik,xjk,xil,xjl,&
         oik,ojk,oil,ojl,qi,qj,qik,qjk,qil,qjl,sqijkl,zconst
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    zconst=1./(16.)
 
    oi=omeg(xi)
    oj=omeg(xj)
    ok=omeg(xk)
    ol=omeg(xl)
 
    xik = vabs(xi,xk,thi,thk)
    xjk = vabs(xj,xk,thj,thk)
    xil = vabs(xi,xl,thi,thl)
    xjl = vabs(xj,xl,thj,thl)
    oik=omeg(xik)
    ojk=omeg(xjk)
    oil=omeg(xil)
    ojl=omeg(xjl)
 
    qi=oi**2/g
    qj=oj**2/g
    qik=oik**2/g
    qjk=ojk**2/g
    qil=oil**2/g
    qjl=ojl**2/g
    sqijkl=sqrt(ok*ol/(oi*oj))
    u = zconst*sqijkl*( 2.*(xi**2*qj+xj**2*qi)-qi*qj*(&
         qik+qjk+qil+qjl) )
    RETURN
  END FUNCTION u
  !
  REAL function w2(xi,xj,xk,xl,thi,thj,thk,thl)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *W2(XI,XJ,XK,XL,THI,THJ,THK,THL)
    !
    !-----------------------------------------------------------------------
    !
    !***  *W2*  DETERMINES THE CONTRIBUTION OF THE DIRECT FOUR-WAVE
    !              INTERACTIONS OF GRAVITY WAVES OF THE TYPE
    !              A_2^*A_3A_4.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV,AND CRAWFORD ET AL)
    !
    !     INTERFACE.
    !     ----------
    !              *W(XI,XJ,XK,XL)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !                      *XL*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    REAL xi,xj,xk,xl,thi,thj,thk,thl
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    w2= u(xi,xj,xk,xl,thi-pi,thj-pi,thk,thl)+&
         u(xk,xl,xi,xj,thk,thl,thi-pi,thj-pi)-&
         u(xk,xj,xi,xl,thk,thj-pi,thi-pi,thl)-&
         u(xi,xk,xj,xl,thi-pi,thk,thj-pi,thl)-&
         u(xi,xl,xk,xj,thi-pi,thl,thk,thj-pi)-&
         u(xl,xj,xk,xi,thl,thj-pi,thk,thi-pi)
    RETURN
  END FUNCTION w2
  !
  REAL function v2(xi,xj,xk,xl,thi,thj,thk,thl)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *V2(XI,XJ,XK,XL,THI,THJ,THK,THL)
    !
    !-----------------------------------------------------------------------
    !
    !***  *V2*  DETERMINES THE CONTRIBUTION OF THE VIRTUAL
    !           FOUR-WAVE INTERACTIONS OF GRAVITY WAVES.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV,AND
    !              CRAWFORD ET AL)
    !
    !     INTERFACE.
    !     ----------
    !              *V2(XI,XJ,XK,XL)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !                      *XL*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    common/precis/doublep
    LOGICAL doublep
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL del1,xi,xj,xk,xl,thi,thj,thk,thl,oi,oj,ok,ol,ri,rj,rk,rl,&
         rij,rik,rli,rjl,rjk,rkl,thij,thik,thli,thjl,thjk,thkl,oij,&
         oik,ojl,ojk,oli,okl,xnik,xnjl,xnjk,xnil,ynil,ynjk,ynjl,ynik,&
         znij,znkl,zpij,zpkl,thlj,thil,thkj,thki,thji,thlk
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    IF (doublep) THEN
      del1=10.**(-5)
    ELSE
      del1=10.**(-2)
    ENDIF
 
 
    ri=xi+del1
    rj=xj+del1/2.
    rk=xk+del1/3.
    rl=xl+del1*(1.+1./2.-1./3.)
 
    oi=omeg(ri)
    oj=omeg(rj)
    ok=omeg(rk)
    ol=omeg(rl)
 
    rij  = vabs(ri,rj,thi,thj)
    thij = vdir(ri,rj,thi,thj)
 
    rik  = vabs(ri,rk,thi,thk-pi)
    thik = vdir(ri,rk,thi,thk-pi)
 
    rli  = vabs(rl,ri,thl,thi-pi)
    thli = vdir(xl,xi,thl,thi-pi)
 
    rjl  = vabs(rj,rl,thj,thl-pi)
    thjl = vdir(rj,rl,thj,thl-pi)
 
    rjk  = vabs(rj,rk,thj,thk-pi)
    thjk = vdir(rj,rk,thj,thk-pi)
 
    rkl  = vabs(rk,rl,thk,thl)
    thkl = vdir(rk,rl,thk,thl)
 
    oij=omeg(rij)
    oik=omeg(rik)
    ojl=omeg(rjl)
    ojk=omeg(rjk)
    oli=omeg(rli)
    okl=omeg(rkl)
 
    xnik = ok+oik-oi
    xnjl = oj+ojl-ol
    xnjk = ok+ojk-oj
    xnil = oi+oli-ol
 
    ynil = ol+oli-oi
    ynjk = oj+ojk-ok
    ynjl = ol+ojl-oj
    ynik = oi+oik-ok
 
    znij = oij-oi-oj
    znkl = okl-ok-ol
    zpij = oij+oi+oj
    zpkl = okl+ok+ol
 
    thlj = thjl-pi
    thil = thli-pi
    thkj = thjk-pi
    thki = thik-pi
    thji = thij-pi
    thlk = thkl-pi
 
    v2= vmin(ri,rk,rik,thi,thk,thik)*vmin(rl,rj,rjl,thl,thj,thlj)*&
         (1./xnik+1./xnjl)&
         +vmin(rj,rk,rjk,thj,thk,thjk)*vmin(rl,ri,rli,thl,thi,thli)*&
         (1./xnjk+1./xnil)&
         +vmin(ri,rl,rli,thi,thl,thil)*vmin(rk,rj,rjk,thk,thj,thkj)*&
         (1./ynil+1./ynjk)&
         +vmin(rj,rl,rjl,thj,thl,thjl)*vmin(rk,ri,rik,thk,thi,thki)*&
         (1./ynjl+1./ynik)&
         +vmin(rij,ri,rj,thij,thi,thj)*vmin(rkl,rk,rl,thkl,thk,thl)*&
         (1./znij+1./znkl)&
         +vplus(rij,ri,rj,thji,thi,thj)*vplus(rkl,rk,rl,thlk,thk,thl)*&
         (1./zpij+1./zpkl)
 
    v2 = -v2
 
    RETURN
  END FUNCTION v2
  !
  REAL function w1(xi,xj,xk,xl,thi,thj,thk,thl)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *W1(XI,XJ,XK,XL,THI,THJ,THK,THL)
    !
    !-----------------------------------------------------------------------
    !
    !***  *W1*  DETERMINES THE NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY WAVES OF THE TYPE
    !              A_2A_3A_4.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV,AND CRAWFORD ET AL)
    !
    !     INTERFACE.
    !     ----------
    !              *W1(XI,XJ,XK,XL)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !                      *XL*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL xi,xj,xk,xl,thi,thj,thk,thl
    !
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    w1= -u(xi,xj,xk,xl,thi-pi,thj,thk,thl)-&
         u(xi,xk,xj,xl,thi-pi,thk,thj,thl)-&
         u(xi,xl,xj,xk,thi-pi,thl,thj,thk)+&
         u(xj,xk,xi,xl,thj,thk,thi-pi,thl)+&
         u(xj,xl,xi,xk,thj,thl,thi-pi,thk)+&
         u(xk,xl,xi,xj,thk,thl,thi-pi,thj)
 
    w1=w1/3.
 
    RETURN
  END FUNCTION w1
  !
  REAL function w4(xi,xj,xk,xl,thi,thj,thk,thl)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *W4(XI,XJ,XK,XL,THI,THJ,THK,THL)
    !
    !-----------------------------------------------------------------------
    !
    !***  *W4*  DETERMINES THE NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY WAVES of the type
    !              A_^*A_3^*A_4^*.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV,AND CRAWFORD ET AL)
    !
    !     INTERFACE.
    !     ----------
    !              *W4(XI,XJ,XK,XL)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !                      *XL*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL xi,xj,xk,xl,thi,thj,thk,thl
    !
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
 
    w4= u(xi,xj,xk,xl,thi,thj,thk,thl)+&
         u(xi,xk,xj,xl,thi,thk,thj,thl)+&
         u(xi,xl,xj,xk,thi,thl,thj,thk)+&
         u(xj,xk,xi,xl,thj,thk,thi,thl)+&
         u(xj,xl,xi,xk,thj,thl,thi,thk)+&
         u(xk,xl,xi,xj,thk,thl,thi,thj)
 
 
    w4=w4/3.
 
    RETURN
  END FUNCTION w4
 
  REAL function b3(xi,xj,xk,xl,thi,thj,thk,thl)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *B3(XI,XJ,XK,XL,THI,THJ,THK,THL)
    !
    !-----------------------------------------------------------------------
    !
    !***  *B3*  WEIGHTS OF THE A_2^*A_3^*A_4 PART OF THE
    !           CANONICAL TRANSFORMATION.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV,AND CRAWFORD ET AL)
    !
    !     INTERFACE.
    !     ----------
    !              *B3(XI,XJ,XK,XL)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !                      *XL*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    common/precis/doublep
    LOGICAL doublep
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL del1,xi,xj,xk,xl,thi,thj,thk,thl,oi,oj,ok,ol,ri,rj,rk,rl,&
         rij,rji,rik,rki,rlj,rjl,rjk,rkj,rli,ril,rlk,rkl,thij,thji,&
         thik,thki,thlj,thjl,thjk,thkj,thli,thil,thlk,thkl,zijkl
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    IF (doublep) THEN
      del1=10.**(-5)
    ELSE
      del1=0.01
    ENDIF
 
    ri=xi
    rj=xj
    rk=xk
    rl=xl
 
    oi=omeg(ri)+del1
    oj=omeg(rj)+del1
    ok=omeg(rk)+del1
    ol=omeg(rl)+del1
 
    rij  = vabs(ri,rj,thi,thj)
    thij = vdir(ri,rj,thi,thj)
 
    rji  = vabs(rj,ri,thj,thi)
    thji = vdir(rj,ri,thj,thi)
 
    rik  = vabs(ri,rk,thi,thk)
    thik = vdir(ri,rk,thi,thk)
 
    rki  = vabs(rk,ri,thk,thi)
    thki = vdir(rk,ri,thk,thi)
 
    rlj  = vabs(rl,rj,thl,thj-pi)
    thlj = vdir(rl,rj,thl,thj-pi)
 
    rjl  = vabs(rj,rl,thj,thl-pi)
    thjl = vdir(rj,rl,thj,thl-pi)
 
    rjk  = vabs(rj,rk,thj,thk)
    thjk = vdir(rj,rk,thj,thk)
 
    rkj  = vabs(rk,rj,thk,thj)
    thkj = vdir(rk,rj,thk,thj)
 
    rli  = vabs(rl,ri,thl,thi-pi)
    thli = vdir(rl,ri,thl,thi-pi)
 
    ril  = vabs(ri,rl,thi,thl-pi)
    thil = vdir(ri,rl,thi,thl-pi)
 
    rlk  = vabs(rl,rk,thl,thk-pi)
    thlk = vdir(rl,rk,thl,thk-pi)
 
    rkl  = vabs(rk,rl,thk,thl-pi)
    thkl = vdir(rk,rl,thk,thl-pi)
 
    zijkl = oi+oj+ok-ol
 
    b3= -1./zijkl*(2.*( &
         vmin(rl,ri,rli,thl,thi,thli)*a1(rjk,rj,rk,thjk,thj,thk)&
         -vmin(rij,ri,rj,thij,thi,thj)*a1(rl,rk,rlk,thl,thk,thlk)&
         -vmin(rik,ri,rk,thik,thi,thk)*a1(rl,rj,rlj,thl,thj,thlj)&
         -vplus(rj,ri,rji,thj,thi,thji-pi)*a1(rk,rl,rkl,thk,thl,thkl)&
         -vplus(rk,ri,rki,thk,thi,thki-pi)*a1(rj,rl,rjl,thj,thl,thjl)&
         +vmin(ri,rl,ril,thi,thl,thil)*a3(rj,rk,rjk,thj,thk,thjk-pi))&
         +3.*w1(rl,rk,rj,ri,thl,thk,thj,thi) )
 
    RETURN
  END FUNCTION b3
  !
  REAL function b4(xi,xj,xk,xl,thi,thj,thk,thl)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *B4(XI,XJ,XK,XL,THI,THJ,THK,THL)
    !
    !-----------------------------------------------------------------------
    !
    !***  *B4*  WEIGHTS OF THE A_2^*A_3^*A_4^* PART OF THE CANONICAL
    !           TRANSFORMATION.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV,AND CRAWFORD ET AL)
    !
    !     INTERFACE.
    !     ----------
    !              *B4(XI,XJ,XK,XL)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !                      *XL*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    common/precis/doublep
    LOGICAL doublep
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL del1,xi,xj,xk,xl,thi,thj,thk,thl,oi,oj,ok,ol,ri,rj,rk,rl,&
         rij,rik,ril,rjl,rjk,rkl,thij,thik,thil,thjl,thjk,thlk,thkl,&
         zijkl
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
 
 
    ri=xi
    rj=xj
    rk=xk
    rl=xl
 
    oi=omeg(ri)
    oj=omeg(rj)
    ok=omeg(rk)
    ol=omeg(rl)
 
 
    rij  = vabs(ri,rj,thi,thj)
    thij = vdir(ri,rj,thi,thj)
 
    rik  = vabs(ri,rk,thi,thk)
    thik = vdir(ri,rk,thi,thk)
 
    ril  = vabs(ri,rl,thi,thl)
    thil = vdir(ri,rl,thi,thl)
 
    rjl  = vabs(rj,rl,thj,thl)
    thjl = vdir(rj,rl,thj,thl)
 
    rjk  = vabs(rj,rk,thj,thk)
    thjk = vdir(rj,rk,thj,thk)
 
    rkl  = vabs(rk,rl,thk,thl)
    thkl = vdir(rk,rl,thk,thl)
 
 
    zijkl = oi+oj+ok+ol
 
    b4= -1./zijkl*(2./3.*(    &
         vplus(rij,ri,rj,thij-pi,thi,thj)*a1(rkl,rk,rl,thkl,thk,thl)&
         +vplus(rik,ri,rk,thik-pi,thi,thk)*a1(rjl,rj,rl,thjl,thj,thl)&
         +vplus(ril,ri,rl,thil-pi,thi,thl)*a1(rjk,rj,rk,thjk,thj,thk)&
         +vmin(rik,ri,rk,thik,thi,thk)*a3(rjl,rj,rl,thjl-pi,thj,thl)&
         +vmin(ril,ri,rl,thil,thi,thl)*a3(rjk,rj,rk,thjk-pi,thj,thk)&
         +vmin(rij,ri,rj,thij,thi,thj)*a3(rkl,rk,rl,thkl-pi,thk,thl) )&
         +w4(ri,rj,rk,rl,thi,thj,thk,thl) )
 
    RETURN
  END FUNCTION b4
  !
  REAL function b1(xi,xj,xk,xl,thi,thj,thk,thl)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *B1(XI,XJ,XK,XL,THI,THJ,THK,THL)
    !
    !-----------------------------------------------------------------------
    !
    !***  *B1*  WEIGHTS OF THE A_2A_3A_4 PART OF THE CANONICAL
    !           TRANSFORMATION.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV,AND CRAWFORD ET AL)
    !
    !     INTERFACE.
    !     ----------
    !              *B1(XI,XJ,XK,XL)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !                      *XL*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    common/precis/doublep
    LOGICAL doublep
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL del1,xi,xj,xk,xl,thi,thj,thk,thl,oi,oj,ok,ol,ri,rj,rk,rl,&
         rij,rji,rik,rki,rjl,rjk,rli,ril,rkl,thij,thji,&
         thik,thki,thjl,thjk,thli,thil,thkl,zijkl
    !
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
 
    ri=xi
    rj=xj
    rk=xk
    rl=xl
 
    oi=omeg(ri)
    oj=omeg(rj)
    ok=omeg(rk)
    ol=omeg(rl)
 
    rij  = vabs(ri,rj,thi,thj-pi)
    thij = vdir(ri,rj,thi,thj-pi)
 
    rji  = vabs(rj,ri,thj,thi-pi)
    thji = vdir(rj,ri,thj,thi-pi)
 
    rik  = vabs(ri,rk,thi,thk-pi)
    thik = vdir(ri,rk,thi,thk-pi)
 
    rki  = vabs(rk,ri,thk,thi-pi)
    thki = vdir(rk,ri,thk,thi-pi)
 
    ril  = vabs(ri,rl,thi,thl-pi)
    thil = vdir(ri,rl,thi,thl-pi)
 
    rli  = vabs(rl,ri,thl,thi-pi)
    thli = vdir(rl,ri,thl,thi-pi)
 
    rjl  = vabs(rj,rl,thj,thl)
    thjl = vdir(rj,rl,thj,thl)
 
    rjk  = vabs(rj,rk,thj,thk)
    thjk = vdir(rj,rk,thj,thk)
 
    rkl  = vabs(rk,rl,thk,thl)
    thkl = vdir(rk,rl,thk,thl)
 
    zijkl = oi-oj-ok-ol
 
    b1= -1./zijkl*(2./3.*(    &
         min(ri,rj,rij,thi,thj,thij)*a1(rkl,rk,rl,thkl,thk,thl)&
         +vmin(ri,rk,rik,thi,thk,thik)*a1(rjl,rj,rl,thjl,thj,thl)&
         +vmin(ri,rl,ril,thi,thl,thil)*a1(rjk,rj,rk,thjk,thj,thk)&
         +vmin(rk,ri,rki,thk,thi,thki)*a3(rjl,rj,rl,thjl-pi,thj,thl)&
         +vmin(rl,ri,rli,thl,thi,thli)*a3(rjk,rj,rk,thjk-pi,thj,thk)&
         +vmin(rj,ri,rji,thj,thi,thji)*a3(rkl,rk,rl,thkl-pi,thk,thl) &
         ) +w1(ri,rj,rk,rl,thi,thj,thk,thl) )
    RETURN
  END FUNCTION b1
  !
  REAL function b2(xi,xj,xk,xl,thi,thj,thk,thl)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *B2(XI,XJ,XK,XL,THI,THJ,THK,THL)
    !
    !-----------------------------------------------------------------------
    !
    !***  *B2*  WEIGHTS OF THE A_2^*A_3A_4 PART OF THE CANONICAL
    !           TRANSFORMATION.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR FOUR
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV,AND CRAWFORD ET AL)
    !
    !     INTERFACE.
    !     ----------
    !              *B2(XI,XJ,XK,XL)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !                      *XL*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    common/precis/doublep
    LOGICAL doublep
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL del1,xi,xj,xk,xl,thi,thj,thk,thl,oi,oj,ok,ol,ri,rj,rk,rl,&
         rij,rik,rki,rjl,rlj,rjk,rkj,rli,ril,rkl,thij,&
         thik,thki,thjl,thlj,thjk,thkj,thli,thil,thkl,zijkl
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
 
    ri=xi
    rj=xj
    rk=xk
    rl=xl
 
    rij  = vabs(ri,rj,thi,thj)
    thij = vdir(ri,rj,thi,thj)
 
    rik  = vabs(ri,rk,thi,thk-pi)
    thik = vdir(ri,rk,thi,thk-pi)
 
    rki  = vabs(rk,ri,thk,thi-pi)
    thki = vdir(rk,ri,thk,thi-pi)
 
    ril  = vabs(ri,rl,thi,thl-pi)
    thil = vdir(ri,rl,thi,thl-pi)
 
    rli  = vabs(rl,ri,thl,thi-pi)
    thli = vdir(rl,ri,thl,thi-pi)
 
    rjl  = vabs(rj,rl,thj,thl-pi)
    thjl = vdir(rj,rl,thj,thl-pi)
 
    rlj  = vabs(rl,rj,thl,thj-pi)
    thlj = vdir(rl,rj,thl,thj-pi)
 
    rjk  = vabs(rj,rk,thj,thk-pi)
    thjk = vdir(rj,rk,thj,thk-pi)
 
    rkj  = vabs(rk,rj,thk,thj-pi)
    thkj = vdir(rk,rj,thk,thj-pi)
 
    rkl  = vabs(rk,rl,thk,thl)
    thkl = vdir(rk,rl,thk,thl)
 
    b2=  a3(ri,rj,rij,thi,thj,thij-pi)*a3(rk,rl,rkl,thk,thl,thkl-pi)&
         +a1(rj,rk,rjk,thj,thk,thjk)*a1(rl,ri,rli,thl,thi,thli)&
         +a1(rj,rl,rjl,thj,thl,thjl)*a1(rk,ri,rki,thk,thi,thki)&
         -a1(rij,ri,rj,thij,thi,thj)*a1(rkl,rk,rl,thkl,thk,thl)&
         -a1(ri,rk,rik,thi,thk,thik)*a1(rl,rj,rlj,thl,thj,thlj)&
         -a1(ri,rl,ril,thi,thl,thil)*a1(rk,rj,rkj,thk,thj,thkj)
 
 
    RETURN
  END FUNCTION b2
  !
  REAL function a1(xi,xj,xk,thi,thj,thk)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *A1(XI,XJ,XK,THI,THJ,THK)
    !
    !-----------------------------------------------------------------------
    !
    !***  *A1*  AUXILIARY SECOND-ORDER COEFFICIENT.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV)
    !
    !     INTERFACE.
    !     ----------
    !              *VMIN(XI,XJ,XK)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    common/precis/doublep
    LOGICAL doublep
    INTEGER mdw
    REAL depth,alpha,gam_j,deptha,depthd
    REAL del1,xi,xj,xk,thi,thj,thk,oi,oj,ok
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    IF (doublep) THEN
      del1 = 10.**(-8)
    ELSE
      del1 = 10.**(-4)
    ENDIF
 
    oi=omeg(xi)+del1
    oj=omeg(xj)+del1
    ok=omeg(xk)+del1
 
    a1 = -vmin(xi,xj,xk,thi,thj,thk)/(oi-oj-ok)
 
    RETURN
  END FUNCTION a1
  !
  REAL function a2(xi,xj,xk,thi,thj,thk)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *A2(XI,XJ,XK,THI,THJ,THK)
    !
    !-----------------------------------------------------------------------
    !
    !***  *A2*  AUXILIARY SECOND-ORDER FUNCTION.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV)
    !
    !     INTERFACE.
    !     ----------
    !              *VMIN(XI,XJ,XK)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    REAL del1,xi,xj,xk,thi,thj,thk
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    a2 = -2.*a1(xk,xj,xi,thk,thj,thi)
    RETURN
  END FUNCTION a2
  !
  REAL function a3(xi,xj,xk,thi,thj,thk)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *A3(XI,XJ,XK,THI,THJ,THK)
    !
    !-----------------------------------------------------------------------
    !
    !***  *A3*  AUXILIARY SECOND-ORDER FUNCTION.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE
    !              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE
    !              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV)
    !
    !     INTERFACE.
    !     ----------
    !              *VMIN(XI,XJ,XK)*
    !                      *XI*  - WAVE NUMBER
    !                      *XJ*  - WAVE NUMBER
    !                      *XK*  - WAVE NUMBER
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/precis/doublep
    LOGICAL doublep
    REAL del1,oi,oj,ok,xi,xj,xk,thi,thj,thk
    !
    !***  1. DETERMINE NONLINEAR TRANSFER.
    !     --------------------------------
    !
    IF (doublep) THEN
      del1 = 10.**(-8)
    ELSE
      del1 = 10.**(-4)
    ENDIF
 
 
    oi=omeg(xi)+del1
    oj=omeg(xj)+del1
    ok=omeg(xk)+del1
 
    a3 = -vplus(xi,xj,xk,thi,thj,thk)/(oi+oj+ok)
    RETURN
  END FUNCTION a3
 
  !
  REAL function omeg(x)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *OMEG(X)*
    !
    !-----------------------------------------------------------------------
    !
    !
    !***  *OMEG*   DETERMINES THE DISPERSION RELATION FOR GRAVITY
    !              WAVES.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES DISPERSION RELATION FOR GRAVITY-
    !              WAVES IN THE IDEAL CASE OF NO CURRENT.
    !
    !     INTERFACE.
    !     ----------
    !              *OMEG(X)*
    !                      *X*  - WAVE NUMBER
    !
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,depthd
    REAL d,xk,x,t
 
    d = depth
    xk = abs(x)
    t = tanh(xk*d)
    omeg=sqrt(g*xk*t)
 
    RETURN
  END FUNCTION omeg
  !
  REAL function vg(x)
    !-----------------------------------------------------------------------
    !
    !***  *REAL FUNCTION* *VG(X)*
    !
    !-----------------------------------------------------------------------
    !
    !***  *VG*   DETERMINES THE GROUP VELOCITY FOR GRAVITY- WAVES.
    !
    !     PETER JANSSEN
    !
    !     PURPOSE.
    !     --------
    !
    !              GIVES GROUP VELOCITY FOR GRAVITY-
    !              WAVES IN THE IDEAL CASE OF NO CURRENT.
    !
    !     INTERFACE.
    !     ----------
    !              *VG(X)*
    !              *X*  - WAVE NUMBER
    !
    !     METHOD.
    !     -------
    !              NONE
    !
    !     EXTERNALS.
    !     ----------
    !              NONE.
    !
    !-----------------------------------------------------------------------
    !
    IMPLICIT NONE
    common/const/depth,alpha,mdw,gam_j,depthd
    INTEGER mdw
    REAL depth,alpha,gam_j,depthd
    REAL d,xk,x,xd
 
    d = depth
    xk = abs(x)
    xd = xk*depth
 
    vg = 0.5*sqrt(g*tanh(xd)/xk)*(1.+2.*xd/sinh(2.*xd))
 
    RETURN
  END FUNCTION vg
  !---------------------------------------------------------------------
  REAL function aki(om,beta)
    ! This function gives the wavenumber ...
    !---------------------------------------------------------------------
    !
    IMPLICIT NONE
    REAL om,beta,g,ebs,akm1,akm2,ao,akp,bo,th,sth
 
    g =9.806
    ebs=0.0001
    akm1=om**2/(4.*g )
    akm2=om/(2.*sqrt(g*beta))
    ao=max(akm1,akm2)
10  CONTINUE
    akp=ao
    bo=beta*ao
    !     IF (BO.GT.10) GO TO 20
    IF (bo.GT.20.) GO TO 20
    th=g*ao*tanh(bo)
    sth=sqrt(th)
    ao=ao+(om-sth)*sth*2./(th/ao+g*bo/cosh(bo)**2)
    IF (abs(akp-ao).GT.ebs*ao) GO TO 10
    aki=ao
    RETURN
20  CONTINUE
    aki=om**2/g
    RETURN
  END FUNCTION aki
  !
  REAL function vabs(xi,xj,thi,thj)
    !
    !---------------------------------------------------------------------
    !
    IMPLICIT NONE
    REAL xi,xj,thi,thj,arg
 
    arg = xi**2+xj**2+2.*xi*xj*cos(thi-thj)
 
    IF (arg.LE.0.) THEN
      vabs = 0.
    ELSE
      vabs = sqrt(arg)
    ENDIF
 
    RETURN
  END FUNCTION vabs
  !
  REAL function vdir(xi,xj,thi,thj)
    !
    !---------------------------------------------------------------------
    !
    IMPLICIT NONE
    REAL xi,xj,thi,thj,eps,y,x
 
    eps = 0.
 
    y = xj*sin(thj-thi)
    x = xi+xj*cos(thj-thi)+eps
    vdir = atan2(y,x)+thi
    IF (x.EQ.0.) vdir = 0.
 
    RETURN
  END FUNCTION vdir
  !/
  !/ End of module W3CANOMD -------------------------------------------- /
  !/
END MODULE w3canomd