WW3/w3sic5md_8F90_source.html

#include "w3macros.h"

!/ ------------------------------------------------------------------- /

MODULE w3sic5md

  !/

  !/                  +-----------------------------------+

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |           Q. Liu                  |

  !/                  |           E. Rogers               |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         19-May-2021 |

  !/                  +-----------------------------------+

  !/

  !/    15-Mar-2016 : Origination.                        ( version 5.10 )

  !/                                                      ( Q. Liu )

  !/    15-Mar-2016 : Started from w3sic1/2/3/4 module    ( Q. Liu )

  !/

  !/    24-Apr-2017 : Adding more filters                 ( Q. Liu )

  !/

  !/    29-Apr-2017 : Introducing CMPLX_TANH2             ( Q. Liu )

  !/

  !/    02-Jun-2017 : Update to version 5.16              ( Q. Liu )

  !/

  !/    17-Jun-2017 : Remove some unnecessary lines       ( Q. Liu )

  !/                 (cg_ice, detla function, complx_tanh etc.)

  !/

  !/    20-Aug-2018 : Ready to be merged to master (v6.06)( Q. Liu)

  !/

  !/    19-May-2021 : Incl. the RP and M2 model           ( Q. Liu)

  !/

  !/ 1. Purpose :

  !     Calculate ice source term S_{ice} according to different ice

  !     models:

  !     * 'FS': the viscoelastic, extended Fox and Squire sea ice model

  !             (Mosig et al. 2015)

  !     * 'RP': the viscoelastic, Robinson and Palmer model (Mosig et al.

  !             2015)

  !     * 'M2': the order 3 power law model proposed by Meylan et al.

  !             (2018)

  !

  !     Reference:

  !     Mosig, J.E.M., F. Montiel, and V. A. Squire (2015):

  !     Comparison of viscoelastic-type models for ocean wave attenuation

  !     in ice-covered seas, J. Geophys. Res. Oceans, 120, 6072–6090,

  !     doi:10.1002/2015JC010881.

  !

  !     Meylan, M.H., L. Bennetts, J. Mosig, W. Rogers, M. Doble, and

  !     M. Peter (2018): Dispersion relations, power laws, and energy loss

  !     for waves in the marginal ice zone. J. Geophys. Res. Oceans, 123,

  !     3322–3335, https://doi.org/10.1002/2018JC013776.

  !

  !     Liu, Q., W. E. Rogers, A. Babanin, J. Li, and C. Guan (2020):

  !     Spectral Modeling of Ice-Induced Wave Decay. J. Phys. Oceanogr.,

  !     50 (6), 1583–1604.

  !

  !  2. Variables and types :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !      KSP       Int.  Private  the kind parameter for single precision

  !                               real variables

  !      KDP       Int.  Private  Same as KSP but for double precision

  !      KSPC      Int.  Private  the kind parameter for single precision

  !                               complex variables

  !      KDPC      Int.  Private  Same as KSPC but for double precision

  !      ERRTOL    Real  Private  A real parameter used for "==" test

  !     ----------------------------------------------------------------

  !

  !  3. Subroutines and functions :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !      W3SIC5    Subr. Public   Ice source term

  !      W3IC5WNCG Subr. Public   Wavenumber and group velocity of ice-

  !                               coupled waves

  !      FSDISP    Subr. Public   Solving the ice-coupled wave dispersion

  !      BALANCING_MATRIX

  !                Subr. Private  Balancing the matrix before we try to

  !                               find its eigenvalues

  !      EIG_HQR   Subr. Private  QR algorithm for real Hessenberg matrix

  !                               (eigenvalues-finding algorithm)

  !      POLYROOTS Subr. Private  Finding roots of a general polynomial

  !      NR_CORR   Func. Private  Get the Newton-Raphson correction term

  !                               for iteration

  !      NR_ROOT   Func. Private  Newton-Raphson algorithm for solving

  !                               the ice-coupled wave dispersion

  !      CMPLX_SINH, CMPLX_COSH, CMPLX_TANH2

  !                Func. Private  sinh, cosh, tanh for complex inputs

  !      INIT_RANDOM_SEED

  !                Subr. Private  Initialize the random seed based on

  !                               the system's time

  !     ----------------------------------------------------------------

  !

  !  4. Subroutines and functions used :

  !      See subroutine documentation

  !

  !  5. Remarks :

  !

  !  6. Switches :

  !      See subroutine documentation

  !

  !  7. Source code:

  !/

  !/ ------------------------------------------------------------------- /

  IMPLICIT NONE

  !/

  PUBLIC  :: w3sic5, w3ic5wncg, fsdisp

  PRIVATE :: balancing_matrix, eig_hqr, polyroots

  PRIVATE :: nr_corr, nr_root

  PRIVATE :: cmplx_sinh, cmplx_cosh, cmplx_tanh2

  PRIVATE :: init_random_seed

  !/

  PRIVATE :: ksp, kdp, kspc, kdpc, errtol

  !/ ------------------------------------------------------------------- /

  !/ Parameter list

  !     Kind for single- and double-precision real type

  INTEGER, PARAMETER :: KSP    = kind(1.0)

  INTEGER, PARAMETER :: KDP    = kind(1.0d0)

  !

  !     Kind for single- and double-precision complex type

  INTEGER, PARAMETER :: KSPC   = kind((1.0, 1.0))

  INTEGER, PARAMETER :: KDPC   = kind((1.0d0, 1.0d0))

  REAL, PARAMETER    :: ERRTOL = 1.e-12

  !/

CONTAINS

  !/ ------------------------------------------------------------------- /

  !/

  SUBROUTINE w3sic5 (A, DEPTH, CG, WN, IX, IY, S, D)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |           E. Rogers               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         19-May-2021 |

    !/                  +-----------------------------------+

    !/

    !/    23-Mar-2016 : Origination                         ( version 5.10 )

    !/                                                      ( Q. Liu)

    !/    23-Mar-2016 : Started from w3sic1/2/3/4 subr.     ( Q. Liu)

    !/    05-Apr-2016 : Options for Cg_{ice} or Cg          ( Q. Liu)

    !/    25-Apr-2017 : Add more filters                    ( Q. Liu)

    !/    20-Aug-2018 : Ready to be merged to master (v6.06)( Q. Liu)

    !/

    !/    Copyright 2009 National Weather Service (NWS),

    !/       National Oceanic and Atmospheric Administration.  All rights

    !/       reserved.  WAVEWATCH III is a trademark of the NWS.

    !/       No unauthorized use without permission.

    !/

    !/ 1. Purpose :

    !     Calculate ice source term S_{ice} according to 3 different sea ice

    !     models (Mosig et al. 2015, Meylan et al. 2018, Liu et al. 2020)

    !

    !  2. Method :

    !     Regarding i/o (general to all Sice modules): S_{ice} source term

    !     is calculated using up to 5 parameters read from input files.

    !     These parameters are allowed to vary in space and time.

    !     The parameters control the exponential decay rate k_i.

    !     Since there are 5 parameters, this permits description of

    !     dependence of k_i on frequency or wavenumber.

    !

    !     Sea ice affects the wavenumber k of wind-generated ocean waves.

    !     The ice-modified wavenumber can be expressed as a complex number

    !     k = k_r + i * k_i, with the real part k_r representing impact of

    !     the sea ice on the physical wavelength and propagation speeds,

    !     producing something analogous to shoaling and refraction by

    !     bathymetry, whereas the imaginary part of the complex

    !     wavenumber, k_i, is an exponential decay coefficient

    !     k_i(x,y,t,sigma) (depending on location, time and frequency,

    !     respectively), representing wave attenuation, and can be

    !     introduced in a wave model such as WW3 as S_ice/E=-2*Cg*k_i,

    !     where S_ice is one of several dissipation mechanisms, along

    !     with whitecapping, for example, S_ds=S_wc+S_ice+⋯. The k_r -

    !     modified by ice would enter the model via the C calculations

    !     on the left-hand side of the governing equation.The fundamentals

    !     are straightforward, e.g. Rogers and Holland (2009 and

    !     subsequent unpublished work) modified a similar model, SWAN

    !     (Booij et al. 1999) to include the effects of a viscous mud

    !     layer using the same approach (k = k_r + i*k_i) previously.

    !

    !     General approach is analogous to Rogers and Holland (2009)

    !         approach for mud.

    !     See text near their eq. 1 :

    !       k        = k_r  +  i * k_i

    !       eta(x,t) = Real( a * exp( i * ( k * x - sigma * t ) ) )

    !       a        = a0 * exp( -k_i * x )

    !       S / E    = -2 * Cg * k_i (see also Komen et al. (1994, pg. 170)

    !

    !     Following W3SBT1 as a guide, equation 1 of W3SBT1 says:

    !         S = D * E

    !     However, the code of W3SBT1 has

    !         S = D * A

    !     This leads me to believe that the calling routine is

    !         expecting "S/sigma" not "S"

    !     Thus we will use D = S/E = -2 * Cg * k_i

    !        (see also documentation of W3SIC1)

    !

    !     Notes regarding numerics:

    !     -------------------------

    !     Experiments with constant k_i values suggest that results may be

    !        dependent on resolution if insufficient resolution is used.

    !        For detailed information, see documentation of W3SIC1.

    !

    !     Note regarding applicability/validity:

    !     --------------------------------------

    !     Similar to the Wang and Shen model used in w3sic3md, the 3 models

    !     used here are empirical medium models as well which treat the sea

    !     ice cover as a continuum and use 1/2 empirical rheological para-

    !     meters, i.e., the effective shear modulus of ice G and the effec-

    !     tive viscosity η to characterize sea ices of various type. Please

    !     see the documentation of w3sic3md for a detailed discussion of

    !     this kind of model.

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       A       R.A.  I   Action density spectrum (1-D).

    !       DEPTH   Real  I   Local water depth.

    !       CG      R.A.  I   Group velocities.

    !       WN      R.A.  I   Wavenumbers

    !       IX,IY   I.S.  I   Grid indices.

    !       S       R.A.  O   Source term (1-D version).

    !       D       R.A.  O   Diagonal term of derivative (1-D version).

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing (!/S switch).

    !      PRT2DS    Subr. W3ARRYMD Print plot output (!/T0 switch).

    !      OUTMAT    Subr. W3ARRYMD Matrix output (!/T1 switch).

    !      W3IC5WNCG Subr. /        Wavenumber and group velocity of ice-

    !                               coupled waves

    !     ----------------------------------------------------------------

    !      * / means this module

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3SRCE    Subr. W3SRCEMD Source term integration.

    !      W3EXPO    Subr.   N/A    ASCII Point output post-processor.

    !      W3EXNC    Subr.   N/A    NetCDF Point output post-processor.

    !      GXEXPO    Subr.   N/A    GrADS point output post-processor.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     None.

    !

    !  7. Remarks :

    !

    !     If ice parameter 1 is zero, no calculations are made.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !     !/T   Enable general test output.

    !     !/T0  2-D print plot of source term.

    !     !/T1  Print arrays.

    !

    ! 10. Source code :

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_T

    USE w3odatmd, ONLY: ndst

#endif

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

#ifdef W3_T0

    USE w3arrymd, ONLY: prt2ds

#endif

#ifdef W3_T1

    USE w3arrymd, ONLY: outmat

#endif

    !/

    USE constants, ONLY: tpi

    USE w3servmd,  ONLY: extcde

    USE w3odatmd,  ONLY: ndse, iaproc, naproc, naperr

    USE w3gdatmd,  ONLY: nk, nth, nspec, sig, mapwn, ic5pars

    USE w3idatmd,  ONLY: inflags2, icep1, icep2, icep3, icep4, icei

    !

    IMPLICIT NONE

    !

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    REAL, INTENT(IN)        :: cg(nk), wn(nk), a(nspec), depth

    REAL, INTENT(OUT)       :: s(nspec), d(nspec)

    INTEGER, INTENT(IN)     :: ix, iy

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

#ifdef W3_T0

    INTEGER                :: ith

    REAL                   :: dout(nk,nth)

#endif

    !/

    REAL                    :: icecoef1, icecoef2, icecoef3, &

         icecoef4, iceconc

    REAL, DIMENSION(NK)     :: d1d, wn_r, wn_i

    !     REAL                    :: TWN_R, TWN_I

    INTEGER                 :: ik, ikth

    LOGICAL                 :: noice

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SIC5')

#endif

    !

    ! 0.  Initializations ------------------------------------------------ /

    d        = 0.

    d1d      = 0.

    wn_r     = 0.

    wn_i     = 0.

    !

    icecoef1 = 0.

    icecoef2 = 0.

    icecoef3 = 0.

    icecoef4 = 0.

    iceconc  = 0.

    !

    ! Set the ice parameters from input

    IF (inflags2(-7)) THEN

      icecoef1 = icep1(ix, iy) ! ice thickness h_i

    ELSE

      IF ( iaproc .EQ. naperr )                       &

           WRITE (ndse,1001) 'ICE PARAMETER 1 (HICE)'

      CALL extcde(2)

    ENDIF

    !

    IF (inflags2(-6)) THEN

      icecoef2 = icep2(ix, iy) ! effective viscosity of ice η

    ELSE

      IF ( iaproc .EQ. naperr )                       &

           WRITE (ndse,1001) 'ICE PARAMETER 2 (VISC)'

      CALL extcde(2)

    ENDIF

    !

    IF (inflags2(-5)) THEN

      icecoef3 = icep3(ix, iy) ! density of ice ρ_i

    ELSE

      IF ( iaproc .EQ. naperr )                       &

           WRITE (ndse,1001) 'ICE PARAMETER 3 (DENS)'

      CALL extcde(2)

    ENDIF

    !

    IF (inflags2(-4)) THEN

      icecoef4 = icep4(ix, iy) ! effective shear modulus of ice G

    ELSE

      IF ( iaproc .EQ. naperr )                       &

           WRITE (ndse,1001) 'ICE PARAMETER 4 (ELAS)'

      CALL extcde(2)

    ENDIF

    !

    IF (inflags2(4)) iceconc = icei(ix, iy) ! ice concentration

    !

    ! 1. No ice --------------------------------------------------------- /

    noice = .false.

    !     Zero ice thickness

    !     Very small ice thickness may cause problems in POLYROOTS because

    !     the first coefficient C1 may be very close to zero. So we regard

    !     cases where hice is less than 0.0001 as no ice.

    !     IF (ICECOEF1 < ERRTOL) NOICE = .TRUE.

    IF (icecoef1 < 0.0001) noice = .true.

    !     zero ice concentration

    IF (inflags2(4) .AND. iceconc < errtol) noice = .true.

    !

    ! Calculate the decay rate k_i

    IF ( noice ) THEN

      d1d = 0.

      !

      ! 2. Ice ------------------------------------------------------------- /

    ELSE

      !         W3IC5WNCG(WN_R, WN_I, CG, HICE, IVISC, RHOI, ISMODG, HWAT)

      CALL w3ic5wncg(wn_r, wn_i, cg, icecoef1, icecoef2, &

           icecoef3, icecoef4, depth)

      ! recall that D=S/E=-2*Cg_{ice}*k_i

      ! In some cases, the FS model yields very large Cg_{ice}, which

      ! subquently may result in numerical failure due to the violation of CFL

      ! conditions, therefore we still use ice-free group velocity to advect

      ! wave packets.

      !

      DO ik = 1, nk

        d1d(ik) = -2.0 * cg(ik) * wn_i(ik)

      END DO

    END IF

    !

    ! 2.1 Fill diagonal matrix

    DO ikth = 1, nspec

      d(ikth) = d1d(mapwn(ikth))

    END DO


    s = d * a

    !

    ! ... Test output of arrays

    !

#ifdef W3_T0

    DO ik=1, nk

      DO ith=1, nth

        dout(ik,ith) = d(ith+(ik-1)*nth)

      END DO

    END DO

    CALL prt2ds (ndst, nk, nk, nth, dout, sig(1:), '  ', 1.,    &

         0.0, 0.001, 'Diag Sice', ' ', 'NONAME')

#endif

    !

#ifdef W3_T1

    CALL outmat (ndst, d, nth, nth, nk, 'diag Sice')

#endif

    !

    ! Formats

    !

1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/   &

         '     ',a,' IS NOT DEFINED IN ww3_shel.inp.')


    !/

    !/ End of W3SIC5------------------------------------------------------ /

    !/

  END SUBROUTINE w3sic5

  !/ ------------------------------------------------------------------- /

  !/

  SUBROUTINE w3ic5wncg(WN_R, WN_I, CG, HICE, IVISC, RHOI, ISMODG, &

       HWAT)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |           E. Rogers               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         25-Apr-2017 |

    !/                  +-----------------------------------+

    !/

    !/    17-Apr-2016 : Origination                         ( version 5.10)

    !/                                                      ( Q. Liu )

    !/    17-Apr-2016 : Start from W3IC3WNCG_CHENG          ( Q. Liu )

    !/

    !/ 1. Purpose:

    !     Calculation of complex wavenumber arrays for ice-coupled waves.

    !

    !     This also allows us to use Cg_ice in the advection part of the

    !     radiative transfer energy equation (RTE). --- abandoned in the end

    !

    !  2. Method:

    !     Using the  Fox-Squire dispersion relations to get (kr, ki) and

    !     then get cg by cg = dσ / dk (here dk uses kr)

    !

    !  3. Parameters:

    !

    !     Parameter list:

    !     ----------------------------------------------------------------

    !     Name     Type   Intent   Description

    !     ----------------------------------------------------------------

    !     WN_R     R.A.   I/O      the real. part of the wave number

    !     WN_I     R.A.   I/O      the imag. part of the wave number

    !     CG       R.A.   I        group velocity (m s^{-1})

    !     HICE     Real.  I        thickness of ice (m)

    !     IVISC    Real.  I        viscosity parameter of ice (m^2 s^{-1})

    !     RHOI     Real.  I        the density of ice (kg m^{-3})

    !     ISMODG   Real.  I        effecitive shear modulus G of ice (Pa)

    !     HWAT     Real.  I        water depth

    !     ----------------------------------------------------------------

    !     * the intent of WN_R/I must be inout

    !     * CG is unchanged but still kept here because some legacy reasons.

    !

    !  4. Subroutines used:

    !

    !     Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !     FSDISP    Subr. /        dispersion relations for ice-coupled waves

    !     CGINICE5  Subr. /        group velocity for given (σ, kr) array

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !     Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !     W3SIC5    Subr. Public   Ice source term

    !     W3WAVE    Subr. W3WAVEMD WW3 integration

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    USE constants, ONLY: tpi

    USE w3gdatmd,  ONLY: nk, sig

    USE w3odatmd,  ONLY: ndse, iaproc, naproc, naperr

    USE w3servmd,  ONLY: extcde

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    REAL, INTENT(INOUT)   :: wn_r(:), wn_i(:)

    REAL, INTENT(IN)      :: cg(:)

    REAL, INTENT(IN)      :: hice, ivisc, rhoi, ismodg, hwat

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    REAL, ALLOCATABLE     :: sigma(:)

    INTEGER               :: kl, ku, ik

    REAL                  :: twn_r, twn_i

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3IC5WNCG')

#endif

    !/

    ! Initialize SIGMA {in w3gdatmd: SIG (0: NK+1)}

    IF (ALLOCATED(sigma)) DEALLOCATE(sigma); ALLOCATE(sigma(SIZE(cg)))

    sigma = 0.


    IF (SIZE(wn_r, 1) .EQ. nk) THEN

      kl    = 1

      ku    = nk

      sigma = sig(1:nk)

    ELSE IF (SIZE(wn_r,1) .EQ. nk+2) THEN

      kl    = 1

      ku    = nk+2

      sigma = sig(0:nk+1)

    ELSE

      IF ( iaproc .EQ. naperr ) WRITE(ndse,900) 'W3IC5WNCG'

      CALL extcde(3)

    END IF

    !

    ! Fox-Squire dispersion

    DO ik = kl, ku

      !         FSDISP(HICE, IVISC, RHOI, ISMODG, HWAT, WT, WNR, WNI)

      CALL fsdisp(hice, ivisc, rhoi, ismodg, hwat, tpi/sigma(ik), &

           twn_r, twn_i)

      wn_r(ik) = twn_r

      wn_i(ik) = twn_i

    END DO

    !

    DEALLOCATE(sigma)

    !

900 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/   &

         '     Subr. ', a, ': Cannot determine bounds of&

         & wavenumber array.'/)

    !/

    !/ End of W3IC5WNCG -------------------------------------------------- /

    !/

  END SUBROUTINE w3ic5wncg

  !/ ------------------------------------------------------------------- /

  !/

  SUBROUTINE fsdisp(HICE, IVISC, RHOI, ISMODG, HWAT, WT, WNR, WNI)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         19-May-2021 |

    !/                  +-----------------------------------+

    !/

    !/    17-Mar-2016 : Origination                         ( version 5.10)

    !/                                                      ( Q. Liu)

    !/    17-Mar-2016 : Start from the Matlab code `FoxSquire.m` (provided

    !/                  by Prof. Vernon Squire from University of Otago)

    !/                                                      ( Q. Liu)

    !/    25-Apr-2017 : Add more filters                    ( Q. Liu)

    !/

    !/    19-May-2021 : Incl. RP and M2 ice models          ( Q. Liu)

    !/

    !  1. Purpose :

    !

    !     Calculate the complex wavenumber for waves in ice according to

    !     three different sea ice models, i.e., FS, RP and M2 (see Liu et

    !     al. 2020)

    !

    !  2. Method :

    !     Mainly solving the dispersion relations of FS and RP models (

    !     Eqs. (20, 24, 25)) in Mosig et al. (2015))

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     Name     Type   Intent   Description

    !     ----------------------------------------------------------------

    !     HICE     Real.  IN       thickness of ice (m)

    !     IVISC    Real.  IN       viscosity parameter of ice (m^2 s^{-1})

    !     RHOI     Real.  IN       the density of ice (kg m^{-3})

    !     ISMODG   Real.  IN       effecitive shear modulus G of ice (Pa)

    !     HWAT     Real.  IN       water depth

    !     WT       Real.  IN       wave period (s; 1/freq)

    !     WNR      Real.  Out      the real. part of the wave number

    !     WNI      Real.  Out      the imag. part of the wave number

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !      POLYROOTS Subr. /        Find the roots of a general polynomial

    !      NR_ROOT   Func. /        Newton-Raphson root finding

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3IC5WNCG Subr. /        Wavenumber and group velocity of ice-

    !                               coupled waves

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     See Format 1000, 1001, 1002

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    USE constants, ONLY: grav, tpi

    USE w3dispmd,  ONLY: wavnu1

    USE w3servmd,  ONLY: extcde

    USE w3odatmd,  ONLY: ndse, iaproc, naproc, naperr

    USE w3gdatmd,  ONLY: ic5pars

    USE w3gsrumd,  ONLY: w3inan

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    REAL, INTENT(IN)      :: hice, ivisc, rhoi, ismodg, hwat, wt

    REAL, INTENT(OUT)     :: wnr, wni

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !

    REAL                  :: ic5minig, ic5minwt, ic5maxkratio, &

         ic5maxki, ic5minhw, ic5vemod

    REAL                  :: tismodg, twt, tratio, thw

    REAL, PARAMETER       :: nu = 0.3, rhow = 1025.

    !     COMPLEX               :: GV, C1

    !     REAL                  :: SIGMA, C2, WNO, CGO, THKH,  &

    COMPLEX               :: gv, c1, c2

    REAL                  :: sigma, wno, cgo, thkh,  &

         rtrl(5), rtim(5), rtang(5)

    INTEGER               :: ireal

    !     COMPLEX(KDPC)         :: GUESS, CROOT, C1D

    !     REAL(KDP)             :: C2D, HWATD

    COMPLEX(KDPC)         :: guess, croot, c1d, c2d

    REAL(kdp)             :: hwatd

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'FSDISP')

#endif

    ! Note, same as W3IC3WNCG_xx in w3sic3md :

    !     HICE   →   ICE1

    !     IVISC  →   ICE2

    !     RHOI   →   ICE3

    !     ISMODG →   ICE4

    ! 0.  Initializations ------------------------------------------------ *

    ! Set limiters

    !

    ! When G = 0, the FS method does not provide a solution. It is not

    ! unexpected because the FS model is originally devised as a

    ! thin elastic plate model in which elasticity is necessary.

    !

    ! The FS algorithm may also have issues for very short wave periods,

    ! shallow waters and low G (e.g., T~3 s, d~10 m, hi~0.5 m, G<10^6 Pa)

    !

    ic5minig     = ic5pars(1) ! Minimum G

    ic5minwt     = ic5pars(2) ! Minimum T

    ic5maxkratio = ic5pars(3) ! Maximum k_{ow}/k_r

    ic5maxki     = ic5pars(4) ! Maximum k_i

    ic5minhw     = ic5pars(5) ! Minimum d

    ic5vemod     = ic5pars(9) ! Model selected 1: EFS, 2: RP, 3: M2

    !

    tismodg  = max(ic5minig, ismodg)

    twt      = max(ic5minwt, wt)

    thw      = max(ic5minhw, hwat)

    !

    ! G <= 0. is not allowed

    IF (abs(tismodg) < errtol) THEN

      IF ( iaproc .EQ. naperr ) WRITE(ndse, 1000) 'FSDISP'

      CALL extcde(1)

    END IF

    !

    ! σ = 2π / T

    sigma = tpi / twt

    !

    IF (abs(ic5vemod - 1.) < errtol) THEN

      ! Complex shear modulus Gv = G - i σ ρ η (EFS model)

      gv = cmplx(tismodg, -1. * sigma * rhoi * ivisc)

      !

      ! -------------------------------------------------------------------- *

      ! Note that Eq. (24) in Mosig et al. (2015) can be written like below:

      ! (c1 * k^5 + c2 * k) * tanh(HWAT*k) - 1 = 0

      ! Most Important part of this module --------------------------------- *

      c1 = gv * hice**3. / (6. * rhow * sigma**2.)

      !

      ! To be divided by (1-NU) or multiplied by (1+NU) ??

      ! Beam model: then multiplied by (1+ν)

      ! Plate model: then divided by (1-ν)

      ! The beam version is more theoretically (J.E.M. Mosig, personal

      ! communication, 2016), although there is only very marginal difference

      ! between this two version as (1+NU = 1.3 and 1/(1-NU) ~ 1.4)

      c1 = c1 * (1+nu)

      !         C1 = C1 / (1-NU)

      !

      ! C2

      !         C2 = GRAV / SIGMA**2. - RHOI * HICE / RHOW

      c2 = cmplx(grav / sigma**2. - rhoi * hice / rhow, 0.)

      !

    ELSE IF (abs(ic5vemod - 2.) < errtol) THEN

      ! See Appendix of Liu et al. (2020) - RP model

      c1 = cmplx(tismodg * hice**3. * (1+nu) / (6. * rhow * sigma**2.), 0.)

      c2 = cmplx(grav/sigma**2. - rhoi * hice / rhow,              &

           -1. * ivisc / (rhow * sigma))

      !

    ELSE IF (abs(ic5vemod - 3.) > errtol) THEN

      WRITE(ndse, 1003) 'FSDISP', ic5vemod

      CALL extcde(4)

    END IF

    ! Use the dispersion in open water to get an approximation of

    ! tanh(HWAT * k). We can also roughly use the dispersion in deep

    ! water case, that is tanh(HWAT*k) ~ 1.

    ! Wavenumber in the open water

    !     WAVNU1(SI, H, K, CG)

    CALL wavnu1(sigma, thw, wno, cgo)

    thkh = tanh(wno * thw)

    !

    IF (abs(ic5vemod - 1.) < errtol .OR. abs(ic5vemod - 2.) < errtol) THEN

      ! Get the first guess of the complex wavenumber

      CALL polyroots(6, &

           (/real(real(c1))*thkh, 0., 0., 0., real(real(c2))*thkh, -1./),&

           rtrl, rtim)

      rtang = atan2(rtim, rtrl)

      !

      ! There should only be one real root in RTRL + i * RTIM because in

      ! this case (ivisc=0) the original viscoelastic-type model reduced to

      ! the thin elastic plate model which has only one real solution.

      ! Find its index ...

      !

      ireal = minloc(abs(rtang), dim=1)

      IF (rtrl(ireal) <= 0. .OR. abs(rtim(ireal)) > errtol) THEN

        IF ( iaproc .EQ. naperr ) WRITE(ndse, 1001) 'FSDISP'

        CALL extcde(2)

      END IF

      !

      ! Get the first guess for iteration

      guess = rtrl(ireal) * exp(cmplx(0., 1e-6))

      !

      ! Newton-Raphson method

      ! Turn c1, c2, hwat to be double

      c1d = c1; c2d = c2; hwatd = thw

      croot = nr_root(c1d, c2d, hwatd, guess)

      wnr = real(real(croot))

      wni = real(aimag(croot))

      !

    ELSE IF (abs(ic5vemod - 3.) < errtol) THEN ! M2

      ! Model with Order 3 Power Law (section 6.2 in Meylan et al. (2018, JGR-Ocean))

      ! Based on my understanding, the wavelength does not change because

      ! the elasticity is not considered in this model

      wnr = wno ! Open-water wavenumber

      ! Eq. (53) in Meylan et al. (2018)

      wni = hice * ivisc * sigma**3. / (rhow * grav**2.)

    END IF

    !

    ! RATIO Check

    ! Using the ratio k0 / kr as a basic check for the reliability of

    ! FSDISP. The FS dispersion relation can give a very different kr from

    ! k0, especially for small wave periods (k0/kr is as high as 100).

    ! From my tests, using IC5MAXKRATIO = 1000. can basically detect most

    ! spurious solutions (although not all of them)

    !

    ! ISNAN Check

    ! Common ways used are:

    ! NAN = SQRT(-1.) or

    ! a /= a then a is NaN or

    ! ISNAN func (supported by gfortran & ifort)

    !        --- ISNAN -> W3INAN because ISNAN is not supported by pgi

    ! For very few cases, we can get nan | negative ki | kr

    !

    ! (N.B.) NaN problem solved by using CMPLX_TANH2

    !

    tratio = wno / wnr

    IF (w3inan(wnr) .OR. w3inan(wni) .OR. wnr <= 0 .OR. wni <= 0. &

         .OR. tratio >= ic5maxkratio) THEN

      IF ( iaproc .EQ. naperr )                       &

           WRITE(ndse, 1002) 'FSDISP', hice, ivisc, tismodg, hwat, twt, &

           wno, wnr, wni

      CALL extcde(3)

    END IF

    !

    ! Filter high ki

    wni = min(ic5maxki, wni)

    !

    ! FORMAT

1000 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/   &

         '     Subr. ', a, ': Zero shear modulus G is not allowed&

         & in the FS viscoelastic model'/)

    !

1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/   &

         '     Subr. ', a, ': get a bad first guess'/)

    !

1002 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/   &

         ' -----------------------------------------------------'/&

         '     Subr. ', a,'   : get NaN/NeG/Huge kr or ki for'   /&

         ' -----------------------------------------------------'/&

         '  Ice thickness     : ', f9.1, ' m'/         &

         '  Ice viscosity     : ', e9.2, ' m2/s'/      &

         '  Ice shear modulus : ', e9.2, ' Pa' /       &

         '  Water depth       : ', f9.1, ' m'/         &

         '  Wave period       : ', f10.2, ' s'/        &

         '  Wave number (Ko)  : ', f11.3, ' rad/m'/    &

         '  Wave number (Kr)  : ', f11.3, ' rad/m'/    &

         '  Attenu. Rate (Ki) : ', e9.2, ' /m'/)

    !

1003 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/   &

         '     Subr. ', a, ': Unknown VE model (', f5.0, ')'/)

    !/

    !/ End of FSDISP ----------------------------------------------------- /

    !/

  END SUBROUTINE fsdisp

  !/ ------------------------------------------------------------------- /

  !/

  SUBROUTINE balancing_matrix(NMAT, MATRIX)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         15-Mar-2016 |

    !/                  +-----------------------------------+

    !/

    !/    15-Mar-2016 : Origination                         ( version 5.10)

    !/                                                      ( Q. Liu )

    !/    15-Mar-2016 : Borrowed from Numerical Recipes in Fortran

    !/                                                      ( Q. Liu )

    !  1. Purpose :

    !     Reducing the sensitivity of eigenvalues to rounding errors during

    !     the execution of some algorithms.

    !

    !  2. Method :

    !     The errors in the eigensystem found by a numerical procedure are

    !     generally proportional to the Euclidean norm of the matrix, that

    !     is, to the square root of the sum of the squares of the elements

    !     (sqrt(sum(a_{i, j} ** 2.)). The idea of balancing is to use

    !     similarity transformations to make corresponding rows and columns

    !     of the matrix have comparable norms, thus reducing the overall

    !     norm of the matrix while leaving the eigenvalues unchanged. Note

    !     that the symmetric matrix is already balanced.

    !

    !     The output is matrix that is balanced in the norm given by

    !     summing the absolute magnitudes of the matrix elements(

    !     sum(abs(a_{i, j})) ). This is more efficient than using the

    !     Euclidean norm, and equally effective: a large reduction in

    !     one norm implies a large reduction in the other.

    !

    !     For the details of this method, please refer to

    !     1) Numerical Recipes in Fortran 77 (Volume 1, 2nd Edition)

    !        [Chapter 11.5 / subroutine balanc]

    !     2) Numerical Recipes in Fortran 90 (Volume 2)

    !        [Chapter B11 / subroutine balanc]

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     Name     Type   Intent   Description

    !     ----------------------------------------------------------------

    !     NMAT     Int.   I        The size of one dimension of MATRIX

    !     MATRIX   R.A.   I/O      A matrix with the shape (NMAT, NMAT)

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing (!/S switch).

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      POLYROOTS Subr. /        Find the roots of polynomials

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     None.

    !

    !  7. Remarks:

    !     Balancing only needs marginal computational efforts but can

    !     substantially improve the accuracy of the eigenvalues computed

    !     for a badly balanced matrix. It is therefore recommended that

    !     you always balance nonsymmetric matrices.

    !

    !     Given a (NMAT, NMAT) MATRIX, this routine replaces it by a

    !     balanced matrix with identical eigenvalues. A symmetric matrix is

    !     already balanced and is unaffected by this procedure.

    !

    !  8. Structure :

    !

    !     See the source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    INTEGER, INTENT(IN)    :: nmat

    REAL, INTENT(INOUT)    :: matrix(nmat, nmat)

    !/ ------------------------------------------------------------------- /

    !/ Local parameter

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    ! the parameter radx is the machine's floating-point radix

    REAL, PARAMETER        :: radx = radix(matrix), &

         sqradx = radx ** 2

    INTEGER                :: i, last

    REAL                   :: c, f, g, r, s

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    CALL strace (ient, 'BALANCING_MATRIX')

#endif

    !

    DO

      last = 1

      DO i = 1, nmat

        ! Calculate row and column norms

        c = sum( abs(matrix(:, i)) ) - matrix(i, i)

        r = sum( abs(matrix(i, :)) ) - matrix(i, i)

        ! If both are non-zero

        IF (c /= 0.0 .AND. r /= 0.0) THEN

          ! Find the integer power of the machine radix that comes closest to

          ! balancing the matrix (get G, F from C, R)

          g = r / radx

          f = 1.0

          s = c + r

          DO

            IF (c >= g) EXIT

            f = f * radx

            c = c * sqradx

          END DO

          !

          g = r * radx

          DO

            IF (c <= g) EXIT

            f = f / radx

            c = c / sqradx

          END DO

          !

          IF ( (c+r)/f < 0.95*s) THEN

            last = 0

            g = 1.0 / f

            ! Apply similarity tranformation

            matrix(i, :) = matrix(i, :) * g

            matrix(:, i) = matrix(:, i) * f

          END IF

        END IF

      END DO

      IF (last /= 0) EXIT

    END DO

    !/

    !/ End of subroutine BALANCING_MATRIX -------------------------------- /

    !/

  END SUBROUTINE balancing_matrix

  !/ ------------------------------------------------------------------- /

  !/

  SUBROUTINE eig_hqr (NMAT, HMAT, EIGR, EIGI)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         17-Mar-2016 |

    !/                  +-----------------------------------+

    !/

    !/    16-Mar-2016 : Origination                         ( version 5.10)

    !/                                                      ( Q. Liu )

    !/    16-Mar-2016 : Borrowed from Numerical Recipes in Fortran

    !/                                                      ( Q. Liu )

    !/    17-Mar-2016 : Update the NR code v2.08 to v2.10   ( Q. Liu )

    !/

    !  1. Purpose :

    !

    !     When we calculate the eigenvalues of a general matrix, we first

    !     reduce the matrix to a simpler form (e.g., Hessenberg form) and

    !     then we perform the iterative procedures.

    !

    !     A upper Hessenberg matrix has zeros everywhere below the diagnal

    !     except for the first subdiagonal row. For example, in the 6x6

    !     case, the non-zero elements are:

    !                   |x x x x x x|

    !                   |x x x x x x|

    !                   |  x x x x x|

    !                   |    x x x x|

    !                   |      x x x|

    !                   |        x x|

    !

    !     This subroutine uses QR algorithm to get the eigenvalues of a

    !     Hessenberg matrix. So make sure the input array HMAT is a

    !     Hessenberg-type matrix.

    !

    !  2. Method :

    !     QR algorithm for real Hessenberg matrices.

    !     (I did not understand this algorithm well, so I could not give

    !      any detailed explanations)

    !

    !     For the details of this HQR method, please refer to

    !     1) Numerical Recipes in Fortran 77 (Volume 1, 2nd Edition)

    !        [Chapter 11.6 / subroutine hqr]

    !     2) Numerical Recipes in Fortran 90 (Volume 2)

    !        [Chapter B11 / subroutine hqr]

    !

    !     Note that there is a bug in the `hqr` subroutine in NR v2.08.

    !     See http://numerical.recipes/latest-known-bugs.html. Please use

    !     the updated code in NR v2.10.

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     Name     Type   Intent   Description

    !     ----------------------------------------------------------------

    !     NMAT     Int.   I        the size of one dimension of HMAT

    !     HMAT     R.A.   I/O      the Hessenberg-type matrix (NMAT, NMAT)

    !     EIGR     R.A.   O        the real part of the N eigenvalues

    !     EIGI     R.A.   O        the imag part of the N eigenvalues

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      POLYROOTS Subr. /        Find the roots of polynomials

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     See Format 1001

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    USE w3servmd, ONLY: extcde

    USE w3odatmd, ONLY: ndse, iaproc, naproc, naperr

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    INTEGER, INTENT(IN)  :: nmat

    REAL, INTENT(INOUT)  :: hmat(nmat, nmat)

    REAL, INTENT(OUT)    :: eigr(nmat), eigi(nmat)

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    !/

    INTEGER              :: i, its, k, l, m, nn, mnnk, idiag

    REAL                 :: anorm, p, q, r, s, t, u, v, w, x, y, z

    REAL                 :: pp(nmat)

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'EIG_HQR')

#endif

    !

    ! Compute matrix norm for possible use in locating single small

    ! subdiagonal element.

    !

    ! Note the speciality of Hessenberg matrix :

    ! Elements below the diagonal are zeros except for the first

    ! subdiagonal row. It might be more accurate if we use a mask array

    ! to mask all zero elments.

    !

    anorm = sum(abs(hmat))

    nn = nmat

    ! Gets changed only by an exceptional shift.

    t = 0.0

    ! Begin search for next eigenvalue: "do while nn >= 1"

    DO

      IF (nn < 1) EXIT

      its=0

      ! Begin iteration

      iterate:DO

        ! Look for single small subdiagonal element.

        small: DO l=nn, 2, -1

          s = abs( hmat(l-1, l-1) ) + abs( hmat(l, l) )

          !                  IF (S == 0.0) S = ANORM

          IF (abs(s) < errtol) s = anorm

          !                  IF ( ABS(HMAT(L, L-1)) + S == S ) THEN

          IF ( abs(hmat(l, l-1)) < errtol ) THEN

            hmat(l, l-1) = 0.0

            EXIT small

          END IF

        END DO small

        x = hmat(nn, nn)

        ! One root found

        IF (l == nn) THEN

          eigr(nn) = x + t

          eigi(nn) = 0.0

          nn=nn-1

          ! Go back for next eigenvalue

          EXIT iterate

        END IF

        y = hmat(nn-1, nn-1)

        w = hmat(nn, nn-1) * hmat(nn-1, nn)

        ! Two roots found . . .

        IF (l == nn-1) THEN

          p = 0.5 * (y - x)

          q = p**2 + w

          z = sqrt( abs(q) )

          x = x + t

          ! . . . A real pair . . .

          IF (q >= 0.0) THEN

            z = p + sign(z, p)

            eigr(nn) = x + z

            eigr(nn-1) = eigr(nn)

            IF (z /= 0.0) eigr(nn) = x - w/z

            eigi(nn) = 0.0

            eigi(nn-1) = 0.0

            ! . . . A complex pair

          ELSE

            eigr(nn) = x + p

            eigr(nn-1) = eigr(nn)

            eigi(nn) = z

            eigi(nn-1) = -z

          END IF

          nn=nn-2

          ! GO BACK FOR NEXT EIGENVALUE.

          EXIT iterate

        END IF

        ! NO ROOTS FOUND. CONTINUE ITERATION.

        IF (its == 30) THEN

          IF ( iaproc .EQ. naperr ) WRITE(ndse, 1001) 'EIG_HQR'

          CALL extcde(2)

        END IF

        ! FORM EXCEPTIONAL SHIFT.

        IF (its == 10 .OR. its == 20) THEN

          t = t + x

          ! Add -X to the diagonal of HMAT

          DO idiag = 1, nn

            hmat(idiag, idiag) = hmat(idiag, idiag) + (-x)

          END DO

          s = abs(hmat(nn, nn-1)) + abs(hmat(nn-1, nn-2))

          x = 0.75 * s

          y = x

          w = -0.4375 * s**2

        END IF

        its = its + 1

        ! Form shift and then look for 2 consecutive small subdiagonal elements.

        DO m = nn-2, l, -1

          z = hmat(m, m)

          r = x - z

          s = y - z

          ! Equation (11.6.23).

          p = (r * s - w) / hmat(m+1, m) + hmat(m, m+1)

          q = hmat(m+1, m+1) - z - r - s

          r = hmat(m+2, m+1)

          ! Scale to prevent overflow or underflow

          s = abs(p) + abs(q) + abs(r)

          p = p / s

          q = q / s

          r = r / s

          IF (m == l) EXIT

          u = abs( hmat(m, m-1) ) * ( abs(q) + abs(r) )

          v = abs(p) * ( abs(hmat(m-1, m-1)) + abs(z) + &

               abs( hmat(m+1, m+1) ))

          ! Equation (11.6.26)

          IF (u+v == v) EXIT

        END DO

        DO i= m+2, nn

          hmat(i, i-2) = 0.0

          IF (i /= m+2) hmat(i, i-3)=0.0

        END DO

        ! Double QR step on rows L to NN and columns M to NN

        DO k=m, nn-1

          IF (k /= m) THEN

            ! Begin setup of householder vector

            p = hmat(k, k-1)

            q = hmat(k+1, k-1)

            r = 0.0

            IF (k /= nn-1) r = hmat(k+2, k-1)

            x = abs(p) + abs(q) + abs(r)

            IF (x /= 0.0) THEN

              ! Scale to prevent overflow or underflow

              p = p / x

              q = q / x

              r = r / x

            END IF

          END IF

          s = sign(sqrt(p**2 + q**2 + r**2), p)

          IF (s /= 0.0) THEN

            IF (k == m) THEN

              IF (l /= m) hmat(k, k-1) = -hmat(k, k-1)

            ELSE

              hmat(k, k-1) = -s * x

            END IF

            ! Equations (11.6.24).

            p = p + s

            x = p / s

            y = q / s

            z = r / s

            q = q / p

            ! READY FOR ROW MODIFICATION.

            r = r / p

            pp(k:nn) = hmat(k, k:nn) + q * hmat(k+1, k:nn)

            IF (k /= nn-1) THEN

              pp(k:nn) = pp(k:nn) + r * hmat(k+2, k:nn)

              hmat(k+2, k:nn) = hmat(k+2, k:nn) - &

                   pp(k:nn)*z

            END IF

            hmat(k+1, k:nn) = hmat(k+1, k:nn) - pp(k:nn) * y

            hmat(k, k:nn) = hmat(k, k:nn) - pp(k:nn) * x

            ! COLUMN MODIFICATION.

            mnnk = min(nn, k+3)

            pp(l:mnnk) = x * hmat(l:mnnk, k) + y * &

                 hmat(l:mnnk, k+1)

            IF (k /= nn-1) THEN

              pp(l:mnnk) = pp(l:mnnk) + z*hmat(l:mnnk, k+2)

              hmat(l:mnnk, k+2) = hmat(l:mnnk, k+2) - &

                   pp(l:mnnk) * r

            END IF

            hmat(l:mnnk, k+1) = hmat(l:mnnk, k+1) - &

                 pp(l:mnnk) * q

            hmat(l:mnnk, k) = hmat(l:mnnk, k) - pp(l:mnnk)

          END IF

        END DO

        ! GO BACK FOR NEXT ITERATION ON CURRENT EIGENEND DO VALUE.

      END DO iterate

    END DO

    !

    ! Formats

1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/  &

         '     Subr. ', a, ': TOO MANY ITERATIONS'/)

    !/ ------------------------------------------------------------------- /

    !/

    !/ End of EIG_HQR ---------------------------------------------------- /

    !/

  END SUBROUTINE eig_hqr

  !/ ------------------------------------------------------------------- /

  !/

  SUBROUTINE polyroots(NPC, PCVEC, RTRL, RTIM)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         16-Mar-2016 |

    !/                  +-----------------------------------+

    !/

    !/    16-Mar-2016 : Origination                         ( version 5.10)

    !/                                                      ( Q. Liu )

    !/    16-Mar-2016 : Started from Numerical Recipes in Fortran

    !/                                                      ( Q. Liu )

    !/

    !  1. Purpose :

    !

    !     Find the roots of arbitrary polynomials through finding the

    !     eigenvalues of companion matrix.

    !

    !  2. Method :

    !     Suppose we have a general polynomial, which reads

    !     P(x) = c_n * x^n + c_{n-1} * x^{n-1} + ... + c_1 * x + c_0

    !

    !     Then finding the roots of P(x) is equivalent to find the eigen-

    !     values of the special n x n companion matrix A

    !         |  -c_{n-1}/c_n   -c_{n-2}/c_n   ...   -c_1/c_n   -c_0/c_n |

    !         |  1              0              ...   0          0        |

    !     A = |  0              1              ...   0          0        |

    !         |  :              :                    :          :        |

    !         |  0              0                    1          0        |

    !

    !     In fact, P(x) is the characteristic polynomial of matrix A, i.e.,

    !     P(x) = det|A-xI| and x is the eigenvalues of A (this is a

    !     Hessenberg matrix).

    !

    !     In this subrountine, we will use the two subroutines above

    !     (BALANCING_MATRIX & EIG_HQR) to get the complex eigenvalues of

    !     an abitrary Hessenberg matrix

    !

    !     For the details of this method, please refer to

    !     1) Numerical Recipes in Fortran 77 (Volume 1, 2nd Edition)

    !        [Chapter 9.5 / Eigenvalue Methods / subroutine zrhqr]

    !     2) Numerical Recipes in Fortran 90 (Volume 2)

    !        [Chapter B9 / subroutine zrhqr]

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     Name     Type   Intent   Description

    !     ----------------------------------------------------------------

    !     NPC      Int.   I        The # of the Polynomial coefficients

    !                              (from c_n to c_0)

    !     PCVEC    R.A.   I        The 1d vector for the Polynomial

    !                              coefficients  [c_n, c_{n-1}, ..., c_0]

    !     RTRL     R.A.   O        The real part of all of the roots

    !                              shape: [NPC-1]

    !     RTIM     R.A.   O        The real part of all of the roots

    !                              shape: [NPC-1]

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !      Name                Type    Module    Description

    !     ------------------------------------- ---------------------------

    !      STRACE              Subr.   W3SERVMD  Subroutine tracing.

    !      BALANCING_MATRIX    Subr.   /         Balancing matrix

    !      EIG_HQR             Subr.   /         Finding eigenvalues

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      FSDISP    Subr. /        Solving the dispersion relations

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     See Format 1001

    !

    !  7. Remarks :

    !     The built-in MATLAB function <roots> uses the same method to

    !     find roots of a general polynomial. But perhaps MATLAB uses

    !     different methods to find eigenvalues of the companion matrix.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    USE w3servmd, ONLY: extcde

    USE w3odatmd, ONLY: ndse, iaproc, naproc, naperr

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    INTEGER, INTENT(IN)   :: npc

    REAL, INTENT(IN)      :: pcvec(npc)

    REAL, INTENT(OUT)     :: rtrl(npc-1), rtim(npc-1)

    !/

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    REAL                  :: hess(npc-1, npc-1)

    INTEGER               :: j

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'POLYROOTS')

#endif

    !

    !

    IF (abs(pcvec(1)) < errtol) THEN

      IF ( iaproc .EQ. naperr ) WRITE(ndse, 1001) 'POLYROOTS'

      CALL extcde(2)

    END IF

    !

    ! Generate the Hessenberg matrix

    hess = 0.

    hess(1, :) = -1 * pcvec(2:) / pcvec(1)

    DO j = 1, npc-2

      hess(j+1, j) = 1.

    END DO


    ! Balancing the matrix HESS

    CALL balancing_matrix(npc-1, hess)

    ! Eigenvalues of the matrix HESS

    CALL eig_hqr(npc-1, hess, rtrl, rtim)


    ! Formats

1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/  &

         '     Subr. ', a, ': the coeff. of x^n must not be 0'/)

    !/

    !/ End of POLYROOTS -------------------------------------------------- /

    !/

  END SUBROUTINE polyroots

  !/ ------------------------------------------------------------------- /

  !/

  FUNCTION nr_corr(K, C1, C2, H)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         19-May-2021 |

    !/                  +-----------------------------------+

    !/

    !/    18-Mar-2016 : Origination.                        ( version 5.10 )

    !/                                                      ( Q. Liu )

    !/    18-Mar-2016 : Start from the Matlab code `FoxSquire.m` (provided

    !/                  by Prof. Vernon Squire from University of Otago)

    !/                                                      ( Q. Liu )

    !/    24-Mar-2016 : Adding the cmplx_sinh/cosh/tanh     ( Q. Liu )

    !/

    !/    19-May-2021 : Change types of few input arguments ( Q. Liu )

    !/

    !  1. Purpose :

    !

    !     Calculate the corrected term in the Newton-Raphson root-finding

    !     method (Must use double precision)

    !

    !  2. Method :

    !     Suppose we want to find the root of f(x) = 0, then according to

    !     the Newton-Raphson method, the root is iteratively updated by the

    !     formula below:

    !

    !         x_{i+1} = x_{i} - f(x_{i}) / f'(x_{i}),

    !

    !     where f'(x) denotes the derivative of f(x). In this function,

    !     our f(x) reads (see also subr. FSDISP)

    !

    !         f(x) = (c1 * k**4 + c2) * k * tanh(kH) -1

    !

    !     we finally will get the Newton-Raphson correted term, i.e.,

    !

    !         dx = f(x_{i}) / f'(x_{i})

    !

    !     For the details of this method, please refer to

    !     1) Numerical Recipes in Fortran 77 (Volume 1, 2nd Edition)

    !        Chapter 9.4

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     Name     Type      Intent   Description

    !     ----------------------------------------------------------------

    !     K        CMPL.(D)  I        complex wave number

    !     C1       CMPL.(D)  I        C1 in FSDISP

    !     C2       Real.(D)  I        C2 in FSDISP

    !     H        Real.(D)  I        water depth

    !     NR_CORR  CMPL.(D)  O        Newton-Raphson corrected term (DK)

    !     ----------------------------------------------------------------

    !     * (D) means double precision

    !

    !  4. Subroutines used :

    !

    !      Name        Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE      Subr. W3SERVMD Subroutine tracing.

    !      CMPLX_SINH  Func. /         sinh for complex var.

    !      CMPLX_COSH  Func. /         cosh for complex var.

    !      CMPLX_TANH2 Func. /         tanh for complex var.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      NR_ROOT   Func. /        Newton-Raphson root finding.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !     COMPLEX(KDPC), INTENT(IN)     :: K, C1

    !     REAL(KDP), INTENT(IN)         :: C2, H

    COMPLEX(KDPC), INTENT(IN)     :: k, c1, c2

    REAL(kdp), INTENT(IN)         :: h

    COMPLEX(KDPC)                 :: nr_corr

    !/

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    ! A rough value to differentiate deep water case from finite water case

    REAL(kdp), PARAMETER          :: kh_lim = 7.5

    COMPLEX(KDPC)                 :: lam, lampr, fv, df, tkh

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'NR_CORR')

#endif

    ! f(k) = (c1 * k**4 + c2) * k * tanh(k*H) - 1

    !      = lam * k * tanh(k*H) - 1

    !

    tkh   = k * h

    lam   = c1 * k**4 + c2

    ! the derivative of (lam * k)

    lampr = 5 * c1 * k**4 + c2

    !

    IF (real(real(tkh)) <= kh_lim) THEN

      !         KH is small enough

      !         FV = LAM * K * SINH(K*H) - COSH(K*H)

      !         DF = LAM * (K*H) * COSH(K*H) + (LAMPR - H) * SINH(K*H)

      fv = lam * k * cmplx_sinh(tkh) - cmplx_cosh(tkh)

      df = lam * tkh * cmplx_cosh(tkh) + (lampr-h) * cmplx_sinh(tkh)

    ELSE

      !         FV = LAM * K * TANH(K*H) - 1

      !         DF = LAM * K * H + (LAMPR - H) * TANH(K*H)

      !         DF = LAMPR * TANH(K*H) + LAM * K * H / (COSH(K*H) **2)

      fv = lam * k * cmplx_tanh2(tkh) - 1

      df = lampr * cmplx_tanh2(tkh) + lam * tkh * &

           (1 - cmplx_tanh2(tkh) ** 2.)

    END IF

    !

    nr_corr = fv / df

    !/

    !/ End of NR_CORR ---------------------------------------------------- /

    !/

  END FUNCTION nr_corr

  !/ ------------------------------------------------------------------- /

  !/

  FUNCTION nr_root(C1, C2, H, GUESS)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         19-May-2021 |

    !/                  +-----------------------------------+

    !/

    !/    18-Mar-2016 : Origination.                        ( version 5.10 )

    !/                                                      ( Q. Liu )

    !/    18-Mar-2016 : Start from the Matlab code `FoxSquire.m` (provided

    !/                  by Prof. Vernon Squire from University of Otago)

    !/                                                      ( Q. Liu )

    !/

    !/    19-May-2021 : Change types of few input arguments ( Q. Liu )

    !/

    !  1. Purpose :

    !

    !     The iterative procedure of the Newton-Raphson method

    !

    !  2. Method :

    !     See the document of Subr. NR_CORR (Must use double precision)

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     Name     Type      Intent  Description

    !     ----------------------------------------------------------------

    !     C1       CMPL.(D)  I       C1 in FS dipsersion relations

    !                                See the doc. of subr. NR_CORR

    !     C2       REAL (D)  I       C2 in FS dipsersion relations

    !     H        REAL (D)  I       water depth

    !     GUESS    CMPL.(D)  I       the first guess obtained from POLYROOTS

    !     NR_ROOT  CMPL.(D)  O       the calculated complex wave number.

    !     ----------------------------------------------------------------

    !     * (D) means double precision

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !      NR_CORR   Func. /        Newton-Raphson correction term

    !      INIT_RANDOM_SEED

    !                Subr. /        Initialize the random seed based on

    !                               the system's time

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      FSDISP    Subr. /        Solve  FS dispersion relations

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    USE w3servmd, ONLY: extcde

    USE w3odatmd, ONLY: ndse, iaproc, naproc, naperr

    USE w3gdatmd, ONLY: ic5pars

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    !     COMPLEX(KDPC), INTENT(IN)     :: C1, GUESS

    !     REAL(KDP), INTENT(IN)         :: C2, H

    COMPLEX(KDPC), INTENT(IN)     :: c1, guess, c2

    REAL(kdp), INTENT(IN)         :: h

    COMPLEX(KDPC)                 :: nr_root

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    COMPLEX(KDPC)                 :: k0, k1, dk

    INTEGER                       :: iter

    REAL                          :: tranval

    REAL                          :: ic5maxiter, ic5rkick, ic5kfilter

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'NR_ROOT')

#endif

    !/ Set parameters

    ic5maxiter = ic5pars(6)

    ic5rkick   = ic5pars(7) ! 0: False, 1: True

    ic5kfilter = ic5pars(8)

    !

    k0 = guess

    dk = nr_corr(k0, c1, c2, h)

    k1 = k0 - dk

    iter = 0


    IF (ic5rkick > 0.5) CALL init_random_seed()

    !

    DO WHILE (abs(dk) > errtol)

      k0 = k1

      dk = nr_corr(k0, c1, c2, h)

      k1 = k0 - dk

      iter = iter + 1

      !

      ! Random kick to avoid converging to evanescent modes

      ! Note: do not use RAND(1) because it alway gives a same random no.

      ! The built in function of RAND is not available in <ifort>, use

      ! random_seed/number instead.

      !

      ! Based on many tests, I found the random kick & the corridor excluded

      ! from imaginary axis are kind of helpful to avoid spurious solutions.

      ! However, it may also lead to no solutions returned, especially for

      ! high G and high T.

      !

      IF (ic5rkick > 0.5 .AND. abs(real(k1)) < ic5kfilter) THEN

        !             K1 = K1 + 2*RAND(0)

        CALL random_number(tranval)

        k1 = k1 + 2 * tranval

      END IF

      !

      IF (iter >= ic5maxiter) THEN

        IF ( iaproc .EQ. naperr ) WRITE(ndse, 1001) 'NR_ROOT'

        CALL extcde(1)

      END IF

      !

    END DO

    !

    nr_root = k1

    !

    ! Formats

1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/  &

         '     Subr. ', a, ': TOO MANY ITERATIONS'/)

    !

    !/

    !/ End of NR_ROOT ---------------------------------------------------- /

    !/

  END FUNCTION nr_root

  !/ ------------------------------------------------------------------- /

  !/

  FUNCTION cmplx_sinh(X)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         24-Mar-2016 |

    !/                  +-----------------------------------+

    !/

    !/    24-Mar-2016 : Origination.                        ( version 5.10 )

    !/                                                      ( Q. Liu )

    !/

    !  1. Purpose :

    !

    !     For a number of compilers, the built-in function sinh, cosh and

    !     tanh do not support the complex inputs. So here I write an

    !     external one.

    !

    !  2. Method :

    !

    !     sinh(x) = (e**x - e**(-x)) / 2 (The built in function exp supports

    !     complex input)

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     Name     Type      Intent  Description

    !     ----------------------------------------------------------------

    !      X       CMPL(D)     I     a double-precision complex var.

    !     ----------------------------------------------------------------

    !      * Note, this subr. will be only called by NR_CORR,

    !      so for simplicity, I only use double-precision complex var.

    !      as input.

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name       Type  Module   Description

    !     ----------------------------------------------------------------

    !      NR_CORR    Subr. /        Newton-Raphson correction term.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    COMPLEX(KDPC), INTENT(IN)  :: x

    COMPLEX(KDPC)              :: cmplx_sinh

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    CALL strace (ient, 'CMPLX_SINH')

#endif

    !/

    cmplx_sinh = (exp(x) - exp(-x)) * 0.5

    !/

    !/ End of CMPLX_SINH ------------------------------------------------- /

    !/

  END FUNCTION cmplx_sinh

  !/ ------------------------------------------------------------------- /

  !/

  FUNCTION cmplx_cosh(X)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         24-Mar-2016 |

    !/                  +-----------------------------------+

    !/

    !/    24-Mar-2016 : Origination.                        ( version 5.10 )

    !/                                                      ( Q. Liu )

    !/

    !  1. Purpose :

    !

    !     For a number of compilers, the built-in function sinh, cosh and

    !     tanh do not support the complex inputs. So here I write an

    !     external one.

    !

    !  2. Method :

    !

    !     cosh(x) = (e**x + e**(-x)) / 2 (The built in function exp supports

    !     complex input)

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     Name     Type      Intent  Description

    !     ----------------------------------------------------------------

    !      X       CMPL(D)     I     a double-precision complex var.

    !     ----------------------------------------------------------------

    !      * Note, this subr. will be only called by NR_CORR,

    !      so for simplicity, I only use double-precision complex var.

    !      as input.

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name       Type  Module   Description

    !     ----------------------------------------------------------------

    !      NR_CORR    Subr. /        Newton-Raphson correction term.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    COMPLEX(KDPC), INTENT(IN)  :: x

    COMPLEX(KDPC)              :: cmplx_cosh

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    CALL strace (ient, 'CMPLX_COSH')

#endif

    !/

    cmplx_cosh = (exp(x) + exp(-x)) * 0.5

    !/

    !/ End of CMPLX_COSH ------------------------------------------------- /

    !/

  END FUNCTION cmplx_cosh

  !/ ------------------------------------------------------------------- /

  !/

  FUNCTION cmplx_tanh2(X)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         24-Mar-2016 |

    !/                  +-----------------------------------+

    !/

    !/    24-Mar-2016 : Origination.                        ( version 5.10 )

    !/                                                      ( Q. Liu )

    !/

    !  1. Purpose :

    !     We may encounter overflow error for the above tanh function as kh

    !     becomes huge. This is another version of tanh function

    !

    !  2. Method :

    !

    !     See https://en.wikipedia.org/wiki/Hyperbolic_function

    !     tanh(x) = (exp(x) - exp(-x)) / (exp(x) + exp(-x))

    !             = (1 - exp(-2x)) / (1 + exp(-2x))

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !     Name     Type      Intent  Description

    !     ----------------------------------------------------------------

    !      X       CMPL(D)     I     a double-precision complex var.

    !     ----------------------------------------------------------------

    !      * Note, this subr. will be only called by NR_CORR, so for

    !        simplicity, I only use double-precision complex var. as input.

    !

    !  4. Subroutines used :


    !      Name        Type  Module   Description

    !     ----------------------------------------------------------------

    !      STRACE      Subr. W3SERVMD Subroutine tracing.

    !      CMPLX_SINH  Func. /        sinh for complex var.

    !      CMPLX_COSH  Func. /        cosh for complex var.

    !     ----------------------------------------------------------------

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      NR_CORR   Subr. /        Newton-Raphson correction term.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !       Calculating tanh in this way may have problems when x ->

    !       -inf. But in our cases x is alway >0.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    COMPLEX(KDPC), INTENT(IN)  :: x

    COMPLEX(KDPC)              :: cmplx_tanh2

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    CALL strace (ient, 'CMPLX_TANH2')

#endif

    !/

    cmplx_tanh2 = (1 - exp(-2*x)) / (1 + exp(-2*x))

    !/

    !/ End of CMPLX_TANH2 ------------------------------------------------ /

    !/

  END FUNCTION cmplx_tanh2

  !/

  !/ ------------------------------------------------------------------- /

  !/

  SUBROUTINE init_random_seed()

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         24-Mar-2016 |

    !/                  +-----------------------------------+

    !/

    !/    24-Mar-2016 : Origination.                        ( version 5.10 )

    !/                                                      ( Q. Liu )

    !/    24-Mar-2016 : Borrowed from Fortran Wiki          ( Q. Liu )

    !

    !  1. Purpose :

    !

    !     Initialize the random seed based on the system's time.

    !

    !  2. Method :

    !

    !     See http://fortranwiki.org/fortran/show/random_seed

    !

    !  3. Parameters :

    !

    !  4. Subroutines used :

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      NR_ROOT   Func. /        Newton-Raphson root finding.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

    !/

    INTEGER                            :: i, n, clock

    INTEGER, DIMENSION(:), ALLOCATABLE :: seed

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    CALL strace (ient, 'INIT_RANDOM_SEED')

#endif

    !/

    CALL random_seed(SIZE = n)

    ALLOCATE(seed(n))

    !

    CALL system_clock(count=clock)

    !

    seed = clock + 37 * (/ (i - 1, i = 1, n) /)

    CALL random_seed(put = seed)

    !

    DEALLOCATE(seed)

    !/

    !/ End of INIT_RANDOM_SEED ------------------------------------------- /

    !/

  END SUBROUTINE init_random_seed

  !/

  !/ ------------------------------------------------------------------- /

  !/

  !/ End of module W3SIC5MD -------------------------------------------- /

  !/

END MODULE w3sic5md

!/ ------------------------------------------------------------------- /