Generic shallow-water Boltzmann integral (FBI or TSA). More...

Functions/Subroutines
subroutine	insnl4
	It reads look-up tables (generated by gridsetr) if file exists otherwise it must generate the look-up tables file. More...

subroutine	w3snl4 (A, CG, WN, DEPTH, S, D)
	Interface module for TSA type nonlinear interactions. More...

subroutine	gridsetr (dep, wka1, cgnrng1)
	This routine sets up the geometric part of the Boltzmann integral. More...

subroutine	shloxr (dep, wk1x, wk1y, wk3x, wk3y)
	General locus solution for input vectors (wk1x,wk1y). More...

subroutine	shlocr (dep, wk1x, wk1y, wk3x, wk3y)
	N/A. More...

subroutine	cplshr (w1x0, w1y0, w2x0, w2y0, w3x0, w3y0, h, csq, irng, krng, kang, ipt)
	Calculates four-wave Boltzmann coupling coefficient in shallow water. More...

subroutine	optsa2 (nrmn, nrmx, npk, fpk, nbins, wka, cga)
	Splits the Action Density into two parts. More...

subroutine	snlr_fbi (pha, ialt)
	For a given Action Density array dens1(k,theta), computes the rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to wave-wave interaction. More...

subroutine	snlr_tsa (pha, ialt)
	For a given Action Density array dens1(k,theta), computes the rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to wave-wave interaction. More...

subroutine	interp2 (X)
	Interpolate bi-linearly to fill in tsa/fbi & diag/diag2 arrays. More...

real function	wkfnc (f, dep)
	Calculate the Wavenumber k (rad/m) as function of frequency 'f' (Hz) and depth 'dep' (m). More...

real function	cgfnc (f, dep, cvel)
	Calculate the Group velocity Cg (m/s) as function of frequency 'f' (Hz), depth 'dep' (m) and phase speed 'cvel' (m/s). More...

Variables
integer, parameter	ismo = 1

integer, parameter	npts = 30

integer, parameter	ndep = 37

integer	nrng

integer	nzz

integer	kzone

integer	nb2fp

integer	nang

integer	na2p1

integer	np2p1

real	dfrq

real	f0

real	ainc

real	twopi

real, dimension(:), allocatable	frqa

real, dimension(:), allocatable	oma

real, dimension(:), allocatable	angl

real, dimension(:), allocatable	sinan

real, dimension(:), allocatable	cosan

real, dimension(:), allocatable	dep_tbl

integer, dimension(:,:,:,:), allocatable	kref2_tbl

integer, dimension(:,:,:,:), allocatable	kref4_tbl

integer, dimension(:,:,:,:), allocatable	jref2_tbl

integer, dimension(:,:,:,:), allocatable	jref4_tbl

real, dimension(:,:,:,:), allocatable	wtk2_tbl

real, dimension(:,:,:,:), allocatable	wtk4_tbl

real, dimension(:,:,:,:), allocatable	wta2_tbl

real, dimension(:,:,:,:), allocatable	wta4_tbl

real, dimension(:,:,:,:), allocatable	tfac2_tbl

real, dimension(:,:,:,:), allocatable	tfac4_tbl

real, dimension(:,:,:,:), allocatable	grad_tbl

real, dimension(:,:), allocatable	pha_tbl

integer, dimension(:,:,:), allocatable	kref2

integer, dimension(:,:,:), allocatable	kref4

integer, dimension(:,:,:), allocatable	jref2

integer, dimension(:,:,:), allocatable	jref4

real, dimension(:,:,:), allocatable	wtk2

real, dimension(:,:,:), allocatable	wtk4

real, dimension(:,:,:), allocatable	wta2

real, dimension(:,:,:), allocatable	wta4

real, dimension(:,:,:), allocatable	tfac2

real, dimension(:,:,:), allocatable	tfac4

real, dimension(:,:,:), allocatable	grad

real, dimension(:), allocatable	wk2x

real, dimension(:), allocatable	wk2y

real, dimension(:), allocatable	wk4x

real, dimension(:), allocatable	wk4y

real, dimension(:), allocatable	ds

real, dimension(:,:), allocatable	ef2

real, dimension(:), allocatable	ef1

real, dimension(:,:), allocatable	dens1

real, dimension(:,:), allocatable	dens2

real, dimension(:,:), allocatable	tsa

real, dimension(:,:), allocatable	diag

real, dimension(:,:), allocatable	fbi

real, dimension(:,:), allocatable	diag2

Detailed Description

Generic shallow-water Boltzmann integral (FBI or TSA).

Interface module for TSA type nonlinear interactions. Based on Resio and Perrie (2008) and Perrie and Resio (2009).

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Function/Subroutine Documentation

◆ cgfnc()

real function w3snl4md::cgfnc	(	real, intent(in)	f,
		real, intent(in)	dep,
		real, intent(in)	cvel
	)

Calculate the Group velocity Cg (m/s) as function of frequency 'f' (Hz), depth 'dep' (m) and phase speed 'cvel' (m/s).

  This routine uses the identity
         sinh(2x) = 2*tanh(x)/(1-tanh(x)**2)
    to avoid extreme sinh(2x) for large x.
    thus,    2kd/sinh(2kd) = kd(1-tanh(kd)**2)/tanh(kd)
    Group velocity Cg (m/s) is returned in "cgfnc".

Parameters

f	Frequency (Hz).
dep	Depth (m).
cvel	Phase speed (m/s).

Returns: cgfnc Group velocity (m/s).

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 5409 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!    it returns:  group velocity (m/s) in cgfnc
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !! 1. Purpose :
     !!
     !!    Calculate the Group velocity Cg (m/s) as function of
     !!    frequency 'f' (Hz), depth 'dep' (m) and phase speed 'cvel' (m/s)
     !!
     !! 2. Method :
     !!
     !!    This routine uses the identity
     !!         sinh(2x) = 2*tanh(x)/(1-tanh(x)**2)
     !!    to avoid extreme sinh(2x) for large x.
     !!    thus,    2kd/sinh(2kd) = kd(1-tanh(kd)**2)/tanh(kd)
     !!    Group velocity Cg (m/s) is returned in "cgfnc"
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    twopi     Real  Public    I   = TPI; WW3 2*pi=8.*atan(1.) (radians)
     !!    ------------------------------------------------------------------
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     IMPLICIT NONE
     !!
     !!    Parameter list
     !!    --------------
     real,    intent(in)  :: f, dep, cvel
     !!
     !!    Local variables
     !!    ---------------
     real                 :: wkd, tkd
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     wkd   = twopi * f*dep/cvel           !* kd
     tkd   = tanh(wkd)                    !* tanh(kd)
     cgfnc = 0.5*cvel*(1.+wkd*(1.-tkd**2)/tkd)
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     RETURN
     !!

References twopi.

Referenced by insnl4().

◆ cplshr()

subroutine w3snl4md::cplshr	(	real, intent(in)	w1x0,
		real, intent(in)	w1y0,
		real, intent(in)	w2x0,
		real, intent(in)	w2y0,
		real, intent(in)	w3x0,
		real, intent(in)	w3y0,
		real, intent(in)	h,
		real, intent(out)	csq,
		integer, intent(in)	irng,
		integer, intent(in)	krng,
		integer, intent(in)	kang,
		integer, intent(in)	ipt
	)

Calculates four-wave Boltzmann coupling coefficient in shallow water.

Calculates four-wave Boltzmann coupling coefficient in shallow water given k1,k2,k3 and following at least Hasselmann (1962) and probably Herterich and Hasselmann (1982). Dimensional wavenumbers are (wnx0,wny0), n = 1,3, h = depth, csq = coupling coefficient. This is the same as Don's cplesh, except within the algorithm, wavenumbers are made dimensionless with h and frequencies with sqrt(h/g), g = gravitational acceleration (the idea is to simplify and speed up the calculations while keeping a reasonable machine resolution of the result). At the end, dimensionless csqhat is redimensioned as csq = csqhat/(h**6) so it is returned as a dimensional entity.

This calculation can be a touchy bird, so we use double precision for internal calculations, using single precision for input and output.

Parameters

[in]	w1x0
[in]	w1y0
[in]	w2x0
[in]	w2y0
[in]	w3x0
[in]	w3y0
[in]	h
[out]	csq
[in]	irng
[in]	krng
[in]	kang
[in]	ipt

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 2876 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!
     !!    it returns: the coupling coefficient  csq
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !! 1. Purpose :
     !!
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    Calculates four-wave Boltzmann coupling coefficient in shallow   !
     !!    water given k1,k2,k3 and following at least Hasselmann (1962)    !
     !!    and probably Herterich and Hasselmann (1982).  Dimensional       !
     !!    wavenumbers are (wnx0,wny0), n = 1,3, h = depth, csq = coupling  !
     !!    coefficient.  This is the same as Don's cplesh, except within    !
     !!    the algorithm, wavenumbers are made dimensionless with h and     !
     !!    frequencies with sqrt(h/g), g = gravitational acceleration (the  !
     !!    idea is to simplify and speed up the calculations while keeping  !
     !!    a reasonable machine resolution of the result).  At the end,     !
     !!    dimensionless csqhat is redimensioned as csq = csqhat/(h**6)     !
     !!    so it is returned as a dimensional entity.                       !
     !!                                                                     !
     !!    This calculation can be a touchy bird, so we use double precision!
     !!    for internal calculations, using single precision for input and  !
     !!    output.                                                          !
     !!                                                                     !
     !!    -----------------------------------------------------------------#
     !!
     !! 2. Method :
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    ------------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!    ------------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!     gridsetr  Subr. W3SNL4MD Setup geometric integration grid
     !!    ------------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!
     !!    !/S  Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     IMPLICIT NONE
     !!
     !!    Parameter list
     !!    --------------
     integer, intent(in)  :: irng,krng, kang,ipt
     real,    intent(in)  :: w1x0,w1y0, w2x0,w2y0, w3x0,w3y0
     real,    intent(in)  :: h          !* depth 'dep'
     real,    intent(out) :: csq
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    Local Parameters & variables
     !!    -----------------------------
     integer              :: ipass
     double precision     :: hh
     double precision     :: s1,    s2,    s3
     double precision     :: k1x,   k2x,   k3x
     double precision     :: k1y,   k2y,   k3y
     double precision     :: k1,    k2,    k3
     double precision     :: om1,   om2,   om3
     double precision     :: som1,  som2,  som3
     double precision     :: om1sq, om2sq, om3sq
     double precision     :: k23,   k23x,  k23y
     double precision     :: dot23, dot123
     double precision     :: omsq23
     !!mpc
     double precision     :: k1sq, k2sq, k3sq, k23sq
     double precision     :: tanh_k1, tanh_k2, tanh_k3, tanh_k23
     !!mpc---
     double precision     :: k1x0,  k2x0,  k3x0,  k1zx
     double precision     :: k1y0,  k2y0,  k3y0,  k1zy
     double precision     :: di,    e
     double precision     :: p1,    p2,    p3,    p4
     double precision     :: t1,    t2,    t3,    t4,    t5
     !!
     double precision     :: csqhatd, csqd
     double precision     :: scple
     double precision     :: pi4
     !!
     !!eps Bash; added +eps to avoid dividing by 0.0. Dividing by 0.0 causes NaN
     !eps  double precision     :: eps
     !!
     !!    Bash; Added domsq23 = denominator of  t1      in cplshr
     !!          and   sumom   = denominator of  csqhatd in cplshr
     double precision     :: domsq23, sumom
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!    initial constants
     !!    ------------------
     hh    = dble(h)                     !* single to dbl precision
     pi4   = 0.785398175d0               !* Set = PI/4 as in CONSTANTS
     !eps  eps   = 1.d-12                      !* set eps to a very small number
     scple = 0.d0                        !* initialize accumulator
     !!ini
     !!    initialize returned variable 'csq'
     !!    ----------------------------------
     csq   = 0.d0
     !!ini---
     !!    ------------------------------------------------------------------
     !!
     do ipass=1,3
       !p1
       if (ipass .eq. 1) then            !* initial pass (+1,+1,-1)
         s1   =  1.d0
         s2   =  1.d0
         s3   = -1.d0
         k1x0 = dble(w1x0) * hh          !* norm. k elements with h
         k1y0 = dble(w1y0) * hh
         k2x0 = dble(w2x0) * hh
         k2y0 = dble(w2y0) * hh
         k3x0 = dble(w3x0) * hh
         k3y0 = dble(w3y0) * hh
         !p1
         !p2
       else if (ipass .eq. 2) then       !* 1st permutation (+1,-1,+1)
         s1   =  1.d0
         s2   = -1.d0
         s3   =  1.d0
         k1zx = k1x0
         k1zy = k1y0
         k1x0 = k2x0
         k1y0 = k2y0
         k2x0 = k3x0
         k2y0 = k3y0
         k3x0 = k1zx
         k3y0 = k1zy
         !p2
         !p3
       else                              !* 2nd permutation (-1,+1,+1)
         s1   = -1.d0
         s2   =  1.d0
         s3   =  1.d0
         k1zx = k1x0
         k1zy = k1y0
         k1x0 = k2x0
         k1y0 = k2y0
         k2x0 = k3x0
         k2y0 = k3y0
         k3x0 = k1zx
         k3y0 = k1zy
         !p3
       end if
       !!k19p1
       !!k19p1 Note: na2p1=nang/2+1   !* this is the angle index opposite to iang=1
       !k19p1  if (krng.ne.irng .and. kang.eq.na2p1 .and. ipt.eq.1 .and.     &
       !k19p1                                        ipass.eq.1) go to 10
       !!k19p1---
       !!
       k1x = s1 * k1x0                    !* sign the norm'ed k parts
       k1y = s1 * k1y0
       k2x = s2 * k2x0
       k2y = s2 * k2y0
       k3x = s3 * k3x0
       k3y = s3 * k3y0
       !!mpc
       !mpc    k1    = dsqrt(k1x**2 + k1y**2)     !* normalized |k|
       !mpc    k2    = dsqrt(k2x**2 + k2y**2)
       !mpc    k3    = dsqrt(k3x**2 + k3y**2)
       !!mpc---
       k1sq  = (k1x*k1x + k1y*k1y)        !* normalized |k| **2
       k2sq  = (k2x*k2x + k2y*k2y)
       k3sq  = (k3x*k3x + k3y*k3y)
       k1    = dsqrt(k1sq)                !* normalized |k|
       k2    = dsqrt(k2sq)
       k3    = dsqrt(k3sq)
       !!mpc---
       !!
       !!mpc
       !mpc    om1   = dsqrt(k1*dtanh(k1))        !* norm. omega (by sqrt(h/g))
       !mpc    om2   = dsqrt(k2*dtanh(k2))
       !mpc    om3   = dsqrt(k3*dtanh(k3))
       !mpc    om1sq = om1**2
       !mpc    om2sq = om2**2
       !mpc    om3sq = om3**2
       !!mpc---
       tanh_k1 = dtanh(k1)
       tanh_k2 = dtanh(k2)
       tanh_k3 = dtanh(k3)
       om1sq   = k1*tanh_k1
       om2sq   = k2*tanh_k2
       om3sq   = k3*tanh_k3
       om1     = dsqrt(om1sq)             !* norm. omega (by sqrt(h/g))
       om2     = dsqrt(om2sq)
       om3     = dsqrt(om3sq)
       !!mpc---
       !!
       som1 = s1 * om1                    !* sign the norm'ed omega's
       som2 = s2 * om2
       som3 = s3 * om3
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       dot23  = k2x*k3x + k2y*k3y        !*  vector k2 dot vector k3
       !!
       k23x   = k2x + k3x                !* (vector k2  +  vector k3)_x
       k23y   = k2y + k3y                !* (vector k2  +  vector k3)_y
       !!
       !!mpc
       !mpc    k23    = dsqrt(k23x**2+k23y**2)   !* |vector k2  +  vector k3|
       !!mpc---
       k23sq  = (k23x*k23x + k23y*k23y)
       k23    = dsqrt(k23sq)             !* |vector k2  +  vector k3|
       !!mpc---
       !!
       !!mpc
       !mpc    omsq23 = k23 * dtanh(k23)         !* norm sq frq of v.k2+v.k3
       !!mpc---
       tanh_k23 = dtanh(k23)
       omsq23   = k23 * tanh_k23         !* norm sq frq of v.k2+v.k3
       !!mpc---
       !!
       dot123 = k1x*k23x + k1y*k23y      !* v.k1 dot (v.k2 + v.k3)
       !!      ----------------------------------------------------------------
       !!
       !!      note: the "i**2" factor from some reference is included in this term
       !!
       !!mpc
       !mpc    di = -(som2+som3)*(om2sq*om3sq-dot23)+0.5d0 *                 &
       !mpc          (som2*(k3**2-om3sq**2)+som3*(k2**2-om2sq**2))
       !!mpc---
       di = -(som2+som3)*(om2sq*om3sq-dot23)+0.5d0 *                 &
            (som2*(k3sq-om3sq*om3sq)+som3*(k2sq-om2sq*om2sq))
       !!mpc---
       !!
       e  = 0.5d0*(dot23-som2*som3*(om2sq+om3sq+som2*som3))
       !!
       p1 = 2.d0 * (som1+som2+som3) * (om1sq*omsq23 - dot123)
       !!
       !!mpc
       !mpc    p2 = -som1 * (k23**2 - omsq23**2)
       !mpc    p3 = -(som2+som3) * (k1**2 - om1sq**2)
       !mpc    p4 = k1**2 - om1sq**2
       !!mpc---
       !!      equation  p2 rewritten to preserve numerical precision
       !!      equations p3, p4 rearranged to avoid recomputations.
       p2 = -som1 * (k23sq*(1 - tanh_k23*tanh_k23))
       p4 = (k1sq*(1-tanh_k1*tanh_k1))
       p3 = -(som2+som3) * p4
       !!mpc---
       !!      ----------------------------------------------------------------
       !!
       !!      Bash; added & used  variable domsq23 = denominator of t1
       domsq23 = omsq23 - ((som2+som3)**2)  !* Bash; needed for test below
       !!      ----------------------------------------------------------------
       !!
       !!cp4   Bash; with !cp4 ON, test if ( domsq23 .eq. 0.d0 )
       !cp4    if ( domsq23 .eq. 0.d0 ) then     !* Bash; this test was needed
       !!                                        !* when !k19p1 & !hv were OFF
       !!        domsq23=0.0 Dividing by 0.0 causes NaN; here we avoid it
       !cp4      t1 = 0.d0
       !eps      t1 = di * (p1+p2+p3) / (domsq23+eps)    !* Add eps to denominator
       !!                                                !* and may be to numerator
       !cp4    endif
       !!cp4---
       !!      Bash; with !cp4 OFF, don't test if ( domsq23 .eq. 0.d0 )
       !!      domsq23 is not = 0.0 (when !k19p1 & !hv were OFF)
       !b      t1 = di * (p1+p2+p3) / (omsq23 - ((som2+som3)**2))
       t1 = di * (p1+p2+p3) / (domsq23)
       !!cp4---
       !!      ----------------------------------------------------------------
       !!
       t2 = -di * som1 * (om1sq+omsq23)
       t3 = e * ((som1**3) * (som2+som3) - dot123 - p4)
       t4 = 0.5d0 * som1 * dot23 *                                   &
            ((som1+som2+som3) * (om2sq+om3sq) + som2*som3*(som2+som3))
       !!
       !!mpc
       !mpc    t5 = -0.5d0 * som1 *                                          &
       !mpc        (om2sq * (k3**2) * (som1+som2 + 2.d0 * som3)  +           &
       !mpc         om3sq * (k2**2) * (som1+som3 + 2.d0 * som2))
       !!mpc---
       t5 = -0.5d0 * som1 *                                          &
            (om2sq * (k3sq) * (som1+som2 + 2.d0 * som3)  +            &
            om3sq * (k2sq) * (som1+som3 + 2.d0 * som2))
       !!mpc---
       !!
       scple = scple + t1 + t2 + t3 + t4 + t5
       !!
     end do  ! do ipass=1,3
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!as  HH did division by 3 after adding 3 terms
     !as   scple = scple/3.d0
     !!as---
     !!
     !!    Bash; Added sumom   = denominator of  csqhatd in cplshr
     sumom = om1*om2*om3*(om2+om3-om1)
     !b    csqhatd = scple*scple*pi4/(om1*om2*om3*(om2+om3-om1))  !* Bash; ok
     csqhatd = scple*scple*pi4/(sumom)
     !!    ------------------------------------------------------------------
     !!
     csqd    = csqhatd / (hh**6)
     csq     = sngl(csqd)                 !* from dbl to single precision
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     RETURN
     !!

Referenced by gridsetr().

◆ gridsetr()

subroutine w3snl4md::gridsetr	(	real, intent(in)	dep,
		real, dimension(nrng), intent(in)	wka1,
		real, intent(in)	cgnrng1
	)

This routine sets up the geometric part of the Boltzmann integral.

This routine sets up the geometric part of the Boltzmann integral based on a grid of wave frequencies and directions, with wave- numbers related to frequency and depth by linear dispersion. It is adapted from Don's original code with changes to modify the indexing so there are fewer unused elements, and a number of algo- rithmic changes that are mathematically equivalent to Don's but take advantage of intrinsic functions to form smooth results with less reliance on if statements.

It calls locus-solving routines shloxr and shlocr and coupling coefficient routine cplshr. If shlocr does not converge, ierr_gr will be something other than 0 and the routine will terminate, returning ierr_gr to the calling program (see shlocr).

It returns array grad(,,), which is an estimate of the product C(k1,k2,k3,k4)*H(k1,k3,k4)*ds/|dW/dn| (where n and the k's are all vectors) as given, for example, by Eq.(7) of 'Nonlinear energy fluxes and the finite depth equilibrium range in wave spectra,' by Resio, Pihl, Tracy and Vincent (2001, JGR, 106(C4), p. 6985), as well as arrays for indexing, interpolating and weighting locus- based wavenumber vectors within the discrete solution grid.

Parameters

[in]	dep	Depth (m)
[in]	wka1	Wavenumber array at one depth dim=(nrng)
[in]	cgnrng1	Group velocity at nrng at one depth.

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 1481 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!    ------------------------------------------------------------------
     !!
     !!    it returns: kref2,kref4, jref2,jref4, wtk2,wtk4, wta2,wta4,
     !!                tfac2,tfac4  and    grad       all dim=(npts,nang,nzz)
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !! 1. Purpose :
     !!
     !!    -------------------------------------------------------------------------#
     !!                                                                             !
     !!    This routine sets up the geometric part of the Boltzmann integral        !
     !!    based on a grid of wave frequencies and directions, with wave-           !
     !!    numbers related to frequency and depth by linear dispersion.  It         !
     !!    is adapted from Don's original code with changes to modify the           !
     !!    indexing so there are fewer unused elements, and a number of algo-       !
     !!    rithmic changes that are mathematically equivalent to Don's but          !
     !!    take advantage of intrinsic functions to form smooth results with        !
     !!    less reliance on if statements.                                          !
     !!                                                                             !
     !!    It calls locus-solving routines shloxr and shlocr and coupling           !
     !!    coefficient routine cplshr.  If shlocr does not converge, ierr_gr        !
     !!    will be something other than 0 and the routine will terminate,           !
     !!    returning ierr_gr to the calling program (see shlocr).                   !
     !!                                                                             !
     !!    It returns array grad(,,), which is an estimate of the product           !
     !!    C(k1,k2,k3,k4)*H(k1,k3,k4)*ds/|dW/dn| (where n and the k's are all       !
     !!    vectors) as given, for example, by Eq.(7) of 'Nonlinear energy           !
     !!    fluxes and the finite depth equilibrium range in wave spectra,'          !
     !!    by Resio, Pihl, Tracy and Vincent (2001, JGR, 106(C4), p. 6985),         !
     !!    as well as arrays for indexing, interpolating and weighting locus-       !
     !!    based wavenumber vectors within the discrete solution grid.              !
     !!    -------------------------------------------------------------------------#
     !!
     !!
     !! 2. Method :
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    nrng      int.  Public    I   # of freq. or rings
     !!    nang      int.  Public    I   # of angles
     !!    npts      int.  Public    I   # of points on the locus
     !!    nzz       int.  Public    I   linear irngxkrng = (NK*(NK+1))/2
     !!    kzone     int.  Public    I   zone of influence = INT(alog(4.0)/alog(dfrq))
     !!    na2p1     int.  Public    I   = nang/2 + 1
     !!    np2p1     int.  Public    I   = npts/2 + 1
     !!    ------------------------------------------------------------------
     !!
     !!    dfrq      Real  Public    I   frequency multiplier for log freq. spacing
     !!    f0        Real  Public    I   = frqa(1); first freq. (Hz)
     !!    twopi     Real  Public    I   = TPI; WW3i 2*pi = 8.*atan(1.) (radians)
     !!    ainc      Real  Public    I   = DTH; WW3 angle increment (radians)
     !!    dep       Real  local     I   = depth (m)
     !!    frqa      R.A.  Public    I   = oma(:)/twopi WW3 frequency array dim=(nrng)
     !!    angl      R.A.  Public    I   = TH(1:NTH);   WW3 angles array       dim=(nang)
     !!    sinan     R.A.  Public    I   = ESIN(1:NTH); WW3 sin(angl(:)) array dim=(nang)
     !!    cosan     R.A.  Public    I   = ECOS(1:NTH); WW3 cos(angl(:)) array dim=(nang)
     !!    wka1      R.A.  local     I   = wavenumber array at one depth dim=(nrng)
     !!    cgnrng1   Real  local     I   = Group Vel. at nrng at one depth
     !!    ------------------------------------------------------------------
     !!
     !!    *** The 11 grid integration geometry arrays at one given depth
     !!    *** from gridsetr.            dim=(npts,nang,nzz,ndep)
     !!    kref2     I.A.  Public    O   Index of reference wavenumber for k2
     !!    kref4     I.A.  Public    O   Idem for k4
     !!    jref2     I.A.  Public    O   Index of reference angle      for k2
     !!    jref4     I.A.  Public    O   Idem for k4
     !!    wtk2      R.A.  Public    O   k2 Interpolation weigth along wavenumbers
     !!    wtk4      R.A.  Public    O   Idem for k4
     !!    wta2      R.A.  Public    O   k2 Interpolation weigth along angles
     !!    wta4      R.A.  Public    O   Idem for k4
     !!    tfac2     R.A.  Public    O   Norm. for interp Action Density at k2
     !!    tfac4     R.A.  Public    O   Idem for k4
     !!    grad      R.A.  Public    O   Coupling and gradient term in integral
     !!                                  grad = C * H * g**2 * ds / |dW/dn|
     !!    ------------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!     shloxr    Subr. W3SERVMD General locus solution for input vectors
     !!                              k1 & k3 when |k1| .eq. |k3|
     !!     shlocr    Subr. W3SERVMD General locus solution for input vectors
     !!                              k1 & k3 when |k1| .ne. |k3|
     !!     cplshr    Subr. W3SERVMD Calculates Boltzmann coupling coefficient
     !!                              in shallow water
     !!    ------------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!     insnl4    Subr. W3SNL4MD initialize the grid geometry
     !!    ------------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!
     !!    !/S  Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     IMPLICIT NONE
     !!
     !!    Parameter list
     !!    --------------
     real,    intent(in)  :: dep
     real,    intent(in)  :: wka1(nrng), cgnrng1  !* Use new names locally
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    Local Parameters & variables
     !!    -----------------------------
     integer              :: irng,krng, iang,kang, ipt
     integer              :: iizz, izz, ir, i
     integer              :: kmax               !* = min(irng+kzone, nrng)
     !!
     real                 :: g, gsq
     real                 :: alf0,aldfrq, wk1x,wk1y, wk3x,wk3y
     !!
     real                 :: wn2,th2, wn2d,tnh2, om2,f2,cg2, tt2,w2
     real                 :: wn4,th4, wn4d,tnh4, om4,f4,cg4, tt4,w4
     real                 :: dWdnsq,dWdn, dif13,dif14, er
     !!
     !!hv  Bash; with !hv ON, move Heaviside section up & don't use var Heaviside
     !hv   real                 :: Heaviside
     !!hv---
     !!
     real                 :: csq
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!    initial constants
     !!    ------------------
     g      = 9.806                              !* set = GRAV as in CONSTANTS
     gsq    = 96.157636                          !* set = GRAV**2
     !!
     alf0   = alog(frqa(1))                      !* ln(f0) for ir calc. below
     aldfrq = alog(dfrq)                         !* ln(dfrq)       "
     !!
     !!ini
     !!    initialize array grad
     !!    ----------------------
     grad(:,:,:)  = 0.0
     kref2(:,:,:) = 0
     kref4(:,:,:) = 0
     jref2(:,:,:) = 0
     jref4(:,:,:) = 0
     wtk2(:,:,:)  = 0.0
     wtk4(:,:,:)  = 0.0
     wta2(:,:,:)  = 0.0
     wta4(:,:,:)  = 0.0
     tfac2(:,:,:) = 0.0
     tfac4(:,:,:) = 0.0
     !!ini---
     !!------------------------------------------------------------------------------
     !!
     !!
     !!    irng and iang are k1 parameters; krng and kang are k3 parameters
     iang = 1                               !* set = 1 and will remain = 1
     !!
     !!20
     do 20 irng=1,nrng
       !!kz
       kmax = min(irng+kzone, nrng)   !* Bash; Sometimes a locus pt is outside nrng
       !kz     kmax = min(irng+kzone, nrng-1) !* Bash; Taking 1 out will not affect kzone, try it
       !!kz---
       !!kz---
       !!
       wk1x   = wka1(irng)
       wk1y   = 0.0                         !* set = 0.0 and will remain = 0.0
       iizz = (nrng-1)*(irng-1)-((irng-2)*(irng-1))/2
       !!30
       !!kz
       do 30 krng=irng,kmax
         !!kz---
         !kz     do 30 krng=irng,nrng
         !!
         !!        Bash; check1 - change this ratio from > 4 to > 3   and
         !!              make it consistent with similar test done in subr. snlr_'s
         !kz       if ( frqa(krng)/frqa(irng) .gt. 2. ) go to 30  !* Bash; use .gt. 2 for speed
         !kz       if ( frqa(krng)/frqa(irng) .gt. 3. ) go to 30  !* original snlr_'s
         !kz       if ( frqa(krng)/frqa(irng) .gt. 4. ) go to 30  !* original gridsetr
         !!kz---
         izz = krng+iizz
         !!40
         do 40 kang=1,nang
           !!
           wk3x = wka1(krng)*cosan(kang)
           wk3y = wka1(krng)*sinan(kang)
           !!
           if ( krng.eq.irng ) then              !* wn3 = wn1
             !!
             !!ba1         Bash; skip k1 but keep the opposite angle to k1 - orig setting
             !!ba1               remember here iang = 1
             if ( kang .eq. 1 ) go to 40         !* th3 = th1
             !!ba1---
             !!            ----------------------------------------------------------
             !!            -- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
             call shloxr ( dep,       wk1x,wk1y,wk3x,wk3y )
             !!            -- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
             !!            it returns: wk2x, wk2y, wk4x, wk4y & ds all dim=(npts)
             !!                        and all are PUBLIC
             !!            -- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
             !!            ----------------------------------------------------------
             !!
           else                                  !* wn3 > wn1
             !!
             !!            ----------------------------------------------------------
             !!            -- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
             call shlocr ( dep,       wk1x,wk1y,wk3x,wk3y )
             !!            -- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
             !!            it returns: wk2x, wk2y, wk4x, wk4y & ds all dim=(npts)
             !!                        and all are PUBLIC
             !!            -- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
             !!            ----------------------------------------------------------
             !!
           end if   !!  if ( krng.eq.irng )
           !!
           !!          set the Heaviside coefficient
           !b          dif13 = (wk1x-wk3x)**2   +  (wk1y-wk3y)**2        !* wk1y = 0.0
           !b          dif13 = (wk1x-wk3x)**2   +       (wk3y)**2
           dif13 = (wk1x-wk3x)*(wk1x-wk3x) + wk3y*wk3y
           !!50
           do 50 ipt=1,npts
             !!
             !!xlc1        Bash; skip k1 but keep the opposite angle to k1 - original setting
             if ( kang.eq.1 ) then                        !* th3=+th1, iang=1
               if (ipt.eq.1 .or. ipt.eq.np2p1) go to 50   !* skip x-axis loci
             end if
             !!xlc1---
             !!            ----------------------------------------------------------
             !!
             !!hv          Bash; with !hv ON, move Heaviside section from below to here
             !!            Bash moved this section here. *** Check first compute after ***
             !!            Skip first then compute only if Heaviside=1, without using it
             !!            i.e. compute only if dif13.le.dif14 with Heaviside=1 omitted.
             !!            Note; with !hv option is ON, you don't need to turn options
             !!            ----  !k19p1 nor !cp4 ON.aYou only need one of the three.
             !!            ----------------------------------------------------------
             !!            set the Heaviside coefficient
             !b            dif14 = (wk1x-wk4x(ipt))**2 + (wk1y-wk4y(ipt))**2 !* wk1y=0.0
             !b            dif14 = (wk1x-wk4x(ipt))**2 +      (wk4y(ipt))**2
             dif14 = (wk1x-wk4x(ipt))*(wk1x-wk4x(ipt)) +             &
                  wk4y(ipt)*wk4y(ipt)
             !!
             if ( dif13 .gt. dif14 ) go to 50    !* skip, don't compute
             !!
             !b            if ( dif13 .gt. dif14 ) then
             !b               Heaviside = 0.                   !* Eq(12) of RPTV
             !b               go to 50
             !b            else
             !b               Heaviside = 1.                   !* Eq(11) of RPTV
             !b            end if
             !!hv---
             !!            ----------------------------------------------------------
             !!
             !!            Set the coupling coefficient for ipt'th locus position
             !!            ----------------------------------------------------------
             !!            -- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
             call cplshr ( wk4x(ipt),wk4y(ipt), wk3x,wk3y,           &
                  wk2x(ipt),wk2y(ipt), dep, csq,            &
                  irng,krng,kang,ipt )
             !!            -- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
             !!            it returns: the coupling coefficient  csq
             !!            -- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
             !!            ----------------------------------------------------------
             !!
             !!wn2         Set parameters related to ipt'th locus wavenumber vector k2
             !!            ----------------------------------------------------------
             !b            wn2  = sqrt(wk2x(ipt)**2 + wk2y(ipt)**2)               !* k2
             wn2  = sqrt(wk2x(ipt)*wk2x(ipt) + wk2y(ipt)*wk2y(ipt)) !* k2
             th2  = atan2(wk2y(ipt),wk2x(ipt))                !* k2 direction
             if ( th2 .lt. 0. ) th2 = th2 + twopi             !* +ve in radians
             wn2d = wn2*dep                                   !* k2*depth
             tnh2 = tanh(wn2d)                                !* tanh(k2*depth)
             om2  = sqrt(g*wn2*tnh2)                          !* omega2 (rad)
             !b            cg2  = 0.5*(om2/wn2)*(1.+wn2d*(1.-tnh2**2)/tnh2) !* group velocity
             cg2  = 0.5*(om2/wn2)*(1.+wn2d*((1./tnh2)-tnh2))  !* group velocity
             f2   = om2/twopi                                 !* f2 (Hz)
             !!            ----------------------------------------------------------
             !!
             !!wn4         Set parameters related to ipt'th locus wavenumber vector k4
             !!            ----------------------------------------------------------
             !b            wn4  = sqrt(wk4x(ipt)**2 + wk4y(ipt)**2)
             wn4  = sqrt(wk4x(ipt)*wk4x(ipt) + wk4y(ipt)*wk4y(ipt))
             th4  = atan2(wk4y(ipt),wk4x(ipt))
             if ( th4 .lt. 0. ) th4 = th4 + twopi
             wn4d = wn4*dep
             tnh4 = tanh(wn4d)
             om4  = sqrt(g*wn4*tnh4)
             !b            cg4  = 0.5*(om4/wn4)*(1.+wn4d*(1.-tnh4**2)/tnh4)
             cg4  = 0.5*(om4/wn4)*(1.+wn4d*((1./tnh4)-tnh4))
             f4   = om4/twopi
             !!            ----------------------------------------------------------
             !!
             !!
             !!hv          Bash; with !hv ON, move Heaviside section up
             !!            Bash moved this section up. Check first compute after.
             !!            ----------------------------------------------------------
             !!            set the Heaviside coefficient
             !b            dif14 = (wk1x-wk4x(ipt))**2 + (wk1y-wk4y(ipt))**2 !* wk1y=0.0
             !b            dif14 = (wk1x-wk4x(ipt))**2 +      (wk4y(ipt))**2
             !hv           dif14 = (wk1x-wk4x(ipt))*(wk1x-wk4x(ipt)) +             &
             !hv                                          wk4y(ipt)*wk4y(ipt)
             !hv           if ( dif13 .gt. dif14 ) then
             !hv              Heaviside = 0.                   !* Eq(12) of RPTV
             !hv           else
             !hv              Heaviside = 1.                   !* Eq(11) of RPTV
             !hv           end if
             !!hv---
             !!            ----------------------------------------------------------
             !!
             !!
             !!            dWdn is the same as sqrt(zzsum) in Don's code, here reduced to a
             !!            simpler but mathematically equivalent form that should vary
             !!            smoothly between deep and intermediate water owing to identities
             !!            using the computer's tanh() function
             !!            ----------------------------------------------------------
             !!
             !!            set grad(,,);
             !!            looks like the g^2 goes with csq (Webb'1978, eq. A2)
             !!            ----------------------------------------------------------
             !!
             !b            dWdnsq = cg2**2  - 2.*cg2*cg4 * cos(th2-th4) + cg4**2
             dwdnsq = cg2*cg2 - 2.*cg2*cg4 * cos(th2-th4) + cg4*cg4
             !!            ----------------------------------------------------------
             !!
             dwdn               = sqrt(dwdnsq)
             !!            ----------------------------------------------------------
             !!
             !!hv          Bash; with !hv ON, don't use var Heaviside (by here it's = 1.0)
             grad(ipt,kang,izz) =           ds(ipt)*csq*gsq/dwdn
             !!hv---
             !hv           grad(ipt,kang,izz) = Heaviside*ds(ipt)*csq*gsq/dWdn
             !!hv---
             !!            ----------------------------------------------------------
             !!            ==========================================================
             !!
             !!
             !!            Set interpolation, indexing and weight parameters for
             !!            computations along wavenumber radials
             !!            ----------------------------------------------------------
             !!
             !!f2          --------------------
             if ( f2 .lt. f0 ) then
               wtk2(ipt,kang,izz)  = 1.
               tfac2(ipt,kang,izz) = 0.
               kref2(ipt,kang,izz) = 1
             else
               ir = 1+int((alog(f2)-alf0)/aldfrq)
               if ( ir+1 .gt. nrng ) then
                 wtk2(ipt,kang,izz)  = 0.
                 er = (wka1(nrng)/wn2)**(2.5)
                 tt2= er*(cg2/cgnrng1)*(frqa(nrng)/f2)*(wka1(nrng)/wn2)
                 tfac2(ipt,kang,izz) = tt2
                 kref2(ipt,kang,izz) = nrng - 1
               else
                 w2 = (f2-frqa(ir))/(frqa(ir+1)-frqa(ir))
                 wtk2(ipt,kang,izz)  = 1. - w2
                 tfac2(ipt,kang,izz) = 1.
                 kref2(ipt,kang,izz) = ir
               end if
             end if
             !!            ----------------------------------------------------------
             !!
             !!f4          --------------------
             if ( f4 .lt. f0 ) then
               wtk4(ipt,kang,izz)  = 1.
               tfac4(ipt,kang,izz) = 0.
               kref4(ipt,kang,izz) = 1
             else
               ir = 1+int((alog(f4)-alf0)/aldfrq)
               if ( ir+1 .gt. nrng ) then
                 wtk4(ipt,kang,izz)  = 0.
                 er = (wka1(nrng)/wn4)**2.5
                 tt4= er*(cg4/cgnrng1)*(frqa(nrng)/f4)*(wka1(nrng)/wn4)
                 tfac4(ipt,kang,izz) = tt4
                 kref4(ipt,kang,izz) = nrng - 1
               else
                 w2 = (f4-frqa(ir))/(frqa(ir+1)-frqa(ir))
                 wtk4(ipt,kang,izz)  = 1. - w2
                 tfac4(ipt,kang,izz) = 1.
                 kref4(ipt,kang,izz) = ir
               end if
             end if
             !!            ----------------------------------------------------------
             !!
             !!
             !!            Set indexing and weight parameters for computations around
             !!            azimuths; it appears that jref2 and jref4 should be bounded
             !!            between 0 and nang-1 so that when iang (=1,nang) is added in
             !!            the integration section, the proper bin index will arise;
             !!            the weights wta2 and wta4 seem to be the fractional bin
             !!            widths between th2 or th4 and the next increasing
             !!            directional bin boundary
             !!            ----------------------------------------------------------
             !!
             i = int(th2/ainc)
             wta2(ipt,kang,izz)  = 1. - abs(th2-i*ainc)/ainc
             if ( i .ge. nang )  i = i - nang
             jref2(ipt,kang,izz) = i
             !mpc          jref2(ipt,kang,izz) = MOD(i,nang) !* is this better that the above two lines?
             !!
             i = int(th4/ainc)
             wta4(ipt,kang,izz)  = 1. - abs(th4-i*ainc)/ainc
             if ( i .ge. nang )  i = i - nang
             jref4(ipt,kang,izz) = i
             !mpc          jref4(ipt,kang,izz) = MOD(i,nang) !* is this better that the above two lines?
             !!
 50        end do                              !* end of ipt loop
           !!
 40      end do                                !* end of kang loop
         !!
 30    end do                                  !* end of krng loop
       !!
 20  end do                                    !* end of irng loop
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     RETURN
     !!

References ainc, cosan, cplshr(), dfrq, ds, f0, frqa, grad, jref2, jref4, kref2, kref4, kzone, nang, np2p1, npts, shlocr(), shloxr(), sinan, tfac2, tfac4, twopi, wk2x, wk2y, wk4x, wk4y, wta2, wta4, wtk2, and wtk4.

Referenced by insnl4().

◆ insnl4()

subroutine w3snl4md::insnl4

It reads look-up tables (generated by gridsetr) if file exists otherwise it must generate the look-up tables file.

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 219 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!
     !!    it returns: 11 look-up tables arrays dim=(npts,nang,nzz,ndep)
     !!                kref2_tbl, kref4_tbl, jref2_tbl, jref4_tbl,
     !!                wtk2_tbl,  wtk4_tbl,  wta2_tbl,  wta4_tbl,
     !!                tfac2_tbl, tfac4_tbl & grad_tbl
     !!                plus       pha_tbl  dim=(nrng,ndep)
     !!                and        dep_tbl  dim=(ndep)
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !! 1. Purpose :
     !!    It reads look-up tables (generated by gridsetr) if file exists
     !!    otherwise it must generate the look-up tables file
     !!
     !! 2. Method :
     !!    See subr gridsetr and subr W3SNL4 (or subr. W3IOGR)
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    nrng      int.  Public    I   # of freq. or rings
     !!    nang      int.  Public    I   # of angles
     !!    npts      int.  Public    I   # of points on the locus
     !!    ndep      int.  Public    I   # of depths in look-up tables
     !!    dfrq      Real  Public    I   frequency multiplier for log freq. spacing
     !!    dep_tbl   R.A.  Public    O   depthes in Look-up tables arrays dim=(ndep)
     !!    grdfname  chr.  Local     -   Look-up tables filename (C*80)
     !!    ------------------------------------------------------------------
     !!
     !!    *** The 11 look-up tables for grid integration geometry arrays
     !!    *** at all selected 'ndep' depths defined in dep_tbl(ndep)' array
     !!    *** from gridsetr.            dim=(npts,nang,nzz,ndep)
     !!    kref2_tbl I.A.  Public    O   Index of reference wavenumber for k2
     !!    kref4_tbl I.A.  Public    O   Idem for k4
     !!    jref2_tbl I.A.  Public    O   Index of reference angle      for k2
     !!    jref4_tbl I.A.  Public    O   Idem for k4
     !!    wtk2_tbl  R.A.  Public    O   k2 Interpolation weigth along wavenumbers
     !!    wtk4_tbl  R.A.  Public    O   Idem for k4
     !!    wta2_tbl  R.A.  Public    O   k2 Interpolation weigth along angles
     !!    wta4_tbl  R.A.  Public    O   Idem for k4
     !!    tfac2_tbl R.A.  Public    O   Norm. for interp Action Density at k2
     !!    tfac4_tbl R.A.  Public    O   Idem for k4
     !!    grad_tbl  R.A.  Public    O   Coupling and gradient term in integral
     !!                                  grad = C * H * g**2 * ds / |dW/dn|
     !!    ------------------------------------------------------------------
     !!
     !!    *** The 11 grid integration geometry arrays at one given depth
     !!    *** from gridsetr.            dim=(npts,nang,nzz,ndep)
     !!    kref2     I.A.  Public    O   Index of reference wavenumber for k2
     !!    kref4     I.A.  Public    O   Idem for k4
     !!    jref2     I.A.  Public    O   Index of reference angle      for k2
     !!    jref4     I.A.  Public    O   Idem for k4
     !!    wtk2      R.A.  Public    O   k2 Interpolation weigth along wavenumbers
     !!    wtk4      R.A.  Public    O   Idem for k4
     !!    wta2      R.A.  Public    O   k2 Interpolation weigth along angles
     !!    wta4      R.A.  Public    O   Idem for k4
     !!    tfac2     R.A.  Public    O   Norm. for interp Action Density at k2
     !!    tfac4     R.A.  Public    O   Idem for k4
     !!    grad      R.A.  Public    O   Coupling and gradient term in integral
     !!                                  grad = C * H * g**2 * ds / |dW/dn|
     !!    ------------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!     gridsetr  Subr. W3SERVMD Calc. the 11 grid geometry arrays for one depth
     !!    ------------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!     STRACE    Subr. W3SERVMD Subroutine tracing.
     !!     W3SNL4    Subr. W3SNL4MD Interface for TSA  nonlinear interactions
     !!     - or -    the option below was not used
     !!     W3IOGR    Subr. W3INITMD Initialization (called by W3SHEL or WMINIT)
     !!    ------------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!    !/S  Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
 #ifdef W3_S
     USE w3servmd, ONLY: strace
 #endif
     USE w3odatmd, ONLY: ndse, ndst, ndso
 #ifdef W3_MPI
     USE wmmdatmd, ONLY: nmpscr, improc
 #endif
     USE constants, ONLY: file_endian
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     IMPLICIT NONE
     !!
     !!    ==================================================================
     !!
     !!    Local variables & Parameters
     !!    ----------------------------
 #ifdef W3_S
     INTEGER, SAVE           :: IENT = 0
 #endif
     !!
     integer              :: irng               !* dummy integer
     integer              :: nd                 !* dummy integer
     !!
     logical              :: unavail = .true.
     logical              :: file_exists
     character            :: grdfname*80
     integer              :: io_unit
     !!
     !!
     !!    Dimension wv# array and
     !!    declare local var. dep2, cvel & cgnrng at nrng & depth 'nd'
     real                 :: wka2(nrng)  !* wv# array at depth nd (local)
     real                 :: pha2(nrng)  !* pha array at depth nd (local)
     real                 :: dep2        !* = dep_tbl(nd) depth at nd
     real                 :: cvel        !* Phase Velocity at (nrng,nd)
     real                 :: cgnrng      !* Group Velocity at (nrng,nd)
     real                 :: dwka        !* dummy storage for dk
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
 #ifdef W3_S
     CALL strace (ient, 'W3SNL4')
 #endif
     !!
     !!    ==================================================================
     !!
     !!
     !!-1  Make-up the filename from the main parameters
     !!    ------------------------------------------------------------------
     !b    example filename; grdfname = 'grd_1.1025_35_36_30_37.dat'
     !!    grdfname = 'grd_dfrq_nr_na_np_nd.dat'
     !!       where    fm = freq. mult. (dfrq) ex. dfrq = 1.1025 (F6.4)
     !!                nr = # of rings  (nrng) ex. nrng = 35     (I2.2)
     !!                na = # of angles (nang) ex. nang = 36     (I2.2)
     !!                np = # of points (npts) ex. npts = 30     (I2.2)
     !!                nd = # of depths (ndep) ex. nd   = 37     (I2.2)
     write(grdfname,'(A,F6.4,A,I2.2,A,I2.2,A,I2.2,A,I2.2,A)')        &
          'grd_', dfrq,'_', nrng,'_', nang,'_',  npts,'_', ndep, '.dat'
     !!    ==================================================================
     !!
     !!
     !!-2  Check if the propre gridsetr Look-up tables file is available.
     !!    if available read it and if not must generate it (by calling gridsetr)
     !!    ------------------------------------------------------------------
     INQUIRE ( file=grdfname, exist=file_exists )
     !!    Assign an unused UNIT number to io_unit.
     !!    Note; It's important to look for an available unused number
     io_unit   = 60
     do while (unavail)
       io_unit = io_unit + 1
       INQUIRE ( io_unit, opened=unavail )
     enddo
     !prt  print *, 'io_unit = ', io_unit
     !!    ==================================================================
     !!
     !!
     !!
     !!
     IF ( file_exists ) THEN
       !!
       !!-3    File exists open it and read it
       !!
 #ifdef W3_MPI
       if ( improc .eq. nmpscr ) then
 #endif
         write ( ndso, 900 ) grdfname
 #ifdef W3_MPI
       end if
 #endif
       !!
       open (unit=io_unit, file=grdfname, status='old',              &
            access='sequential', action='read', form='unformatted', convert=file_endian)
       read (io_unit)  kref2_tbl, kref4_tbl, jref2_tbl, jref4_tbl,   &
            wtk2_tbl,  wtk4_tbl,  wta2_tbl,  wta4_tbl,    &
            tfac2_tbl, tfac4_tbl, grad_tbl,               &
            pha_tbl,   dep_tbl
       close (io_unit)
       !!      ----------------------------------------------------------------
       !!
     ELSE      !* ELSE IF ( file_exists )
       !!
       !!
       !!-4    File does not exist, create it here
       !!
 #ifdef W3_MPI
       if ( improc .eq. nmpscr ) then
 #endif
         write ( ndso, 901 ) grdfname
 #ifdef W3_MPI
       end if
 #endif
       !!      ----------------------------------------------------------------
       !!
       !!-4a   Define Look-up tables depth array 'dep_tbl(ndep)' for ndep=37
       !!      with depths are +ve values
       !!      ----------------------------------------------------------------
       dep_tbl(1:ndep) =                                             &
            (/  2.,  4.,  6.,  8., 10., 12., 14., 16., 18., 20.,   &
            25., 30., 35., 40., 45., 50., 55., 60., 65., 70.,   &
            80., 90.,100.,110.,120.,130.,140.,150.,160.,170.,   &
            220.,270.,320.,370.,420.,470.,520.  /)
       !prt    print *, ' ndep = ', ndep
       !prt    print *, ' dep_tbl(1:ndep) = ', dep_tbl
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       do nd = 1,ndep
         !!
         !!
         !!-4b     For given new depth dep2 = dep_tbl(nd) calculate
         !!        a new array wka2(:) & new cgnrng corresp. to this depth
         dep2 = dep_tbl(nd)
         do irng=1,nrng
           wka2(irng) = wkfnc(frqa(irng),dep2)
         end do
         cvel    = oma(nrng)/wka2(nrng)         !* Phase Vel. at (nrng,nd)
         cgnrng  = cgfnc(frqa(nrng),dep2,cvel)  !* Group Vel. at (nrng,nd)
         !!        --------------------------------------------------------------
         !!        ==============================================================
         !!
         !!
         !!-4c     Call gridsetr for this depth at nd
         !!        --------------------------------------------------------------
         !!        -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         call gridsetr ( dep2, wka2, cgnrng )
         !!        -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         !!        it returns: 11 gridsetr arrays which are declared PUBLIC
         !!                    kref2,kref4, jref2,jref4, wtk2,wtk4, wta2,wta4,
         !!                    tfac2,tfac4  and   grad   all dim=(npts,nang,nzz)
         !!        -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         !!        --------------------------------------------------------------
         !!        ==============================================================
         !!
         !!-4d     Store in Look-up tables arrays at depth bin # 'nd'
         kref2_tbl(:,:,:,nd) = kref2(:,:,:)
         kref4_tbl(:,:,:,nd) = kref4(:,:,:)
         jref2_tbl(:,:,:,nd) = jref2(:,:,:)
         jref4_tbl(:,:,:,nd) = jref4(:,:,:)
         wtk2_tbl(:,:,:,nd)  = wtk2(:,:,:)
         wtk4_tbl(:,:,:,nd)  = wtk4(:,:,:)
         wta2_tbl(:,:,:,nd)  = wta2(:,:,:)
         wta4_tbl(:,:,:,nd)  = wta4(:,:,:)
         tfac2_tbl(:,:,:,nd) = tfac2(:,:,:)
         tfac4_tbl(:,:,:,nd) = tfac4(:,:,:)
         grad_tbl(:,:,:,nd)  = grad(:,:,:)
         !!        --------------------------------------------------------------
         !!        ==============================================================
         !!
         !!
         !!-4e     Calculate pha2(:) at nd depth and store it in pha_tbl(:,:)
         !!        pha2()=k*dk*dtheta, the base area at a grid intersection
         !!        for use in integration of 2-D Density functions.
         !!        --------------------------------------------------------------
         !!        Below: variable dwka = dk centered at ring 1 (between 0 & 2)
         !!        and computed    pha2(1) = k*dk*dtheta at ring 1
         !!        with wkfnc(frqa(1)/dfrq,dep2) is like wka2(0)
         !!        --assuming frqa(1)/dfrq is like frqa(0)
         dwka    = ( wka2(2) - wkfnc(frqa(1)/dfrq,dep2) ) / 2.
         pha2(1) = wka2(1)*dwka*ainc
         !!
         do irng=2,nrng-1
           !!          Below: variable dwka = dk centered at irng (between irng-1 & irng+1)
           !!          and computed    pha2(irng) = k*dk*dtheta at irng
           dwka       = ( wka2(irng+1) - wka2(irng-1) ) / 2.
           pha2(irng) = wka2(irng)*dwka*ainc
         end do
         !!
         !!        Below: variable dwka = dk centered at nrng (between nrng-1 & nrng+1)
         !!        and computed    pha2(nrng) = k*dk*dtheta at nrng
         !!        with wkfnc(dfrq*frqa(nrng),dep2) is like wka2(nrng+1)
         !!        --assuming dfrq*frqa(nrng) is like frqa(nrng+1)
         dwka       = ( wkfnc(dfrq*frqa(nrng),dep2) - wka2(nrng-1) ) / 2.
         pha2(nrng) = wka2(nrng)*dwka*ainc
         !!        --------------------------------------------------------------
         !!        ==============================================================
         !!
         !!
         !!-4f     Store pha2(:) at nd in  pha_tbl(:,nd) to be added to Look-up tables
         pha_tbl(1:nrng, nd) = pha2(1:nrng)
         !!        --------------------------------------------------------------
         !!        ==============================================================
         !!
         !!
       end do ! nd = 1,ndep
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       !!-5    Ounce the Look-up tables arrays are full write it out to 'io_unit'
       !!
 #ifdef W3_MPI
       if ( improc .eq. nmpscr ) then
 #endif
         write( ndso,902 )
         open (unit=io_unit, file=grdfname, status='new',              &
              access='sequential', action='write', form='unformatted', convert=file_endian)
         write (io_unit) kref2_tbl, kref4_tbl, jref2_tbl, jref4_tbl,   &
              wtk2_tbl,  wtk4_tbl,  wta2_tbl,  wta4_tbl,    &
              tfac2_tbl, tfac4_tbl, grad_tbl,               &
              pha_tbl,   dep_tbl
         close (io_unit)
         write( ndso,903 ) grdfname
 #ifdef W3_MPI
       end if
 #endif
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
     ENDIF     !* End  IF ( file_exists )
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     RETURN
     !!
 900 format ( '  grdfname does exist = ',a/                          &
          '  open, read & close file ' )
     !!
 901 format ( '  grdfname does not exist = ',a/                      &
          '  Generate look-up table arrays ' )
     !!
 902 format ( '  Done generating look-up table arrays ----------- ' )
     !!
 903 format ( '  Done writing & closing grdfname ', a )
     !!

References ainc, cgfnc(), dep_tbl, dfrq, file(), constants::file_endian, frqa, grad, grad_tbl, gridsetr(), wmmdatmd::improc, jref2, jref2_tbl, jref4, jref4_tbl, kref2, kref2_tbl, kref4, kref4_tbl, nang, ndep, w3odatmd::ndse, w3odatmd::ndso, w3odatmd::ndst, wmmdatmd::nmpscr, npts, oma, pha_tbl, w3servmd::strace(), tfac2, tfac2_tbl, tfac4, tfac4_tbl, wkfnc(), wta2, wta2_tbl, wta4, wta4_tbl, wtk2, wtk2_tbl, wtk4, and wtk4_tbl.

Referenced by w3snl4().

◆ interp2()

subroutine w3snl4md::interp2 ( real, dimension(nrng,nang), intent(inout) X )

Interpolate bi-linearly to fill in tsa/fbi & diag/diag2 arrays.

Optionally smooth the interior and the corners after alternating the irng, iang, krng & kang loops in snlr's.

Parameters

[in,out] X Array to be interpolated and smoothed.

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 5018 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!
     !! 1. Purpose :
     !!
     !!    Interpolate bi-linearly to fill in tsa/fbi & diag/diag2 arrays
     !!    and then (optional) smooth the interior and the corners
     !!    after alternating the irng, iang, krng & kang loops in snlr's
     !!
     !! 2. Method :
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    nrng      int.  Public    I   # of freq. or rings
     !!    nang      int.  Public    I   # of angles
     !!    ismo      int.  Local     I   switch;  ismo=0 skip smoothing
     !!                                           ismo.ne.0 do smoothing
     !!    X         R.A.  Local    I/O  Array to be ineterp. & smoothing
     !!                                  its returned to snlr_tsa as tsa or diag
     !!                                           and to snlr_fbi as fbi or diag2
     !!                                           dim=(nrng,nang)
     !!    ------------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!    ----------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!     snlr_tsa  Subr. W3SNL4MD Computes dN(k,theta)/dt for TSA
     !!     --------                 due to wave-wave inter. (set itsa = 1)
     !!     snlr_fbi  Subr. W3SNL4MD Computes dN(k,theta)/dt for FBI
     !!     --------                 due to wave-wave inter. (set itsa = 0)
     !!    ----------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!
     !!    !/S  Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     IMPLICIT NONE
     !!
     !!
     !!    Parameter list
     !!    --------------
     REAL,    INTENT(INOUT)  :: X(nrng,nang)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    Local Parameters
     !!    ----------------
     integer              :: irng, iang
     real                 :: Y(nrng,nang)  !* dummy array used in smoothing
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!-0  Initial Y(:,:) array before it's computed
     !!ini
     y(:,:) = 0.0
     !!ini---
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !!-1  Interpolate using simple 2 point averaging to fill in X
     !!    Remeber: nang must be an even number ==> nang-2 is an even number
     !!    and      nrng must be an odd  number ==> nrng-1 is an even number
     !!    Example numbers used here are for nrng=35 & nang=36
     !!    ------------------------------------------------------------------
     !!
     !!-1a For every calculated iang (1,3,5,..,nang-1=35)
     !!    fill in missing irng's    (2,4,6,..,nrng-1=34)
     do iang=1,nang-1,2      !* = 1,3,5,...,nang-1=35
       do  irng=2,nrng-1,2      !* = 2,4,6,...,nrng-1=34
         x(irng,iang) = 0.5 * ( x(irng-1,iang) + x(irng+1,iang) )
       end do
     end do
     !!    ------------------------------------------------------------------
     !!
     !!-1b Now, for every irng (1,2,3,..,nrng  =35)
     !!    fill missing iang's (2,4,6,..,nang-2=34)
     do irng=1,nrng          !* 1,2,3,..,nrng  =35
       do iang=2,nang-2,2      !* 2,4,6,..,nang-2=34
         x(irng,iang) = 0.5 * ( x(irng,iang-1) + x(irng,iang+1) )
       end do
     end do
     !!    ------------------------------------------------------------------
     !!
     !!-1c for iang = nang (special case since nang is an even number)
     do irng=1,nrng
       x(irng,nang) = 0.5 * ( x(irng,nang-1) + x(irng,1) )
     end do
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !!
     !!    Skip smoothing only if ismo = 0
     !!
     !!
     if ( ismo.eq.0 ) goto 99
     !!
     !!
     !!
     !!-2  Smoothing the 2D array X into array Y
     !!
     !!-2a Smoothing the interior [2;nrng-1] x [2:nang-1]
     !!-   Using 9 points averaged with equal weights.
     !!-   Here use the dummy array so we don't spoil the original array.
     do irng=2,nrng-1
       do iang=2,nang-1
         y(irng,iang)=(x(irng-1,iang-1)+x(irng-1,iang)+x(irng-1,iang+1) + &
              x(irng,  iang-1)+x(irng,  iang)+x(irng,  iang+1) + &
              x(irng+1,iang-1)+x(irng+1,iang)+x(irng+1,iang+1))/9.
       end do
     end do
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !!-3  Smooth first & last line at iang=1 & iang=nang  (special cases)
     !!
     !!-3a Smooth line at iang = 1     (special case)
     !!-   Using 9 points averaged with equal weights.
     do irng=2,nrng-1
       y(irng, 1) = (x(irng-1,nang) + x(irng-1, 1) + x(irng-1, 2) +  &
            x(irng,  nang) + x(irng,   1) + x(irng,   2) +  &
            x(irng+1,nang) + x(irng+1, 1) + x(irng+1, 2) )/9.
     end do
     !!    ------------------------------------------------------------------
     !!
     !!-3b Smooth line at iang = nang  (special case)
     !!-   Using 9 points averaged with equal weights.
     do irng=2,nrng-1
       y(irng,nang)=(x(irng-1,nang-1) +x(irng-1,nang) +x(irng-1,1) + &
            x(irng,  nang-1) +x(irng,  nang) +x(irng,  1) + &
            x(irng+1,nang-1) +x(irng+1,nang) +x(irng+1,1))/9.
     end do
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !!-4  Smooth first & last col. at irng=1 & irng=nrng  (special cases)
     !!
     !!-4a Smooth col. at irng = 1     (low frq. can be skipped)
     !!-   Using 6 points averaged with equal weights.
     do iang=2,nang-1
       y(1,iang)   = (x(1,iang-1) + x(1,iang) + x(1,iang+1) +        &
            x(2,iang-1) + x(2,iang) + x(2,iang+1) )/6.
     end do
     !!    ------------------------------------------------------------------
     !!
     !!-4b Smooth col. at irng = nrng  (high frq. can be skipped)
     !!-   Using 6 points averaged with equal weights.
     do iang=2,nang-1
       y(nrng,iang)=(x(nrng-1,iang-1)+x(nrng-1,iang)+x(nrng-1,iang+1)+ &
            x(nrng, iang-1)+x(nrng, iang)+x(nrng, iang+1) )/6.
     end do
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !!-5  Smooth the 4 corners (optional):  <== Skip no sig. effect
     !!-   Using 6 points averaged with equal weights
     !!
     !!-5a Corner (1, 1)
     y(1, 1)      =( x(1,nang) + x(1, 1) + x(1, 2) +                 &
          x(2,nang) + x(2, 1) + x(2, 2) )/6.0
     !!    ------------------------------------------------------------------
     !!
     !!-5b Corner (nrng,1)
     y(nrng,1)    =( x(nrng-1,nang) + x(nrng-1,1) + x(nrng-1,2) +    &
          x(nrng,  nang) + x(nrng,  1) + x(nrng,  2) )/6.0
     !!    ------------------------------------------------------------------
     !!
     !!-5c Corner (1,nang)
     y(1,nang)    =( x(1,nang-1) + x(1,nang) + x(1, 1) +             &
          x(2,nang-1) + x(2,nang) + x(2, 1) ) / 6.
     !!    ------------------------------------------------------------------
     !!
     !!-5d Corner (nrng,nang)
     y(nrng,nang) =( x(nrng-1,nang-1) +x(nrng-1,nang) +x(nrng-1,1) + &
          x(nrng,  nang-1) +x(nrng,  nang) +x(nrng,  1) )/6.
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !!-6  Final, dump smoothed array Y(:,:) into X(:,:) to be returned
     !!
     !!-6a Done with X(:,:) re-initial before it's replaced by Y(:,:)
     !!ini
     x(:,:) = 0.0
     !!ini---
     !!
     !!-6b Dump smoothed array Y(:,:) into X(:,:) to be returned
     do iang=1,nang
       do irng=1,nrng
         x(irng,iang) = y(irng,iang)
       end do
     end do
     !!    Bash; can simplify in one line
     !b    X(1:nrng, 1:nang) = Y(1:nrng, 1:nang)
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
 99  continue
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     RETURN
     !!

References ismo.

Referenced by snlr_fbi(), and snlr_tsa().

◆ optsa2()

subroutine w3snl4md::optsa2	(	integer, intent(in)	nrmn,
		integer, intent(in)	nrmx,
		integer, intent(in)	npk,
		real, intent(in)	fpk,
		integer, intent(in)	nbins,
		real, dimension(nrng), intent(in)	wka,
		real, dimension(nrng), intent(in)	cga
	)

Splits the Action Density into two parts.

    (1) large-scale part dens1(nrng,nang)  and
    (2) small-scale part dens2(nrng,nang)
    dens1 & dens2 in Polar Action Density (k,theta) space Norm. (in k)

Parameters

[in]	nrmn	Number of first frequency bin in [1,nrng-1].
[in]	nrmx	Number of last frequency bin in [2,nrng].
[in]	npk	Number of peak frequency in [2,nrng-1].
[in]	fpk	Peak frequency [Hz] of initial frequency spectrum.
[in]	nbins	number of bins (see more below).
[in]	wka	Wavenumbers array [1/m].
[in]	cga	Group velocities array [m/s].

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 3259 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!    ------------------------------------------------------------------
     !!
     !!    It returns variables dens1(nrng,nang) and dens2(nrng,nang)
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !! 1. Purpose :
     !!
     !!    Splits the Action Density into two parts:
     !!    (1) large-scale part dens1(nrng,nang)  and
     !!    (2) small-scale part dens2(nrng,nang)
     !!    dens1 & dens2 in Polar Action Density (k,theta) space Norm. (in k)
     !!
     !! 2. Method :
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!op2
     !!    nrmn      int.  Local     I   number of first freq. bin in [1,nrng-1]
     !!    nrmx      int.  Local     I   number of last  freq. bin in [2,nrng]
     !!    npk       int.  Local     I   number of peak frequency  in [2,nrng-1]
     !!    nbins     int.  Local     I   actual # of bins > npk  (incl. nfs)  or
     !!                                  actual # of bins > npk2 (incl. nrng)
     !!                                  to guarantee a min 1 bin in equi. range
     !!                                  (see subr. W3SNL4)
     !!    ------------------------------------------------------------ !!op2
     !!
     !!    nrng      int.  Public    I   # of freq. or rings
     !!    nang      int.  Public    I   # of angles
     !!
     !!    dfrq      Real  Public    I   freq mult. for log freq spacing
     !!    fpk       Real  Public    I   peak freq. [Hz] of initial freq spectrum
     !!    oma       R.A.  Public    I   rel. freq. array (rad*Hz) ----- dim=(nrng)
     !!    frqa      R.A.  Public    I   radian frequencies (Hz) ------- dim=(nrng)
     !!
     !!    ainc      Real  Public    I   angle increment (radians)
     !!    angl      R.A.  Public    I   dir. array (rad) (full circle); dim=(nrng)
     !!    cosan     R.A.  Public    I   cosine angles array ----------- dim=(nang)
     !!    sinan     R.A.  Public    I   sine   angles array ----------- dim=(nang)
     !!    ------------------------------------------------------------------
     !!
     !!    wka       R.A.  Local     I   wavenumbers array [1/m] ------- dim=(nrng)
     !!    cga       R.A.  Local     I   group velocities array [m/s] -- dim=(nrng)
     !!                                  wka & cga arrays are corrsp. to depth 'dep'
     !!    ------------------------------------------------------------------
     !!
     !!    ef2       R.A.  Public    I   2D Energy Density spectrum ef2(theta,f)
     !!                                  = A(theta,k)*2*pi*oma(f)/cga(f) dim=(nrng,nang)
     !!    ef1       R.A.  Public    I   1D Energy Density spectrum ef1(f) dim=(nrng)
     !!    ------------------------------------------------------------------
     !!
     !!    dens1     R.A.  Public    O   large-scale Action Density (k,theta)
     !!                                                     dim=(nrng,nang)
     !!    dens2     R.A.  Public    O   Small-scale Action Density (k,theta)
     !!                                                     dim=(nrng,nang)
     !!    ------------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!    ------------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!     gridsetr  Subr. W3SNL4MD Setup geometric integration grid
     !!    ------------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!
     !!    !/S  Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     !!
     IMPLICIT NONE
     !!
     !!
     !!
     !!    Parameter list
     !!    --------------
     !!op2 Bash; new for optsa2
     integer, intent(in)  :: nrmn, nrmx, nbins
     !!    ------------------------------------------------------------ !!op2
     !!
     integer, intent(in)  :: npk
     real,    intent(in)  :: fpk
     real,    intent(in)  :: wka(nrng),   cga(nrng)
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    Local Parameters & variables
     !!    -----------------------------
     integer              :: irng, iang
     !!
     !!p2
     !!    Bash; Uses of original "psi2(:)", it was very bad (see below)
     !p2   integer              :: n1,  n2,  m,   mm
     !p2   integer              :: nn1, nn2, ii, idif
     !p2   real                 :: q(16)
     !p2   real                 :: emax
     !p2   real                 :: y, qmin, adif
     !!p2---
     !!
     !!p3
     !!    Bash; This is an attempt to replace the original psi2(:)
     !!          with a distr. based on sin()**mm  with 'newmaxang'
     !!          - not good enough (see below)
     !!          !!p4 is an override of !!p3 with mm=4
     !p3   integer              :: n1,  n2,  m,   mm
     !p3   real                 :: q(16)
     !p3   real                 :: y, qmin, adif
     !!    The var. below are needed to find 'newmaxang' used in !p3 & !p4
     integer              :: maxang,    newmaxang
     integer              :: maxangshift
     integer              :: halfangl,  halfangu
     real                 :: ef2maxrow(nang)
     real                 :: ef2shift(nang)
     real                 :: halfmax
     !!p3---
     !!
     !!p4
     !!    Bash; !!p4 is an override of !!p3 with mm=4
     integer              :: n1,  n2
     real                 :: q4
     !!p4---
     !!
     !!eq
     !!    Bash; Use variable equi. range suitable to TSA min condition
     !!          simplifed to one point equi. range nearest to 2.*fp
     !eq   integer              :: neq
     !eq   real                 :: fovfp
     !!eq---
     !!
     integer              :: igam
     real                 :: sum1, fac
     real                 :: beta, gam
     real                 :: fdenp, fr, ratio, z, ddd
     real                 :: sigz               !* sigz  = 0.109
     real                 :: fk(nrng),    fknrm(nrng)
     real                 :: bscl1(nrng), fkscl1(nrng)
     real                 :: psi2(nang)
     real                 :: act2d(nrng,nang)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!ini
     !!    Bash; initialize psi2() array here
     psi2(:)    = 0.0
     !!
     !!    Initialize all the 1d & 2d arrays that are being used
     !!    and especially those that are being returned
     fk(:)      = 0.0
     fknrm(:)   = 0.0
     fkscl1(:)  = 0.0
     bscl1(:)   = 0.0
     act2d(:,:) = 0.0
     dens1(:,:) = 0.0
     dens2(:,:) = 0.0
     !!ini---
     !!    ------------------------------------------------------------------
     !!
     !!
     !!*   Convert 2D Energy Density ef2(f,theta)
     !!         to 2D Polar Action Density act2d(k,theta) Norm. (in k)
     do irng=nrmn,nrmx
       fac   = cga(irng)/(twopi*oma(irng)*wka(irng))
       do iang=1,nang
         act2d(irng,iang) = ef2(irng,iang) * fac
       end do
       !!
       !!*     Convert ef1(f) to fk(k); both are 1d Energy Density
       fk(irng)    = cga(irng)*ef1(irng)/twopi  !* fk(k) energy
       !!
       !!*     Normalize the 1d wavenumber Energy Density fk(k) to give fknrm(k)
       fknrm(irng) = fk(irng)*wka(irng)**2.5    !* fknrm(k) = Norm. fk(k)
     end do
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    Fit parameters to spectrum
     !!    --------------------------
     !!eq
     !eq   sum1 = 0.
     !eq   neq  = 0
     !eq   do 26 irng=nrmn,nrmx
     !eq     fovfp = frqa(irng)/fpk
     !!      Bash; check2 test equilibrium range
     !b      if ( fovfp.ge.1.55.and.fovfp.le.2.45 ) then !* orig   equi range
     !b      if ( fovfp.ge.1.20.and.fovfp.le.2.20 ) then !* wide   equi range
     !b      if ( fovfp.ge.1.90.and.fovfp.le.2.20 ) then !* narrow equi range
     !!      --------------------------------------------------------------!*  <<<<<
     !!      Bash; select variable equi. range suitable to TSA min condition
     !eq     if ( fovfp.ge.(dfrq**(nbins))-0.005 .and.                     &
     !eq          fovfp.le.(dfrq**(nbins))+0.005 ) then  !* narrow equi range  <<<<<
     !!      --------------------------------------------------------------!*  <<<<<
     !eq       sum1 = sum1 + fknrm(irng)
     !eq       neq  = neq + 1
     !eq     endif
     ! 26  end do
     !eq   beta = sum1 / neq
     !eq   gam  = fknrm(npk) / beta
     !!eq---
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    Simplify
     beta = fknrm(npk+nbins)
     gam  = fknrm(npk) / beta
     !!eq---
     !!    ------------------------------------------------------------------
     !!
     do irng=nrmn,nrmx
       fknrm(irng) = fknrm(irng) / beta
     end do
     !!    ==================================================================
     !!
     !!
     !!p2
     !!    Construct Directional Distribution "psi2(:)" - original option
     !!    ------------------------------------------------------------------
     !!
     !!    Solve for Normalizing Coefficient for Integral [1.0/(cos**m)]
     !!    Note: n1, n2 spans half circle (from -pi/2 to +pi/2 going through 0.)
     !p2   n1 = -nang/4 + 1
     !p2   n2 =  nang/4 + 1
     !p2   do m=1,16
     !p2     sum1 = 0.
     !p2     do iang=n1,n2
     !p2       ii = iang
     !p2       if ( iang .lt. 1 ) ii = iang + nang
     !p2       sum1 = sum1 + cosan(ii)**m
     !       end do
     !p2     q(m) = 1./(sum1*ainc)
     !     end do
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    Find peak direction "maxang" in ef2() at "npk" the peak in ef1()
     !!    needed to define the energy spreading factor y=ef2(npk,maxang)/ef1(npk)
     !!    Bash; Note; This original "psi2(:)" was simply very bad because the drift
     !!          in "maxang" location causing the 2D Snl to lose symmetry
     !p2   emax   = 0.
     !p2   maxang = 0
     !p2   do iang=1,nang
     !p2     if ( ef2(npk,iang).gt.emax ) then
     !p2       emax   = ef2(npk,iang)
     !p2       maxang = iang                          !* in [1,nang]
     !p2     endif
     !     end do
     !p2   y  = ef2(npk,maxang)/ef1(npk)              !* Bash; Energy Spread
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    Compare value of peak with q-array for closest fit to cos()**m at peak
     !p2   mm   = 1
     !p2   qmin = abs(q(1)-y)
     !p2   do m=2,16
     !p2     adif = abs(q(m)-y)
     !p2     if ( adif.lt.qmin ) then
     !p2       qmin = adif
     !p2       mm   = m
     !p2     endif
     !     end do
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !p2   nn1 = maxang - nang/4           !* nn1 in [-8, 27], -ve/+ve (incl. 0)
     !p2   nn2 = maxang + nang/4           !* nn2 in [10, 45], all +ve  (no  0)
     !p2   do iang=nn1,nn2              !* Bash; nn1 -> nn2 covers half circle
     !p2     ii = iang                                !* ii always in range [1,nang]
     !p2     if ( ii .lt.    1 ) ii = ii + nang       !* ""
     !p2     if ( ii .gt. nang ) ii = ii - nang       !* ""
     !p2     idif = iabs(maxang-iang) + 1             !* =10,9,..,2,1,2,..,9,10
     !p2     psi2(ii) = q(mm) * cos(angl(idif))**mm   !* Normalized psi2 distr.
     !     end do
     !!p2---
     !!    ==================================================================
     !!
     !!
     !!p3
     !!    Construct New Directional Distribution "psi2(:)"
     !!    In an attempt to replace the original psi2(:) with
     !!    a distribution based on sin()**mm with 'newmaxang' - not good enough
     !!    ------------------------------------------------------------------
     !!
     !!    Solve for Normalizing Coefficient for Integral [1.0/(sin()**m)]
     !!    Note: n1, n2 spans half circle (from 0 to +pi)
     !p3   n1 =  1
     !p3   n2 =  nang/2 + 1
     !p3   do m=1,16
     !p3     sum1 = 0.
     !p3     do iang=n1,n2
     !p3       sum1 = sum1 + sinan(iang)**m
     !       end do
     !p3     q(m) = 1./(sum1*ainc)
     !     end do
     !!p3---
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!p3
     ef2maxrow(:) = ef2(npk,:)
     maxang       = maxloc(ef2maxrow,1)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    Shift the row so that the max is at location of 90 degress
     !!          Negative shift is to the right, Postive to the left
     !!          halfangl - lower angular limit of the half maximum
     !!          halfangu - upper angular limit of the half maximum
     ef2shift(:) = cshift( ef2maxrow(:), (maxang-1-nang/4) )
     halfangu    = nang/4+2
     halfangl    = nang/4
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     halfmax = 0.5 * ef2(npk,maxang)
     do while((ef2shift(halfangu).gt.halfmax).and.(halfangu.lt.nang/2))
       halfangu = halfangu + 1
     enddo
     do while((ef2shift(halfangl).gt.halfmax).and.(halfangl.gt.1))
       halfangl = halfangl - 1
     enddo
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    Convert angles indices with respect to peak
     !!       e.g. halfangl should go to halfangl - (nang/4+1)
     !!            halfangu should go to halfangu - (nang/4+1)
     halfangl = halfangl - (nang/4+1)
     halfangu = halfangu - (nang/4+1)
     !!
     !!    Now average the positions, round to nearest integer.
     !!    -ve result means the centre is one greater than it should be.
     maxangshift = nint( 0.5 * (halfangl + halfangu) )
     newmaxang   = maxang + maxangshift
     if (newmaxang .lt.    1) newmaxang = newmaxang + nang
     if (newmaxang .gt. nang) newmaxang = newmaxang - nang
     !!p3---
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!p3
     !!    Bash; need this section if you want to try sin()**mm with 'newmaxang'
     !p3   y = ef2(npk,newmaxang) / ef1(npk)          !* New Energy Spread
     !!
     !!    Compare value of peak with q-array for closest fit to sin()**m at peak
     !!    This !p3 section is needed for use with  sin()**mm
     !!    Bash; Note; This new "psi2(:)" although better than original "psi2(:)"
     !!          it was still not good enough:  the 2D Energy was OK but
     !!          the 2D Snl, now with better symmetry, didn't always have the side lobes.
     !p3   mm   = 1
     !p3   qmin = abs(q(1)-y)
     !p3   do m=2,16
     !p3     adif = abs(q(m)-y)
     !p3     if ( adif.lt.qmin ) then
     !p3       qmin = adif
     !p3       mm   = m
     !p3     endif
     !     end do
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    Final step, use 'mm' for sin()**mm
     !p3   psi2(n1:n2) = (sinan(n1:n2))**mm           !* Un-norm. psi2 distr.
     !p3   psi2(n1:n2) = q(mm) * psi2(n1:n2)          !* Norm. psi2 distr.
     !!    Rotate peak to correct angle
     !p3   psi2(:) = CSHIFT( psi2(:), newmaxang-1+nang/4 )
     !!p3---
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !!p4
     !!    !!p4 is an override of !!p3 with mm=4, so go straight to Final step
     !!    Note; all you need from !!p3 is the  "newmaxang"
     !!    So it's a sin()**4 distr. shifted to "newmaxang" - worked very well
     !!    ------------------------------------------------------------------
     !!
     !!    Solve for Normalizing Coefficient for Integral [1.0/(sin()**4)]
     !!    Note: n1, n2 spans half circle (from 0 to +pi)
     n1 =  1
     n2 =  nang/2 + 1
     !p4   sum1 = 0.
     !p4   do iang=n1,n2
     !p4     sum1 = sum1 + sinan(iang)**4
     !     end do
     !p4   q4 = 1.0/(sum1*ainc)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    Change the angles that aren't zero  (0 deg to +180 deg)
     psi2(n1:n2) = (sinan(n1:n2))**4                 !* Un-norm. psi2 distr.
     q4          = 1.0/(sum(psi2(n1:n2))*ainc)
     psi2(n1:n2) = q4 * psi2(n1:n2)                  !* Norm. psi2 distr.
     !!
     !!    Rotate peak to correct angle
     psi2(:) = cshift( psi2(:), newmaxang-1+nang/4 )
     !!p4---
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !!
     !!    Estimate parametric spectrum and deviation from parametric spectrum
     !!    ------------------------------------------------------------------
     igam  = (gam-0.4)*10 + 0.5
     sigz  = 0.109
     gam   = igam/10. + 0.4
     fdenp = gam * beta / wka(npk)**2.5
     !!
     !!
     do irng=nrmn,nrmx
       fr = frqa(irng) / fpk
       if ( fr.le.1.0001 ) then
         if ( fr.ge.0.85 ) then
           ratio = 1.-(1.-fr)*0.7/0.15
         else
           ratio = 0.3*exp(-17.3*(0.85-fr))
         endif
         fkscl1(irng) = fdenp*ratio
         bscl1(irng)  = fkscl1(irng)/oma(irng)
       else
         z = 0.5*((fr-1.)/sigz)**1.2
         if ( z.gt.6. ) z = 6.
         ratio = 1.+exp(-z)*(gam-1.)
         fkscl1(irng) = beta*ratio/wka(irng)**2.5
         bscl1(irng)  = fkscl1(irng)/oma(irng)
       endif
       !!
       do iang=1,nang
         ddd = bscl1(irng) * psi2(iang) / wka(irng)  !* large-scale
         dens1(irng,iang) = ddd                      !* large-scale
         dens2(irng,iang) = act2d(irng,iang) - ddd   !* small-scale
       end do
     end do ! do irng=nrmn,nrmx
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     RETURN
     !!

References ainc, dens1, dens2, ef1, ef2, frqa, oma, sinan, and twopi.

Referenced by w3snl4().

◆ shlocr()

subroutine w3snl4md::shlocr	(	real, intent(in)	dep,
		real, intent(in)	wk1x,
		real, intent(in)	wk1y,
		real, intent(in)	wk3x,
		real, intent(in)	wk3y
	)

N/A.

  With wavenumber vector n identified by (wknx,wkny), its magnitude
  given by wkn = sqrt(wknx**2+wkny**2) and its associated radian
  frequency given by sign = sqrt[g*wkn*tanh(wkn*dep)], where g is
  gravitational acceleration and dep is water depth, the four-wave
  resonance condition is satisfied along a locus of pts defined by

  [1]  (wk1x,wk1y) + (wk2x,wk2y) - (wk3x,wk3y) - (wk4x,wk4y) = 0

  [2]  sig1 + sig2 - sig3 - sig4 = 0

  Because of the influence of depth, it is convenient to define new
  vectors (wnx,wny) = (wknx*dep,wkny*dep) with magnitudes wn =
  sqrt(wnx**2+wny**2) = wkn*dep such that a dimensionless frequency
  is sign*sqrt(dep/g) = sqrt[wkn*dep*tanh(wkn*dep)]
                      = sqrt[wn*tanh(wn)].
  With these definitions and vectors (wk1x,wk1y) and (wk3x,wk3y)
  given as input, we can write (with some rearrangement
  of [1] and [2])

  [3]  w3x - w1x = px = w2x - w4x
  [4]  w3y - w1y = py = w2y - w4y
  [5]  sqrt[w3*tanh(w3)] - sqrt[w1*tanh(w1)] = q
       = sqrt[w2*tanh(w2)] - sqrt[w4*tanh(w4)]

  With dimensionless vector (px,py) = (w3x-w1x,w3y-w1y) [magnitude
  p = sqrt(px**2+py**2), direction atan2(py,px)] and dimensionless
  frequency difference q = sqrt(w3*tanh(w3)] - sqrt(w1*tanh(w1)]
  defined by input parameters, we see from [3] and [4] that
  (w4x,w4y) = (w2x-px,w2y-py) [magnitude w4 = sqrt((w2x-px)**2 +
  (w2y-py)**2)] and thus from [5] we must basically find elements
  w2x and w2y that satisfy

  [6] sqrt[sqrt(w2x**2+w2y**2)*tanh(sqrt(w2x**2+w2y**2))] -
      sqrt[sqrt((w2x-px)**2+(w2y-py)**2) *
           tanh(sqrt((w2x-px)**2+(w2y-py)**2] = q

  The locus curve defined by the set of pts (w2x,w2y) crosses the
  p-vector axis at two points; one with magnitude w2=rmin*p with
  0.5 < rmin < 1.0 and one with magnitude w2=rmax*p with rmax > 1.
  We first isolate rmin, rmax using various iterative algorithms,
  and then find locus pts that are not on the p-vector axis with
  another iterative scheme.  At the end, we un-normalize the w2
  and w4 vectors to find the wk2 and wk4 vectors.

Parameters

[in]	dep
[in]	wk1x
[in]	wk1y
[in]	wk3x
[in]	wk3y

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 2306 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!    ------------------------------------------------------------------
     !!
     !!    it returns: wk2x, wk2y, wk4x, wk4y & ds all dim=(npts)
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !! 1. Purpose :
     !!
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    With wavenumber vector n identified by (wknx,wkny), its magnitude!
     !!    given by wkn = sqrt(wknx**2+wkny**2) and its associated radian   !
     !!    frequency given by sign = sqrt[g*wkn*tanh(wkn*dep)], where g is  !
     !!    gravitational acceleration and dep is water depth, the four-wave !
     !!    resonance condition is satisfied along a locus of pts defined by !
     !!                                                                     !
     !!    [1]  (wk1x,wk1y) + (wk2x,wk2y) - (wk3x,wk3y) - (wk4x,wk4y) = 0   !
     !!                                                                     !
     !!    [2]  sig1 + sig2 - sig3 - sig4 = 0                               !
     !!                                                                     !
     !!    Because of the influence of depth, it is convenient to define new!
     !!    vectors (wnx,wny) = (wknx*dep,wkny*dep) with magnitudes wn =     !
     !!    sqrt(wnx**2+wny**2) = wkn*dep such that a dimensionless frequency!
     !!    is sign*sqrt(dep/g) = sqrt[wkn*dep*tanh(wkn*dep)]                !
     !!                        = sqrt[wn*tanh(wn)].                         !
     !!    With these definitions and vectors (wk1x,wk1y) and (wk3x,wk3y)   !
     !!    given as input, we can write (with some rearrangement            !
     !!    of [1] and [2])                                                  !
     !!                                                                     !
     !!    [3]  w3x - w1x = px = w2x - w4x                                  !
     !!    [4]  w3y - w1y = py = w2y - w4y                                  !
     !!    [5]  sqrt[w3*tanh(w3)] - sqrt[w1*tanh(w1)] = q                   !
     !!         = sqrt[w2*tanh(w2)] - sqrt[w4*tanh(w4)]                     !
     !!                                                                     !
     !!    With dimensionless vector (px,py) = (w3x-w1x,w3y-w1y) [magnitude !
     !!    p = sqrt(px**2+py**2), direction atan2(py,px)] and dimensionless !
     !!    frequency difference q = sqrt(w3*tanh(w3)] - sqrt(w1*tanh(w1)]   !
     !!    defined by input parameters, we see from [3] and [4] that        !
     !!    (w4x,w4y) = (w2x-px,w2y-py) [magnitude w4 = sqrt((w2x-px)**2 +   !
     !!    (w2y-py)**2)] and thus from [5] we must basically find elements  !
     !!    w2x and w2y that satisfy                                         !
     !!                                                                     !
     !!    [6] sqrt[sqrt(w2x**2+w2y**2)*tanh(sqrt(w2x**2+w2y**2))] -        !
     !!        sqrt[sqrt((w2x-px)**2+(w2y-py)**2) *                         !
     !!             tanh(sqrt((w2x-px)**2+(w2y-py)**2] = q                  !
     !!                                                                     !
     !!    The locus curve defined by the set of pts (w2x,w2y) crosses the  !
     !!    p-vector axis at two points; one with magnitude w2=rmin*p with   !
     !!    0.5 < rmin < 1.0 and one with magnitude w2=rmax*p with rmax > 1. !
     !!    We first isolate rmin, rmax using various iterative algorithms,  !
     !!    and then find locus pts that are not on the p-vector axis with   !
     !!    another iterative scheme.  At the end, we un-normalize the w2    !
     !!    and w4 vectors to find the wk2 and wk4 vectors.                  !
     !!    -----------------------------------------------------------------#
     !!
     !! 2. Method :
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    npts      int.  Public    I   # of points on the locus
     !!    np2p1     int.  Public    I   = npts/2 + 1
     !!    ----------------------------------------------------------------
     !!
     !!    *** arrays wk2x,wk2y, wk4x,wk4y & ds are related to locus solutioni
     !!        for given vectors K1 & k3        all have dim=(npts)
     !!    wk2x      R.A.  Public    O   x_cmp of vector k2 solution on the locus
     !!    wk2y      R.A.  Public    O   y_cmp of vector k2 solution on the locus
     !!    wk4x      R.A.  Public    O   x_cmp of vector k4 solution on the locus
     !!    wk4y      R.A.  Public    O   y_cmp of vector k4 solution on the locus
     !!    ds        R.A.  Public    O   element length along the locus
     !!                                  (npts*ds circles the whole locus)
     !!    ----------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!    ----------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!     gridsetr  Subr. W3SNL4MD Setup geometric integration grid
     !!    ----------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!
     !!    !/S  Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     USE w3servmd, ONLY: extcde
     USE w3odatmd, ONLY: ndse
     !!
     !!
     IMPLICIT NONE
     !!
     !!    Parameter list
     !!    --------------
     real,    intent(in)  :: dep
     real,    intent(in)  :: wk1x,wk1y, wk3x,wk3y
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    Local Parameters & variables
     !!    -----------------------------
     integer              :: n, np, nnp, nplace
     integer              :: ierr_gr
     !!
     real                 :: p,  px,  py,  q,    qrtp,  qsqp,        &
          dr, dth, thp, dphi, cphi,               &
          w1, w1x, w1y, wk1,  w3, w3x, w3y, wk3,  &
          rold, rold1, rold2, rnew, rnew1, rnew2, &
          pxod, pyod,  zpod,                      &
          t, t1, t2, t3, tm, tp, ds1, ds2,        &
          rmin, rmax, rcenter, rradius
     !!
     double precision     :: dbt3,dbt4,dbt5,dbt6, dbz, dbp, dbqrtp,  &
          cdthold, cdthnew, wate1, wate2
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!ini
     !!    initial all returned arrays before they are computed
     !!    -----------------------------------------------------
     wk2x(:) = 0.0
     wk2y(:) = 0.0
     wk4x(:) = 0.0
     wk4y(:) = 0.0
     ds(:)   = 0.0
     !!ini---
     !!    ------------------------------------------------------------------
     !!
     !!
     !b    wk1   = sqrt(wk1x**2 + wk1y**2)         !* k1=wk1x since wk1y=0.0
     wk1   = wk1x
     !b    wk3   = sqrt(wk3x**2 + wk3y**2)
     wk3   = sqrt(wk3x*wk3x + wk3y*wk3y)
     !!
     w1    = wk1  * dep
     w1x   = wk1x * dep
     !b    w1y   = wk1y * dep
     w1y   = wk1y                            !* wk1y=0.0
     !!
     w3    = wk3  * dep
     w3x   = wk3x * dep
     w3y   = wk3y * dep
     !!
     px    = w3x - w1x
     py    = w3y - w1y
     !b    p     = sqrt(px**2 + py**2)             !* argument never = 0.0
     p     = sqrt(px*px + py*py)             !* argument never = 0.0
  
     thp   = atan2(py,px)
     q     = sqrt(w3*tanh(w3)) - sqrt(w1*tanh(w1))
     qrtp  = q / sqrt(p)
     qsqp  = qrtp*qrtp
     !!
     !!
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    for (w2x,w2y) = rmin*(px,py) (locus crossing the p-vector axis   !
     !!    nearest the origin), we have (w4x,w4y) = (w2x,w2y) - (px,py)     !
     !!    = rmin*(px,py) - (px,py) = (rmin - 1)*(px,py); note that because !
     !!    rmin < 1, the length of (w4x,w4y) is w4 = (1 - rmin)*p; then [6] !
     !!    takes the simpler form                                           !
     !!                                                                     !
     !!    [7] sqrt(rmin*p*tanh(rmin*p)) - sqrt[(1-rmin)*p*tanh((1-rmin)*p)]!
     !!        = q                                                          !
     !!                                                                     !
     !!    assuming the tanh() functions are slowly varying and can be      !
     !!    treated as separate entities, [7] can be written as a quadratic  !
     !!    in sqrt(rmin), i.e.,                                             !
     !!                                                                     !
     !!               2*qrtp*sqrt(rmin*p)                                   !
     !!    [8] rmin - ------------------- sqrt(rmin) +                      !
     !!                        T                                            !
     !!                                                                     !
     !!                                   qsqp-p*tanh((1-rmin)*p)           !
     !!                                   -----------------------  =  0,    !
     !!                                              T                      !
     !!                                                                     !
     !!    where   T    = tanh(rmin*p) + tanh((1-rmin)*p),                  !
     !!            qrtp = q/sqrt(p)  and qsqp = qrtp**2                     !
     !!                                                                     !
     !!    the square of the most positive root of [8] can (with some       !
     !!    algebra) be written                                              !
     !!                                                                     !
     !!    [9] rmin =                                                       !
     !!   (1/T**2)*{qsqp*[tanh(rmin*p)-tanh((1-rmin)*p)]+T*tanh((1-rmin)*p) !
     !!           +2*qrtp*sqrt[tanh(rmin*p)*tanh((1-rmin)*p)]*sqrt(T-qsqp)} !
     !!                                                                     !
     !!    setting rnew=rmin on the LHS and rold=rmin on the RHS in all     !
     !!    instances allows the creation of an iterative algorithm for rmin;!
     !!    convergence can be slow in general, so a coarse search for the   !
     !!    crossing of [9] with the rnew=rold line is conducted first, then !
     !!    a weighted iterative replacement loop is executed until the      !
     !!    desired accuracy is achieved;                                    !
     !!    Note that if p is sufficiently large, all tanh() -> 1            !
     !!    and [9] becomes the analytic expression                          !
     !!            rmin = 0.5*[1 + qrtp*sqrt(2-qsqp)]                       !
     !!                                                                     !
     !!    following is the coarse search using rold1 = 0.5,0.1,1.0,        !
     !!                                         rold2 = rold1 + 0.1:        !
     !!                                                                     !
     !!    -----------------------------------------------------------------#
     !!
     !!
     ierr_gr = 0
     !!
     rold1 = 0.5
     tp    = tanh(rold1 * p)
     tm    = tanh((1.-rold1) * p)
     t     = tp + tm
     t1    = qsqp * (tp-tm)
     t2    = t  * tm
     t3    = 2. * qrtp * sqrt(tp*tm) * sqrt(t-qsqp)
     rnew1 = (t1 + t2 + t3) / (t*t)
     !!
     !!
     do n=1,4
       rold2 = rold1 + 0.1
       tp    = tanh(rold2 * p)
       tm    = tanh((1.-rold2) * p)
       t     = tp + tm
       t1    = qsqp * (tp-tm)
       t2    = t  * tm
       t3    = 2. * qrtp * sqrt(tp*tm) * sqrt(t-qsqp)
       rnew2 = (t1 + t2 + t3) / (t*t)
       if ( rnew2 .lt. rold2 ) then
         rold = (rold2*rnew1-rold1*rnew2)/(rold2-rold1-rnew2+rnew1)
         go to 11
       end if
       rold1 = rold2
       rnew1 = rnew2
     end do  ! do n=1,4
     rold = 0.9                       !* default if not otherwise found
 11  continue
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    iterative replacement search for rmin
     do n=1,50
       tp   = tanh(rold * p)
       tm   = tanh((1.-rold) * p)
       t    = tp + tm
       t1   = qsqp * (tp-tm)
       t2   = t  * tm
       t3   = 2. * qrtp*sqrt(tp*tm)*sqrt(t-qsqp)
       rnew = (t1 + t2 + t3) / (t*t)
       if ( abs(rnew-rold) .lt. 0.00001 ) then
         rmin = rnew
         go to 21
       end if
       rold = 0.5 * (rold + rnew)
     end do
     ierr_gr = ierr_gr + 1  !* set 1's flag in ierr_gr if no convergence
     rmin = rnew
 21  continue
     !!    ------------------------------------------------------------------
     !!
     !!    set (dimensional) wavenumber components for this point on locus
     wk2x(1) =  rmin * px / dep
     wk2y(1) =  rmin * py / dep
     wk4x(1) = (rmin-1.) * px / dep
     wk4y(1) = (rmin-1.) * py / dep
     !!
     !!
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    for (w2x,w2y) = rmax*(px,py) (locus crossing the p-vector axis   !
     !!    farthest from the origin), we have (w4x,w4y)=(w2x,w2y) - (px,py) !
     !!    = rmax*(px,py) - (px,py) = (rmax - 1)*(px,py);                   !
     !!    here, because rmax > 1, the length of (w4x,w4y) is               !
     !!    w4 = (rmax - 1)*p; then [6] takes the form                       !
     !!                                                                     !
     !!    [10] sqrt(rmax*p*tanh(rmax*p))-sqrt[(rmax-1)*p*tanh((rmax-1)*p)] !
     !!         = q                                                         !
     !!                                                                     !
     !!    rearranging terms, squaring both sides and again rearranging     !
     !!    terms yields                                                     !
     !!                                                                     !
     !!    [11] rmax*p*[tanh(rmax*p) - tanh((rmax-1)*p)]                    !
     !!         = 2*q*sqrt(tanh(rmax*p))*sqrt(rmax*p) +                     !
     !!           q**2 + p*tanh((rmax-1)*p)                                 !
     !!                                                                     !
     !!    because the difference of the two tanh()'s on the LHS tend to    !
     !!    make the whole term small, we solve for rmax from the rapidly    !
     !!    varying part of the first term on the RHS, i.e.,                 !
     !!                                                                     !
     !!                         [tanh((rmax-1)*p)+qsqp + rmax*T]**2         !
     !!    [12]          rmax = ----------------------------------- ,       !
     !!                                 4*qsqp*tanh(rmax*p)                 !
     !!                                                                     !
     !!    where, in this algorithm, T = tanh(rmax*p) - tanh((rmax-1)*p);   !
     !!    as for rmin in [9], setting rnew=rmax on the LHS and rold=rmax   !
     !!    in all instances on the RHS allows the formation of an iterative !
     !!    algorithm; initially, we only know rmax > 1 so we do a coarse    !
     !!    search in the 10's place out to some reasonably big number to    !
     !!    try to find the place where [12] crosses the rnew = rold line    !
     !!    (if this fails, we set an error flag); in refining the estimate, !
     !!    it appears that [12] can get a little squirrely, so we do a      !
     !!    brute force successive decimation search to nplace decimal places!
     !!    to home in on the answer; note that if p is big enough for the   !
     !!    tanh()'s to reach unity, [12] becomes exact and                  !
     !!    rmax = [(1 + qsqp)**2]/(4*qsqp)                                  !
     !!    following is the coarse search with rold = 1,10,2001:            !
     !!                                                                     !
     !!    -----------------------------------------------------------------#
     !!
     rold = 1.0
     do n=1,200
       rold = rold + 10.
       tp   = tanh(rold * p)
       tm   = tanh((rold-1.) * p)
       t    = tp - tm
       t1   = tm + qsqp
       t2   = 4. * tp * qsqp
       rnew = ((t1+rold*t)**2) / t2
       if ( rnew .lt. rold ) then
         rold = rold - 10.
         go to 31
       end if
     end do
     ierr_gr = ierr_gr + 10    !* set 10's place in ierr_gr if no sol'n
 31  continue
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    successive decimation search to refine rmax
     dr = 10.
     do nplace=1,6
       dr = dr/10.
       do n=1,10
         rold = rold + dr
         tp   = tanh(rold * p)
         tm   = tanh((rold-1.) * p)
         t    = tp - tm
         t1   = tm + qsqp
         t2   = 4. * tp * qsqp
         rnew = ((t1+rold*t)**2) / t2
         if ( rnew .lt. rold ) then
           rold = rold - dr
           go to 51
         end if
       end do
 51    continue
     end do
     !!
     rmax = rold
     !!
     !!    set (dimensional) wavenumber components for this locus point
     !!
     wk2x(np2p1) =  rmax * px / dep
     wk2y(np2p1) =  rmax * py / dep
     wk4x(np2p1) = (rmax-1.) * px / dep
     wk4y(np2p1) = (rmax-1.) * py / dep
     !!
     !!
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    search for cos(dth) for off-p-vector solutions; use a circle     !
     !!    centered on the p-vector axis at a distance                      !
     !!    rcenter = 0.5*(rmax+rmin)  from the  origin with a               !
     !!    radius  = 0.5*(rmax-rmin); radii from the center of the circle   !
     !!    at successive angle increments np*dphi intersect the circle at   !
     !!    distances r*p from the origin of the p vector such that          !
     !!                                                                     !
     !!    [13]  r**2 = rradius**2 + rcenter**2 -                           !
     !!                 2*rcenter*rradius*cos(np*dphi)                      !
     !!                                                                     !
     !!    and makes an angle dth with the p vector satisfying              !
     !!                                                                     !
     !!    [14]   cdth = cos(dth) = (rcenter/r) - (rradius/r)*cos(np*dphi)  !
     !!                                                                     !
     !!    we then rotate this vector, holding its length=r*p constant and  !
     !!    successively estimating cdth (using the above equation as an     !
     !!    initial guess) until it intersects the locus curve; some algebra !
     !!    yields the estimation equation as                                !
     !!                                                                     !
     !!                  r**2 + 1      [sqrt(r*tanh(rp)) - q/sqrt(p)]**4    !
     !!    [15] cdthnew= ------- - ---------------------------------------- !
     !!                    2*r     2*r*[tanh(p*sqrt(r**2-2*r*cdthold+1))]**2!
     !!                                                                     !
     !!    we use a weighted new estimate of cdthold with the weights based !
     !!    on the argument of the tanh() function in the denominator        !
     !!    (if the argument is big, tanh() -> 1 and cdthnew is found in one !
     !!    pass; for small arguments, convergence is faster with equal      !
     !!    weighting of old and new estimates; all this is empirical to try !
     !!    to increase speed); double precision is used to gain enough      !
     !!    accuracy when the arccos is taken; note that if p is big enough  !
     !!    for all tanh()'s -> 1, [15] is exact and                         !
     !!    cdthnew = cdth = [r**2 + 1 - (sqrt(r) - qrtp)**4]/(2*r)          !
     !!                                                                     !
     !!    -----------------------------------------------------------------#
     !!
     !!
     rcenter = 0.5 * (rmax + rmin)
     rradius = 0.5 * (rmax - rmin)
     t1      = rradius**2 + rcenter**2
     t2      = 2. * rradius * rcenter
     dphi    = 6.283185308 / float(npts)
     pxod    = px / dep
     pyod    = py / dep
     !!
     dbp     = dble(p)
     dbqrtp  = dble(qrtp)
     !!
     !!
     do np=2,npts/2                        !* np = 2 --> 15
       !!
       cphi    = cos(float(np-1)*dphi)
       dbz     = dsqrt(dble(t1-t2*cphi))
       cdthold = dble(rcenter-rradius*cphi) / dbz
       dbt3    = (dbz*dbz) + 1.d0
       dbt4    = dbt3 / (2.d0*dbz)
       dbt5    = ((dsqrt(dbz*dtanh(dbz*dbp))-dbqrtp)**4)/(2.d0*dbz)
       dbt6    = dbp * dsqrt(dbt3-2.d0*dbz*cdthold)
       !!
       if ( dbt6 .gt. 0.55d0 ) then
         wate1 = dtanh(dbt6)
         wate2 = 1.d0 - wate1
       else
         wate1 = 0.5d0
         wate2 = 0.5d0
       end if
       !!
       do  n=1,25
         cdthnew = dbt4 - dbt5 / ((dtanh(dbt6))**2)
         if ( dabs(cdthnew-cdthold) .lt. 0.0000001d0 ) go to 71
         cdthold = wate1 * cdthnew + wate2 * cdthold
         dbt6    = dbp * dsqrt(dbt3-2.d0*dbz*cdthold)
       end do
       ierr_gr = ierr_gr + 100   !* add to 100's place for every failure
 71    continue
       !!
       dth  = sngl(dacos(cdthnew))
       zpod = sngl(dbz) * p / dep
       !!
       wk2x(np) = zpod * cos(thp+dth)
       wk2y(np) = zpod * sin(thp+dth)
       wk4x(np) = wk2x(np) - pxod
       wk4y(np) = wk2y(np) - pyod
       !!
       nnp = npts-np+2                        !* for npts = 30
       !!                                             !* nnp = 30 --> 17
       !!
       wk2x(nnp) = zpod * cos(thp-dth)
       wk2y(nnp) = zpod * sin(thp-dth)
       wk4x(nnp) = wk2x(nnp) - pxod
       wk4y(nnp) = wk2y(nnp) - pyod
       !!
     end do ! do np=2,npts/2
     !!    ------------------------------------------------------------------
     !!
     !!
     if ( ierr_gr .ne. 0 ) then
       write ( ndse,1000 ) ierr_gr
       CALL extcde ( 60 )
     endif
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    set arc length ds as the sum of half the segment lengths on either
     !!    side of a given point
     !!
     ds1   = sqrt((wk2x(2)-wk2x(1))**2+(wk2y(2)-wk2y(1))**2)
     ds(1) = ds1
     do np=3,npts/2+1
       ds2   = sqrt((wk2x(np)-wk2x(np-1))**2+(wk2y(np)-wk2y(np-1))**2)
       ds(np-1)      = 0.5*(ds1+ds2)
       ds(npts-np+3) = ds(np-1)
       ds1           = ds2
     end do
     ds(npts/2+1)    = ds2
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     RETURN
     !!
 1000 format ( ' W3SNL4 Error : In shlocr. Error from gridset ',i10)
     !!

References ds, w3servmd::extcde(), w3odatmd::ndse, np2p1, npts, wk2x, wk2y, wk4x, and wk4y.

Referenced by gridsetr().

◆ shloxr()

subroutine w3snl4md::shloxr	(	real, intent(in)	dep,
		real, intent(in)	wk1x,
		real, intent(in)	wk1y,
		real, intent(in)	wk3x,
		real, intent(in)	wk3y
	)

General locus solution for input vectors (wk1x,wk1y).

    General locus solution for input vectors (wk1x,wk1y) and
    (wk3x,wk3y) of the same magnitude but NOT in the same direction
    (or singularness will occur), output vectors (wk2x,wk2y) and
    (wk4x,wk4y), and element length ds along locus curve:

    With wavenumber vector n identified by (wknx,wkny), its magnitude
    given by wkn = sqrt(wknx**2+wkny**2) and its associated radian
    frequency given by sign = sqrt[g*wkn*tanh(wkn*dep)], where g is
    gravitational acceleration and dep is water depth, the four-wave
    resonance condition is satisfied along a locus of pts defined by

    [1]  (wk1x,wk1y) + (wk2x,wk2y) - (wk3x,wk3y) - (wk4x,wk4y) = 0

    [2]  sig1 + sig2 - sig3 - sig4 = 0

    In the case where k1 [= sqrt(wk1x**2+wk1y**2)] is equal to k3
    [= sqrt(wk3x**2+wk3y**2)], we have by dispersion,

    [3]  sig1 = sqrt[g*k1*tanh(k1*h)] = sqrt[g*k3*tanh(k3*h)] = sig3

    so sig1 - sig3 = 0 and [2] becomes sig2 = sig4, where, again by
    dispersion,

    [4]  sig2 = sqrt(g*k2*tanh(k2*h)] = sig4 = sqrt(g*k4*tanh(k4*h)]

    and consequently k2 = k4.  This simplifies the locus solution
    considerably, and it can be shown that the (wk2x,wk2y) locus is
    along the perpendicular bisector of the (px,py) vector given by

    [5]  (px,py) = (wk3x-wk1x,wk3y-wk1y)

    and thereby from [1]
    [6]  (wk4x,wk4y) = (wk2x,wk2y) - (px,py)

    We note that these loci are independent of depth, although depth
    is used to set the length of the locus line by requiring that its
    range on either side of the p vector correspond to a wave with a
    wkx freq four times that of the k1 vector (the locus line is made
    up of npts segments of length ds; the outer edges of the terminal
    segments satisfy the length constraint; vectors k2 and k4 extend
    to segment centers and will sufficiently approximate the length
    constraint).  As compared to srshlocr.f, we can do all
    calculations here in dimensional space.

Parameters

[in]	dep
[in]	wk1x
[in]	wk1y
[in]	wk3x
[in]	wk3y

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 2011 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!    ------------------------------------------------------------------
     !!
     !!    it returns: wk2x, wk2y, wk4x, wk4y & ds all dim=(npts)
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !! 1. Purpose :
     !!
     !!    -------------------------------------------------------------------------#
     !!                                                                             !
     !!    General locus solution for input vectors (wk1x,wk1y) and                 !
     !!    (wk3x,wk3y) of the same magnitude but NOT in the same direction          !
     !!    (or singularness will occur), output vectors (wk2x,wk2y) and             !
     !!    (wk4x,wk4y), and element length ds along locus curve:                    !
     !!                                                                             !
     !!    With wavenumber vector n identified by (wknx,wkny), its magnitude        !
     !!    given by wkn = sqrt(wknx**2+wkny**2) and its associated radian           !
     !!    frequency given by sign = sqrt[g*wkn*tanh(wkn*dep)], where g is          !
     !!    gravitational acceleration and dep is water depth, the four-wave         !
     !!    resonance condition is satisfied along a locus of pts defined by         !
     !!                                                                             !
     !!    [1]  (wk1x,wk1y) + (wk2x,wk2y) - (wk3x,wk3y) - (wk4x,wk4y) = 0           !
     !!                                                                             !
     !!    [2]  sig1 + sig2 - sig3 - sig4 = 0                                       !
     !!                                                                             !
     !!    In the case where k1 [= sqrt(wk1x**2+wk1y**2)] is equal to k3            !
     !!    [= sqrt(wk3x**2+wk3y**2)], we have by dispersion,                        !
     !!                                                                             !
     !!    [3]  sig1 = sqrt[g*k1*tanh(k1*h)] = sqrt[g*k3*tanh(k3*h)] = sig3         !
     !!                                                                             !
     !!    so sig1 - sig3 = 0 and [2] becomes sig2 = sig4, where, again by          !
     !!    dispersion,                                                              !
     !!                                                                             !
     !!    [4]  sig2 = sqrt(g*k2*tanh(k2*h)] = sig4 = sqrt(g*k4*tanh(k4*h)]         !
     !!                                                                             !
     !!    and consequently k2 = k4.  This simplifies the locus solution            !
     !!    considerably, and it can be shown that the (wk2x,wk2y) locus is          !
     !!    along the perpendicular bisector of the (px,py) vector given by          !
     !!                                                                             !
     !!    [5]  (px,py) = (wk3x-wk1x,wk3y-wk1y)                                     !
     !!                                                                             !
     !!    and thereby from [1]                                                     !
     !!    [6]  (wk4x,wk4y) = (wk2x,wk2y) - (px,py)                                 !
     !!                                                                             !
     !!    We note that these loci are independent of depth, although depth         !
     !!    is used to set the length of the locus line by requiring that its        !
     !!    range on either side of the p vector correspond to a wave with a         !
     !!    wkx freq four times that of the k1 vector (the locus line is made        !
     !!    up of npts segments of length ds; the outer edges of the terminal        !
     !!    segments satisfy the length constraint; vectors k2 and k4 extend         !
     !!    to segment centers and will sufficiently approximate the length          !
     !!    constraint).  As compared to srshlocr.f, we can do all                   !
     !!    calculations here in dimensional space.                                  !
     !!    -------------------------------------------------------------------------#
     !!
     !! 2. Method :
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    npts      int.  Public    I   # of points on the locus
     !!    np2p1     int.  Public    I   = npts/2 + 1
     !!    ----------------------------------------------------------------
     !!
     !!    *** arrays wk2x,wk2y, wk4x,wk4y & ds are related to locus solutioni
     !!        for given vectors K1 & k3        all have dim=(npts)
     !!    wk2x      R.A.  Public    O   x_cmp of vector k2 solution on the locus
     !!    wk2y      R.A.  Public    O   y_cmp of vector k2 solution on the locus
     !!    wk4x      R.A.  Public    O   x_cmp of vector k4 solution on the locus
     !!    wk4y      R.A.  Public    O   y_cmp of vector k4 solution on the locus
     !!    ds        R.A.  Public    O   element length along the locus
     !!                                  (npts*ds circles the whole locus)
     !!    ----------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!    ----------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!     gridsetr  Subr. W3SNL4MD Setup geometric integration grid
     !!    ----------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!
     !!    !/S  Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     !!wvn
     !wvn  USE W3DISPMD, ONLY: WAVNU2
     !!wvn---
     !!
     IMPLICIT NONE
     !!
     !!    Parameter list
     !!    --------------
     real,    intent(in)  :: dep
     real,    intent(in)  :: wk1x,wk1y, wk3x,wk3y
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    Local Parameters & variables
     !!    -----------------------------
     integer              :: m, n
     real                 :: g
     real                 :: wk1, f1, fx
     real                 :: wkx, db, px, py, p, thp, halfp
     real                 :: a, b, dth, a_halfp
     !a    real                 :: wkfnc      !* real function or use WAVNU2
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!    initial constants
     !!    ------------------
     g     = 9.806                            !* set = GRAV as in CONSTANTS
     !!
     !!ini
     !!    initial all returned arrays before they are computed
     !!    -----------------------------------------------------
     wk2x(:) = 0.0
     wk2y(:) = 0.0
     wk4x(:) = 0.0
     wk4y(:) = 0.0
     ds(:)   = 0.0
     !!ini---
     !!    ------------------------------------------------------------------
     !!
     !!wkx Bash; Try use wkx = wka(nrng) instead of  wkx = wkfnc(4.*f1,dep)
     !b    wk1   = sqrt(wk1x**2+wk1y**2)           !* k1=wk1x since wk1y=0.0
     wk1   = wk1x
     f1    = sqrt(g*wk1*tanh(wk1*dep))/twopi !* f1=f(k1,dep)
     fx    = 4. * f1                         !* fx = 4*f1
     !!
     !!wvn Bash; Try use subr WAVNU2 instead of function wkfnc
     !wvn  call WAVNU2(twopi*fx,dep,wkx,cgx)       !* => wkx=k(4*f1,dep) & cgx=Cg(4*f1,dep)
     wkx   = wkfnc(fx,dep)                   !* +> wkx=k(4*f1,dep)  (cgx not needed)
     !!wvn---
     !!    ------------------------------------------------------------------
     !!
     db    = wkx/float(np2p1-1)              !* locus length increment
     !!
     !!
     px    = wk3x - wk1x
     py    = wk3y - wk1y                     !* wk1y=0.0 can be omitted
     !b    p     = sqrt(px**2 + py**2)             !* argument never = 0.0
     p     = sqrt(px*px + py*py)
     thp   = atan2(py,px)
     halfp = 0.5*p
     !!
     !!
     do n=np2p1,npts                      !* for npts = 30
       !!                                            !* n = 16 --> 30
       !!
       !b       b   = 0.5 * db  + float(n-np2p1) * db
       b   = db * (0.5 + float(n-np2p1))
       !b       a   = sqrt(1. + (2.*b/p)**2)
       a   = sqrt(1. + (2.*b/p)*(2.*b/p))
       dth = acos(1./a)
       a_halfp = a*halfp
       !!
       !b       wk2x(n) = a*halfp * cos(thp+dth)
       !b       wk2y(n) = a*halfp * sin(thp+dth)
       wk2x(n) = a_halfp * cos(thp+dth)
       wk2y(n) = a_halfp * sin(thp+dth)
       wk4x(n) = wk2x(n) - px
       wk4y(n) = wk2y(n) - py
       ds(n)   = db
       !!
       m       = npts - n + 1                !* m = 15 --> 1
       !b       wk2x(m) = a*halfp * cos(thp-dth)
       !b       wk2y(m) = a*halfp * sin(thp-dth)
       wk2x(m) = a_halfp * cos(thp-dth)
       wk2y(m) = a_halfp * sin(thp-dth)
       wk4x(m) = wk2x(m) - px
       wk4y(m) = wk2y(m) - py
       ds(m)   = db
       !!
     end do !  do n=np2p1,npts
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     RETURN
     !!

References ds, np2p1, npts, twopi, wk2x, wk2y, wk4x, wk4y, and wkfnc().

Referenced by gridsetr().

◆ snlr_fbi()

subroutine w3snl4md::snlr_fbi	(	real, dimension(nrng), intent(in)	pha,
		integer, intent(in)	ialt
	)

For a given Action Density array dens1(k,theta), computes the rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to wave-wave interaction.

For a given Action Density array dens1(k,theta), computes the rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to wave-wave interaction, as well as some ancillary arrays relating to positive and negative fluxes and their integrals.

Parameters

[in]	pha	Pha = kdkdtheta.
[in]	ialt	Integer switch.

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 3758 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!    ------------------------------------------------------------------
     !!
     !!    it returns: fbi & diag2 all dim=(nrng,nang)
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !! 1. Purpose :
     !!
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    For a given Action Density array dens1(k,theta), computes the    !
     !!    rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to   !
     !!    wave-wave interaction, as well as some ancillary arrays          !
     !!    relating to positive and negative fluxes and their integrals.    !
     !!                                                                     !
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    Compute:                                                         !
     !!    --------                                                         !
     !!    for both -tsa and -fbi                                           !
     !!    + sumint    contains scale 1 contibution for Snl -tsa & Snl -fbi !
     !!                                                                     !
     !!    for -tsa                                                         !
     !!    + sumintsa  contains tsa approximation to Snl -tsa               !
     !!                                                                     !
     !!    for -fbi                                                         !
     !!    + sumintp   contains scale 2 contribution to Snl -fbi            !
     !!    + sumintx   contains cross interactions between scales 1 and 2   !
     !!    -----------------------------------------------------------------#
     !!
     !! 2. Method :
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    nrng      int.  Public    I   # of freq. or rings
     !!    nang      int.  Public    I   # of angles
     !!    npts      int.  Public    I   # of points on the locus
     !!    nzz       int.  Public    I   linear irng x krng = (NK*(NK+1))/2
     !!    ialt      int.  Public    I   integer switch ialt=2; do alternate
     !!                                                 ialt=1; do not alternate
     !!    kzone     int.  Public    I   zone of influence = INT(alog(4.0)/alog(dfrq))
     !!    na2p1     int.  Public    I   = nang/2 + 1
     !!    np2p1     int.  Public    I   = npts/2 + 1
     !!    dfrq      real  Public    I   frequency multiplier for log freq. spacing
     !!    frqa      R.A.  Public    I   radian frequencies (Hz);   dim=(nrng)
     !!    pha       R.A.  local     I   pha = k*dk*dtheta      ;   dim=(nrng)
     !!    ------------------------------------------------------------------
     !!
     !!    *** The 11 grid integration geometry arrays at one given depth
     !!    *** from gridsetr.            dim=(npts,nang,nzz,ndep)
     !!    kref2     I.A.  Public    I   Index of reference wavenumber for k2
     !!    kref4     I.A.  Public    I   Idem for k4
     !!    jref2     I.A.  Public    I   Index of reference angle      for k2
     !!    jref4     I.A.  Public    I   Idem for k4
     !!    wtk2      R.A.  Public    I   k2 Interpolation weigth along wavenumbers
     !!    wtk4      R.A.  Public    I   Idem for k4
     !!    wta2      R.A.  Public    I   k2 Interpolation weigth along angles
     !!    wta4      R.A.  Public    I   Idem for k4
     !!    tfac2     R.A.  Public    I   Norm. for interp Action Density at k2
     !!    tfac4     R.A.  Public    I   Idem for k4
     !!    grad      R.A.  Public    I   Coupling and gradient term in integral
     !!                                  grad = C * H * g**2 * ds / |dW/dn|
     !!    ------------------------------------------------------------------
     !!
     !!    *** large & small scale Action Density from optsa dim=(nrng,nang)
     !!    dens1     R.A.  Public    I   lrg-scl Action Density (k,theta);
     !!    dens2     R.A.  Public    I   Sml-scl Action Density (k,theta);
     !!    ------------------------------------------------------------------
     !!
     !!    for both -tsa and -fbi
     !!    sumint    R.A.  local     O   contains scale 1 contribution to Snl
     !!                                                       dim=(nrng,nang)
     !!    for -tsa
     !!    sumintsa  R.A.  local     O   contains tsa approximation to Snl -tsa
     !!                                                       dim=(nrng,nang)
     !!    for -fbi
     !!    sumintp   R.A.  local     O   contains scale 2 contribution to Snl -fbi
     !!                                                       dim=(nrng,nang)
     !!    sumintx   R.A.  local     O   contains cross interactions " "   "  -fbi
     !!                                                       dim=(nrng,nang)
     !!    ------------------------------------------------------------------
     !!
     !!    for -tsa; The 2 returned arrays tsa & diag dim=(nrng,nang)
     !!    tsa       R.A.  Public    O   Snl-tsa = sumint + sumintsa
     !!    diag      R.A.  Public    O   Snl-tsa diagonal term = [dN/dn1]
     !!    ------------------------------------------------------------------
     !!
     !!    for -fbi; The 2 returned arrays fbi & diag2 dim=(nrng,nang)
     !!    fbi       R.A.  Public    O   Snl-fbi = sumint + sumintp  + sumintx
     !!    diag2     R.A.  Public    O   Snl-fbi diagonal term = [dN/dn1]
     !!    ------------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!    ----------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!     gridsetr  Subr. W3SNL4MD Setup geometric integration grid
     !!    ----------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!
     !!    !/S  Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     IMPLICIT NONE
     !!
     !!    Parameter list
     !!    --------------
     real,    intent(in)  :: pha(nrng)
     integer, intent(in)  :: ialt
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    Local Parameters & variables
     !!    -----------------------------
     !!    for both -tsa and -fbi
     integer              :: irng,krng, iang,kang
     integer              :: ipt, iizz, izz
     integer              :: kmax
     integer              :: ia2, ia2p, k2, k2p
     integer              :: ia4, ia4p, k4, k4p
     integer              :: nref
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for -tsa
     !tsa  integer              :: nklimit, nalimit
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for both -tsa and -fbi
     real                 :: d1,   d3,   d2,   d4
     real                 :: dp1,  dp3
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for both -tsa and -fbi
     !!    but for -tsa they are being calc. inside  if/endif if test is successful
     !!    and for -fbi they are being calc. outside if/endif always
     real                 :: dz4,  dz5
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for both -tsa and -fbi
     real                 :: dx13, ds13, dxp13, dsp13
     real                 :: dgm,  t31,  tr31
     real                 :: w2,   w2p,  wa2,  wa2p,  d2a,  d2b,  tt2
     real                 :: w4,   w4p,  wa4,  wa4p,  d4a,  d4b,  tt4
     real                 :: sumint(nrng,nang)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for -tsa
     !tsa  real                 :: dz2a, dz3a,  ttsa,   trtsa
     !tsa  real                 :: ddn1, ddn3,  diagk1, diagk3
     !tsa  real                 :: sumintsa(nrng,nang)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for -fbi
     real                 :: dp2,   dp4
     real                 :: d2pa,  d4pa
     real                 :: d2pb,  d4pb
     real                 :: dz1,   dz2,   dz3,   dz6,   dz7,   dz8
     real                 :: dgmp,  tp31,  trp31, dzsum, txp31, trx31
     !!
     !!    for -fbi; Bash added 4 new terms for a full expression of diag2 term
     real                 :: ddp1,  ddp2,  ddp3,  ddp4   !* ddpi=di+dpi for i=1,4
     real                 :: dd2n1, dd2n3, diag2k1, diag2k3
     real                 :: sumintp(nrng,nang)
     real                 :: sumintx(nrng,nang)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!    for -tsa
     !!    Bash; hardwire these two parameters
     !tsa  nklimit = 6
     !tsa  nalimit = 6
     !!
     !!
     !!ini
     !!    Bash; move initialization of all returned arrays from below to here
     !!    ------------------------------------------------------------------
     !!    for both -tsa and -fbi
     !!    sumint is now initialized here instead of below!
     sumint(:,:)   = 0.0
     !!
     !!    for -tsa
     !!    sumintsa are now initialized here instead of below!
     !tsa  sumintsa(:,:) = 0.0
     !tsa  tsa(:,:)      = 0.0
     !tsa  diag(:,:)     = 0.0
     !!
     !!    for -fbi
     !!    sumintp and sumintx are now initialized here instead of below!
     sumintp(:,:)  = 0.0
     sumintx(:,:)  = 0.0
     fbi(:,:)      = 0.0
     diag2(:,:)    = 0.0
     !!ini---
     !!    ------------------------------------------------------------------
     !!    ------------------------------------------------------------------
     !!    ##################################################################
     !!
     !!
     !!    for -tsa
     !tsa  ddn1    = 0.0                           !* for -tsa diag  [dN/dn1]
     !tsa  ddn3    = 0.0                           !* for -tsa diag  [dN/dn3]
     !!
     !!    for -fbi
     dd2n1   = 0.0                           !* for -fbi diag2 [dN/dn1]
     dd2n3   = 0.0                           !* for -fbi diag2 [dN/dn3]
     !!
     !!
     !!
     !!50
     do 50 irng=1,nrng,ialt
       !!kz
       kmax = min(irng+kzone, nrng)   !* Bash; Sometimes a locus pt is outside nrng
       !kz     kmax = min(irng+kzone, nrng-1) !* Bash; Taking 1 out will not affect kzone, try it
       !!kz---
       !!
       iizz = (nrng-1)*(irng-1) - ((irng-2)*(irng-1))/2
       !!      ----------------------------------------------------------------
       !!
       !!
       !!60
       do 60 iang=1,nang,ialt
         !!
         !!        for both -tsa and -fbi
         d1   = dens1(irng,iang)
         dp1  = dens2(irng,iang)
         !!
         !!        for -fbi
         ddp1 = d1+dp1             !! for full expression of diag2 term
         !!
         !!70
         !!kz
         !kz       do 70 krng=irng,nrng
         do 70 krng=irng,kmax,ialt
           !!
           !!          for both -tsa and -fbi
           !!          Bash; check5  be consistent with gridsetr
           !!          moved here from below (was after do 80 kang=1,nang)
           !!          and changed go to 80 into go to 70 (i.e. go to next krng)
           !kz         if ( frqa(krng)/frqa(irng) .gt. 4. ) go to 70  !* original gridsetr
           !kz         if ( frqa(krng)/frqa(irng) .gt. 3. ) go to 70  !* original snlr_'s
           !kz         if ( frqa(krng)/frqa(irng) .gt. 2. ) go to 70  !* Bash; use .gt. 2
           !!kz---
           !!
           izz = krng + iizz
           !!          ------------------------------------------------------------
           !!
           !!80
           do 80 kang=1,nang,ialt
             !!
             !!            for both -tsa and -fbi
             !!ba1         Bash; Remove self interaction
             !!                  skip k1 but keep the opposite angle to k1 - original setting
             if ( krng.eq.irng ) then              !* wn3 = wn1
               if ( kang.eq.iang ) go to 80        !* th3 = th1
             endif
             !!ba1---
             !!            ----------------------------------------------------------
             !!
             !!            for both -tsa and -fbi
             d3   = dens1(krng,kang)
             dp3  = dens2(krng,kang)
             !!
             !!            for -fbi
             ddp3 = d3+dp3         !! for full expression of diag2 term
             !!
             !!
             !!            for both -tsa and -fbi
             nref = kang - iang + 1
             if ( nref .lt. 1 ) nref = nref + nang
             !!
             !!
             !!            for both -tsa and -fbi
             !!            Bash; check5  be consistent with gridsetr
             !!                  and move this test above right after do 70 krng=irng,nrng
             !x            if ( frqa(krng)/frqa(irng) .gt. 4. ) go to 80  !* gridsetr
             !b            if ( frqa(krng)/frqa(irng) .gt. 3. ) go to 80  !* original
             !!
             !!
             !!            for both -tsa and -fbi
             t31     = 0.0             !* must be reset to 0.0
             !!
             !!            for -tsa
             !tsa          ttsa    = 0.0             !* must be reset to 0.0
             !tsa          diagk1  = 0.0             !* must be reset to 0.0
             !tsa          diagk3  = 0.0             !* must be reset to 0.0
             !!
             !!            for -fbi
             tp31    = 0.0             !* must be reset to 0.0
             txp31   = 0.0             !* must be reset to 0.0
             diag2k1 = 0.0             !* must be reset to 0.0
             diag2k3 = 0.0             !* must be reset to 0.0
             !!
             !!            for both -tsa and -fbi
             dx13  = d1*d3
             ds13  = d3-d1
             dxp13 = dp1*dp3
             dsp13 = dp3-dp1
             !!            ----------------------------------------------------------
             !!
             !!90
             do 90 ipt=1,npts
               !!
               !!              for both -tsa and -fbi
               !!              save time by skipping insignificant contributions
               !!e-30
               !e-30           if ( grad(ipt,nref,izz) .lt. 1.e-30 ) go to 90
               !!e-30---
               if ( grad(ipt,nref,izz) .lt. 1.e-15 ) go to 90
               !!e-30---
               !!              --------------------------------------------------------
               !!
               !!xlc1          Bash; skip k1 but keep the opposite angle to k1 - original setting
               !xlc1           if ( kang.eq.iang ) then                     !* th3=+th1
               !xlc1             if (ipt.eq.1 .or. ipt.eq.np2p1) go to 90   !* skip x-axis loci
               !xlc1           end if
               !!xlc1---
               !!              --------------------------------------------------------
               !!
               !!
               !!2             Estimation of Density for wave #2
               !!
               !!              for both -tsa and -fbi
               k2  = kref2(ipt,nref,izz)
               k2p = k2 + 1
               w2  = wtk2(ipt,nref,izz)
               w2p = 1. - w2
               !!
               !!              for both -tsa and -fbi
               ia2 = iang + jref2(ipt,nref,izz)
               if ( ia2 .gt. nang )  ia2  = ia2  - nang
               !!
               !!              for both -tsa and -fbi
               ia2p = ia2 + 1
               if ( ia2p .gt. nang ) ia2p = ia2p - nang
               !!
               !!              for both -tsa and -fbi
               wa2  = wta2(ipt,nref,izz)
               wa2p = 1. - wa2
               d2a  = w2 * dens1(k2,ia2)  + w2p * dens1(k2p,ia2)
               d2b  = w2 * dens1(k2,ia2p) + w2p * dens1(k2p,ia2p)
               tt2  = tfac2(ipt,nref,izz)
               d2   = (wa2*d2a  + wa2p*d2b)  * tt2
               !!
               !!              for -fbi
               d2pa = w2 * dens2(k2,ia2)  + w2p * dens2(k2p,ia2)
               d2pb = w2 * dens2(k2,ia2p) + w2p * dens2(k2p,ia2p)
               !!
               !!              for -fbi
               dp2  = (wa2*d2pa + wa2p*d2pb) * tt2   !* for -fbi
               ddp2 = d2+dp2       !! for full expression of diag2 term
               !!              ========================================================
               !!
               !!
               !!4             Estimation of Density for wave #4
               !!
               !!              for both -tsa and -fbi
               k4  = kref4(ipt,nref,izz)
               k4p = k4 + 1
               w4  = wtk4(ipt,nref,izz)
               w4p = 1. - w4
               !!
               !!              for both -tsa and -fbi
               ia4 = iang + jref4(ipt,nref,izz)
               if ( ia4 .gt. nang )  ia4  = ia4  - nang
               !!
               !!              for both -tsa and -fbi
               ia4p= ia4 + 1
               if ( ia4p .gt. nang ) ia4p = ia4p - nang
               !!
               !!              for both -tsa and -fbi
               wa4  = wta4(ipt,nref,izz)
               wa4p = 1. - wa4
               d4a  = w4*dens1(k4,ia4)  + w4p*dens1(k4p,ia4)
               d4b  = w4*dens1(k4,ia4p) + w4p*dens1(k4p,ia4p)
               tt4  = tfac4(ipt,nref,izz)
               d4   = (wa4*d4a  + wa4p*d4b)  * tt4
               !!
               !!              for -fbi
               d4pa = w4*dens2(k4,ia4)  + w4p*dens2(k4p,ia4)
               d4pb = w4*dens2(k4,ia4p) + w4p*dens2(k4p,ia4p)
               !!
               !!              for -fbi
               dp4  = (wa4*d4pa + wa4p*d4pb) * tt4   !* for -fbi
               ddp4 = d4+dp4       !! for full expression of diag2 term
               !!              ========================================================
               !!
               !!
               !!              for both -tsa and -fbi
               dgm  = dx13*(d4-d2) + ds13*d4*d2        !* dgm=B of R&P'08 eqn(8)
               !!                                         !* represents Broad Scale interactions
               t31  = t31  + dgm  * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!
               !!              for -fbi
               dgmp = dxp13*(dp4-dp2) + dsp13*dp4*dp2  !* dgmp=L of R&P'08 eqn(8)
               !!                                          !* represents Local Scale interactions
               tp31 = tp31 + dgmp * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!              ========================================================
               !!
               !!
               !!              for -tsa : -diag
               !!              use this expression for the diagonal term
               !!              whose derivation neglect "dp2" & "dp4"
               !tsa            ddn1   = (d3+dp3)*(d4-d2) - d4*d2              !* dN/dn1
               !tsa            ddn3   = (d1+dp1)*(d4-d2) + d4*d2              !* dN/dn3
               !tsa            diagk1 = diagk1 + ddn1 * grad(ipt,nref,izz)
               !tsa            diagk3 = diagk3 + ddn3 * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!
               !!              for -fbi : -diag2
               !!              use the full expression for the diagonal terms
               !!              whose derivation keeps all large + small scale
               dd2n1   = ddp3*(ddp4-ddp2) - ddp4*ddp2         !* dN/dn1
               dd2n3   = ddp1*(ddp4-ddp2) + ddp4*ddp2         !* dN/dn3
               diag2k1 = diag2k1 + dd2n1 * grad(ipt,nref,izz)
               diag2k3 = diag2k3 + dd2n3 * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!              ========================================================
               !!
               !!
               !!              for -fbi
               dz1  = dx13    * (dp4-dp2)
               dz2  = d1*dp3  * ((d4-d2)+(dp4-dp2))
               dz3  = d3*dp1  * ((d4-d2)+(dp4-dp2))
               !!
               !!              for -fbi (calc. dz4 & dz5 here)
               dz4  = dxp13   * (d4-d2)
               dz5  = d2*d4   *  dsp13
               !!
               !!              for -tsa
               !!              Cross-interactions between parametric and perturbation
               !!              that occur only when k3 is close enough to k1
               !!              Bash; added an extra check on (nang-nalimit)
               !b              if ( iabs(irng-krng).lt.nklimit .and.                 &
               !b                   iabs(iang-kang).lt.nalimit )    then        !* original
               !!
               !tsa            if (     (krng-irng).lt.nklimit .and.                 &
               !tsa               ( iabs(kang-iang).lt.nalimit .or.                  &
               !tsa                 iabs(kang-iang).gt.(nang-nalimit) ) )  then !* Bash
               !!
               !!                for -tsa (calc. dz4 & dz5 here)
               !tsa              dz4  = dxp13   * (d4-d2)
               !tsa              dz5  = d2*d4   *  dsp13
               !tsa              dz2a = d1*dp3 * (d4-d2)
               !tsa              dz3a = d3*dp1 * (d4-d2)
               !!
               !tsa              ttsa = ttsa + (dz4+dz5+dz2a+dz3a)*grad(ipt,nref,izz)
               !!
               !tsa            endif
               !!              --------------------------------------------------------
               !!
               !!              for -fbi
               dz6   = d2*dp4  * (ds13+dsp13)
               dz7   = d4*dp2  * (ds13+dsp13)
               dz8   = dp2*dp4 *  ds13
               dzsum = dz1 + dz2 + dz3 + dz4 + dz5 + dz6 + dz7 + dz8
               txp31 = txp31 + dzsum * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!              ========================================================
               !!
               !!
 90          end do                        !* end of ipt (locus) loop
             !!            ----------------------------------------------------------
             !!
             !!
             !!            multiply the following components by factor 2. in here
             !!
             !!            for both -tsa and -fbi
             tr31  = 2. * t31
             !!
             !!            for -tsa
             !tsa          trtsa = 2. * ttsa
             !!
             !!            for -fbi
             trp31 = 2. * tp31
             trx31 = 2. * txp31
             !!
             !!            for -tsa : -diag
             !tsa          diagk1  = 2. * diagk1
             !tsa          diagk3  = 2. * diagk3
             !!
             !!            for -fbi : -diag2
             diag2k1 = 2. * diag2k1
             diag2k3 = 2. * diag2k3
             !!            ----------------------------------------------------------
             !!
             !!            for both -tsa and -fbi
             sumint(irng,iang)  = sumint(irng,iang)  + tr31*pha(krng)
             sumint(krng,kang)  = sumint(krng,kang)  - tr31*pha(irng)
             !!            ----------------------------------------------------------
             !!
             !!            for -tsa
             !tsa          sumintsa(irng,iang)= sumintsa(irng,iang)+ trtsa*pha(krng)
             !tsa          sumintsa(krng,kang)= sumintsa(krng,kang)- trtsa*pha(irng)
             !!            ----------------------------------------------------------
             !!
             !!            for -fbi
             sumintp(irng,iang) = sumintp(irng,iang) + trp31*pha(krng)
             sumintp(krng,kang) = sumintp(krng,kang) - trp31*pha(irng)
             !!
             !!            for -fbi
             sumintx(irng,iang) = sumintx(irng,iang) + trx31*pha(krng)
             sumintx(krng,kang) = sumintx(krng,kang) - trx31*pha(irng)
             !!            ----------------------------------------------------------
             !!
             !!            for -tsa : -diag
             !tsa          diag(irng,iang) = diag(irng,iang)  + diagk1*pha(krng)
             !tsa          diag(krng,kang) = diag(krng,kang)  - diagk3*pha(irng)
             !!            ----------------------------------------------------------
             !!
             !!            for -fbi : -diag2
             diag2(irng,iang) = diag2(irng,iang) + diag2k1*pha(krng)
             diag2(krng,kang) = diag2(krng,kang) - diag2k3*pha(irng)
             !!            ----------------------------------------------------------
             !!
 80        end do                              !* end of kang loop
           !!
 70      end do                                !* end of krng loop
         !!
 60    end do                                  !* end of iang loop
       !!
 50  end do                                    !* end of irng loop
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!    Final sum-up to get Snl and diag. term to be returned
     !!
     !!    for -tsa
     !tsa  tsa(:,:)   = sumint(:,:) + sumintsa(:,:)
     !b    diag(:,:)  = diag(:,:)   !* is Ok, already summed up
     !!
     !!    for -fbi
     fbi(:,:)   = sumint(:,:) + sumintp(:,:) + sumintx(:,:)
     !b    diag2(:,:) = diag2(:,:)  !* is Ok, already summed up
     !!    --------------------------------------------------------------------------
     !!    ==========================================================================
     !!
     !!
     !!alt Call interp2 only if ialt=2,
     !!    Interpolate bi-linearly to fill in tsa/fbi & diag/diag2 arrays
     !!    after alternating the irng, iang, krng & kang loops above
     !!    ------------------------------------------------------------------
     if ( ialt.eq.2 ) then
       !!      for -tsa
       !tsa    call interp2 ( tsa )
       !tsa    call interp2 ( diag )
       !!
       !!      for -fbi
       call interp2 ( fbi )
       call interp2 ( diag2 )
     endif
     !!alt---
     !!    --------------------------------------------------------------------------
     !!    ==========================================================================
     !!
     RETURN
     !!

References dens1, dens2, diag2, fbi, grad, interp2(), jref2, jref4, kref2, kref4, kzone, npts, tfac2, tfac4, wta2, wta4, wtk2, and wtk4.

Referenced by w3snl4().

◆ snlr_tsa()

subroutine w3snl4md::snlr_tsa	(	real, dimension(nrng), intent(in)	pha,
		integer, intent(in)	ialt
	)

For a given Action Density array dens1(k,theta), computes the rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to wave-wave interaction.

For a given Action Density array dens1(k,theta), computes the rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to wave-wave interaction, as well as some ancillary arrays relating to positive and negative fluxes and their integrals.

Parameters

[in]	pha	Pha = kdkdtheta.
[in]	ialt	Integer switch.

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 4391 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!    ------------------------------------------------------------------
     !!
     !!    it returns: tsa & diag  all dim=(nrng,nang)
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !!
     !! 1. Purpose :
     !!
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    For a given Action Density array dens1(k,theta), computes the    !
     !!    rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to   !
     !!    wave-wave interaction, as well as some ancillary arrays          !
     !!    relating to positive and negative fluxes and their integrals.    !
     !!                                                                     !
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    Compute:                                                         !
     !!    --------                                                         !
     !!    for both -tsa and -fbi                                           !
     !!    + sumint    contains scale 1 contibution for Snl -tsa & Snl -fbi !
     !!                                                                     !
     !!    for -tsa                                                         !
     !!    + sumintsa  contains tsa approximation to Snl -tsa               !
     !!                                                                     !
     !!    for -fbi                                                         !
     !!    + sumintp   contains scale 2 contribution to Snl -fbi            !
     !!    + sumintx   contains cross interactions between scales 1 and 2   !
     !!    -----------------------------------------------------------------#
     !!
     !! 2. Method :
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    nrng      int.  Public    I   # of freq. or rings
     !!    nang      int.  Public    I   # of angles
     !!    npts      int.  Public    I   # of points on the locus
     !!    nzz       int.  Public    I   linear irng x krng = (NK*(NK+1))/2
     !!    ialt      int.  Public    I   integer switch ialt=2; do alternate
     !!                                                 ialt=1; do not alternate
     !!    kzone     int.  Public    I   zone of influence = INT(alog(4.0)/alog(dfrq))
     !!    na2p1     int.  Public    I   = nang/2 + 1
     !!    np2p1     int.  Public    I   = npts/2 + 1
     !!    dfrq      real  Public    I   frequency multiplier for log freq. spacing
     !!    frqa      R.A.  Public    I   radian frequencies (Hz);   dim=(nrng)
     !!    pha       R.A.  local     I   pha = k*dk*dtheta      ;   dim=(nrng)
     !!    ------------------------------------------------------------------
     !!
     !!    *** The 11 grid integration geometry arrays at one given depth
     !!    *** from gridsetr.            dim=(npts,nang,nzz,ndep)
     !!    kref2     I.A.  Public    I   Index of reference wavenumber for k2
     !!    kref4     I.A.  Public    I   Idem for k4
     !!    jref2     I.A.  Public    I   Index of reference angle      for k2
     !!    jref4     I.A.  Public    I   Idem for k4
     !!    wtk2      R.A.  Public    I   k2 Interpolation weigth along wavenumbers
     !!    wtk4      R.A.  Public    I   Idem for k4
     !!    wta2      R.A.  Public    I   k2 Interpolation weigth along angles
     !!    wta4      R.A.  Public    I   Idem for k4
     !!    tfac2     R.A.  Public    I   Norm. for interp Action Density at k2
     !!    tfac4     R.A.  Public    I   Idem for k4
     !!    grad      R.A.  Public    I   Coupling and gradient term in integral
     !!                                  grad = C * H * g**2 * ds / |dW/dn|
     !!    ------------------------------------------------------------------
     !!
     !!    *** large & small scale Action Density from optsa dim=(nrng,nang)
     !!    dens1     R.A.  Public    I   lrg-scl Action Density (k,theta);
     !!    dens2     R.A.  Public    I   Sml-scl Action Density (k,theta);
     !!    ------------------------------------------------------------------
     !!
     !!    for both -tsa and -fbi
     !!    sumint    R.A.  local     O   contains scale 1 contribution to Snl
     !!                                                       dim=(nrng,nang)
     !!    for -tsa
     !!    sumintsa  R.A.  local     O   contains tsa approximation to Snl -tsa
     !!                                                       dim=(nrng,nang)
     !!    for -fbi
     !!    sumintp   R.A.  local     O   contains scale 2 contribution to Snl -fbi
     !!                                                       dim=(nrng,nang)
     !!    sumintx   R.A.  local     O   contains cross interactions " "   "  -fbi
     !!                                                       dim=(nrng,nang)
     !!    ------------------------------------------------------------------
     !!
     !!    for -tsa; The 2 returned arrays tsa & diag dim=(nrng,nang)
     !!    tsa       R.A.  Public    O   Snl-tsa = sumint + sumintsa
     !!    diag      R.A.  Public    O   Snl-tsa diagonal term = [dN/dn1]
     !!    ------------------------------------------------------------------
     !!
     !!    for -fbi; The 2 returned arrays fbi & diag2 dim=(nrng,nang)
     !!    fbi       R.A.  Public    O   Snl-fbi = sumint + sumintp  + sumintx
     !!    diag2     R.A.  Public    O   Snl-fbi diagonal term = [dN/dn1]
     !!    ------------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!    ----------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ----------------------------------------------------------------
     !!     gridsetr  Subr. W3SNL4MD Setup geometric integration grid
     !!    ----------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!
     !!    !/S  Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     IMPLICIT NONE
     !!
     !!    Parameter list
     !!    --------------
     real,    intent(in)  :: pha(nrng)
     integer, intent(in)  :: ialt
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!
     !!
     !!    Local Parameters & variables
     !!    -----------------------------
     !!    for both -tsa and -fbi
     integer              :: irng,krng, iang,kang
     integer              :: ipt, iizz, izz
     integer              :: kmax
     integer              :: ia2, ia2p, k2, k2p
     integer              :: ia4, ia4p, k4, k4p
     integer              :: nref
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for -tsa
     integer              :: nklimit, nalimit
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for both -tsa and -fbi
     real                 :: d1,   d3,   d2,   d4
     real                 :: dp1,  dp3
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for both -tsa and -fbi
     !!    but for -tsa they are being calc. inside  if/endif if test is successful
     !!    and for -fbi they are being calc. outside if/endif always
     real                 :: dz4,  dz5
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for both -tsa and -fbi
     real                 :: dx13, ds13, dxp13, dsp13
     real                 :: dgm,  t31,  tr31
     real                 :: w2,   w2p,  wa2,  wa2p,  d2a,  d2b,  tt2
     real                 :: w4,   w4p,  wa4,  wa4p,  d4a,  d4b,  tt4
     real                 :: sumint(nrng,nang)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for -tsa
     real                 :: dz2a, dz3a,  ttsa,   trtsa
     real                 :: ddn1, ddn3,  diagk1, diagk3
     real                 :: sumintsa(nrng,nang)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     !!    for -fbi
     !fbi  real                 :: dp2,   dp4
     !fbi  real                 :: d2pa,  d4pa
     !fbi  real                 :: d2pb,  d4pb
     !fbi  real                 :: dz1,   dz2,   dz3,   dz6,   dz7,   dz8
     !fbi  real                 :: dgmp,  tp31,  trp31, dzsum, txp31, trx31
     !!
     !!    for -fbi; Bash added 4 new terms for a full expression of diag2 term
     !fbi  real                 :: ddp1,  ddp2,  ddp3,  ddp4   !* ddpi=di+dpi for i=1,4
     !fbi  real                 :: dd2n1, dd2n3, diag2k1, diag2k3
     !fbi  real                 :: sumintp(nrng,nang)
     !fbi  real                 :: sumintx(nrng,nang)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!    for -tsa
     !!    Bash; hardwire these two parameters
     nklimit = 6
     nalimit = 6
     !!
     !!
     !!ini
     !!    Bash; move initialization of all returned arrays from below to here
     !!    ------------------------------------------------------------------
     !!    for both -tsa and -fbi
     !!    sumint is now initialized here instead of below!
     sumint(:,:)   = 0.0
     !!
     !!    for -tsa
     !!    sumintsa are now initialized here instead of below!
     sumintsa(:,:) = 0.0
     tsa(:,:)      = 0.0
     diag(:,:)     = 0.0
     !!
     !!    for -fbi
     !!    sumintp and sumintx are now initialized here instead of below!
     !fbi  sumintp(:,:)  = 0.0
     !fbi  sumintx(:,:)  = 0.0
     !fbi  fbi(:,:)      = 0.0
     !fbi  diag2(:,:)    = 0.0
     !!ini---
     !!    ------------------------------------------------------------------
     !!    ------------------------------------------------------------------
     !!    ##################################################################
     !!
     !!
     !!    for -tsa
     ddn1    = 0.0                           !* for -tsa diag  [dN/dn1]
     ddn3    = 0.0                           !* for -tsa diag  [dN/dn3]
     !!
     !!    for -fbi
     !fbi  dd2n1   = 0.0                           !* for -fbi diag2 [dN/dn1]
     !fbi  dd2n3   = 0.0                           !* for -fbi diag2 [dN/dn3]
     !!
     !!
     !!
     !!50
     do 50 irng=1,nrng,ialt
       !!kz
       kmax = min(irng+kzone, nrng)   !* Bash; Sometimes a locus pt is outside nrng
       !kz     kmax = min(irng+kzone, nrng-1) !* Bash; Taking 1 out will not affect kzone, try it
       !!kz---
       !!
       iizz = (nrng-1)*(irng-1) - ((irng-2)*(irng-1))/2
       !!      ----------------------------------------------------------------
       !!
       !!
       !!60
       do 60 iang=1,nang,ialt
         !!
         !!        for both -tsa and -fbi
         d1   = dens1(irng,iang)
         dp1  = dens2(irng,iang)
         !!
         !!        for -fbi
         !fbi      ddp1 = d1+dp1             !! for full expression of diag2 term
         !!
         !!70
         !!kz
         !kz       do 70 krng=irng,nrng
         do 70 krng=irng,kmax,ialt
           !!
           !!          for both -tsa and -fbi
           !!          Bash; check5  be consistent with gridsetr
           !!          moved here from below (was after do 80 kang=1,nang)
           !!          and changed go to 80 into go to 70 (i.e. go to next krng)
           !kz         if ( frqa(krng)/frqa(irng) .gt. 4. ) go to 70  !* original gridsetr
           !kz         if ( frqa(krng)/frqa(irng) .gt. 3. ) go to 70  !* original snlr_'s
           !kz         if ( frqa(krng)/frqa(irng) .gt. 2. ) go to 70  !* Bash; use .gt. 2
           !!kz---
           !!
           izz = krng + iizz
           !!          ------------------------------------------------------------
           !!
           !!80
           do 80 kang=1,nang,ialt
             !!
             !!            for both -tsa and -fbi
             !!ba1         Bash; Remove self interaction
             !!                  skip k1 but keep the opposite angle to k1 - original setting
             if ( krng.eq.irng ) then              !* wn3 = wn1
               if ( kang.eq.iang ) go to 80        !* th3 = th1
             endif
             !!ba1---
             !!            ----------------------------------------------------------
             !!
             !!            for both -tsa and -fbi
             d3   = dens1(krng,kang)
             dp3  = dens2(krng,kang)
             !!
             !!            for -fbi
             !fbi          ddp3 = d3+dp3         !! for full expression of diag2 term
             !!
             !!
             !!            for both -tsa and -fbi
             nref = kang - iang + 1
             if ( nref .lt. 1 ) nref = nref + nang
             !!
             !!
             !!            for both -tsa and -fbi
             !!            Bash; check5  be consistent with gridsetr
             !!                  and move this test above right after do 70 krng=irng,nrng
             !x            if ( frqa(krng)/frqa(irng) .gt. 4. ) go to 80  !* gridsetr
             !b            if ( frqa(krng)/frqa(irng) .gt. 3. ) go to 80  !* original
             !!
             !!
             !!            for both -tsa and -fbi
             t31     = 0.0             !* must be reset to 0.0
             !!
             !!            for -tsa
             ttsa    = 0.0             !* must be reset to 0.0
             diagk1  = 0.0             !* must be reset to 0.0
             diagk3  = 0.0             !* must be reset to 0.0
             !!
             !!            for -fbi
             !fbi          tp31    = 0.0             !* must be reset to 0.0
             !fbi          txp31   = 0.0             !* must be reset to 0.0
             !fbi          diag2k1 = 0.0             !* must be reset to 0.0
             !fbi          diag2k3 = 0.0             !* must be reset to 0.0
             !!
             !!            for both -tsa and -fbi
             dx13  = d1*d3
             ds13  = d3-d1
             dxp13 = dp1*dp3
             dsp13 = dp3-dp1
             !!            ----------------------------------------------------------
             !!
             !!90
             do 90 ipt=1,npts
               !!
               !!              for both -tsa and -fbi
               !!              save time by skipping insignificant contributions
               !!e-30
               !e-30           if ( grad(ipt,nref,izz) .lt. 1.e-30 ) go to 90
               !!e-30---
               if ( grad(ipt,nref,izz) .lt. 1.e-15 ) go to 90
               !!e-30---
               !!              --------------------------------------------------------
               !!
               !!xlc1          Bash; skip k1 but keep the opposite angle to k1 - original setting
               !xlc1           if ( kang.eq.iang ) then                     !* th3=+th1
               !xlc1             if (ipt.eq.1 .or. ipt.eq.np2p1) go to 90   !* skip x-axis loci
               !xlc1           end if
               !!xlc1---
               !!              --------------------------------------------------------
               !!
               !!
               !!2             Estimation of Density for wave #2
               !!
               !!              for both -tsa and -fbi
               k2  = kref2(ipt,nref,izz)
               k2p = k2 + 1
               w2  = wtk2(ipt,nref,izz)
               w2p = 1. - w2
               !!
               !!              for both -tsa and -fbi
               ia2 = iang + jref2(ipt,nref,izz)
               if ( ia2 .gt. nang )  ia2  = ia2  - nang
               !!
               !!              for both -tsa and -fbi
               ia2p = ia2 + 1
               if ( ia2p .gt. nang ) ia2p = ia2p - nang
               !!
               !!              for both -tsa and -fbi
               wa2  = wta2(ipt,nref,izz)
               wa2p = 1. - wa2
               d2a  = w2 * dens1(k2,ia2)  + w2p * dens1(k2p,ia2)
               d2b  = w2 * dens1(k2,ia2p) + w2p * dens1(k2p,ia2p)
               tt2  = tfac2(ipt,nref,izz)
               d2   = (wa2*d2a  + wa2p*d2b)  * tt2
               !!
               !!              for -fbi
               !fbi            d2pa = w2 * dens2(k2,ia2)  + w2p * dens2(k2p,ia2)
               !fbi            d2pb = w2 * dens2(k2,ia2p) + w2p * dens2(k2p,ia2p)
               !!
               !!              for -fbi
               !fbi            dp2  = (wa2*d2pa + wa2p*d2pb) * tt2   !* for -fbi
               !fbi            ddp2 = d2+dp2       !! for full expression of diag2 term
               !!              ========================================================
               !!
               !!
               !!4             Estimation of Density for wave #4
               !!
               !!              for both -tsa and -fbi
               k4  = kref4(ipt,nref,izz)
               k4p = k4 + 1
               w4  = wtk4(ipt,nref,izz)
               w4p = 1. - w4
               !!
               !!              for both -tsa and -fbi
               ia4 = iang + jref4(ipt,nref,izz)
               if ( ia4 .gt. nang )  ia4  = ia4  - nang
               !!
               !!              for both -tsa and -fbi
               ia4p= ia4 + 1
               if ( ia4p .gt. nang ) ia4p = ia4p - nang
               !!
               !!              for both -tsa and -fbi
               wa4  = wta4(ipt,nref,izz)
               wa4p = 1. - wa4
               d4a  = w4*dens1(k4,ia4)  + w4p*dens1(k4p,ia4)
               d4b  = w4*dens1(k4,ia4p) + w4p*dens1(k4p,ia4p)
               tt4  = tfac4(ipt,nref,izz)
               d4   = (wa4*d4a  + wa4p*d4b)  * tt4
               !!
               !!              for -fbi
               !fbi            d4pa = w4*dens2(k4,ia4)  + w4p*dens2(k4p,ia4)
               !fbi            d4pb = w4*dens2(k4,ia4p) + w4p*dens2(k4p,ia4p)
               !!
               !!              for -fbi
               !fbi            dp4  = (wa4*d4pa + wa4p*d4pb) * tt4   !* for -fbi
               !fbi            ddp4 = d4+dp4       !! for full expression of diag2 term
               !!              ========================================================
               !!
               !!
               !!              for both -tsa and -fbi
               dgm  = dx13*(d4-d2) + ds13*d4*d2        !* dgm=B of R&P'08 eqn(8)
               !!                                         !* represents Broad Scale interactions
               t31  = t31  + dgm  * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!
               !!              for -fbi
               !fbi            dgmp = dxp13*(dp4-dp2) + dsp13*dp4*dp2  !* dgmp=L of R&P'08 eqn(8)
               !!                                          !* represents Local Scale interactions
               !fbi            tp31 = tp31 + dgmp * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!              ========================================================
               !!
               !!
               !!              for -tsa : -diag
               !!              use this expression for the diagonal term
               !!              whose derivation neglect "dp2" & "dp4"
               ddn1   = (d3+dp3)*(d4-d2) - d4*d2              !* dN/dn1
               ddn3   = (d1+dp1)*(d4-d2) + d4*d2              !* dN/dn3
               diagk1 = diagk1 + ddn1 * grad(ipt,nref,izz)
               diagk3 = diagk3 + ddn3 * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!
               !!              for -fbi : -diag2
               !!              use the full expression for the diagonal terms
               !!              whose derivation keeps all large + small scale
               !fbi            dd2n1   = ddp3*(ddp4-ddp2) - ddp4*ddp2         !* dN/dn1
               !fbi            dd2n3   = ddp1*(ddp4-ddp2) + ddp4*ddp2         !* dN/dn3
               !fbi            diag2k1 = diag2k1 + dd2n1 * grad(ipt,nref,izz)
               !fbi            diag2k3 = diag2k3 + dd2n3 * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!              ========================================================
               !!
               !!
               !!              for -fbi
               !fbi            dz1  = dx13    * (dp4-dp2)
               !fbi            dz2  = d1*dp3  * ((d4-d2)+(dp4-dp2))
               !fbi            dz3  = d3*dp1  * ((d4-d2)+(dp4-dp2))
               !!
               !!              for -fbi (calc. dz4 & dz5 here)
               !fbi            dz4  = dxp13   * (d4-d2)
               !fbi            dz5  = d2*d4   *  dsp13
               !!
               !!              for -tsa
               !!              Cross-interactions between parametric and perturbation
               !!              that occur only when k3 is close enough to k1
               !!              Bash; added an extra check on (nang-nalimit)
               !b              if ( iabs(irng-krng).lt.nklimit .and.                 &
               !b                   iabs(iang-kang).lt.nalimit )    then        !* original
               !!
               if (     (krng-irng).lt.nklimit .and.                 &
                    ( iabs(kang-iang).lt.nalimit .or.                  &
                    iabs(kang-iang).gt.(nang-nalimit) ) )  then !* Bash
                 !!
                 !!                for -tsa (calc. dz4 & dz5 here)
                 dz4  = dxp13   * (d4-d2)
                 dz5  = d2*d4   *  dsp13
                 dz2a = d1*dp3 * (d4-d2)
                 dz3a = d3*dp1 * (d4-d2)
                 !!
                 ttsa = ttsa + (dz4+dz5+dz2a+dz3a)*grad(ipt,nref,izz)
                 !!
               endif
               !!              --------------------------------------------------------
               !!
               !!              for -fbi
               !fbi            dz6   = d2*dp4  * (ds13+dsp13)
               !fbi            dz7   = d4*dp2  * (ds13+dsp13)
               !fbi            dz8   = dp2*dp4 *  ds13
               !fbi            dzsum = dz1 + dz2 + dz3 + dz4 + dz5 + dz6 + dz7 + dz8
               !fbi            txp31 = txp31 + dzsum * grad(ipt,nref,izz)
               !!              --------------------------------------------------------
               !!              ========================================================
               !!
               !!
 90          end do                        !* end of ipt (locus) loop
             !!            ----------------------------------------------------------
             !!
             !!
             !!            multiply the following components by factor 2. in here
             !!
             !!            for both -tsa and -fbi
             tr31  = 2. * t31
             !!
             !!            for -tsa
             trtsa = 2. * ttsa
             !!
             !!            for -fbi
             !fbi          trp31 = 2. * tp31
             !fbi          trx31 = 2. * txp31
             !!
             !!            for -tsa : -diag
             diagk1  = 2. * diagk1
             diagk3  = 2. * diagk3
             !!
             !!            for -fbi : -diag2
             !fbi          diag2k1 = 2. * diag2k1
             !fbi          diag2k3 = 2. * diag2k3
             !!            ----------------------------------------------------------
             !!
             !!            for both -tsa and -fbi
             sumint(irng,iang)  = sumint(irng,iang)  + tr31*pha(krng)
             sumint(krng,kang)  = sumint(krng,kang)  - tr31*pha(irng)
             !!            ----------------------------------------------------------
             !!
             !!            for -tsa
             sumintsa(irng,iang)= sumintsa(irng,iang)+ trtsa*pha(krng)
             sumintsa(krng,kang)= sumintsa(krng,kang)- trtsa*pha(irng)
             !!            ----------------------------------------------------------
             !!
             !!            for -fbi
             !fbi          sumintp(irng,iang) = sumintp(irng,iang) + trp31*pha(krng)
             !fbi          sumintp(krng,kang) = sumintp(krng,kang) - trp31*pha(irng)
             !!
             !!            for -fbi
             !fbi          sumintx(irng,iang) = sumintx(irng,iang) + trx31*pha(krng)
             !fbi          sumintx(krng,kang) = sumintx(krng,kang) - trx31*pha(irng)
             !!            ----------------------------------------------------------
             !!
             !!            for -tsa : -diag
             diag(irng,iang) = diag(irng,iang)  + diagk1*pha(krng)
             diag(krng,kang) = diag(krng,kang)  - diagk3*pha(irng)
             !!            ----------------------------------------------------------
             !!
             !!            for -fbi : -diag2
             !fbi          diag2(irng,iang) = diag2(irng,iang) + diag2k1*pha(krng)
             !fbi          diag2(krng,kang) = diag2(krng,kang) - diag2k3*pha(irng)
             !!            ----------------------------------------------------------
             !!
 80        end do                              !* end of kang loop
           !!
 70      end do                                !* end of krng loop
         !!
 60    end do                                  !* end of iang loop
       !!
 50  end do                                    !* end of irng loop
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!    Final sum-up to get Snl and diag. term to be returned
     !!
     !!    for -tsa
     tsa(:,:)   = sumint(:,:) + sumintsa(:,:)
     !b    diag(:,:)  = diag(:,:)   !* is Ok, already summed up
     !!
     !!    for -fbi
     !fbi  fbi(:,:)   = sumint(:,:) + sumintp(:,:) + sumintx(:,:)
     !b    diag2(:,:) = diag2(:,:)  !* is Ok, already summed up
     !!    --------------------------------------------------------------------------
     !!    ==========================================================================
     !!
     !!
     !!alt Call interp2 only if ialt=2,
     !!    Interpolate bi-linearly to fill in tsa/fbi & diag/diag2 arrays
     !!    after alternating the irng, iang, krng & kang loops above
     !!    ------------------------------------------------------------------
     if ( ialt.eq.2 ) then
       !!      for -tsa
       call interp2 ( tsa )
       call interp2 ( diag )
       !!
       !!      for -fbi
       !fbi    call interp2 ( fbi )
       !fbi    call interp2 ( diag2 )
     endif
     !!alt---
     !!    --------------------------------------------------------------------------
     !!    ==========================================================================
     !!
     RETURN
     !!

References dens1, dens2, diag, grad, interp2(), jref2, jref4, kref2, kref4, kzone, npts, tfac2, tfac4, tsa, wta2, wta4, wtk2, and wtk4.

Referenced by w3snl4().

◆ w3snl4()

subroutine w3snl4md::w3snl4	(	real, dimension(nth,nk), intent(in)	A,
		real, dimension(nk), intent(in)	CG,
		real, dimension(nk), intent(in)	WN,
		real, intent(in)	DEPTH,
		real, dimension(nth,nk), intent(out)	S,
		real, dimension(nth,nk), intent(out)	D
	)

Interface module for TSA type nonlinear interactions.

Based on Resio and Perrie (2008) and Perrie and Resio (2009).

Parameters

[in,out]	A	2D Action Density A(NTH,NK) as function of direction (rad) and wavenumber (theta,k).
[in,out]	CG	Group velocities.
[in,out]	WN	Wavenumbers.
[in,out]	DEPTH	Water depth (m).
[in,out]	S	Source term.
[in,out]	D	Diagonal term of derivative.

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 611 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!    ------------------------------------------------------------------
     !!
     !!    it returns: S & D  dim = (NTH,NK) = (nang,nrng)
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !! 1. Purpose :
     !!
     !!    Interface module for TSA type nonlinear interactions.
     !!    Based on Resio and Perrie (2008) and Perrie and Resio (2009)
     !!
     !! 2. Method :
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    A         R.A.            I   2D Action Density A(NTH,NK) as function of
     !!                                  direction (rad) and wavenumber (theta,k)
     !!    CG        R.A.            I   Group velocities             dim=NK
     !!    WN        R.A.            I   Wavenumbers                  dim=NK
     !!    DEPTH     Real            I   Water depth (m)
     !!    S         R.A.            O   Source term.                 dim=(NTH,NK)
     !!    D         R.A.            O   Diagonal term of derivative. dim=(NTH,NK)
     !!    ------------------------------------------------------------------
     !!
     !!    nrng      int.  Public    O   # of freq. or rings
     !!    nang      int.  Public    O   # of angles
     !!    npts      int.  Public    I   # of points on the locus
     !!    ndep      int.  Public    I   # of depths in look-up tables
     !!    nzz       int.  Public    O   linear irngxkrng = (NK*(NK+1))/2
     !!    kzone     int.  Public    O   zone of influence = INT(alog(4.0)/alog(dfrq))
     !!    nb2fp     int.  Public    O   # of bins over fp => dfrq**nb2fp)*fp ~ 2.*fp
     !!                                                    = INT(alog(2.0)/alog(dfrq))
     !!    na2p1     int.  Public    O   = nang/2 + 1
     !!    np2p1     int.  Public    O   = npts/2 + 1
     !!    dfrq      Real  Public    O   frequency multiplier for log freq. spacing
     !!    f0        Real  Public    O   = frqa(1); first freq. (Hz)
     !!    ainc      Real  Public    O   = DTH; WW3 angle increment (radians)
     !!    twopi     Real  Public    O   = TPI; WW3i 2*pi = 8.*atan(1.) (radians)
     !!    oma       R.A.  Public    O   = SIG(1:NK)  WW3 waveumber array   dim=(nrng)
     !!    frqa      R.A.  Public    O   = oma(:)/twopi WW3 frequency array dim=(nrng)
     !!    angl      R.A.  Public    O   = TH(1:NTH);   WW3 angles array       dim=(nang)
     !!    sinan     R.A.  Public    O   = ESIN(1:NTH); WW3 sin(angl(:)) array dim=(nang)
     !!    cosan     R.A.  Public    O   = ECOS(1:NTH); WW3 cos(angl(:)) array dim=(nang)
     !!    dep_tbl   R.A.  Public    I   depthes in Look-up tables arrays dim=(ndep)
     !!    ------------------------------------------------------------------
     !!
     !!    *** The 11 look-up tables for grid integration geometry arrays
     !!    *** at all selected 'ndep' depths defined in dep_tbl(ndep)' array
     !!    *** from gridsetr.            dim=(npts,nang,nzz,ndep)
     !!    kref2_tbl I.A.  Public    I   Index of reference wavenumber for k2
     !!    kref4_tbl I.A.  Public    I   Idem for k4
     !!    jref2_tbl I.A.  Public    I   Index of reference angle      for k2
     !!    jref4_tbl I.A.  Public    I   Idem for k4
     !!    wtk2_tbl  R.A.  Public    I   k2 Interpolation weigth along wavenumbers
     !!    wtk4_tbl  R.A.  Public    I   Idem for k4
     !!    wta2_tbl  R.A.  Public    I   k2 Interpolation weigth along angles
     !!    wta4_tbl  R.A.  Public    I   Idem for k4
     !!    tfac2_tbl R.A.  Public    I   Norm. for interp Action Density at k2
     !!    tfac4_tbl R.A.  Public    I   Idem for k4
     !!    grad_tbl  R.A.  Public    I   Coupling and gradient term in integral
     !!                                  grad = C * H * g**2 * ds / |dW/dn|
     !!    ------------------------------------------------------------------
     !!
     !!    *** The 11 grid integration geometry arrays at one given depth
     !!    *** from gridsetr.            dim=(npts,nang,nzz,ndep)
     !!    kref2     I.A.  Public    O   Index of reference wavenumber for k2
     !!    kref4     I.A.  Public    O   Idem for k4
     !!    jref2     I.A.  Public    O   Index of reference angle      for k2
     !!    jref4     I.A.  Public    O   Idem for k4
     !!    wtk2      R.A.  Public    O   k2 Interpolation weigth along wavenumbers
     !!    wtk4      R.A.  Public    O   Idem for k4
     !!    wta2      R.A.  Public    O   k2 Interpolation weigth along angles
     !!    wta4      R.A.  Public    O   Idem for k4
     !!    tfac2     R.A.  Public    O   Norm. for interp Action Density at k2
     !!    tfac4     R.A.  Public    O   Idem for k4
     !!    grad      R.A.  Public    O   Coupling and gradient term in integral
     !!                                  grad = C * H * g**2 * ds / |dW/dn|
     !!    ------------------------------------------------------------------
     !!
     !!    ef2       R.A.  Public    O   2D Energy Density spectrum ef2(theta,f)
     !!                                  = A(theta,k) * 2*pi*oma(f)/cga(f)
     !!                                                     dim=(nrng,nang)
     !!    ef1       R.A.  Public    O   1D Energy Density spectrum ef1(f)
     !!                                                     dim=(nrng)
     !!
     !!    dens1     R.A.  Public    O   large-scale Action Density (k,theta)
     !!                                                     dim=(nrng,nang)
     !!    dens2     R.A.  Public    O   Small-scale Action Density (k,theta)
     !!                                                     dim=(nrng,nang)
     !!    ------------------------------------------------------------------
     !!
     !!    for -tsa; The 2 returned arrays tsa & diag dim=(nrng,nang)
     !!    tsa       R.A.  Public    O   Snl-tsa = sumint + sumintsa
     !!    diag      R.A.  Public    O   Snl-tsa diagonal term = [dN/dn1]
     !!    ------------------------------------------------------------------
     !!
     !!    for -fbi; The 2 returned arrays fbi & diag2 dim=(nrng,nang)
     !!    fbi       R.A.  Public    O   Snl-fbi = sumint + sumintp  + sumintx
     !!    diag2     R.A.  Public    O   Snl-fbi diagonal term = [dN/dn1]
     !!    ------------------------------------------------------------------
     !!
     !! 4. Subroutines used :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!     STRACE    Subr. W3SERVMD Subroutine tracing.
     !!     optsa2    Subr. W3SERVMD Converts the 2D Energy Density (f,theta)
     !!     ------                   to Polar Action Density (k,theta) Norm. (in k)
     !!                              then splits it into large and small scale
     !!     snlr_tsa  Subr. W3SERVMD Computes dN(k,theta)/dt for TSA
     !!     --------                 due to wave-wave inter. (set itsa = 1)
     !!     snlr_fbi  Subr. W3SERVMD Computes dN(k,theta)/dt for FBI
     !!     --------                 due to wave-wave inter. (set itsa = 0)
     !!    ------------------------------------------------------------------
     !!
     !! 5. Called by :
     !!
     !!     Name      Type  Module   Description
     !!    ------------------------------------------------------------------
     !!     W3SRCE    Subr. w3srcemd Source term integration.
     !!     W3EXPO    Subr. ww3_outp Point output post-processor.
     !!     W3EXNC    Subr. ww3_ounp NetCDF Point output post-processor.
     !!     GXEXPO    Subr. gx_outp  GrADS  Point output post-processor.
     !!    ------------------------------------------------------------------
     !!
     !! 6. Error messages :
     !!
     !!      None.
     !!
     !! 7. Remarks :
     !!
     !! 8. Structure :
     !!
     !!    See source code.
     !!
     !! 9. Switches :
     !!
     !!    !/S   Enable subroutine tracing.
     !!
     !!10. Source code :
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     USE constants, ONLY: tpi
     USE w3gdatmd,  ONLY: nk,  nth,  xfr, dth,  sig, th, ecos, esin, &
          itsa, ialt
     !!    dimension: SIG(0:NK+1),TH(NTH), ECOS(NSPEC+NTH), ESIN(NSPEC+NTH)
     !!
     USE w3servmd, ONLY: extcde
     USE w3odatmd, ONLY: ndse, ndst, ndso
 #ifdef W3_S
     USE w3servmd, ONLY: strace
 #endif
     !!    ==================================================================
     !!
     IMPLICIT NONE
     !!
     !!    Parameter list
     !!    --------------
     REAL,    INTENT(IN)  :: A(NTH,NK), CG(NK), WN(NK), DEPTH
     REAL,    INTENT(OUT) :: S(NTH,NK), D(NTH,NK)
     !!
     LOGICAL, SAVE        :: FIRST_TSA = .true.
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!
     !!    Local Parameters & variables
     !!    -----------------------------
 #ifdef W3_S
     INTEGER, SAVE           :: IENT = 0
 #endif
     integer              :: irng, iang
     integer              :: nd3         !* bin # corresp. to ww3 dep
     real                 :: dep         !* depth (m), get it from WW3 DEPTH
     real                 :: wka(NK)     !* from WW3 WN(1:NK) corresp. to "DEPTH"
     real                 :: cga(NK)     !* from WW3 CG(1:NK) corresp. to "DEPTH"
     real                 :: pha(NK)     !* k*dk*dtheta array corresp. to "DEPTH"
     real                 :: fac         !* twopi*oma()/cga()
     real                 :: sum1        !* dummy variable
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     integer              :: npk         !* bin# at peak frequency fpk
     integer              :: npk2        !* bin# of second peak frequency
     integer              :: npk0        !* dummy int. used in the shuffle of npk's
     integer              :: nsep        !* min # of bins that separates npk & npk2
     !* set nsep = 2
     integer              :: npeaks      !* # of peaks (=0, 1, or 2)
     integer              :: nfs         !* bin# of freq. separation
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     integer              :: nbins       !* actual # of bins > npk  (incl. nfs) or
     !!                                        !* actual # of bins > npk2 (incl. nrng)
     !!                                        !* to guarantee a min 1 bin in equi. range
     real                 :: fpk         !* peak frequency (Hz)
     real                 :: fpk2        !* second peak frequency (Hz)
     real                 :: e1max       !* 1D energy at 1st peak 'fpk'
     real                 :: e1max2      !* 1D energy at 2nd peak 'fpk2'
     real                 :: sumd1       !* sum dens1+dens2 at nfs
     real                 :: sumd2       !* sum dens1+dens2 at nfs+1
     real                 :: densat1     !* averaged dens1  at nfs
     real                 :: densat2     !* averaged dens1  at nfs+1
     !!    --------------------------------------------------------------- &
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
 #ifdef W3_S
     CALL strace (ient, 'W3SNL4')
 #endif
     !!
     !!ini
     !!    Initialization of the output arrays
     !!    before calling TSA subroutines.
     !!    -----------------------------------
     s(:,:) = 0.0
     d(:,:) = 0.0
     !!ini---
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!
     IF ( first_tsa ) THEN
       !!
       !!
       !!-0    Set parameters & constants
       !!      ---------------------------
       nrng  = nk                     !* nrng = NK  must be odd   <---
       nzz   = (nk * (nk+1)) / 2      !* linear irng, krng
       nang  = nth                    !* nang = NTH must be even  <---
       na2p1 = nang/2 + 1             !* mid-angle or angle opposite to 1
       np2p1 = npts/2 + 1             !* mid-index of locus array
       twopi = tpi                    !* twopi = 8.*atan(1.)
       !* get it from WW3 TPI
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       !!-1    Allocate freq & angle related array declared as PUBLIC
       if ( allocated (frqa) )    deallocate (frqa)
       if ( allocated (oma) )     deallocate (oma)
       if ( allocated (angl) )    deallocate (angl)
       if ( allocated (sinan) )   deallocate (sinan)
       if ( allocated (cosan) )   deallocate (cosan)
       if ( allocated (dep_tbl) ) deallocate (dep_tbl)
       allocate(frqa(nrng))
       allocate(oma(nrng))
       allocate(angl(nang))
       allocate(sinan(nang))
       allocate(cosan(nang))
       allocate(dep_tbl(ndep))
       !!      ----------------------------------------------------------------
       !!
       !!-1a   Initialize frequency arrays and related parameters
       !!      --------------------------------------------------
       oma(:)   = sig(1:nk)           !* get it from WW3 SIG(1:NK)
       frqa(:)  = oma(:) / twopi
       f0       = frqa(1)
       dfrq     = xfr                 !* WW3 freq mult. for log freq
       !* get it from WW3 XFR
       !!      ----------------------------------------------------------------
       !!
       !!-1b   Initialize direction arrays and related parameters
       !!      --------------------------------------------------
       angl(:)  = th(1:nth)           !* get it from WW3 TH(1:NTH)
       cosan(:) = ecos(1:nth)         !* get it from WW3 ECOS(1:NTH)
       sinan(:) = esin(1:nth)         !* get it from WW3 ESIN(1:NTH)
       ainc     = dth                 !* WW3 angle increment (radians)
       !* get it from WW3 DTH
       !!      ----------------------------------------------------------------
       !!
       !!-1c   Define kzone & nb2fp
       !!kz
       !!      kzone = zone of freq influence, function of dfrq
       !!      for different values of x = 2,3,4 & 5
       !!      So,    kzone(x) = INT( alog(x)/alog(dfrq) )
       !!      +--------+----------+----------+----------+----------+
       !!      | dfrq   | kzone(2) | kzone(3) | kzone(4) | kzone(5) |
       !!      +--------+----------+----------+----------+----------+
       !!      | 1.05   |    14    |    22    |    28    |    33    |
       !!      +--------+----------+----------+----------+----------+
       !!      | 1.07   |    10    |    16    |    20    |    24    |
       !!      +--------+----------+----------+----------+----------+
       !!      | 1.10   |     7    |    11    |    14    |    17    |
       !!      +--------+----------+----------+----------+----------+
       kzone = int( alog(2.0)/alog(dfrq) )  !* Bash; faster without loss of accuracy
       !kz     kzone = INT( alog(3.0)/alog(dfrq) )  !* as in gridsetr & snlr_'s
       !kz     kzone = INT( alog(4.0)/alog(dfrq) )  !* as in gridsetr & snlr_'s
       !kz     kzone = INT( alog(5.0)/alog(dfrq) )  !* as in gridsetr & snlr_'s
       !!kz---
       !!
       !!op2
       !!      nb2fp = # of bins over fp (not incl. fp) - this depends on dfrq
       !!              so that (dfrq**nb2fp)*fp ~ 2.*fp  (like kzone(2))
       !!              used in 1 bin equi. range
       nb2fp = int( alog(2.0)/alog(dfrq) )  !* for equi. range near 2*fp
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - !!op2
       !!      ================================================================
       !!
       !!
       !!-2    Allocate gridsetr 11 look-up tables arrays
       !!      plus     pha_tbl array dim=(nrng,ndep) declared as PUBLIC
       if ( allocated (kref2_tbl) ) deallocate (kref2_tbl)
       if ( allocated (kref4_tbl) ) deallocate (kref4_tbl)
       allocate(kref2_tbl(npts,nang,nzz,ndep))
       allocate(kref4_tbl(npts,nang,nzz,ndep))
       !!
       if ( allocated (jref2_tbl) ) deallocate (jref2_tbl)
       if ( allocated (jref4_tbl) ) deallocate (jref4_tbl)
       allocate(jref2_tbl(npts,nang,nzz,ndep))
       allocate(jref4_tbl(npts,nang,nzz,ndep))
       !!
       if ( allocated (wtk2_tbl) ) deallocate (wtk2_tbl)
       if ( allocated (wtk4_tbl) ) deallocate (wtk4_tbl)
       allocate(wtk2_tbl(npts,nang,nzz,ndep))
       allocate(wtk4_tbl(npts,nang,nzz,ndep))
       !!
       if ( allocated (wta2_tbl) ) deallocate (wta2_tbl)
       if ( allocated (wta4_tbl) ) deallocate (wta4_tbl)
       allocate(wta2_tbl(npts,nang,nzz,ndep))
       allocate(wta4_tbl(npts,nang,nzz,ndep))
       !!
       if ( allocated (tfac2_tbl) ) deallocate (tfac2_tbl)
       if ( allocated (tfac4_tbl) ) deallocate (tfac4_tbl)
       allocate(tfac2_tbl(npts,nang,nzz,ndep))
       allocate(tfac4_tbl(npts,nang,nzz,ndep))
       !!
       if ( allocated (grad_tbl) ) deallocate (grad_tbl)
       allocate(grad_tbl(npts,nang,nzz,ndep))
       !!
       if ( allocated (pha_tbl) ) deallocate (pha_tbl)
       allocate(pha_tbl(nrng,ndep))
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       !!-3    Allocate gridsetr 11 returned arrays declared as PUBLIC
       if ( allocated (kref2) ) deallocate (kref2)
       if ( allocated (kref4) ) deallocate (kref4)
       allocate(kref2(npts,nang,nzz))
       allocate(kref4(npts,nang,nzz))
       !!
       if ( allocated (jref2) ) deallocate (jref2)
       if ( allocated (jref4) ) deallocate (jref4)
       allocate(jref2(npts,nang,nzz))
       allocate(jref4(npts,nang,nzz))
       !!
       if ( allocated (wtk2) ) deallocate (wtk2)
       if ( allocated (wtk4) ) deallocate (wtk4)
       allocate(wtk2(npts,nang,nzz))
       allocate(wtk4(npts,nang,nzz))
       !!
       if ( allocated (wta2) ) deallocate (wta2)
       if ( allocated (wta4) ) deallocate (wta4)
       allocate(wta2(npts,nang,nzz))
       allocate(wta4(npts,nang,nzz))
       !!
       if ( allocated (tfac2) ) deallocate (tfac2)
       if ( allocated (tfac4) ) deallocate (tfac4)
       allocate(tfac2(npts,nang,nzz))
       allocate(tfac4(npts,nang,nzz))
       !!
       if ( allocated (grad) ) deallocate (grad)
       allocate(grad(npts,nang,nzz))
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       !!-4    Allocate shloxr/shlocr 5 returned arrays declared as PUBLIC
       if ( allocated (wk2x) ) deallocate (wk2x)
       if ( allocated (wk2y) ) deallocate (wk2y)
       allocate(wk2x(npts))
       allocate(wk2y(npts))
       !!
       if ( allocated (wk4x) ) deallocate (wk4x)
       if ( allocated (wk4y) ) deallocate (wk4y)
       allocate(wk4x(npts))
       allocate(wk4y(npts))
       !!
       if ( allocated (ds) ) deallocate (ds)
       allocate(ds(npts))
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       !!-5    Allocate w3snlx/optsa2 2 shared arrays declared as PUBLIC
       if ( allocated (ef2) ) deallocate (ef2)
       if ( allocated (ef1) ) deallocate (ef1)
       allocate(ef2(nrng,nang))
       allocate(ef1(nrng))
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       !!-6    Allocate optsa2 2 returned arrays declared as PUBLIC
       if ( allocated (dens1) ) deallocate (dens1)
       if ( allocated (dens2) ) deallocate (dens2)
       allocate(dens1(nrng,nang))
       allocate(dens2(nrng,nang))
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       !!-7    Allocate snlr_??? 2 returned arrays declared as PUBLIC
       if ( itsa .eq. 1) then
         !!        allocate tsa, diag  used for -tsa
         if ( allocated (tsa) ) deallocate (tsa)
         if ( allocated (diag) ) deallocate (diag)
         allocate(tsa(nrng,nang))
         allocate(diag(nrng,nang))
       elseif ( itsa .eq. 0) then
         !!        allocate fbi, diag2 used for -fbi
         if ( allocated (fbi) ) deallocate (fbi)
         if ( allocated (diag2) ) deallocate (diag2)
         allocate(fbi(nrng,nang))
         allocate(diag2(nrng,nang))
       else
         write ( ndse,1000 ) itsa
         CALL extcde ( 115 )
       endif
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       !!-8    Get the 11 look-up table arrays by calling INSNL4
       !!      ----------------------------------------------------------------
       !!
       !!      ----------------------------------------------------------------
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       call insnl4
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!
       !!      it returns: 11 look-up tables arrays dim=(npts,nang,nzz,ndep)
       !!                  kref2_tbl, kref4_tbl, jref2_tbl, jref4_tbl,
       !!                  wtk2_tbl,  wtk4_tbl,  wta2_tbl,  wta4_tbl,
       !!                  tfac2_tbl, tfac4_tbl & grad_tbl
       !!                  plus       pha_tbl  dim=(nrng,ndep)
       !!                  and        dep_tbl  dim=(ndep)
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!
       first_tsa  = .false.
       !!
       !!
     ENDIF    !! IF ( FIRST_TSA ) THEN
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!    ##################################################################
     !!
     !!
     !!
     !!*i1 Map input ww3 "DEPTH" to "dep" and find corresp. depth bin # "nd3"
     !!    ------------------------------------------------------------------
     dep = depth                !* ww3 depth at a given time & loc.
     nd3 = minloc( abs(dep - dep_tbl(:)), dim=1 )
     !prt  print *, 'DEPTH, corresp depth bin # (nd3)  = ', DEPTH, nd3
     !!    ------------------------------------------------------------------
     !!
     !!
     !!*i2 Map from Look-up tables the 11 gridsetr arrays corresp. to "nd3"
     !!    kref2(:,:,:) -> grad(:,:,:) are used in subrs. "snlr_*"
     !!    ------------------------------------------------------------------
     kref2(:,:,:) = kref2_tbl(:,:,:,nd3)
     kref4(:,:,:) = kref4_tbl(:,:,:,nd3)
     jref2(:,:,:) = jref2_tbl(:,:,:,nd3)
     jref4(:,:,:) = jref4_tbl(:,:,:,nd3)
     wtk2(:,:,:)  = wtk2_tbl(:,:,:,nd3)
     wtk4(:,:,:)  = wtk4_tbl(:,:,:,nd3)
     wta2(:,:,:)  = wta2_tbl(:,:,:,nd3)
     wta4(:,:,:)  = wta4_tbl(:,:,:,nd3)
     tfac2(:,:,:) = tfac2_tbl(:,:,:,nd3)
     tfac4(:,:,:) = tfac4_tbl(:,:,:,nd3)
     grad(:,:,:)  = grad_tbl(:,:,:,nd3)
     pha(:)       = pha_tbl(:,nd3)
     !!    ------------------------------------------------------------------
     !!
     !!
     !!*i3 Map input ww3 arrays "WN(:)" & "CG(:)" to "wka(:)" & "cga(:)"
     !!    Note; Arrays wka(:) & cga(:) corresp to ww3 "DEPTH" & to be used in "optsa2"
     !!    ------------------------------------------------------------------
     wka(:)  = wn(1:nk)               !* Wavenumber array at ww3 "DEPTH"
     cga(:)  = cg(1:nk)               !* Group velocity array at ww3 "DEPTH"
     !!    ------------------------------------------------------------------
     !!
     !!
     !!*i4 Convert input WW3 2D Action Density spectrum "A(theta,k)"
     !!    to 2D Energy Density spectrum "ef2(theta,f)" & reverse indices
     !!    ==>  ef2(f,theta) = A(theta,k) * 2*pi*oma(f)/cga(f)
     !!    ------------------------------------------------------------------
     do irng=1,nrng
       fac = twopi*oma(irng)/cga(irng)
       do  iang=1,nang
         ef2(irng,iang) = a(iang,irng) * fac
       end do
     end do
     !!    ------------------------------------------------------------------
     !!
     !!
     !!*i5 Calculte the 1D Energy Density "ef1(f)"
     !!    ------------------------------------------------------------------
     do irng=1,nrng
       sum1 = 0.0
       do iang=1,nang
         sum1 = sum1 + ef2(irng,iang)
       end do
       ef1(irng) = sum1 * ainc
     end do
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     !!op2
     !!*   Bash;
     !!*   Find 1 or 2 peaks that satisfy TSA min condition (below) ------- *
     !!*   before calling TSA subrs. otherwise bailout (return) ----------- *
     !!*   Bailout & return with init. values of S & D = 0.0 -------------- *
     !!*   nsep   = min # of bins that separates between npk & npk2 (set=2) *
     !!*   nbins  = actual # of bins > npk  (incl. nfs)  -- or --           *
     !!*            actual # of bins > npk2 (incl. nrng)                    *
     !!*            to guarantee a min 1 bin in equi. range                 *
     !!*                                                                    *
     !!*   ===>  In case of just 1 peak the TSA min condition ------------- *  *****
     !!*   ===>  is relative to nrng and is satisfied when ---------------- *  <<<<<
     !!*   ===>  npk.le.nrng-1, to guarantee min 1 bin (incl nrng) > npk -- *  <<<<<
     !!*   ===>  we only need 1 bin in optsa2 to be in the equi. range ---- *  <<<<<
     !!*   ===>  skip if condition is not met ie if  npk.gt.nrng-1  ------- *  <<<<<
     !!*   ---------------------------------------------------------------- *
     !!*                                                                    *
     !!*   ===>  In case of 2 peaks the TSA min condition is applied twice: *  *****
     !!*   ===>   *1) at low freq peak (npk), *2) at high freq peak (npk2)  *  *****
     !!*   ===> *1) TSA min condition for the low freq peak (npk) --------- *  *****
     !!*   ===>  is relative to nfs and is satisfied when ----------------- *  <<<<<
     !!*   ===>  npk.le.nfs-1, to guarantee min 1 bin (incl nfs) > npk2 --- *  <<<<<
     !!*   ===>  we only need 1 bin in optsa2 to be in the equi. range ---- *  <<<<<
     !!*   ===>  skip if condition is not met ie if  npk.gt.nfs-1  -------- *  <<<<<
     !!*                                                                    *
     !!*   ===> *2) TSA min condition for the high freq peak (npk2) ------- *  *****
     !!*   ===>  is relative to nrng and is satisfied when ---------------- *  <<<<<
     !!*   ===>  npk2.le.nrng-1 to guarantee min 1 bin (incl nrng) > npk2   *  <<<<<
     !!*   ===>  we only need 1 bin in optsa2 to be in the equi. range ---- *  <<<<<
     !!*   ===>  skip if condition is not met ie if  npk2.gt.nrng-1 ------- *  <<<<<
     !!*   ---------------------------------------------------------------- *
     !!    ------------------------------------------------------------ !!op2
     !!    ==================================================================
     !!
     !!op2 ctd
     !!    First find the overall peak in ef1(:) with e1max must be > 0.000001
     !!    Starting from low freq. find the Energy max "e1max" and
     !!    corresp. peak freq. "fpk" and its freq. number "npk".
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     npk    = 0
     fpk    = 0.0
     e1max  = 0.0
     npeaks = 0
     !!    Look in the freq range that works for TSA call (see condition below)
     do irng=2,nrng-1             !* last peak loc. is at nrng-1      <<<<<
       !!      Pick the 1st local abs. max in [2,nrng-1] using (ef1(irng).gt.e1max)
       !!      so that if 2 equal adj. peaks are found it will pick the 1st e1max
       !!      encountered (i.e. the lower freq. one)
       !!      --------------------------------------------------------------!*  <<<<<
       if ( ef1(irng).gt.ef1(irng-1) .and. ef1(irng).gt.ef1(irng+1)  &
            .and. ef1(irng).gt.e1max ) then
         !!      --------------------------------------------------------------!*  <<<<<
         npk    = irng               !* update npk
         fpk    = frqa(npk)          !* update fpk
         e1max  = ef1(npk)           !* update e1max
         npeaks = 1
       endif
     end do
     !!    ------------------------------------------------------------------
     !!
     !!B   if a 1st peak is not found (npeaks=0 & e1max=0.0 < eps) or
     !!B   if a 1st peak is found with a tiny peak energy (e1max < eps) or
     !!B   if TSA min condition is not met rel. to nrng (npk.gt.nrng-1)   <<<<<
     !!B   this spectrum is Not suitable for tsa, so don't call tsa
     !!B   just return (don't stop) with init. values of S(:,:) and D(:,:)=0.0
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     if ( e1max.lt.0.000001 ) return
     !!    ------------------------------------------------------------ !!op2
     !!    ==================================================================
     !!
     !!
     !!op2 ctd
     !!    Bash; if we are here (i.e. we did not return) then we must
     !!    have found the 1st good peak (= overall peak) with e1max > eps
     !!
     !!    Now we look for a new 2nd peak that is at least 'nsep' bins away from
     !!    the 1st peak (nsep=2) (i.e. iabs(irng-npk).gt.nsep) before
     !!    calling new "optsa2" (the 2nd peak will have e1max2 < e1max)
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     nsep   = 2
     npk2   = 0
     fpk2   = 0.0
     e1max2 = 0.0
     !!    Again look in the freq range that is in line with TSA min condition
     !!    and find the 2nd highest peak with  eps < e1max2 < e1max
     do irng=2,nrng-1           !* last peak loc. is at nrng-1      <<<<<
       !!      Pick the 2nd local abs. max in [2,nrng-1] that is at least 'nsep'
       !!      bins away from the 1st peak using (ef1(irng).ge.e1max2) so that
       !!      if 2 equal adj. peaks are found it will pick the 2nd e1max2
       !!      encountered (i.e. the higher freq. one)
       !!      --------------------------------------------------------------!*  <<<<<
       if ( ef1(irng).gt.ef1(irng-1) .and. ef1(irng).gt.ef1(irng+1)  &
            .and. ef1(irng).ge.e1max2 .and. iabs(irng-npk).gt.nsep )  then
         !!      --------------------------------------------------------------!*  <<<<<
         npk2   = irng             !* update npk2
         fpk2   = frqa(npk2)       !* update fpk2
         e1max2 = ef1(npk2)        !* update e1max2
         npeaks = 2
       endif
     end do
     !!    ------------------------------------------------------------------
     !!
     !!B   if a 2nd peak is not found (npeaks=1 & e1max2=0.0 < eps)
     !!B   if a 2nd peak is found with a tiny peak energy (e1max2 < eps) or
     !!B   if TSA min condition is not met rel. to nrng (npk2.gt.nrng-1)   <<<<<
     !!B   This 2nd peak is not suitable for tsa, drop it and stay with just 1st peak.
     if ( e1max2.lt.0.000001 ) then
       npeaks = 1
       goto 200           !* skip the remaings tests goto 200
     endif
     !!    ------------------------------------------------------------ !!op2
     !!    ==================================================================
     !!
     !!
     !!
     !!
     !!op2 ctd
     if ( npeaks.eq.2 ) then
       !!-1    Shuffle the 2 peaks (if necessary) to keep npk to be always < npk2
       !!      This says nothing about which peak is the dominant peak
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       if ( npk2.lt.npk ) then
         npk0   = npk2
         npk2   = npk
         npk    = npk0                 !*  this way  npk < npk2  always
         fpk    = frqa(npk)
         fpk2   = frqa(npk2)
       endif
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!
       !!-2    here we have 2 peaks (npeaks=2) with  npk < npk2
       !!      find the freq. separation "nfs" (that divide the freq. regime into 2)
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       nfs = int( (npk+npk2)   / 2.0 )  !* take the lower  bin # to be nfs
       !b      nfs = INT ( (npk+npk2+1) / 2.0 )  !* take the higher bin # to be nfs
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     endif   !! if ( npeaks.eq.2 )
     !!
 200 continue
     !!    ------------------------------------------------------------ !!op2
     !!    ==================================================================
     !!
     !!
     !!    Bash; With the new "optsa2" you are allowed one call (if 1 peak)
     !!          or 2 calls (if 2 peaks) o account for spectra with double peaks.
     !!    Note; when nrmn=1 & nrmx=nrng ==> optsa2 = the old optsa
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!
     if ( npeaks.eq.1 ) then
       !!-1    one call to optsa2 for the whole freq. regime ( 1 --> nrng )
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       nbins = nrng - npk    !* # of bins in (npk, nrng] not incl. npk
       if ( nbins.gt.nb2fp ) nbins=nb2fp !* limit equi. range to ~2.0*fp
       !!      ----------------------------------------------------------------
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       call optsa2 ( 1,nrng,       npk, fpk,  nbins, wka, cga )
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      It returns variables dens1(nrng,nang) and dens2(nrng,nang)
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      ----------------------------------------------------------------
     endif   !! if ( npeaks.eq.1 )
     !!    ==================================================================
     !!
     !!
     if ( npeaks.eq.2 ) then
       !!
       !!-2    Now make two calls to new "optsa2" one for each freq regime.
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       nbins = nfs - npk     !* # of bins in (npk, nfs] not incl. npk
       if ( nbins.gt.nb2fp ) nbins=nb2fp  !* limit equi. range to ~2.0*fp
       !!      ---------------------------------------------------------------
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       call optsa2 ( 1,nfs,        npk, fpk,  nbins, wka, cga )
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      It returns variables dens1(nrng,nang) and dens2(nrng,nang)
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      ----------------------------------------------------------------
       !!
       nbins = nrng - npk2   !* # of bins in (npk2, nrng] not incl. npk2
       if ( nbins.gt.nb2fp ) nbins=nb2fp  !* limit equi. range to ~2.0*fp
       !!      ----------------------------------------------------------------
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       call optsa2 ( nfs+1,nrng,   npk2,fpk2, nbins, wka, cga )
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      It returns variables dens1(nrng,nang) and dens2(nrng,nang)
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      ----------------------------------------------------------------
       !!      ================================================================
       !!
       !!-3    Remove the step like jump (if exists) in dens1() between nfs & nfs+1
       do iang=1,nang
         sumd1 = dens1(nfs,iang)   + dens2(nfs,iang)   !* sum at nfs
         sumd2 = dens1(nfs+1,iang) + dens2(nfs+1,iang) !* sum at nfs+1
         !!
         !!        do 3 bin average for dens1() at nfs   and store in densat1
         densat1 = ( dens1(nfs-1,iang) + dens1(nfs,iang) +           &
              dens1(nfs+1,iang) ) / 3.
         !!        do 3 bin average for dens1() at nfs+1 and store in densat2
         densat2 = ( dens1(nfs,iang)   + dens1(nfs+1,iang) +         &
              dens1(nfs+2,iang) ) / 3.
         !!
         !!        subtitute back into dens1(nfs,iang) & dens1(nfs+1,iang)
         dens1(nfs,iang)   = densat1             ! dens1 at nfs
         dens1(nfs+1,iang) = densat2             ! dens1 at nfs+1
         !!
         !!        recalculate dens2(nfs,iang) & dens2(nfs+1,iang)
         dens2(nfs,iang)   = sumd1 - densat1     ! dens2 at nfs
         dens2(nfs+1,iang) = sumd2 - densat2     ! dens2 at nfs+1
       end do
       !!
     endif   !! if ( npeaks.eq.2 )
     !!    ==================================================================
     !!
     !!
     !!
     !!    -----------------------------------------------------------------#
     !!                                                                     !
     !!    Get Snl source term and its diagonal term  from "snlr"           !
     !!    for -tsa  only       use "snlr_tsa"       itsa = 1               !
     !!    for -fbi  only       use "snlr_fbi"       itsa = 0               !
     !!                                                                     !
     !!    -----------------------------------------------------------------#
     !!
     !!
     if ( itsa .eq. 1) then
       !!
       !!      ----------------------------------------------------------------
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       call snlr_tsa ( pha, ialt )
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      It returns tsa(nrng,nang) & diag(nrng,nang)
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      ----------------------------------------------------------------
       !!
       !!      Pack results in proper format ---------------------------------- *
       !!      S() & D() arrays are to be returned to WW3 in (k,theta) space
       do irng=1,nrng
         do iang=1,nang
           !!        Convert the Norm. (in k) Polar tsa(k,theta) to Polar S(theta,k)
           !!        and  reverse indices back to (iang,irng) as in WW3
           s(iang,irng) = tsa(irng,iang) * wka(irng)   !* <=============
           d(iang,irng) = diag(irng,iang)
           !!        ---------------------------
         end do
       end do
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!
       !!
     elseif ( itsa .eq. 0) then
       !!
       !!
       !!      ----------------------------------------------------------------
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       call snlr_fbi ( pha, ialt )
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      It returns fbi(nrng,nang) & diag2(nrng,nang)
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!      ----------------------------------------------------------------
       !!
       !!      Pack results in proper format ---------------------------------- *
       !!      S() & D() arrays are to be returned to WW3 in (k,theta) space
       do irng=1,nrng
         do iang=1,nang
           !!        Convert the Norm. (in k) Polar fbi(k,theta) to Polar S(theta,k)
           !!        and  reverse indices back to (iang,irng) as in WW3
           s(iang,irng) = fbi(irng,iang) * wka(irng)   !* <=============
           d(iang,irng) = diag2(irng,iang)
           !!        --------------------------------
         end do
       end do
       !!      -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!
     else
       !!
       write( ndse,1000 ) itsa
       CALL extcde ( 130 )
       !!       --- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !!
     endif
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     RETURN
     !!
 1000 format ( ' W3SNL4 Error : Bad itsa value ',i4)
     !!

References ainc, angl, cosan, dens1, dens2, dep_tbl, dfrq, diag, diag2, ds, w3gdatmd::dth, w3gdatmd::ecos, ef1, ef2, w3gdatmd::esin, w3servmd::extcde(), f0, fbi, frqa, grad, grad_tbl, w3gdatmd::ialt, insnl4(), w3gdatmd::itsa, jref2, jref2_tbl, jref4, jref4_tbl, kref2, kref2_tbl, kref4, kref4_tbl, kzone, na2p1, nang, nb2fp, ndep, w3odatmd::ndse, w3odatmd::ndso, w3odatmd::ndst, w3gdatmd::nk, np2p1, npts, nrng, w3gdatmd::nth, nzz, oma, optsa2(), pha_tbl, w3gdatmd::sig, sinan, snlr_fbi(), snlr_tsa(), w3servmd::strace(), tfac2, tfac2_tbl, tfac4, tfac4_tbl, w3gdatmd::th, constants::tpi, tsa, twopi, wk2x, wk2y, wk4x, wk4y, wta2, wta2_tbl, wta4, wta4_tbl, wtk2, wtk2_tbl, wtk4, wtk4_tbl, and w3gdatmd::xfr.

Referenced by gxexpo(), w3exnc(), w3expo(), and w3srcemd::w3srce().

◆ wkfnc()

real function w3snl4md::wkfnc	(	real, intent(in)	f,
		real, intent(in)	dep
	)

Calculate the Wavenumber k (rad/m) as function of frequency 'f' (Hz) and depth 'dep' (m).

  Using what looks like a "Pade approximation" of an inversion
  of the linear wave dispersion relation.

      sigma^2 = gk*tanh(kd),  sigma = 2*pi*f
      Wavenumber k (rad/m) is returned in "wkfnc"

Parameters

f	Frequency (Hz).
dep	Depth (m).

Returns: wkfnc Wavenumber 'k'

Author: Bash Toulany; Michael Casey; William Perrie

Date: 12-Apr-2016

Definition at line 5299 of file w3snl4md.F90.

     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III                 BIO |
     !/                  |           Bash Toulany            |
     !/                  |           Michael Casey           |
     !/                  |           William Perrie          |
     !/                  |                        FORTRAN 90 |
     !/                  | Last update :         12-Apr-2016 |
     !/                  +-----------------------------------+
     !/
     !/    01-Mar-2016 : Origination.                        ( version 5.13 )
     !/
     !!
     !!    it returns: wavenumber 'k' in wkfnc
     !!
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     !! 1. Purpose :
     !!
     !!    Calculate the Wavenumber k (rad/m) as function of
     !!    frequency 'f' (Hz) and depth 'dep' (m).
     !!
     !! 2. Method :
     !!
     !!    Using what looks like a "Pade approximation" of an inversion
     !!    of the linear wave dispersion relation.
     !!        sigma^2 = gk*tanh(kd),  sigma = 2*pi*f
     !!        Wavenumber k (rad/m) is returned in "wkfnc"
     !!
     !! 3. Parameters :
     !!
     !!    Parameter list
     !!    ------------------------------------------------------------------
     !!    Name     Type   Scope    I/O  Description
     !!    ------------------------------------------------------------------
     !!    twopi     Real  Public    I   = TPI; WW3 2*pi=8.*atan(1.) (radians)
     !!    ------------------------------------------------------------------
     !!
     !!    --------------------------------------------------------------- &
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ----------------------------------------------------------------72
     !!    ==================================================================
     !!
     !!
     IMPLICIT NONE
     !!
     !!    Parameter list
     !!    --------------
     real,    intent(in)  :: f, dep
     !!
     !!    Local variables
     !!    ---------------
     real(KIND=8)         :: g, y, x
     !!    -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     !!    ---------------------::-----------------------------------------72
     !!    ##################################################################
     !!------------------------------------------------------------------------------
     !!==============================================================================
     !!
     !!
     g     = 9.806                   !* set = GRAV as in CONSTANTS
     !!
     y     = ( (twopi*f)**2 ) * dep / g   !* sigma^2 d/g
     !!
     !!    --------------------------------------------------------------- &
     x     = y * ( y +                                               &
          1./(1.00000+y*(0.66667+y*(0.35550+y*(0.16084+y*(0.06320   &
          +y*(0.02174+y*(0.00654+y*(0.00171+y*(0.00039+y*0.00011)   &
          )))))))))
     !!    --------------------------------------------------------------- &
     !!
     x     = sqrt(x)                 !* kd
     !!
     wkfnc = x / dep                 !* k
     !!    ------------------------------------------------------------------
     !!    ==================================================================
     !!
     RETURN
     !!

References twopi.

Referenced by insnl4(), and shloxr().

Variable Documentation

◆ ainc

real w3snl4md::ainc

Definition at line 142 of file w3snl4md.F90.

142 real :: ainc, twopi

Referenced by gridsetr(), insnl4(), optsa2(), and w3snl4().

◆ angl

real, dimension(:), allocatable w3snl4md::angl

Definition at line 144 of file w3snl4md.F90.

144 real, allocatable, dimension(:) :: angl, sinan, cosan

Referenced by w3snl4().

◆ cosan

real, dimension(:), allocatable w3snl4md::cosan

Definition at line 144 of file w3snl4md.F90.

Referenced by gridsetr(), and w3snl4().

◆ dens1

real, dimension(:,:), allocatable w3snl4md::dens1

Definition at line 188 of file w3snl4md.F90.

188 real, allocatable, dimension(:,:) :: dens1, dens2

Referenced by optsa2(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ dens2

real, dimension(:,:), allocatable w3snl4md::dens2

Definition at line 188 of file w3snl4md.F90.

Referenced by optsa2(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ dep_tbl

real, dimension(:), allocatable w3snl4md::dep_tbl

Definition at line 145 of file w3snl4md.F90.

145 real, allocatable, dimension(:) :: dep_tbl

Referenced by insnl4(), and w3snl4().

◆ dfrq

real w3snl4md::dfrq

Definition at line 141 of file w3snl4md.F90.

141 real :: dfrq, f0

Referenced by gridsetr(), insnl4(), and w3snl4().

◆ diag

real, dimension(:,:), allocatable w3snl4md::diag

Definition at line 195 of file w3snl4md.F90.

Referenced by snlr_tsa(), and w3snl4().

◆ diag2

real, dimension(:,:), allocatable w3snl4md::diag2

Definition at line 196 of file w3snl4md.F90.

Referenced by snlr_fbi(), and w3snl4().

◆ ds

real, dimension(:), allocatable w3snl4md::ds

Definition at line 173 of file w3snl4md.F90.

Referenced by gridsetr(), shlocr(), shloxr(), and w3snl4().

◆ ef1

real, dimension(:), allocatable w3snl4md::ef1

Definition at line 181 of file w3snl4md.F90.

181 real, allocatable, dimension(:) :: ef1

Referenced by optsa2(), and w3snl4().

◆ ef2

real, dimension(:,:), allocatable w3snl4md::ef2

Definition at line 180 of file w3snl4md.F90.

180 real, allocatable, dimension(:,:) :: ef2

Referenced by optsa2(), and w3snl4().

◆ f0

real w3snl4md::f0

Definition at line 141 of file w3snl4md.F90.

Referenced by gridsetr(), and w3snl4().

◆ fbi

real, dimension(:,:), allocatable w3snl4md::fbi

Definition at line 196 of file w3snl4md.F90.

196 real, allocatable, dimension(:,:) :: fbi, diag2

Referenced by snlr_fbi(), and w3snl4().

◆ frqa

real, dimension(:), allocatable w3snl4md::frqa

Definition at line 143 of file w3snl4md.F90.

143 real, allocatable, dimension(:) :: frqa, oma

Referenced by gridsetr(), insnl4(), optsa2(), and w3snl4().

◆ grad

real, dimension(:,:,:), allocatable w3snl4md::grad

Definition at line 167 of file w3snl4md.F90.

167 real, allocatable, dimension(:,:,:) :: grad

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ grad_tbl

real, dimension(:,:,:,:), allocatable w3snl4md::grad_tbl

Definition at line 156 of file w3snl4md.F90.

156 real, allocatable, dimension(:,:,:,:) :: grad_tbl

Referenced by insnl4(), and w3snl4().

◆ ismo

integer, parameter w3snl4md::ismo = 1

Definition at line 126 of file w3snl4md.F90.

126 integer, parameter :: ismo = 1 !* = 1 do smooth in interp2

Referenced by interp2().

◆ jref2

integer, dimension(:,:,:), allocatable w3snl4md::jref2

Definition at line 163 of file w3snl4md.F90.

163 integer, allocatable, dimension(:,:,:) :: jref2, jref4

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ jref2_tbl

integer, dimension(:,:,:,:), allocatable w3snl4md::jref2_tbl

Definition at line 152 of file w3snl4md.F90.

152 integer, allocatable, dimension(:,:,:,:) :: jref2_tbl, jref4_tbl

Referenced by insnl4(), and w3snl4().

◆ jref4

integer, dimension(:,:,:), allocatable w3snl4md::jref4

Definition at line 163 of file w3snl4md.F90.

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ jref4_tbl

integer, dimension(:,:,:,:), allocatable w3snl4md::jref4_tbl

Definition at line 152 of file w3snl4md.F90.

Referenced by insnl4(), and w3snl4().

◆ kref2

integer, dimension(:,:,:), allocatable w3snl4md::kref2

Definition at line 162 of file w3snl4md.F90.

162 integer, allocatable, dimension(:,:,:) :: kref2, kref4

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ kref2_tbl

integer, dimension(:,:,:,:), allocatable w3snl4md::kref2_tbl

Definition at line 151 of file w3snl4md.F90.

151 integer, allocatable, dimension(:,:,:,:) :: kref2_tbl, kref4_tbl

Referenced by insnl4(), and w3snl4().

◆ kref4

integer, dimension(:,:,:), allocatable w3snl4md::kref4

Definition at line 162 of file w3snl4md.F90.

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ kref4_tbl

integer, dimension(:,:,:,:), allocatable w3snl4md::kref4_tbl

Definition at line 151 of file w3snl4md.F90.

Referenced by insnl4(), and w3snl4().

◆ kzone

integer w3snl4md::kzone

Definition at line 138 of file w3snl4md.F90.

Referenced by gridsetr(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ na2p1

integer w3snl4md::na2p1

Definition at line 139 of file w3snl4md.F90.

Referenced by w3snl4().

◆ nang

integer w3snl4md::nang

Definition at line 139 of file w3snl4md.F90.

139 integer :: nang, na2p1

Referenced by gridsetr(), insnl4(), and w3snl4().

◆ nb2fp

integer w3snl4md::nb2fp

Definition at line 138 of file w3snl4md.F90.

Referenced by w3snl4().

◆ ndep

integer, parameter w3snl4md::ndep = 37

Definition at line 131 of file w3snl4md.F90.

131 integer, parameter :: ndep = 37 !* # of depths in look-up tables

Referenced by insnl4(), and w3snl4().

◆ np2p1

integer w3snl4md::np2p1

Definition at line 140 of file w3snl4md.F90.

140 integer :: np2p1

Referenced by gridsetr(), shlocr(), shloxr(), and w3snl4().

◆ npts

integer, parameter w3snl4md::npts = 30

Definition at line 129 of file w3snl4md.F90.

129 integer, parameter :: npts = 30 !* # of points on the locus

Referenced by gridsetr(), insnl4(), shlocr(), shloxr(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ nrng

integer w3snl4md::nrng

Definition at line 138 of file w3snl4md.F90.

138 integer :: nrng, nzz, kzone, nb2fp

Referenced by w3snl4().

◆ nzz

integer w3snl4md::nzz

Definition at line 138 of file w3snl4md.F90.

Referenced by w3snl4().

◆ oma

real, dimension(:), allocatable w3snl4md::oma

Definition at line 143 of file w3snl4md.F90.

Referenced by insnl4(), optsa2(), and w3snl4().

◆ pha_tbl

real, dimension(:,:), allocatable w3snl4md::pha_tbl

Definition at line 157 of file w3snl4md.F90.

157 real, allocatable, dimension(:,:) :: pha_tbl

Referenced by insnl4(), and w3snl4().

◆ sinan

real, dimension(:), allocatable w3snl4md::sinan

Definition at line 144 of file w3snl4md.F90.

Referenced by gridsetr(), optsa2(), and w3snl4().

◆ tfac2

real, dimension(:,:,:), allocatable w3snl4md::tfac2

Definition at line 166 of file w3snl4md.F90.

166 real, allocatable, dimension(:,:,:) :: tfac2, tfac4

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ tfac2_tbl

real, dimension(:,:,:,:), allocatable w3snl4md::tfac2_tbl

Definition at line 155 of file w3snl4md.F90.

155 real, allocatable, dimension(:,:,:,:) :: tfac2_tbl, tfac4_tbl

Referenced by insnl4(), and w3snl4().

◆ tfac4

real, dimension(:,:,:), allocatable w3snl4md::tfac4

Definition at line 166 of file w3snl4md.F90.

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ tfac4_tbl

real, dimension(:,:,:,:), allocatable w3snl4md::tfac4_tbl

Definition at line 155 of file w3snl4md.F90.

Referenced by insnl4(), and w3snl4().

◆ tsa

real, dimension(:,:), allocatable w3snl4md::tsa

Definition at line 195 of file w3snl4md.F90.

195 real, allocatable, dimension(:,:) :: tsa, diag

Referenced by snlr_tsa(), and w3snl4().

◆ twopi

real w3snl4md::twopi

Definition at line 142 of file w3snl4md.F90.

Referenced by cgfnc(), gridsetr(), optsa2(), shloxr(), w3snl4(), and wkfnc().

◆ wk2x

real, dimension(:), allocatable w3snl4md::wk2x

Definition at line 172 of file w3snl4md.F90.

172 real, allocatable, dimension(:) :: wk2x, wk2y

Referenced by gridsetr(), shlocr(), shloxr(), and w3snl4().

◆ wk2y

real, dimension(:), allocatable w3snl4md::wk2y

Definition at line 172 of file w3snl4md.F90.

Referenced by gridsetr(), shlocr(), shloxr(), and w3snl4().

◆ wk4x

real, dimension(:), allocatable w3snl4md::wk4x

Definition at line 173 of file w3snl4md.F90.

173 real, allocatable, dimension(:) :: wk4x, wk4y, ds

Referenced by gridsetr(), shlocr(), shloxr(), and w3snl4().

◆ wk4y

real, dimension(:), allocatable w3snl4md::wk4y

Definition at line 173 of file w3snl4md.F90.

Referenced by gridsetr(), shlocr(), shloxr(), and w3snl4().

◆ wta2

real, dimension(:,:,:), allocatable w3snl4md::wta2

Definition at line 165 of file w3snl4md.F90.

165 real, allocatable, dimension(:,:,:) :: wta2, wta4

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ wta2_tbl

real, dimension(:,:,:,:), allocatable w3snl4md::wta2_tbl

Definition at line 154 of file w3snl4md.F90.

154 real, allocatable, dimension(:,:,:,:) :: wta2_tbl, wta4_tbl

Referenced by insnl4(), and w3snl4().

◆ wta4

real, dimension(:,:,:), allocatable w3snl4md::wta4

Definition at line 165 of file w3snl4md.F90.

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ wta4_tbl

real, dimension(:,:,:,:), allocatable w3snl4md::wta4_tbl

Definition at line 154 of file w3snl4md.F90.

Referenced by insnl4(), and w3snl4().

◆ wtk2

real, dimension(:,:,:), allocatable w3snl4md::wtk2

Definition at line 164 of file w3snl4md.F90.

164 real, allocatable, dimension(:,:,:) :: wtk2, wtk4

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ wtk2_tbl

real, dimension(:,:,:,:), allocatable w3snl4md::wtk2_tbl

Definition at line 153 of file w3snl4md.F90.

153 real, allocatable, dimension(:,:,:,:) :: wtk2_tbl, wtk4_tbl

Referenced by insnl4(), and w3snl4().

◆ wtk4

real, dimension(:,:,:), allocatable w3snl4md::wtk4

Definition at line 164 of file w3snl4md.F90.

Referenced by gridsetr(), insnl4(), snlr_fbi(), snlr_tsa(), and w3snl4().

◆ wtk4_tbl

real, dimension(:,:,:,:), allocatable w3snl4md::wtk4_tbl

Definition at line 153 of file w3snl4md.F90.

Referenced by insnl4(), and w3snl4().

Functions/Subroutines

Variables

Detailed Description

Function/Subroutine Documentation

◆ cgfnc()

◆ cplshr()

◆ gridsetr()

◆ insnl4()

◆ interp2()

◆ optsa2()

◆ shlocr()

◆ shloxr()

◆ snlr_fbi()

◆ snlr_tsa()

◆ w3snl4()

◆ wkfnc()

Variable Documentation

◆ ainc

◆ angl

◆ cosan

◆ dens1

◆ dens2

◆ dep_tbl

◆ dfrq

◆ diag

◆ diag2

◆ ds

◆ ef1

◆ ef2

◆ f0

◆ fbi

◆ frqa

◆ grad

◆ grad_tbl

◆ ismo

◆ jref2

◆ jref2_tbl

◆ jref4

◆ jref4_tbl

◆ kref2

◆ kref2_tbl

◆ kref4

◆ kref4_tbl

◆ kzone

◆ na2p1

◆ nang

◆ nb2fp

◆ ndep

◆ np2p1

◆ npts

◆ nrng

◆ nzz

◆ oma

◆ pha_tbl

◆ sinan

◆ tfac2

◆ tfac2_tbl

◆ tfac4

◆ tfac4_tbl

◆ tsa

◆ twopi

◆ wk2x

◆ wk2y

◆ wk4x

◆ wk4y

◆ wta2

◆ wta2_tbl

◆ wta4

◆ wta4_tbl

◆ wtk2

◆ wtk2_tbl

◆ wtk4

◆ wtk4_tbl