WW3/w3sic3md_8F90_source.html

#include "w3macros.h"

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

MODULE w3sic3md

  !/

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

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |           E. Rogers               |

  !/                  |           S. Zieger               |

  !/                  |           X. Zhao                 |

  !/                  |           S. Cheng                |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         04-Jan-2016 |

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

  !/

  !/    Updates:

  !/    29-May-2014 : Generalization with turbulent BL

  !/           (F.A. method imported from IC2 by E.R.)    ( version 5.01 )

  !/    04-Jan-2016 : Importing code provided by S. Cheng

  !/           (improved solution methods for Wang and Shen model)

  !/

  !  1. Purpose :

  !

  !     Calculate ice source term S_{ice} according to a viscoelastic sea

  !     ice model (Wang and Shen 2010).

  !

  !     Reference: Wang, R., and H. H. Shen (2010), Gravity waves

  !     propagating into an ice‐covered ocean: A viscoelastic model, J.

  !     Geophys. Res., 115, C06024, doi:10.1029/2009JC005591 .

  !

  !  2. Variables and types :

  !      Name              Type  Scope    Description

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

  !      IC3TABLE_CHENG    Int.  Public   Table of wave number k_r,

  !                                       attenuation k_i and group

  !                                       velocity cg

  !      IC3_DITK          R.A.  Private  Ice thickness increment

  !      IC3_MAXITK        R.A.  Private  Maximum ice thickness, the code

  !                                       may fail for situation with ice

  !                                       thickness larger than this value

  !

  !  3. Subroutines and functions :

  !

  !      Name              Type  Scope    Description

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

  !      W3SIC3            Subr. Public   Ice source term.

  !      BSDET             Func. Private  Calculate the determinant for

  !                                       the dispersion relation.

  !      WN_CMPLX_V1       Func. Private  Calculate complex wavenumber in

  !                                       ice

  !      WN_PRECALC_CHENG  Subr. Private  Calculate complex wavenumber in

  !                                       ice

  !      WN_CMPLX_HF       Func. Private  Like above, but for h-f waves

  !      CMPLX_ROOT_MULLER_CHENG Func. Private  Find root for complex

  !                                             numbers

  !      FUN_ZHAO          Func. Private  Wrapper function for FUNC0/FUNC1

  !      FUNC0_ZHAO        Func. Private

  !      FUNC1_ZHAO        Func. Private

  !      W3IC3WNCG         Subr. Public   Calculate kr,ki and cg for all

  !                                       frequency at each grid point

  !      IC3PRECALC_CHENG  Subr. Private  Calculate kr,ki and cg table for

  !                                       all frequencies and ice

  !                                       thickness from 0~IC3_MAXITK

  !      CGinIC3_CHENG           func. Private  Calculate group velocity

  !                                       related to IC3 model

  !      F_ZHAO            Func. Private  Wrapper function for double/

  !                                       quadruple precision

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

  !

  !  4. Subroutines and functions used :

  !

  !     See subroutine documentation.

  !

  !  5. Remarks :

  !

  !  6. Switches :

  !

  !     See subroutine documentation.

  !

  !  7. Source code :

  !/

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

  !/

  PUBLIC  ::  w3sic3, w3ic3wncg_v1, w3ic3wncg_cheng

  PRIVATE ::  wn_cmplx_v1, wn_cmplx_hf

  PRIVATE ::  cmplx_root_muller_v1, cmplx_root_muller_cheng

  PRIVATE ::  f_zhao_v1, f_zhao_cheng

  PRIVATE ::  func1_zhao, func0_zhao, bsdet

  INTEGER,SAVE ::  calledic3table = 0

  REAL,PRIVATE,PARAMETER    ::  ic3_ditk  = 0.01, ic3_maxitk = 3.

  PUBLIC  ::  ic3table_cheng

  PRIVATE ::  ic3precalc_cheng, cginic3_cheng


CONTAINS

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

  !/

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

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         11-Oct-2013 |

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

    !/

    !/    06-May-2013 : Origination (copied from SICE1)     ( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/

    !/        FIXME   : Move field input to W3SRCE and provide

    !/     (S.Zieger)   input parameter to W3SIC1 to make the subroutine

    !/                : versatile for point output processors ww3_outp

    !/                  and ww3_ounp.

    !/

    !/    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 a viscoelastic sea

    !     ice model (Wang and Shen 2010).

    !

    !     Reference: Wang, R., and H. H. Shen (2010), Gravity waves

    !     propagating into an ice‐covered ocean: A viscoelastic model, J.

    !     Geophys. Res., 115, C06024, doi:10.1029/2009JC005591 .

    !

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

    !

    !  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

    !

    !     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:

    !

    !     The Wang and Shen model is intended as a generalized model for

    !     various types of ice cover. It is a "continuum" model for

    !     which the same model is used from the ice edge to the ice

    !     interior. Though the ice types are expected to be very different

    !     from the edge to the interior, this is accomodated by the relative

    !     importance of the "effective viscosity" and the "modulus of

    !     elasticity". At the ice edge, where one finds frazil ice, pancake

    !     ice, or ice floes much smaller than the wave length, the "viscous"

    !     component of the model is believed to be most appropriate. At the

    !     interior, where one finds a continuous ice sheet, the "elastic

    !     model" component of the generalized visco-elastic model is

    !     expected to be appropriate. In addition to the case of continuous

    !     ice, Wang and Shen argue that the elastic model is also applicable

    !     to ice floes when the floe sizes are large relative to the

    !     wavelength. So to summarize,

    !     * frazil ice, pancake ice, and floes smaller than wavelength :

    !          viscosity dominates

    !     * continuous ice, and floes larger than wavelength :

    !          elasticity dominates

    !     * intermediate conditions: neither dominates

    !     All this is accomodated in WW3 by using non-uniform specification

    !     of viscosity and elasticity.

    !

    !     In the case where a user wishes to utilize only the "viscous

    !     model" aspect of Wang and Shen, and use an alternative scheme for

    !     continous ice and large ice floes, we allow this through the use

    !     of a user-defined namelist parameter "IC3MAXTHK". Floe size is

    !     not (at time of writing) an output available from ice model CICE,

    !     so we use the ice thickness as a way to anticipate the floe size.

    !     When ice thickness exceeds IC3MAXTHK, WW3 will use another model

    !     in place of the Wang and Shen formulation :

    !

    !     S_ice by F.A., an estimation of dissipation by turbulence

    !     at the ice-water interface. It uses only namelists for input, and

    !     no space/time varying input (though of course ice concentration is

    !     space/time varying). Unlike Liu et al. (IC2), it does not use

    !     ice thickness and does not yield a new C|Cg|k (i.e. it is non-

    !     dispersive), but it has the very nice feature of not requiring

    !     an eddy viscosity, which is a major drawback of the Liu et al.

    !     model. That is why we use it here, vs. Liu et al.

    !     (S_ice by Liu et al. and S_ice by F.A. are the two options

    !     available in IC2, i.e. w3sic2md.ftn)

    !

    !     At time of writing (March 23 2015), there may be some problems

    !     with the root-selection of IC3. For example with these settings:

    !     hice=1   ; rho_ice=917.0 ; emod=4.9e+12 ; visc=5e+7 ;

    !     h_water=deep (without using the IC3MAXTHK feature) the solution

    !     is rather irregular:

    !     T=11.11    k_i = 215.0e-5

    !     T=10       k_i = 266.0e-5

    !     T=9        k_i =   1.4e-5

    !     T=8.1      k_i =   1.7e-5

    !     With hice=0.1, ki is monotonically increasing in that range:

    !     T=11.11    k_i = 0.97e-5

    !     T=10       k_i = 2.64e-5

    !     T=9        k_i = 3.90e-5

    !     T=8.1      k_i = 4.47e-5

    !     Of course, when using IC3MAXTHK=0.1, the first example (hice=1)

    !     would switch to the "S_ice by F.A." model, and so this problem

    !     is circumvented.

    !

    !  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 (!/T1 switch).

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

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

    !

    !  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.

    !

    !     Code by S. Cheng sets NOICE=.TRUE. if ISNAN(ICECOEF1).

    !     Comments are "maps may not be compatible"

    !     This feature is not understood by me (ER) and so omitted.

    !

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

    !     On array size, S. Cheng says:

    !     Upon checking the origin of CG, I would say CG = CG_IC3, this is

    !     the topic ‘call W3IC3WNCG twice’. Recently I find they are not

    !     exactly the same due to calculation and smoothing of cg using

    !     several neighbor points. For different input array size, results

    !     are slightly different. In subr. w3wave, the size of input arrays

    !     is 0:NK+1, while in w3sic3, the size of all input arrays is 1:NK.

    !     This array size  difference is reflected in the resulting cg. The

    !     small difference in cg between calling IC3 twice or just calling

    !     once produces small difference in SWH. To eliminate this small

    !     difference, I suggest to keep CG instead of CG_IC3, as well as WN,

    !     WN_I, because other source terms use CG. I confirmed this change

    !     would make results of ICE the same whether calling twice or once

    !     by defining dimension WN_R, WN_I, CG_IC3 as 0:NK+1 instead of

    !     1:NK. Then CG_IC3 = CG.

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

    !     On optimization, S. Cheng says:

    !     For Wang and Shen’s model, D does not change in the loop

    !     corresponding to NSTEPS in subr. W3SRCE. I find the most efficient

    !     and easy way to speed up is that Add D and NSTEPS as inputs of

    !     W3SIC3

    !     If NSTEPS==1

    !          Current Wang and Shen’s model code above.

    !     ELSE

    !          S = D*A

    !     Endif

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

    !

    !  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 :

    !

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

    USE constants, ONLY: tpi, dwat, abmin, delab, sizefwtable,      &

         fwtable, grav

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

    !     USE WMMDATMD, ONLY: IMPROC, NMPERR ! WMMDATMD unavailable to outp

    USE w3servmd, ONLY: extcde

    USE w3gdatmd, ONLY: nk, nth, nspec, sig, mapwn, ic3pars, dden,  &

         flagll, ygrd, gtype, rlgtype

    USE w3idatmd, ONLY: icep1, icep2, icep3, icep4, icep5, icei,    &

         inflags2

#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

    !

    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

    INTEGER                    :: ith

#ifdef W3_T0

    REAL                   :: dout(nk,nth)

#endif

    INTEGER                    :: ikth, ik

    REAL                       :: icecoef1, icecoef2, icecoef3, &

         icecoef4, icecoef5, iceconc

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

    LOGICAL                    :: noice

    REAL                       :: viscm=1.83e-6

    REAL                       :: freq

    !  ............VISCM=1.83E-6 :  molecular viscosity of water at freezing

    REAL                    :: pturb, pvisc, dturb, dvisc,   &

         smooth, re, uorb, aorb, eb,   &

         deli1, deli2, fw, xi, fturb,  &

         maxthk, maxcnc, use_cheng,    &

         use_cgice, fixedhice,   &

         fixedvisc,fixeddens,fixedelas

    INTEGER                 :: ind, is, numin

    !

    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SIC3')

#endif

    !

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

    !

    !

    d        = 0.0

    d1d      = 0.0

    !

    wn_r     = wn

    wn_i     = 0.0

    cg_ic3   = 0.0

    cg_tmp   = 0.0

    !

    icecoef1 = 0.0

    icecoef2 = 0.0

    icecoef3 = 0.0

    icecoef4 = 0.0

    icecoef5 = 0.0

    iceconc  = 0.0

    !

    !     Rename variables to make code easier to read.

    maxthk=ic3pars(1)

    maxcnc=ic3pars(8)

    use_cheng=ic3pars(9)

    use_cgice=ic3pars(12)

    fixedhice=ic3pars(13)

    fixedvisc=ic3pars(14)

    fixeddens=ic3pars(15)

    fixedelas=ic3pars(16)


    !     --- Error checking for input ----------------------------------- /

    !     --- Allow one and only one input option for each variable ------ /

    numin=0

    IF (inflags2(-7)) numin=numin+1

    IF (fixedhice.GE.0.0) numin=numin+1

    IF (numin.NE.1) THEN

      IF ( iaproc .EQ. naperr )                       &

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

      CALL extcde(2)

    ENDIF


    numin=0

    IF (inflags2(-6)) numin=numin+1

    IF (fixedvisc.GE.0.0) numin=numin+1

    IF (numin.NE.1) THEN

      IF ( iaproc .EQ. naperr )                       &

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

      CALL extcde(2)

    ENDIF


    numin=0

    IF (inflags2(-5)) numin=numin+1

    IF (fixeddens.GE.0.0) numin=numin+1

    IF (numin.NE.1) THEN

      IF ( iaproc .EQ. naperr )                       &

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

      CALL extcde(2)

    ENDIF


    numin=0

    IF (inflags2(-4)) numin=numin+1

    IF (fixedelas.GE.0.0) numin=numin+1

    IF (numin.NE.1) THEN

      IF ( iaproc .EQ. naperr )                       &

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

      CALL extcde(2)

    ENDIF


    !     --- Set local value to be used subsequently (ICEPx variables

    !         are not used beyond this point). --------------------------- /

    IF (inflags2(-7)) THEN

      icecoef1 = icep1(ix,iy) ! ice thickness

    ELSE

      icecoef1 = fixedhice

    ENDIF


    IF (inflags2(-6)) THEN

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

    ELSE

      icecoef2 = fixedvisc

    ENDIF


    IF (inflags2(-5)) THEN

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

    ELSE

      icecoef3 = fixeddens

    ENDIF


    IF (inflags2(-4)) THEN

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

    ELSE

      icecoef4 = fixedelas

    ENDIF


    !     ICECOEF5 = ICEP5(IX,IY) ! ICEP5 is inactive in W3SIC3


    IF (inflags2(4))  iceconc = icei(ix,iy)


    !

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

    !

    noice=.false.

    IF (icecoef1==0.0) noice=.true.

    IF (inflags2(4).AND.(iceconc==0.0)) noice=.true.


    IF ( noice ) THEN


      d1d=0.0

      !

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

    ELSEIF ( use_cheng==1.0 .AND.  &

         ((icecoef1.LE.maxthk).OR.(iceconc.LE.maxcnc)) ) THEN


      ! 2.a Write test output ---------------------------------------------- /

#ifdef W3_T38

      WRITE (ndst,9000) depth,icecoef1,icecoef2,icecoef3,icecoef4

#endif


      ! 2.b Make calculations using Cheng routines ------------------------- /


      !     --- Input to routine (part 1): 6 ice parameters from single

      !         precision variables. ---------------------------------------


      CALL w3ic3wncg_cheng(wn_r, wn_i, cg_ic3, icecoef1, icecoef2, &

           icecoef3, icecoef4, depth)

      !

      !    --- calculate source term --------------------------------------- /

      !    --- see Remarks section re: array size -------------------------- /

      IF ( use_cgice==1.0 ) THEN

        cg_tmp=cg_ic3

      ELSE

        cg_tmp=cg

      ENDIF

      DO ik=1, nk

        !            recall that D=S/E=-2*Cg*k_i

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


      END DO


    ELSEIF ( (icecoef1 .LE. maxthk) .OR. (iceconc .LE. maxcnc) ) THEN

      !.......... e.g. if ice thickness is .le. 10 cm

      !...............or concentration is .le. 1.0

      !

      ! 2.a Write test output ---------------------------------------------- /

#ifdef W3_T38

      WRITE (ndst,9000) depth,icecoef1,icecoef2,icecoef3,icecoef4

#endif

      !

      ! 2.b Make calculations using original routines ---------------------- /

      !     --- Input to routine (part 1): 6 ice parameters from single

      !         precision variables. ---------------------------------------


      CALL w3ic3wncg_v1(wn_r, wn_i, cg_ic3, icecoef1, icecoef2, &

           icecoef3, icecoef4, depth     )

      !

      !    --- calculate source term --------------------------------------- /

      IF ( use_cgice==1.0 ) THEN

        cg_tmp=cg_ic3

      ELSE

        cg_tmp=cg

      ENDIF

      DO ik=1, nk

        !            recall that D=S/E=-2*Cg*k_i

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

      END DO

      !

    ELSE ! .. e.g. if ice thickness is .gt. 10 cm

      ! Alternative by F.A., see Remarks section.

      IF (ic3pars(2).GT.0.) THEN

        uorb=0.

        aorb=0.

        fturb = ic3pars(2)

        IF (ic3pars(7).GT.0) THEN

          IF (ygrd(iy,ix).LT.0.AND.gtype.EQ.rlgtype.AND.flagll) &

               fturb = ic3pars(7)

        END IF

        DO ik=1, nk

          eb  = 0.

          DO ith=1, nth

            is=ith+(ik-1)*nth

            eb  = eb  + a(is)

          END DO

          !

          !  UORB and AORB are the variances of the orbital

          !  velocity and surface elevation

          !

          uorb = uorb + eb *sig(ik)**2 * dden(ik) / cg(ik)

          aorb = aorb + eb             * dden(ik) / cg(ik)

          !deep water only

        END DO

        !

        aorb = 2*sqrt(aorb)  ! significant amplitude

        uorb = 2*sqrt(uorb)  ! significant amplitude


        re = uorb*aorb / viscm

        smooth = 0.5*tanh((re-ic3pars(4))/ic3pars(5))

        pturb=(0.5+smooth)

        pvisc=(0.5-smooth)


        xi=(alog10(max(aorb/ic3pars(3),3.))-abmin)/delab

        ind  = min(sizefwtable-1, int(xi))

        deli1= min(1. ,xi-float(ind))

        deli2= 1. - deli1

        fw =fwtable(ind)*deli2+fwtable(ind+1)*deli1

        dturb=-1.* fturb*fw*uorb/grav

      ELSE ! so case of IC3PARS(2).LE.0.

        dturb = 0.

      END IF ! IF (IC3PARS(2).GT.0.)


      DO ik=1, nk

        dvisc = -1. *ic3pars(6) * wn(ik) * sqrt(viscm* sig(ik) / 2.)

        d1d(ik) = pturb*dturb*sig(ik)**2 +  pvisc*dvisc

      END DO


    END IF !   IF ( NOICE ) THEN


    ! 2.c Fill diagional matrix ------------------------------------------ /

    !

    DO ikth=1, nspec

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

    END DO


    !

    ! sign convention (example):

    !        S is from -10e-3 to 0

    !        A is from 0 to 10

    ! See Remarks section re: optimization

    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 W3SIC3 : '/               &

         '     ',a,' REQUIRED ONCE, BUT WAS PROVIDED BY USER '/ &

         '     ',i4,' TIMES.'/)

    !

#ifdef W3_T

9000 FORMAT (' TEST W3SIC3 : depth and 4 ice coef. : ',5e10.3)

#endif

    !/

    !/ End of W3SIC3 ----------------------------------------------------- /

    !/

  END SUBROUTINE w3sic3

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

  !/

  SUBROUTINE w3ic3wncg_v1(WN_R,WN_I,CG,ICE1,ICE2,ICE3,ICE4,DPT)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         25-Oct-2013 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Calculation of complex wavenumber for waves in ice. Outsourced

    !     from W3SIC3 to allow update on wavenumbers and  group

    !     velocities at each time step an ice parameter is updated.

    !

    !  2. Method :

    !

    !     <text here>

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !       WN_R   R. A.  I/O  Wave number (real part)

    !       WN_I   R. A.  I/O  Wave number (imag. part=wave attenuation)

    !       CG     R. A.  I/O  Group velocity

    !       ICE1   REAL    I   Thickness of ice                 [in m]

    !       ICE2   REAL    I   Effective viscosity of ice       [in m2/s]

    !       ICE3   REAL    I   Density of ice                  [in kg/m3]

    !       ICE4   REAL    I   Effective shear modulus of ice   [in Pa]

    !       DPT    REAL    I   Water depth                      [in m]

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

    !

    !  4. Subroutines used :

    !

    !      Name    Type  Module    Description

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

    !     WAVNU1   Subr. W3DISPMD   Wavenumber for waves in open water.

    !     WN_CMPLX Func. W3SIC3MD   Complex wavenumber for waves in ice.

    !     WN_CMPLX_HF Func. W3SIC3MD   Like WN_CMPLX, but for h-f

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

    !

    !  5. Called by :

    !

    !      Name    Type  Module   Description

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

    !      W3SIC3  Subr. W3SIC3MD Ice source term.

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

    !

    !  6. Error messages :

    !

    !     See FORMAT 900.

    !

    !  7. Remarks :

    !

    !     Optional: Cap WN_I at 2.0E-4, since in simple tests with "normal

    !     resolution" (not finer than 1 km), WW3 has trouble resolving the

    !     dissipation if k_i>2e-4. Also, very large values of dissipation

    !     (e.g. k_i=100e-4) with IC3 occurs more in the higher frequencies

    !     which makes WW3 slow down quite a bit. This is done via ICEKILIM

    !     in the namelists.

    !

    !     This function does not get used in update by S. Cheng.

    !     It should be removed if/when "V1" routines are removed.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Source code :

    !/

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

    !/

    USE w3gdatmd, ONLY: nk, sig, ic3pars

    USE w3dispmd, ONLY: wavnu1

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

    USE constants, ONLY: tpi

    !/

    IMPLICIT NONE

    !/

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

    REAL, INTENT(IN)   :: ice1, ice2, ice3, ice4, dpt


    INTEGER            :: ik, kl,ku

    REAL, ALLOCATABLE  :: sigma(:),cg_ic3(:)

    REAL               :: k_ocean, cg_ocean

    DOUBLE PRECISION   :: kh, k_noice, hwat, hice, nu, dice, es_mod

    DOUBLE PRECISION,PARAMETER :: khmax = 18.0d0 ! 18=OK, 19=fails

    DOUBLE COMPLEX     :: wncomplex,wncomplex_old

    REAL               :: stensec

    REAL               :: ic3hilim,ic3kilim

    !

    ALLOCATE( cg_ic3( SIZE(cg) ) )

    ALLOCATE(  sigma( SIZE(cg) ) )

    cg_ic3  = 0.

    sigma   = 0.

    stensec=tpi/10.0 ! sigma for T=10 sec


    ic3hilim=ic3pars(10)

    ic3kilim=ic3pars(11)

    !

    !     --- Input to routine (part 1): set 6 double precision variables

    !        using single precision variables. -------------------------- /

    hwat    = dble(dpt)     ! water depth

    hice    = dble(ice1)    ! ice thickness

    nu      = dble(ice2)    ! "effective viscosity" parameter

    dice    = dble(ice3)    ! density of ice

    es_mod  = dble(ice4)    ! effective shear modulus of ice


    ! Optional: limit ice thickness

    hice=min(dble(ic3hilim),hice)


    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

      WRITE(ndse,900)

      CALL extcde(3)

    END IF

    !

    wncomplex_old=cmplx(0.0d0,0.0d0)


    DO ik = kl,ku

      !     --- Input to routine (part 2): set 2 double precision variables

      !         using single precision variable. --------------------------- /

      CALL wavnu1(sigma(ik),dpt,k_ocean,cg_ocean)

      k_noice = dble(k_ocean)

      !

      !     --- Muller Method fails for deep water: workaround follows ----- /

      kh        = k_noice * hwat ! kh w/out ice

      IF (kh.GT.khmax) THEN

        hwat   = khmax / k_noice

      ENDIF

      !     --- Calculate complex wavenumber ------------------------------- /


      IF((ik.GT.kl).AND.(sigma(ik).GT.stensec))THEN


        wncomplex = wn_cmplx_hf(dble(sigma(ik)),k_noice,es_mod,nu, &

             dice,hice,hwat,dble(sigma(ik-1)),wncomplex_old)

        wncomplex_old=wncomplex


      ELSE


        wncomplex = wn_cmplx_v1(dble(sigma(ik)),k_noice,es_mod,nu, &

             dice,hice,hwat)

        wncomplex_old=wncomplex


      ENDIF


      !     --- Output from function is type of DOUBLE COMPLEX. Set

      !         precision of imaginary to single precision array element --- /

      wn_i(ik)  = real(aimag(wncomplex))


      ! Optional : limit ki

      wn_i(ik)  = min(wn_i(ik),ic3kilim) ! see Remarks above

      wn_r(ik)  = real(wncomplex)


    END DO

    !     --- Update group velocitiy ----

    cg_ic3 = delta(sigma) / delta(wn_r)

    cg     = cg_ic3


    DEALLOCATE(cg_ic3)

    !

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

         '     CANNOT DETERMINE BOUNDS OF WAVENUMBER ARRAY.')

    !

    !/

  END SUBROUTINE w3ic3wncg_v1

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

  !/

  FUNCTION wn_cmplx_v1(SIGMA,WN_O,ES,NU,DICE,HICE,DEPTH)    RESULT(WN)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         30-Oct-2013 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Calculate complex wavenumber for waves in ice.

    !

    !  2. Method :

    !

    !     Wang and Shen (JGR 2010)

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !       WN     CMPLX DBL O  Wave number (imag. part=wave attenuation)

    !       SIGMA  REAL DBL  I  Wave angular frequency           [in rad]

    !       WN_O   REAL DBL  I  Wave number (open water)

    !       ES     REAL DBL  I  Effective shear modulus of ice   [in Pa]

    !       NU     REAL DBL  I  Effective viscosity of ice       [in m2/s]

    !       DICE   REAL DBL  I  Density of ice                  [in kg/m3]

    !       HICE   REAL DBL  I  Thickness of ice                 [in m]

    !       DEPTH  REAL DBL  I  Water depth                      [in m]

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

    !

    !  4. Subroutines used :

    !

    !      Name               Type  Module   Description

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

    !      CMPLX_ROOT_MULLER_CHENG  Func. W3SIC3MD Find root for complex

    !                                         wavenumbers for waves in ice.

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

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

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

    !     W3IC3WNCG_V1    Subr. W3SIC3MD Ice source term.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !   Hayley Shen says,

    !     We have determined that it may not be necessary to use curve

    !     fitting or lookup tables to get the group velocity and the

    !     attenuation coefficient. Attached is a short report with some

    !     sample numerical solutions. To implement the viscoelastic model,

    !     there are 4 fortran programs. According to Xin Zhao, the graduate

    !     student, it is very fast to find roots. I suggest that perhaps you

    !     try the pure viscous case by setting G=0 to start with. nu can be

    !     set at 0.05*ice concentration (m^2/s) to begin with, because for

    !     grease ice Newyear's data showed nu to be about 0.02-0.03 m^2/s.

    !     By setting G=0 in you get exactly the Keller model for pure

    !     viscous layer.

    !

    !   This routine provides the initial guess according to the parameters

    !   of the present case. T>10s use open water, T<10s cases, calculate

    !   T=10s first using open water as the initial guess.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

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

    USE constants, ONLY: tpi

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    DOUBLE PRECISION, INTENT(IN) :: sigma,wn_o,es,nu,dice,hice,depth

    DOUBLE COMPLEX               :: wn                      ! RESULT

    !/

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

    !/ Local parameters

    !/

    INTEGER            :: i, nsub

    DOUBLE PRECISION   :: tt, ts, t

    DOUBLE COMPLEX     :: x0, x1, x2, wn0

    !/

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

    t    = dble(tpi) / sigma

    ts   = 10.

    nsub = int((ts-t) * 10.)

    !/

    IF (hice<0.001) THEN

      wn  = cmplx(wn_o,0.)

    ELSE IF (t.LT.ts) THEN

      x0  = 0.01

      x1  = 0.1

      x2  = 1.0

      wn0 = cmplx_root_muller_v1(x0,x1,x2,0,dble(tpi)/ts, &

           es,nu,dice,hice,depth    )

      x0  = 0.90 * wn0

      x1  = wn0

      x2  = 1.1*wn0

      wn  = cmplx_root_muller_v1(x0,x1,x2,1,dble(tpi)/ts, &

           es,nu,dice,hice,depth    )

      DO i=1,nsub

        x0 = 0.90 * wn

        x1 = wn

        x2 = 1.1 * wn

        tt = ts - (ts-t) / real(nsub) * real(i)

        wn = cmplx_root_muller_v1(x0,x1,x2,1,dble(tpi)/tt, &

             es,nu,dice,hice,depth    )

      ENDDO

    ELSE

      x0  = 0.01

      x1  = 0.1

      x2  = 1.0

      wn0 = cmplx_root_muller_v1(x0,x1,x2,0,sigma,        &

           es,nu,dice,hice,depth    )

      x0  = 0.8 * wn0

      x1  = wn0

      x2  = 1.2 * wn0

      wn  = cmplx_root_muller_v1(x0,x1,x2,1,sigma,        &

           es,nu,dice,hice,depth    )

    ENDIF

    !/

  END FUNCTION wn_cmplx_v1

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

  !/

  FUNCTION wn_cmplx_hf(SIGMA,WN_O,ES,NU,DICE,HICE,DEPTH,SIGMA_LAST, &

       WN_LAST)     RESULT(WN)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. Shen                 |

    !/                  |           E. Rogers               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         17-Apr-2014 |

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

    !/

    !/    15-Jan-2014 : Origination (from WN_CMPLXA.f90)      (H. Shen)

    !/    17-Apr-2014 : Import to WW3                         (E. Rogers)

    !/

    !  1. Purpose :

    !

    !     Calculate complex wavenumber for waves in ice.

    !

    !  2. Method :

    !

    !     Wang and Shen (JGR 2010)

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !       WN     CMPLX DBL O  Wave number (imag. part=wave attenuation)

    !       SIGMA  REAL DBL  I  Wave angular frequency           [in rad]

    !       WN_O   REAL DBL  I  Wave number (open water)

    !       ES     REAL DBL  I  Effective shear modulus of ice   [in Pa]

    !       NU     REAL DBL  I  Effective viscosity of ice       [in m2/s]

    !       DICE   REAL DBL  I  Density of ice                  [in kg/m3]

    !       HICE   REAL DBL  I  Thickness of ice                 [in m]

    !       DEPTH  REAL DBL  I  Water depth                      [in m]

    !       SIGMA_LAST REAL DBL I : Like SIGMA, but of last IK

    !       WN_LAST    REAL DBL I : WN_O of last IK

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

    !

    !  4. Subroutines used :

    !

    !      Name               Type  Module   Description

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

    !      CMPLX_ROOT_MULLER_CHENG  Func. W3SIC3MD Find root for complex

    !                                         wavenumbers for waves in ice.

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

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

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

    !      W3IC3WNCG_V1  Subr. W3SIC3MD Ice source term.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Zhao and Shen.

    !     See notes in FUNCTION WN_CMPLX, not repeated here.

    !     New in this function, Hayley Shen says (Jan 15 2014) :

    !     "To speed up the computation, we need to add a new function

    !     WN_CMPLXA (attached) into the earlier version of the MODULE

    !     W3SIC3MD. When wave period T>=10s, we call old function WN_CMPLX

    !     directly. When T<10s, call the new function WN_CMPLXA with last

    !     calculation step's information: last complex wave number, last

    !     angular wave frequency. The calculation should be from large T

    !     to small T."

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

    USE constants, ONLY: tpi

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    DOUBLE PRECISION, INTENT(IN) :: sigma,wn_o,es,nu,dice,hice,depth

    DOUBLE PRECISION, INTENT(IN) :: sigma_last

    DOUBLE COMPLEX,   INTENT(IN) :: wn_last

    DOUBLE COMPLEX               :: wn                      ! RESULT

    !/

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

    !/ Local parameters

    !/

    INTEGER            :: i, nsub

    DOUBLE PRECISION   :: tt, ts, t

    DOUBLE COMPLEX     :: x0, x1, x2, wn0

    !/

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

    t    = dble(tpi) / sigma

    ts   = dble(tpi) / sigma_last

    nsub = int((ts-t) * 10.)

    !/

    IF (hice<0.001) THEN

      wn  = cmplx(wn_o,0.)

    ELSE

      x0  = 0.90 * wn_last

      x1  = wn_last

      x2  = 1.1 * wn_last

      wn  = cmplx_root_muller_v1(x0,x1,x2,1,dble(tpi)/ts, &

           es,nu,dice,hice,depth    )

      DO i=1,nsub

        x0 = 0.90 * wn

        x1 = wn

        x2 = 1.1 * wn

        tt = ts - (ts-t) / real(nsub) * real(i)

        wn = cmplx_root_muller_v1(x0,x1,x2,1,dble(tpi)/tt, &

             es,nu,dice,hice,depth    )

      ENDDO

    ENDIF

    !/

  END FUNCTION wn_cmplx_hf

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

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

  !/

  FUNCTION cmplx_root_muller_v1(X0, X1, X2, JUDGE, SIGMA, ES, NU, &

       DICE, HICE, DEPTH)           RESULT(P3)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         30-Oct-2013 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Find root.

    !

    !  2. Method :

    !

    !     Muller method for complex equations is a recursive approximation

    !     with initial guess X0, X1, and X2. To the initial guesses a

    !     quadratic parabola is fitted.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !       P3    CMPLX DBL O  Approximation for the root problem

    !       X0    CMPLX DBL I  Initial guess variable

    !       X1    CMPLX DBL I  Initial guess variable

    !       X2    CMPLX DBL I  Initial guess variable

    !       JUDGE INTEGER   I  "switch variable" for F_ZHAO

    !       SIGMA DOUBLE    I  Wave angular frequency

    !       ES    DOUBLE    I  Effective shear modulus of ice

    !       NU    DOUBLE    I  Effective viscosity of ice       [in m2/s]

    !       DICE  DOUBLE    I  Density of ice                   [in kg/m3]

    !       HICE  DOUBLE    I  Thickness of ice                 [in m]

    !       DEPTH DOUBLE    I  Water depth                      [in m]

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

    !

    !  4. Subroutines used :

    !      Name      Type  Module   Description

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

    !      F_ZHAO    Func. W3SIC3MD Wrapper function for root finding.

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

    !

    !  5. Called by :

    !

    !      Name               Type  Module   Description

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

    !      WN_CMPLX_V1             Find root for complex wave-

    !      WN_CMPLX_HF             numbers for waves in ice.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

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

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    DOUBLE COMPLEX                :: p3                     ! RESULT

    DOUBLE COMPLEX, INTENT(IN)    :: x0,x1,x2

    DOUBLE PRECISION, INTENT(IN)  :: sigma,es,nu,dice,hice,depth

    INTEGER, INTENT(IN)           :: judge

    !/

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

    !/ Local parameters

    !/

    INTEGER                   :: i

    INTEGER, PARAMETER        :: imax = 1000

    DOUBLE PRECISION          :: dlta,epsi

    DOUBLE COMPLEX            :: p0,p1,p2

    DOUBLE COMPLEX            :: y0,y1,y2,y3

    DOUBLE COMPLEX            :: a,b,c,q,disc,den1,den2

    !/

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

    p0      = x0

    p1      = x1

    p2      = x2

    p3      = 0.0

    !

    i       = 0

    epsi    = 1.e-5

    dlta    = 1.e-5

    y0      = f_zhao_v1(p0,judge,sigma,es,nu,dice,hice,depth)

    y1      = f_zhao_v1(p1,judge,sigma,es,nu,dice,hice,depth)

    y2      = f_zhao_v1(p2,judge,sigma,es,nu,dice,hice,depth)

    !

    DO i = 1,imax

      q = (p2 - p1) / (p1 - p0)

      a = q * y2 - q * (1.+q) * y1 + q**2. * y0

      b = (2. * q + 1.) * y2 - (1 + q)**2. * y1 + q**2. * y0

      c = (1. + q) * y2

      !

      IF ( abs(a).NE.0. ) THEN


        disc = b**2. - 4 * a * c;

        !

        den1 = ( b + sqrt( disc ) )

        den2 = ( b - sqrt( disc ) )

        !

        IF ( abs( den1 ) .LT. abs( den2 ) )THEN

          p3 = p2 - (p2 - p1) * (2 * c / den2)

        ELSE

          p3 = p2 - (p2 - p1) * (2 * c / den1)

        ENDIF

        !

      ELSE

        !

        IF ( abs(b) .NE. 0. )THEN

          p3 = p2 - (p2 - p1) * (c / b)

        ELSE

          WRITE(ndse,800)

          WRITE(ndse,801)x0,x1,x2

          WRITE(ndse,802)sigma,es,nu,dice,hice,depth

          WRITE(ndse,803)judge

          CALL extcde(2)

        ENDIF

      ENDIF


      y3 = f_zhao_v1(p3,judge,sigma,es,nu,dice,hice,depth);


      IF ( abs(p3-p2).LT.dlta .OR. abs(y3).LT.epsi ) THEN

        RETURN

      ENDIF


      p0 = p1

      p1 = p2

      p2 = p3


      y0 = y1

      y1 = y2

      y2 = y3


    ENDDO

    !

    WRITE(ndse,800)

    WRITE(ndse,801)x0,x1,x2

    WRITE(ndse,802)sigma,es,nu,dice,hice,depth

    WRITE(ndse,803)judge

    CALL extcde(2)

    !

800 FORMAT (/' *** WAVEWATCH III ERROR IN W3SIC3_CMPLX_ROOT_MULLER'/&

         ' :      MULLER METHOD FAILED TO FIND ROOT.'            )

801 FORMAT (/'X0,X1,X2 = ',3(1x,'(',f10.5,',',f10.5,')'))

802 FORMAT (/'SIGMA,ES,NU,DICE,HICE,DEPTH = ',6(1x,f10.5))

803 FORMAT (/'JUDGE = ',i5)

    !/

  END FUNCTION cmplx_root_muller_v1

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

  !/

  FUNCTION f_zhao_v1(X,JUDGE,SIGMA,ES,NU,DICE,HICE,DEPTH) &

       result(fzhao)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         30-Oct-2013 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Decide whether to call sub-function.

    !

    !  2. Method :

    !

    !     Decide based on value of integer "JUDGE"

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !      FZHAO COMPL8  O    Result (double complex)

    !      X     CMPLX8  I    Approximate result (double complex)

    !      JUDGE INTEGR  I    Switch variable

    !      SIGMA DOUBLE  I    Wave angular frequency

    !      ES    DOUBLE  I    Effective shear modulus

    !      NU    DOUBLE  I    Effective viscosity parameter

    !      DICE  DOUBLE  I    Density of ice

    !      HICE  DOUBLE  I    Thickness of ice

    !      DEPTH DOUBLE  I    Water depth

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

    !

    !  4. Subroutines used :

    !

    !      Name        Type  Module   Description

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

    !      FUNC0_ZHAO  Func. W3SIC3MD Function to find root.

    !      FUNC1_ZHAO  Func. W3SIC3MD Function to find root.

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

    !

    !  5. Called by :

    !

    !      Name               Type  Module   Description

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

    !      CMPLX_ROOT_MULLER_V1  Func. W3SIC3MD Find root for complex wave-

    !                                        numbers for waves in ice.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !/

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

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    INTEGER, INTENT(IN)          :: judge

    DOUBLE PRECISION, INTENT(IN) :: sigma,es,nu,dice,hice,depth

    DOUBLE COMPLEX, INTENT(IN)   :: x

    DOUBLE COMPLEX               :: fzhao  ! RESULT

    !/

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

    IF (judge.EQ.0) THEN

      fzhao = func0_zhao(x,sigma,depth)

    ELSE

      fzhao = func1_zhao(x,sigma,es,nu,dice,hice,depth)

    ENDIF

    !

  END FUNCTION f_zhao_v1

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

  !/

  FUNCTION func0_zhao(WN, SIGMA, DEPTH)                RESULT(FUNC0)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         30-Oct-2013 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Calculate the difference between the left and right side

    !     of the dispersion relation. It is called by the Muller method.

    !

    !  2. Method :

    !

    !     <text here>

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !      FUNC0  COMPL DBL  O    Result (double complex)

    !      WN     CMPLX DBL  I    Complex wavenumber

    !      SIGMA  DOUBLE     I    Wave angular frequency

    !      DEPTH  DOUBLE     I    Water depth [in m]

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

    !

    !  4. Subroutines used :

    !

    !  5. Called by :

    !

    !      Name    Type  Module   Description

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

    !      F_ZHAO_V1  Func. W3SIC3MD Function for computation of complex

    !                             wavenumbers for waves in ice.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !   This function does not get used in update by S. Cheng.

    !     It should be removed if/when "V1" routines are removed.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

    !/

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

    USE constants, ONLY: grav

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    DOUBLE COMPLEX, INTENT(IN)    :: wn

    DOUBLE PRECISION, INTENT(IN)  :: sigma, depth

    DOUBLE COMPLEX                :: func0                 ! RESULT

    !/

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

    !/ Local parameters

    !/

    DOUBLE COMPLEX   :: th

    !/

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

    IF (real(wn*depth).LE.4.) THEN

      th = (exp(wn*depth)-exp(-wn*depth)) &

           / (exp(wn*depth)+exp(-wn*depth))

      func0 = sigma**2. - th * wn * dble(grav)

    ELSE

      func0 = sigma**2. -      wn * dble(grav)

    END IF

    !/

  END FUNCTION func0_zhao

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

  !/

  FUNCTION func1_zhao(WN,SIGMA,ES,NU,DICE,HICE,DEPTH)  RESULT(FUNC1)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         11-Oct-2013 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     <text here>

    !

    !  2. Method :

    !

    !     <text here>

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !      FUNC1 CMPLX DBL O  Result (double complex)

    !      WN    CMPLX DBL I  Wavenumber (double complex)

    !      W     REAL  DBL I  Wave angular frequency

    !      ES    REAL  DBL I  Effective shear modulus on ice

    !      NU    REAL  DBL I  Effective viscosity

    !      DICE  REAL  DBL I  Density of ice

    !      HICE  REAL  DBL I  Thickness of ice

    !      DEPTH REAL  DBL I  Water depth

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

    !

    !  4. Subroutines used :

    !

    !      Name    Type  Module    Description

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

    !      BSDET   Func. W3SIC3MD  Calculates the determinant for the

    !                              dispersion relation.

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

    !

    !  5. Called by :

    !

    !      Name    Type  Module   Description

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

    !      F_ZHAO_V1  Func. W3SIC3MD Function for computation of complex

    !                             wavenumbers for waves in ice.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !   This function does not get used in update by S. Cheng.

    !     It should be removed if/when "V1" routines are removed.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !/

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

    USE constants, ONLY: grav, dwat

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    DOUBLE COMPLEX, INTENT(IN)   :: wn

    DOUBLE PRECISION, INTENT(IN) :: sigma, es, nu, dice, hice, depth

    DOUBLE COMPLEX               :: func1                   ! RESULT

    !/

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

    !/ Local parameters

    !/

    DOUBLE COMPLEX   :: ve,alpha,n,m,l,sk,ck,sa,ca,th,thh

    DOUBLE COMPLEX   :: aa(4,4)

    !/

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

    ve    = cmplx( nu, es/dice/sigma )

    alpha = sqrt( wn**2. - sigma/ve * cmplx(0.,1.) )

    n     = sigma + 2. * ve * wn**2. * cmplx(0.,1.)

    l     = 2 * wn * alpha * sigma * ve

    sk    = (exp(wn*hice)-exp(-wn*hice))/2.

    ck    = (exp(wn*hice)+exp(-wn*hice))/2.

    sa    = (exp(alpha*hice)-exp(-alpha*hice))/2.

    ca    = (exp(alpha*hice)+exp(-alpha*hice))/2.

    !

    IF (real(wn*depth).LE.4.) THEN

      th  = (exp(wn*depth)-exp(-wn*depth))                        &

           / (exp(wn*depth)+exp(-wn*depth))

      thh = ( exp(wn*(depth-hice)) - exp(-wn*(depth-hice)) )      &

           / ( exp(wn*(depth-hice)) + exp(-wn*(depth-hice)) )

    ELSE

      th  = 1.0

      thh = 1.0

    END IF

    !

    m     = (dble(dwat)/dice - 1) * dble(grav) * wn               &

         - dble(dwat) / dice * sigma**2 / th

    !

    IF (es.GT.1.e7) THEN

      aa(1,1) = 0.

      aa(1,2) = 2 * cmplx(0.,1.) * wn**2.

      aa(1,3) = alpha**2. + wn**2.

      aa(1,4) = 0.

      !

      aa(2,1) = n * sigma

      aa(2,2) = -wn * dble(grav)

      aa(2,3) = cmplx(0.,1.) * wn * dble(grav)

      aa(2,4) = l

      !

      aa(3,1) = -2. * cmplx(0.,1.) * wn**2. * sk

      aa(3,2) =  2. * cmplx(0.,1.) * wn**2. * ck

      aa(3,3) =  (alpha**2. + wn**2.) * ca

      aa(3,4) = -(alpha**2. + wn**2.) * sa

      !

      aa(4,1) =   n * sigma * ck - m * sk

      aa(4,2) = - n * sigma * sk + m * ck

      aa(4,3) = -cmplx(0.,1.) * m * ca - l * sa

      aa(4,4) =  cmplx(0.,1.) * m * sa + l * ca

      !

      func1   = bsdet(aa,4)

    ELSE

      func1 = sigma**2. - th*wn*dble(grav) - th*dice/dble(dwat)* &

           (wn**2.*dble(grav)**2.*sk*sa - (n**4. + 16.* &

           ve**4.*wn**6.*alpha**2.)*sk*sa - 8. &

           *wn**3.*alpha*ve**2.*n**2.*(ck*ca-1.))/(4.*wn**3. &

           *alpha*ve**2.*sk*ca+n**2.*sa*ck-dble(grav)*wn*sk*sa)

    ENDIF

    !/

  END FUNCTION func1_zhao

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

  !/

  FUNCTION bsdet(AA, N)                                  RESULT(DET)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         11-Oct-2013 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     This subroutine calculates the determinant for the

    !     dispersion relation.

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !      Parameter list

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

    !      AA    R.A.       I/O  Square array type of REAL

    !      N     INT        I    Size of array (number of rows/cols)

    !      DET   CMPLX DBLE I/O  Determinant (double complex)

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

    !

    !  4. Subroutines used :

    !

    !  5. Called by :

    !

    !      Name        Type  Module    Description

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

    !      FUNC1_ZHAO  Func. W3SIC3MD  Function for computation of complex

    !                                  wavenumbers for waves in ice.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !   This function does not get used in update by S. Cheng.

    !     It should be removed if/when "V1" routines are removed.

    !

    !  8. Source code :

    !

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

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    INTEGER, INTENT(IN)        :: n

    DOUBLE COMPLEX, INTENT(IN) :: aa(n,n)

    DOUBLE COMPLEX             :: det                       ! RESULT

    !/

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

    !/ Local parameters

    !/

    INTEGER                    :: k, i, j, is, js

    DOUBLE COMPLEX             :: f, d, mat(n,n)

    DOUBLE PRECISION           :: q

    !/

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

    mat = aa

    f = 1.0

    det = 1.0

    loop100: DO k = 1,n-1

      q = 0.0

      loop10a: DO i = k,n

        loop10b: DO j = k,n

          IF (abs(mat(i,j)).GT.q) THEN

            q = abs(mat(i,j))

            is = i

            js = j

          END IF

        END DO loop10b

      END DO loop10a

      IF (q+1.0.EQ.1.0) THEN

        det = 0.0

        RETURN

      END IF

      IF (is.NE.k) THEN

        f = -f

        loop20: DO j = k,n

          d         = mat(k,j)

          mat(k,j)  = mat(is,j)

          mat(is,j) = d

        END DO loop20

      END IF

      IF (js.NE.k) THEN

        f = -f

        loop30: DO i = k,n

          d         = mat(i,js)

          mat(i,js) = mat(i,k)

          mat(i,k)  = d

        END DO loop30

      END IF

      det = det * mat(k,k)

      loop50: DO i = k+1,n

        d = mat(i,k) / mat(k,k)

        loop40: DO j = k+1,n

          mat(i,j) = mat(i,j) - d * mat(k,j)

        END DO loop40

      END DO loop50

    END DO loop100

    !/

    det =  f * det * mat(n,n)

    !/

    !/ End of BSDET ------------------------------------------------------ /

    !/

  END FUNCTION bsdet

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

  FUNCTION delta(X) RESULT(DX)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         22-Oct-2013 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.12 )

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     This function calculates bin withs for any discretized function.

    !     May be used for numerical integration and differentiation.

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !      Parameter list

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

    !      X     R.A.      I  Array type of REAL

    !      DX    R.A       O  Bin widths if X

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

    !

    !  4. Remarks :

    !

    !     This function does not get used in update by S. Cheng.

    !     It should be removed if/when "V1" routines are removed.

    !

    !  5. Called by :

    !     W3IC3WNCG_V1

    !

    !  6. Source code :

    !/

    IMPLICIT NONE

    !/

    REAL, INTENT(IN)  :: x(:)

    REAL, ALLOCATABLE :: dx(:)

    INTEGER           :: ix, nx

    !/

    nx = SIZE(x,1)

    ALLOCATE(dx(nx))

    dx = 0.

    !/

    DO ix = 1,nx

      IF (ix==1) THEN

        dx(ix) = (x(ix+1)-x(ix  ))

      ELSE IF (ix==nx) THEN

        dx(ix) = (x(ix  )-x(ix-1))

      ELSE

        dx(ix) = (x(ix+1)-x(ix-1)) / 2.

      END IF

    END DO

    !/

  END FUNCTION delta


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

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

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

  ! Start of new codes (or new variants) provided by S. Cheng

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

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

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


  SUBROUTINE w3ic3wncg_cheng(WN_R,WN_I,CG,ICE1,ICE2,ICE3,ICE4,DPT)

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |           X. Zhao                 |

    !/                  |           S. Cheng                |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Jan-2016 |

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

    !/

    !  1. Purpose :

    !

    !     Calculation of complex wavenumber for waves in ice. Outsourced

    !     from W3SIC3 to allow update on wavenumbers and  group

    !     velocities at each time step an ice parameter is updated.

    !

    !  2. Method :

    !

    !     <text here>

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !       WN_R   R. A.  I/O  Wave number (real part)

    !       WN_I   R. A.  I/O  Wave number (imag. part=wave attenuation)

    !       CG     R. A.  I/O  Group velocity

    !       ICE1   REAL    I   Thickness of ice                 [in m]

    !       ICE2   REAL    I   Effective viscosity of ice       [in m2/s]

    !       ICE3   REAL    I   Density of ice                  [in kg/m3]

    !       ICE4   REAL    I   Effective shear modulus of ice   [in Pa]

    !       DPT    REAL    I   Water depth                      [in m]

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

    !

    !  4. Subroutines used :

    !

    !      Name    Type  Module    Description

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

    !     WAVNU1   Subr. W3DISPMD   Wavenumber for waves in open water.

    !     WN_CMPLX Func. W3SIC3MD   Complex wavenumber for waves in ice.

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

    !

    !  5. Called by :

    !

    !      Name    Type  Module   Description

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

    !      W3SIC3  Subr. W3SIC3MD Ice source term.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !     On pre-calculation table, S. Cheng says:

    !     Instead of interpolation, WN_R of arbitrary HICE is approximated

    !     by IC3WN_R in the look up table. IC3WN_R is related to an ice

    !     thickness in the table, which is closest to HICE and less than

    !     HICE

    !

    !     Fix submitted by Sukun March 2017:

    !         replace :

    !            I = MIN(INT(HICE/IC3_DITK)+1,ITKNUM)

    !         with :

    !            I = MIN(NINT(HICE/IC3_DITK),ITKNUM)

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Source code :

    !/

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

    !/

    USE w3gdatmd, ONLY: nk,sig, ic3pars

    USE w3adatmd, ONLY: ic3wn_r, ic3wn_i, ic3cg

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

    USE w3dispmd, ONLY: wavnu1

    !/

    IMPLICIT NONE

    !/

    REAL, INTENT(IN)   :: ice1, ice2, ice3, ice4, dpt

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

    REAL, ALLOCATABLE  :: sigma(:)

    !

    INTEGER     :: i, i1, i2, ik, kl,ku, itknum

    COMPLEX(8)  :: wncomplex, x0,x1,x2, wnr, wnl

    REAL(8)     :: depth, hice, nu, dice, es_mod, rr, k_ocean, &

         cg_ocean

    REAL        :: ic3hilim,ic3kilim


    ic3hilim=ic3pars(10)

    ic3kilim=ic3pars(11)


    ALLOCATE(  sigma( SIZE(cg) ) )

    sigma = 0.

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

      kl    = 1

      ku    = nk

      i1    = 1

      i2    = nk

      sigma = sig(1:nk)

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

      kl    = 1

      ku    = nk+2

      i1    = 0

      i2    = nk+1

      sigma = sig(0:nk+1)

    ELSE

      WRITE(ndse,900)

      CALL extcde(3)

    END IF

    depth   =   dpt     ! water depth

    hice    =   ice1    ! ice thickness


    ! Optional: limit ice thickness

    hice=min(dble(ic3hilim),hice)

    nu      =   ice2    ! "effective viscosity" parameter

    dice    =   ice3    ! density of ice

    es_mod  =   ice4    ! effective shear modulus of ice

    !

    itknum  = ceiling(ic3_maxitk/ic3_ditk)

    i       = min(nint(hice/ic3_ditk),itknum)


    ! Find values in pre-calculated look-up table

    ! See Remarks section.

    wn_r    =   ic3wn_r(i1:i2, i)

    wn_i    =   ic3wn_i(i1:i2, i)

    cg      =   ic3cg(i1:i2, i)

    rr      =   0.01

    !/ --- CHECK If it's shallow water situation, then it needs recalculate

    !              kr,ki,cg

    DO ik   =   kl,ku

      IF (wn_r(ik)*depth>4.0)THEN ! exit do-loop

        EXIT  ! assume kr is proportional to frequency

      ELSE

        x1  = cmplx(wn_r(ik),wn_i(ik))

        x0  = x1*(1-rr)

        x2  = x1*(1+rr)

        wncomplex = cmplx_root_muller_cheng(x0,x1,x2,1, &

             dble(sigma(ik)),es_mod,nu,dice,hice,depth)

        wn_i(ik)  = real(aimag(wncomplex))  !  ki

        wn_r(ik)  = real(wncomplex)         !  kr

      ENDIF

    ENDDO


    CALL smooth_k(wn_r,wn_i,sigma,ku-kl+ 1,0)

    ! Optional : limit ki

    DO ik   =   kl,ku

      wn_i(ik)  = min(wn_i(ik),ic3kilim)

    ENDDO


!!!     --- Update group velocitiy ----

    CALL cginic3_cheng(cg,sigma,wn_r,ku-kl+1)


    DEALLOCATE(sigma)


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

         '     CANNOT DETERMINE BOUNDS OF WAVENUMBER ARRAY.')

  END SUBROUTINE w3ic3wncg_cheng

  !

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

  SUBROUTINE ic3table_cheng(ICE2,ICE3,ICE4)

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |           X. Zhao                 |

    !/                  |           S. Cheng                |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Jan-2016 |

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

    !/

    !  1. Purpose :

    !

    !     It's a Preprocess part to create a table of wave nubmer,

    !        attenuation and group velocity

    !     for all ice thickness in deep water situation for main

    !        computation.

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !  4. Subroutines used :

    !

    !      Name               Type  Module   Description

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

    !      CMPLX_ROOT_MULLER_CHENG  Func. W3SIC3MD Find root for complex

    !                                         wavenumbers for waves in ice.

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

    !

    !  5. Called by :

    !

    !      Name              Type   Module   Description

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

    !      W3WAVE             Subr. W3WAVEMD

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Cheng and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on Aug 25 2015

    !

    !   **UNRESOLVED BUG** This routine should be called again if ice

    !     rheology (visc., elast.) changes (either time or space).

    !     It doesn't. We need to set CALLEDIC3TABLE=0 if either parameter is

    !     changed.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

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

    USE w3gdatmd, ONLY: nk,sig

    USE w3adatmd, ONLY: ic3wn_r, ic3wn_i, ic3cg

    USE w3idatmd, ONLY: inflags2

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

    USE constants, ONLY: grav

    !

    IMPLICIT NONE

    REAL       :: ice1, ice2, ice3, ice4, dpt

    INTEGER ::  i1, i2, jitk, itknum, ik


    dpt    = 999.


    itknum = ceiling(ic3_maxitk/ic3_ditk)

    ic3wn_r(:,0) = sig**2/grav

    ic3wn_i(:,0) = 0

    DO jitk = 1,itknum

      ice1 = jitk*ic3_ditk  !HICE

      CALL ic3precalc_cheng(ic3wn_r(:,jitk), ic3wn_i(:,jitk), &

           ic3cg(:,jitk), ice1, ice2, ice3, ice4, dpt)

    ENDDO


    RETURN

  END SUBROUTINE ic3table_cheng


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

  SUBROUTINE ic3precalc_cheng(WN_R,WN_I,CG,ICE1,ICE2,ICE3,ICE4,DPT)

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |           X. Zhao                 |

    !/                  |           S. Cheng                |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Jan-2016 |

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

    !/

    !  1. Purpose :

    !

    !     Preprocess part to create a table of wave nubmer, attenuation and

    !     group velocity for all ice thickness in deep water situation for

    !     main computation.

    !

    !  2. Method :

    !

    !     Calculate them use Muller's method

    !

    !  3. Parameters :

    !

    !  4. Subroutines used :

    !

    !      Name               Type  Module   Description

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

    !      CMPLX_ROOT_MULLER_CHENG  Func. W3SIC3MD Find root for complex

    !                                       wavenumbers for waves in ice.

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

    !

    !  5. Called by :

    !

    !      Name              Type  Module   Description

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

    !      IC3TABLE_CHENG      Subr. W3SIC3MD  Create a table of kr,ki,cg

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Cheng and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on Aug 25 2015

    !

    !   Sukun Cheng in reference to MIN(WN_I(IK),--):

    !   "This artificial limitation reduces IC3 model’s effect,

    !       though it saves some time.

    !   ki > 2.e-4  is a common truth for angular frequency larger

    !   than the value around 2 or 3 depending on other parameters."

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

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

    USE w3gdatmd, ONLY: nk, sig

    USE w3dispmd, ONLY: wavnu1

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

    !/

    IMPLICIT NONE

    !/

    REAL, INTENT(INOUT):: wn_r(0:nk+1),wn_i(0:nk+1), cg(0:nk+1)

    REAL, INTENT(IN)   :: ice1, ice2, ice3, ice4, dpt


    INTEGER            :: ik, kl,ku,ix,num,switchid

    REAL               :: k_ocean,cg_ocean

    REAL(8)            :: kh, k_noice, depth, hice, nu, dice, &

         es_mod,sigmam1

    COMPLEX(8)         :: wncomplex, x0,x1,x2,wnm1,wnm2,wn_o

    !

    !     --- Input to routine

    depth   =   dpt     ! water depth

    hice    =   ice1    ! ice thickness

    nu      =   ice2    ! "effective viscosity" parameter

    dice    =   ice3    ! density of ice

    es_mod  =   ice4    ! effective shear modulus of ice


    kl    = 0

    ku    = nk+1

    switchid = 0


    DO ik = kl,ku

      CALL wavnu1(sig(ik),dpt,k_ocean,cg_ocean)

      k_noice = k_ocean


      !     --- Calculate complex wavenumber ------------------------------- /

      IF(ik<kl+2)THEN

        wnm1     = cmplx(0.,0.)

        wnm2     = cmplx(0.,0.)

        sigmam1  = 0.0

      ELSE

        wnm1     = cmplx(wn_r(ik-1),wn_i(ik-1))

        wnm2     = cmplx(wn_r(ik-2),wn_i(ik-2))

        sigmam1  = sig(ik-1)

      ENDIF


      wn_o = cmplx(k_noice,0.0)

      CALL wn_precalc_cheng(wncomplex,dble(sig(ik)),dble(sigmam1),wn_o, &

           wnm1,wnm2,es_mod,nu,dice,hice,depth,switchid,ik,ku)

      !     --- Output to routine

      wn_r(ik)  = real(wncomplex)         !  kr

      wn_i(ik)  = imag(wncomplex)         !  ki

    ENDDO

    CALL smooth_k(wn_r,wn_i,sig,ku-kl+1,switchid)

    CALL cginic3_cheng(cg,sig,wn_r,ku-kl+1)

    !/

  END SUBROUTINE ic3precalc_cheng


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


  SUBROUTINE wn_precalc_cheng(WN,SIGMA,SIGMAM1,WN_O,WNM1,WNM2,ES_MOD, &

       NU,DICE,HICE,DEPTH,SWITCHID,IK,KU)

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           S. Cheng                |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |           X. Zhao                 |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Jan-2016 |

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

    !/

    !  1. Purpose :

    !

    !     Calculate complex wavenumber for waves in ice.

    !

    !  2. Method :

    !

    !     Wang and Shen (JGR 2010)

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !       WN     CMPLX  O  Wave number (imag. part=wave attenuation)

    !       SIGMA  REAL   I  Wave angular frequency           [in rad]

    !       WN_O   REAL   I  Wave number (open water)

    !       ES     REAL   I  Effective shear modulus of ice   [in Pa]

    !       NU     REAL   I  Effective viscosity of ice       [in m2/s]

    !       DICE   REAL   I  Density of ice                  [in kg/m3]

    !       HICE   REAL   I  Thickness of ice                 [in m]

    !       DEPTH  REAL   I  Water depth                      [in m]

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

    !

    !  4. Subroutines used :

    !

    !      Name               Type  Module   Description

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

    !      CMPLX_ROOT_MULLER_CHENG  Func. W3SIC3MD Find root for complex

    !                                         wavenumbers for waves in ice.

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

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

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

    !      IC3PRECALC_CHENG    xxx  xxx  xxxx

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Updated authors: Cheng and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on Aug 25 2015

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !   Hayley Shen says,

    !     We have determined that it may not be necessary to use curve

    !     fitting or lookup tables to get the group velocity and the

    !     attenuation coefficient. Attached is a short report with some

    !     sample numerical solutions. To implement the viscoelastic model,

    !     there are 4 fortran programs. According to Xin Zhao, the graduate

    !     student, it is very fast to find roots. I suggest that perhaps you

    !     try the pure viscous case by setting G=0 to start with. nu can be

    !     set at 0.05*ice concentration (m^2/s) to begin with, because for

    !     grease ice Newyear’s data showed nu to be about 0.02-0.03 m^2/s.

    !     By setting G=0 in you get exactly the Keller model for pure

    !     viscous layer.

    !

    !   This routine provides the initial guess according to the parameters

    !   of the present case. T>10s use open water, T<10s cases, calculate

    !   T=10s first using open water as the initial guess.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

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

    USE constants, ONLY: tpi

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    REAL(8)       :: SIGMA,SIGMAM1,ES_MOD,NU,DICE,HICE,DEPTH

    COMPLEX(8)    :: WN_O, WNM1, WNM2, WN0, WN1,WN2

    COMPLEX(8),INTENT(OUT) :: WN ! RESULT

    !/

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

    !/ Local parameters

    !/

    INTEGER       :: I,IX,NUM,SWITCHID,IK,KU

    REAL(8)       :: RR, R2, EPS, SIGMA0, KAPPA, DIS

    COMPLEX(8)    :: X0, X1, X2, kp, ks, Gv

    !/

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

    !   compute shear wave and pressure wave modes

    gv = cmplx(es_mod,-sigma*nu*dice)

    kp = sigma/sqrt(4*gv/dice)

    ks = 2*kp

    !   RR,R2 are empirical coefficients

    rr  = 0.2

    r2  = 4

    eps = 1.d-10   ! assuming variable =0, if it is less than this value

    !   compute root 1

    !   initial guesses from wave number of open water

    !

    x1  = wn_o    ! initial guess

    x0  = x1*(1-rr)

    x2  = x1*(1+rr)


    wn1 = cmplx_root_muller_cheng(x0,x1,x2,1,sigma, &

         es_mod,nu,dice,hice,depth)

    !   if root finder failed, or found shear wave mode,

    !   redo searching using simplified dispersion relation form, index 2

    IF (real(wn1)<0.or.abs(wn1-ks)/abs(ks)<0.03 .or. &

         abs(wn1-kp)/abs(kp)<0.03.and.es_mod>1.e7*nu)THEN

      wn1 = cmplx_root_muller_cheng(x0,x1,x2,2,sigma, &

           es_mod,nu,dice,hice,depth)

    ENDIF

    !   similar as wn1, but search with opposite order of inital guesses

    wn2 = cmplx_root_muller_cheng(x2,x1,x0,1,sigma, &

         es_mod,nu,dice,hice,depth)

    IF (real(wn2)<0.or.abs(wn2-ks)/abs(ks)<0.03)THEN

      wn2 = cmplx_root_muller_cheng(x2,x1,x0,2,sigma, &

           es_mod,nu,dice,hice,depth)

    ENDIF

    IF(abs(real(wn1)-real(wn_o))<abs(real(wn2)-real(wn_o)))THEN

      wn0 = wn1

    ELSE

      wn0 = wn2

    ENDIF


    IF(sigmam1==0.)THEN

      wn = wn0

    ELSE

      !      compute root 2

      !      Calculate the other wave number based on last frequency

      !      in the frequency array

      r2    = max(abs(wnm1-wnm2)/abs(wnm2), abs((sigma-sigmam1)/sigma))

      DO i  = 1,7,2

        IF(i<3.AND.sigma>4)THEN

          cycle

        ENDIF

        num = 2**(i-1)

        rr  = min(max(r2/real(num,8),0.001d0),0.5d0)

        x1  = wnm1

        DO ix  = 1,num

          x0 = x1*(1-rr)

          x2 = x1*(1+rr)

          sigma0 = sigmam1 + (1.0*ix)/num*(sigma-sigmam1)

          wn = cmplx_root_muller_cheng(x0,x1,x2,1,sigma0, &

               es_mod,nu,dice,hice,depth)

          IF(real(wn)<0)THEN   ! try another searching direction

            wn = cmplx_root_muller_cheng(x2,x1,x0,1,sigma0, &

                 es_mod,nu,dice,hice,depth)

          ENDIF

          IF(real(wn)<0)THEN   ! try another dispersion relation form

            wn = cmplx_root_muller_cheng(x0,x1,x2,3,sigma0, &

                 es_mod,nu,dice,hice,depth)

          ENDIF

          kp = sigma0/sqrt(4*gv/dice)

          ks = 2*kp

          !               set 3 means simple dispersion relation form

          IF(abs(wn-ks)/abs(ks)<0.03.OR.real(wn)<0)THEN

            wn = cmplx_root_muller_cheng(x0,x1,x2,2,sigma0, &

                 es_mod,nu,dice,hice,depth)

          ENDIF

          x1 = wn

          IF( real(wn)<0.99*real(wnm1).OR. abs(wn-wn0)>eps .AND. &

               (abs(x1-wn)/abs(wn)>0.3 .OR.                  &

               imag(wn)/(imag(x1)+eps)>10.OR. real(wn)<0))THEN

            EXIT  ! redo with smaller intervals

          ENDIF

          x0 = x1

          x1 = wn

        ENDDO !            DO IX  = 1,NUM


        IF(ix==num+1)THEN

          EXIT

        ENDIF

        wn = x1  !/ --- if exit of inner loop: give an approximate wn

      ENDDO !         DO I  = 1,7,2


      ! part 2

      ! assume found two roots, choose from the two condidates.

      kp = sigma/sqrt(4*gv/dice)

      IF(real(wn0)>0.AND.imag(wn0)>=0.AND.abs(wn - wn0)>eps)THEN

        !do switch at last 3 points is not worth numerically.

        !For v>5.e-2, it is low chance to be wrong at the last 3 points

        ! we suppose to use two mode in the future

        IF(nu>0.AND.ik>=ku-1.OR.   &

             nu>0.AND.switchid/=0)THEN

          ! assume one mode switch for general viscoelastic model

        ELSE

          dis   = abs(real(wn)-real(wn_o))/abs(real(wn0)-real(wn_o))

          kappa = (imag(wn)+ eps)/(imag(wn0) + eps)

          ! wn0 has smaller attenuation and closer to k0

          IF ((dis >= 1 .AND. kappa>=1 .AND.  &

               imag(wn0)>=0.1*imag(wnm1).AND. &

               abs(wn-kp)<abs(wn0-kp)) .OR.  &

               ! wn0 has larger attenuation and closer to k0

               ( dis>=1 .AND. kappa<1 .AND. &

               ((kappa> 0.2 .AND. imag(wn0)/real(wn0)<0.5).or. &

               abs(real(wn)-real(kp))<abs(real(wn0)-real(kp)))).OR.&

               ! wn0 has smaller attenuation and farther to k0,

               ! wn0 could be wrong at high G

               ! (kappa>1 .and. dis<1 .and. dis>  0.8 ))then

               ( kappa>1 .AND. dis<1 .AND. &

               abs(real(wn)-real(wnm1))>   &

               abs(real(wn0)-real(wnm1)) ))THEN

            wn = wn0

            switchid = ik

            ! wn0 has lager attenuation and farther to k0

          ELSEIF(dis<1 .AND. kappa<1) THEN

            ! keep wn without change

          ENDIF !  IF ((DIS >= 1 .AND. KAPPA>=1 .AND.  &

        ENDIF !   IF(NU>0.AND.IK>=KU-2.OR.   &


        ! choose dominant mode is farther than pressure wave.

        !but it doens't work for high viscosity.

        if (abs(wn-kp)/abs(kp)<0.03) then

          wn = wn0

          switchid = ik

        endif

        !if (real(wn)<=real(wnm1))then

        !    wn = wn0

        !    switchid = ik

        !endif


        IF (real(wn0)>real(wnm1).AND.real(wn0)<real(wn).and. &

             abs(imag(wn0)-imag(wnm1))<abs(imag(wn)-imag(wnm1)))THEN

          ! mainly for pure elastic case has multiple branchs

          wn = wn0

          switchid = ik

        ENDIF

      ENDIF ! IF(REAL(WN0)>0.AND.IMAG(WN0)>=0.AND.ABS(WN - WN0)....

    ENDIF   !      IF(SIGMAM1==0.)THEN


    IF(real(wn)<0)THEN

      print*, "MULLER METHOD FAILED, ES_MOD,NU,HICE:",es_mod,nu,hice

    ENDIF

    RETURN


  END SUBROUTINE wn_precalc_cheng

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

  !/

  SUBROUTINE cginic3_cheng(CG,SIGMA,WN_R,N)

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |           X. Zhao                 |

    !/                  |           S. Cheng                |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Jan-2016 |

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

    !  1. Purpose :

    !

    !     Calculate group velocity in ic3 model

    !

    !  2. Method :

    !

    !     finite differece

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !       CG      R. A.  I/O  Group velocity

    !       SIMGA   R. A.  angular frequency

    !       WN_R    R. A.  wave number

    !       N       Int.   array size

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

    !

    !  4. Subroutines used :

    !

    !      None.

    !

    !  5. Called by :

    !

    !      Name             Type  Module   Description

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

    !      IC3PRECALC_CHENG    Subr. W3SIC3MD Create table of kr,ki,cg for

    !                 deep water

    !      W3IC3WNCG_CHENG     Subr. W3SIC3MD Calculate kr,ki,cg

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !      Smooth function is used due to jump problem when wave modes

    !             switch

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Source code :

    !/

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

    !/

    REAL, INTENT(IN)  :: SIGMA(1:N),WN_R(1:N)

    REAL, INTENT(OUT) :: CG(1:N)

    INTEGER           :: N

    !/ LOCAL variables

    INTEGER           :: IK, M

    REAL              :: CG1,CG2,CG3,CG0(1:N)

    !/

    ik = 1

    cg(ik)=(sigma(ik+1)-sigma(ik))/(wn_r(ik+1)-wn_r(ik))

    ik = n

    cg(ik)=(sigma(ik)-sigma(ik-1))/(wn_r(ik)-wn_r(ik-1))

    DO ik = 2,n-1

      cg1 = (sigma(ik)-sigma(ik-1))/(wn_r(ik)-wn_r(ik-1))

      cg2 = (sigma(ik+1)-sigma(ik))/(wn_r(ik+1)-wn_r(ik))

      cg(ik) = 2.0/(1./cg1 + 1./cg2)

    END DO


    RETURN

  END SUBROUTINE cginic3_cheng


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

  ! numerically smooth WN_R and WN_I by linear interpolation

  SUBROUTINE smooth_k(WN_R,WN_I,SIGMA,N,SWITCHID)

    REAL, INTENT(IN)  :: SIGMA(N)

    REAL              :: WN_R(N), WN_I(N),DIFF(N),REMOVEID(N)

    INTEGER           :: N,I,J,SWITCHID

    !

    diff = 0

    ! remove kr in mode switch zone,

    ! if it is a local extremum or suddenly increasing

    DO j = 1,3 ! 3 times to guarantee wavenumber increases monotonically

      removeid = 0

      DO i = 2,n

        diff(i) = wn_r(i) - wn_r(i-1)

      ENDDO

      DO i = 3,n

        IF(diff(i)<=0.OR.diff(i)>3*diff(i-1).OR.switchid==i)THEN

          removeid(i) = 1

          removeid(i-1) = 1

        ENDIF

      ENDDO

      ! fill removed location with kr,ki by interpolation

      DO i = 2,n-1

        IF (removeid(i) ==1) THEN

          wn_r(i) = wn_r(i-1) + (wn_r(i+1)-wn_r(i-1))/    &

               (sigma(i+1)-sigma(i-1))*(sigma(i)-sigma(i-1))

          wn_i(i) = wn_i(i-1) + (wn_i(i+1)-wn_i(i-1))/    &

               (sigma(i+1)-sigma(i-1))*(sigma(i)-sigma(i-1))

        ENDIF

      ENDDO

    ENDDO

    ! mode switch upward at the last frequencies

    IF (diff(n)>3*diff(n-1))THEN

      i = n

      wn_r(i) = wn_r(i-1) + (wn_r(i-1)-wn_r(i-2))/   &

           (sigma(i-1)-sigma(i-2))*(sigma(i)-sigma(i-1))

      wn_i(i) = wn_i(i-1) + (wn_i(i-1)-wn_i(i-2))/   &

           (sigma(i-1)-sigma(i-2))*(sigma(i)-sigma(i-1))

    ENDIF

    RETURN

  END SUBROUTINE smooth_k


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


  FUNCTION cmplx_root_muller_cheng(X0, X1, X2, JUDGE,  &

       SIGMA,ES,NU,DICE,HICE,DEPTH )           RESULT(P3)

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |           X. Zhao                 |

    !/                  |           S. Cheng                |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Jan-2016 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Find root.

    !

    !  2. Method :

    !

    !     Muller method for complex equations is a recursive approximation

    !     with initial guess X0, X1, and X2. To the initial guesses a

    !     quadratic parabola is fitted.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !       P3    CMPLX DBL O  Approximation for the root problem

    !       X0    CMPLX DBL I  Initial guess variable

    !       X1    CMPLX DBL I  Initial guess variable

    !       X2    CMPLX DBL I  Initial guess variable

    !       SIGMA DOUBLE    I  Wave angular frequency

    !       ES    DOUBLE    I  Effective shear modulus of ice

    !       NU    DOUBLE    I  Effective viscosity of ice       [in m2/s]

    !       DICE  DOUBLE    I  Density of ice                   [in kg/m3]

    !       HICE  DOUBLE    I  Thickness of ice                 [in m]

    !       DEPTH DOUBLE    I  Water depth                      [in m]

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

    !

    !  4. Subroutines used :

    !

    !      Name               Type  Module   Description

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

    !      F_ZHAO_CHENG       xxxxx xxxx     xxxx

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

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

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

    !     WN_PRECALC_CHENG xxxxx    xxxx

    !     W3IC3WNCG_CHENG  xxxxx    xxxx

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

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

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    COMPLEX(8)            :: P3                     ! RESULT

    COMPLEX(8), INTENT(IN):: X0,X1,X2

    REAL(8), INTENT(IN)   :: SIGMA,ES,NU,DICE,HICE,DEPTH

    INTEGER, INTENT(IN)   :: JUDGE

    !/

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

    !/ Local parameters

    !/


    INTEGER               :: I

    INTEGER, PARAMETER    :: IMAX = 200

    REAL(8)               :: DLTA,EPSI

    COMPLEX(8)            :: P0,P1,P2

    COMPLEX(8)            :: Y0,Y1,Y2,Y3

    COMPLEX(8)            :: A,B,C,Q,DISC,DEN1,DEN2

    !/

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


    p0      =   x0

    p1      =   x1

    p2      =   x2

    p3      = 0.d0


    epsi    = 1.d-8

    dlta    = 1.d-8   ! relax it may cause error, too rigorous is not good neither.

    y0      = f_zhao_cheng(judge,p0,sigma,es,nu,dice,hice,depth)

    y1      = f_zhao_cheng(judge,p1,sigma,es,nu,dice,hice,depth)

    y2      = f_zhao_cheng(judge,p2,sigma,es,nu,dice,hice,depth)


    DO i = 1,imax

      q = (p2 - p1) / (p1 - p0)

      a = q * y2 - q * (1.d0+q) * y1 + q**2.d0 * y0

      b = (2.d0 * q + 1.d0) * y2 - (1.d0 + q)**2.d0 * y1 &

           + q**2.d0 * y0

      c = (1.d0 + q) * y2


      IF ( abs(a).NE.0.d0 ) THEN


        disc = b**2.d0 - 4.d0 * a * c;


        den1 = ( b + sqrt( disc ) )

        den2 = ( b - sqrt( disc ) )


        IF ( abs( den1 ) .LT. abs( den2 ) )THEN

          p3 = p2 - (p2 - p1) * (2.d0 * c / den2)

        ELSE

          p3 = p2 - (p2 - p1) * (2.d0 * c / den1)

        ENDIF


      ELSE


        IF ( abs(b) .NE. 0.d0 )THEN

          p3 = p2 - (p2 - p1) * (c / b)

        ELSE

          p3 = p2

          RETURN

        ENDIF

      ENDIF


      IF (imag(p3).LT.0)THEN

        p3 = real(p3) - cmplx(0.,1.)*imag(p3)

      ENDIF

      IF (real(p3).LT.0)THEN

        p3 = -real(p3) + cmplx(0.,1.)*imag(p3)

      ENDIF


      IF(nu==0)THEN

        p3 = cmplx(real(p3),0)

      ENDIF

      y3 = f_zhao_cheng(judge,p3,sigma,es,nu,dice,hice,depth);

      IF ( abs(p3-p2).LT.dlta .AND. abs(y3).LT.epsi ) THEN

        ! exit before finding a true root,Result may not be accurate

        RETURN

      ENDIF


      p0 = p1

      p1 = p2

      p2 = p3


      y0 = y1

      y1 = y2

      y2 = y3

    ENDDO


    p3 = cmplx(-100.,0)

    RETURN


  END FUNCTION cmplx_root_muller_cheng

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

  !/

  FUNCTION f_zhao_cheng(JUDGE,X,SIGMA,ES,NU,DICE,HICE,DEPTH) &

       result(fzhao)

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |           X. Zhao                 |

    !/                  |           S. Cheng                |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Jan-2016 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Decide whether to call sub-function.

    !

    !  2. Method :

    !

    !     Decide based on value of integer "JUDGE"

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !      FZHAO COMPL8  O    Result (double complex)

    !      X     CMPLX8  I    Approximate result (double complex)

    !      JUDGE INTEGR  I    Switch variable

    !      SIGMA DOUBLE  I    Wave angular frequency

    !      ES    DOUBLE  I    Effective shear modulus

    !      NU    DOUBLE  I    Effective viscosity parameter

    !      DICE  DOUBLE  I    Density of ice

    !      HICE  DOUBLE  I    Thickness of ice

    !      DEPTH DOUBLE  I    Water depth

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

    !

    !  4. Subroutines used :

    !

    !      Name        Type  Module   Description

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

    !      DRFUN_dble  Func. W3SIC3MD Function to find root with double

    !      precision.

    !      DRFUN_quad  Func. W3SIC3MD Function to find root with quadruple

    !      precision.

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

    !

    !  5. Called by :

    !

    !      Name               Type  Module   Description

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

    !      CMPLX_ROOT_MULLER_CHENG  Func. W3SIC3MD Find root for complex

    !                                         wavenumbers for waves in ice.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Updated authors: Cheng and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on Aug 25 2015

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !/

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

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    INTEGER,INTENT(IN)           :: JUDGE

    DOUBLE PRECISION, INTENT(IN) :: SIGMA,ES,NU,DICE,HICE,DEPTH

    DOUBLE COMPLEX, INTENT(IN)   :: X

    DOUBLE COMPLEX               :: FZHAO  ! RESULT

    !/

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

    IF (judge >0) THEN

      fzhao = drfun_dble_cheng(x,sigma,es,nu,dice,hice,depth,judge)

    ELSEIF(judge ==0)THEN

      fzhao = drfun_quad_cheng(x,sigma,es,nu,dice,hice,depth)

    ENDIF

    !

  END FUNCTION f_zhao_cheng

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

  !/

  FUNCTION drfun_dble_cheng(WN,SIGMA,ES,NU,DICE,HICE,DEPTH,JUDGE) &

       result(func1)

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |           X. Zhao                 |

    !/                  |           S. Cheng                |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Jan-2016 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Return dispersion relation function value for root finding

    !

    !  2. Method :

    !

    !     function based on dispersion relation derived by Wang and Shen

    !     2010

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !      FUNC1 CMPLX DBL O  Result (COMPLEX(8))

    !      WN    CMPLX DBL I  Wavenumber (COMPLEX(8))

    !      W     REAL  DBL I  Wave angular frequency

    !      ES    REAL  DBL I  Effective shear modulus on ice

    !      NU    REAL  DBL I  Effective viscosity

    !      DICE  REAL  DBL I  Density of ice

    !      HICE  REAL  DBL I  Thickness of ice

    !      DEPTH REAL  DBL I  Water depth

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

    !

    !  4. Subroutines used :

    !

    !      Name    Type  Module    Description

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

    !       None.

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

    !

    !  5. Called by :

    !

    !      Name    Type  Module   Description

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

    !      F_ZHAO_CHENG   xxx    xxxx      xxxxx

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Updated authors: Cheng and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on Aug 25 2015

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !/

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

    USE constants, ONLY: grav, dwat

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    COMPLEX(8), INTENT(IN)   :: wn

    REAL(8), INTENT(IN)      :: sigma, es, nu, dice, hice, depth

    COMPLEX(8)               :: func1,aa(4,4)                   ! RESULT

    INTEGER                  :: judge

    !/

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

    !/ Local parameters

    !/

    COMPLEX(8)   :: ve,alpha,n,m,l,sk,ck,sa,ca,th,thh,temp,j1,j2

    !/

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

    ve    = cmplx( nu, es/dice/sigma )

    alpha = sqrt( wn**2. - sigma/ve * cmplx(0.,1.d0) )

    n     = sigma + 2. * ve * wn**2. * cmplx(0.,1.d0)


    temp  = exp(wn*hice)

    sk    = (temp - 1.d0/temp)/2.d0

    ck    = (temp + 1.d0/temp)/2.d0


    temp  = exp(alpha*hice)

    sa    = (temp - 1.d0/temp)/2.d0

    ca    = (temp + 1.d0/temp)/2.d0

    !

    temp  = (wn*depth)

    IF ( real(temp).LT.18.d0 ) THEN

      temp = exp(temp)

      th = (temp - 1./temp)/(temp + 1./temp)

    ELSE

      th = 1.d0

    ENDIF

    !      JUDGE==3 is not used yet

    IF ((es>=1.e5.AND.judge/=2).or.judge==3) THEN

      l     = 2 * wn * alpha * sigma * ve

      m     = (dble(dwat)/dice - 1) * dble(grav) * wn             &

           - dble(dwat) / dice * sigma**2 / th

      aa(1,1) = 0.

      aa(1,2) = 2 * cmplx(0.,1.) * wn**2.

      aa(1,3) = alpha**2. + wn**2.

      aa(1,4) = 0.

      !

      aa(2,1) = n * sigma

      aa(2,2) = -wn * dble(grav)

      aa(2,3) = cmplx(0.,1.) * wn * dble(grav)

      aa(2,4) = l

      !

      aa(3,1) = -2. * cmplx(0.,1.) * wn**2. * sk

      aa(3,2) =  2. * cmplx(0.,1.) * wn**2. * ck

      aa(3,3) =  (alpha**2. + wn**2.) * ca

      aa(3,4) = -(alpha**2. + wn**2.) * sa

      !

      aa(4,1) =   n * sigma * ck - m * sk

      aa(4,2) = - n * sigma * sk + m * ck

      aa(4,3) = -cmplx(0.,1.) * m * ca - l * sa

      aa(4,4) =  cmplx(0.,1.) * m * sa + l * ca

      !

      func1   = bsdet(aa,4)

    ELSE

      j1    = dice/dble(dwat)*(wn**2.*dble(grav)**2.*sk*sa - (n**4. &

           + 16.* ve**4.*wn**6.*alpha**2.)*sk*sa - 8. &

           *wn**3.*alpha*ve**2.*n**2.*(ck*ca-1.))

      j2    = (4.*wn**3.*alpha*ve**2.*sk*ca+n**2.*sa*ck &

           -dble(grav)*wn*sk*sa)

      IF (judge==2)THEN

        func1 = (sigma**2. - th*wn*dble(grav)) - th*j1/(j2+1.e-20)

      ELSEIF (judge==1)THEN

        func1 = (sigma**2. - th*wn*dble(grav))*j2 - th*j1

      ENDIF

    ENDIF

    !/

  END FUNCTION drfun_dble_cheng

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


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

  !/

  FUNCTION drfun_quad_cheng(WN,SIGMA,ES,NU,DICE,HICE,DEPTH) &

       result(func1)

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           E. Rogers               |

    !/                  |           S. Zieger               |

    !/                  |           X. Zhao                 |

    !/                  |           S. Cheng                |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         13-Jan-2016 |

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

    !/

    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )

    !/                                                        (E. Rogers)

    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)

    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Return dispersion relation function value for root finding

    !

    !  2. Method :

    !

    !     Use quadruple precision for computation

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !      FUNC1 CMPLX DBL O  Result (COMPLEX(8))

    !      WN    CMPLX DBL I  Wavenumber (COMPLEX(8))

    !      W     REAL  DBL I  Wave angular frequency

    !      ES    REAL  DBL I  Effective shear modulus on ice

    !      NU    REAL  DBL I  Effective viscosity

    !      DICE  REAL  DBL I  Density of ice

    !      HICE  REAL  DBL I  Thickness of ice

    !      DEPTH REAL  DBL I  Water depth

    !

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

    !

    !  4. Subroutines used :

    !

    !      Name    Type  Module    Description

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

    !      None.

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

    !

    !  5. Called by :

    !

    !      Name    Type  Module   Description

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

    !      F_ZHAO_CHENG   xxx    xxxx      xxxxx

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !   Updated authors: Cheng and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on Aug 25 2015

    !   Original authors: Zhao and Shen.

    !   This code is based on Fortran code provided by Hayley Shen (Clarkson

    !     University) to Erick Rogers (NRL) on April 19 2013.

    !   ER: S. Cheng had "COMPLEX(16)" for the local parameters. This is not

    !     supported by my compiler (gfortran on linux machine), so I changed

    !     it to "COMPLEX(8)"

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !/

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

    USE constants, ONLY: grav, dwat

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    COMPLEX(8), INTENT(IN)   :: wn

    REAL(8), INTENT(IN)      :: sigma, es, nu, dice, hice, depth

    COMPLEX(8)               :: func1                   ! RESULT

    !/

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

    !/ Local parameters

    !/

    COMPLEX(8)   :: ve,alpha,n,m,l,sk,ck,sa,ca,th,thh,temp,j1,j2

    !/

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

    ve    = cmplx( nu, es/dice/sigma )

    alpha = sqrt( wn**2. - sigma/ve * cmplx(0.,1.d0) )

    n     = sigma + 2. * ve * wn**2. * cmplx(0.,1.d0)


    temp  = exp(wn*hice)

    sk    = (temp - 1.d0/temp)/2.d0

    ck    = (temp + 1.d0/temp)/2.d0


    temp  = exp(alpha*hice)

    sa    = (temp - 1.d0/temp)/2.d0

    ca    = (temp + 1.d0/temp)/2.d0

    !

    temp  = (wn*depth)

    IF ( real(temp).LT.18.d0 ) THEN

      temp = exp(temp)

      th = (temp - 1./temp)/(temp + 1./temp)

    ELSE

      th = 1.d0

    ENDIF

    !

    j1    = dice/dble(dwat)*(wn**2.*dble(grav)**2.*sk*sa - (n**4. &

         + 16.* ve**4.*wn**6.*alpha**2.)*sk*sa - 8. &

         *wn**3.*alpha*ve**2.*n**2.*(ck*ca-1.))

    j2    = (4.*wn**3.*alpha*ve**2.*sk*ca+n**2.*sa*ck &

         -dble(grav)*wn*sk*sa)

    func1 = (sigma**2. - th*wn*dble(grav))*j2 - th*j1

    !/

  END FUNCTION drfun_quad_cheng

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

  ! End of new codes (or new variants) provided by S. Cheng

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

  !/

  !/ End of module W3SIC3MD -------------------------------------------- /

  !/

END MODULE w3sic3md