WW3/w3sic4md_8F90_source.html

#include "w3macros.h"

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

MODULE w3sic4md

  !/

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

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |           C. Collins              |

  !/                  |           E. Rogers               |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         21-Jan-2015 |

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

  !/

  !/    For updates see W3SIC4 documentation.

  !/

  !  1. Purpose :

  !

  !     Calculate ice source term S_{ice} according to simple methods.

  !          Attenuation is a function of frequency and specified directly

  !          by the user. Example: a function is based on an exponential fit to

  !          the empirical data of Wadhams et al. (1988).

  !

  !  2. Variables and types :

  !

  !  3. Subroutines and functions :

  !

  !      Name      Type  Scope    Description

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

  !      W3SIC4    Subr. Public   ice source term.

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

  !

  !  4. Subroutines and functions used :

  !

  !     See subroutine documentation.

  !

  !  5. Remarks :

  !

  !     Documentation of IC4:

  !         1) Collins and Rogers, NRL Memorandum report 2017

  !         ---> "A Source Term for Wave Attenuation by Sea

  !               Ice in WAVEWATCH III® : IC4"

  !         ---> describes original IC4 methods, 1 to 6

  !         2) Rogers et al., NRL Memorandum report 2018a

  !         ---> "Forecasting and hindcasting waves in and near the

  !              marginal ice zone: wave modeling and the ONR “Sea

  !              State” Field Experiment"

  !         ---> IC4 method 7 added

  !         2) Rogers et al., NRL Memorandum report 2018b

  !         ---> "Frequency Distribution of Dissipation of Energy of

  !               Ocean Waves by Sea Ice Using Data from Wave Array 3 of

  !               the ONR “Sea State” Field Experiment"

  !         ---> New recommendations for IC4 Method 2 (polynomial fit)

  !              and IC4 Method 6 (step function via namelist)

  !

  !     Other source material :

  !        *** Wadhams et al. JGR 1988

  !        *** Meylan et al. GRL 2014

  !        *** Kohout & Meylan JGR 2008 in Horvat & Tziperman Cryo. 2015

  !        *** Kohout et al. Nature 2014

  !        *** Doble et al. GRL 2015

  !        *** Rogers et al. JGR 2016

  !        *** Meylan et al. JGR 2018

  !        *** Yu et al. JGR 2019

  !        *** Liu et al. JPO 2020

  !        *** Rogers et al. CRST 2021 (RMK2021)

  !        *** Rogers et al. tech. rep. 2021 (RYW2021)

  !        *** Yu et al. CRST 2022

  !        *** Yu JMSE 2022

  !

  !  6. Switches :

  !

  !     See subroutine documentation.

  !

  !  7. Source code :

  !/

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

  !/

  PUBLIC :: w3sic4

  !/

CONTAINS

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

  SUBROUTINE w3sic4 (A, DEPTH, CG, IX, IY, S, D)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           C. Collins              |

    !/                  |           E. Rogers               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         24-Feb-2017 |

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

    !/

    !/    03-Dec-2015 : Origination                         ( version 5.09 )

    !/                     (starting from IC1)                (C. Collins)

    !/    03-Dec-2015 : W3SIC4 created, Methods 1,2,3,4       (C. Collins)

    !/    21-Jan-2016 : IC4 added to NCEP repository          (E. Rogers)

    !/    27-Jan-2016 : Method 5 added (step function)        (E. Rogers)

    !/    08-Apr-2016 : Method 6 added (namelist step funct.) (E. Rogers)

    !/    24-Feb-2017 : Corrections to Methods 1,2,3,4        (E. Rogers)

    !/    13-Apr-2017 : Method 7 added (Doble et al. 2015)    (E. Rogers)

    !/    11-Jan-2024 : Method 8 added (Meylan et al. 2018)   (E. Rogers)

    !/    11-Jan-2024 : Method 9 added (Rogers et al., 2021)

    !/                                      denoted "RYW2021" (E. Rogers)

    !/

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

    !

    !     S_{ice} source term using 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.

    !

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

    !

    !  2. Method :

    !

    !     Apply parametric/empirical functions, e.g. from the literature.

    !     1) Exponential fit to Wadhams et al. 1988, Table 2

    !     2) Polynomial fit, Eq. 3 from Meylan et al. 2014

    !     3) Quadratic fit to Kohout & Meylan'08 in Horvat & Tziperman'15

    !        Here, note that their eqn is given as ln(alpha)=blah, so we

    !        have alpha=exp(blah).

    !        Note from ER:

    !        This implementation has two things to keep in mind:

    !          1) This is a scattering model, applied as dissipation,

    !             which is not correct.

    !          2) This is not actually HT15! The alpha of HT15 has

    !             different meaning from alpha of CR17, as follows:

    !             HT15: decay is exp(-alpha*Lambda) where Lambda

    !                   is the number of floes encountered.

    !             CR17: decay is exp(-alpha*x)

    !             Thus, CR17's implementation of HT15 is equivalent to

    !             the actual HT15 only if one assumes one floe encountered

    !             per meter. This is very strong attenuation, as shown in

    !             Figure 3 of CR17! This problem might be fixed by computing

    !             an encounter interval length scale from an a_ice and d_ice

    !             provided by the user...or a length scale provided by the

    !             user.

    !             See also: page 3 of Rogers et al. (RYW2021).

    !     4) Eq. 1 from Kohout et al. 2014

    !

    !     5) Simple step function for ki as a function of frequency

    !          with up to 4 "steps". Controlling parameters KIx and FCx are

    !          read in as input fields, so they may be nonstationary and

    !          non-uniform in the same manner that ice concentration and

    !          water levels may be nonstationary and non-uniform.

    !                                          444444444444

    !                               33333333333

    !                   222222222222

    !      1111111111111

    !                  ^            ^          ^

    !                  |            |          |

    !                  5            6          7

    !      Here, 1 indicates ki=KI1=ICECOEF1 (ICEP1)

    !            2 indicates ki=KI2=ICECOEF2 (ICEP2)

    !            3 indicates ki=KI3=ICECOEF3 (ICEP3)

    !            4 indicates ki=KI4=ICECOEF4 (ICEP4)

    !            5 indicates freq cutoff #1 =FC5=ICECOEF5 (ICEP5)

    !            6 indicates freq cutoff #2 =FC6=ICECOEF6 (MUDD)

    !            7 indicates freq cutoff #3 =FC7=ICECOEF7 (MUDT)

    !     freq cutoff is given in Hz, freq=1/T (not sigma)

    !     Examples using hindcast, inversion with uniform ki:

    !        5.1) Beaufort Sea, AWAC mooring, 2012, Aug 17 to 20

    !             0.0418 Hz to 0.15 Hz : ki=10e-6

    !             0.15 Hz to 0.175 Hz : ki=11e-6

    !             0.175 Hz to 0.25 Hz : ki=15e-6

    !             0.25 Hz to 0.5 Hz : ki=25e-6

    !        5.2) Beaufort Sea, AWAC mooring, 2012, Oct 27 to 30

    !             0.0418 Hz to 0.1 Hz : ki=5e-6

    !             0.1 Hz to 0.12 Hz : ki=7e-6

    !             0.12 Hz to 0.16 Hz : ki=15e-6

    !             0.16 Hz to 0.5 Hz : ki=100e-6

    !             ICEP1=KI1=5.0e-6

    !             ICEP2=KI2=7.0e-6

    !             ICEP3=KI3=15.0e-6

    !             ICEP4=KI4=100.0e-6

    !             ICEP5=FC5=0.10

    !             MUDD=FC6=0.12

    !             MUDT=FC7=0.16

    !             In terms of the 3-character IDs for "Homogeneous field

    !             data" in ww3_shel.inp, these are, respectively, IC1, IC2,

    !             IC3, IC4, IC5, MDN, MTH, and so this might look like:

    !                'IC1' 19680606 000000   5.0e-6

    !                'IC2' 19680606 000000   7.0e-6

    !                'IC3' 19680606 000000   15.0e-6

    !                'IC4' 19680606 000000   100.0e-6

    !                'IC5' 19680606 000000   0.10

    !                'MDN' 19680606 000000   0.12

    !                'MTH' 19680606 000000   0.16

    !

    !     6) Simple step function for ki as a function of frequency

    !          with up to 16 "steps". Controlling parameters KIx and FCx are

    !          read in as namelist parameters, so they are stationary and

    !          uniform. (If 16 steps is not enough, the number of steps can be

    !          increased at compile time by changing NIC4 in w3gdatmd.ftn.)

    !          The last non-zero FCx value should be a large number, e.g. 99 Hz

    !

    !                                          4444444444  <--- ki=ic4_ki(4)

    !                               3333333333             <--- ki=ic4_ki(3)

    !                   2222222222                         <--- ki=ic4_ki(2)

    !      11111111111                                     <--- ki=ic4_ki(1)

    !                 ^           ^           ^         ^

    !                 |           |           |         |

    !          ic4_fc(1)   ic4_fc(2)   ic4_fc(3) ic4_fc(4)=large number

    !       Example: Beaufort Sea, AWAC mooring, 2012, Oct 27 to 30

    !           &SIC4  IC4METHOD = 6,

    !                  IC4KI =    0.50E-05,   0.70E-05,   0.15E-04,

    !                             0.10E+00,   0.00E+00,   0.00E+00,

    !                             0.00E+00,   0.00E+00,   0.00E+00,

    !                             0.00E+00,

    !                  IC4FC =    0.100,      0.120,      0.160,

    !                             99.00,      0.000,      0.000,

    !                             0.000,      0.000,      0.000,

    !                             0.000

    !                             /

    !

    !     7) Doble et al. (GRL 2015), eq. 3. This is a function of ice

    !        thickness and wave period.

    !        ALPHA  = 0.2*(T^(-2.13)*HICE or

    !        ALPHA  = 0.2*(FREQ^2.13)*HICE

    !

    !     8) Meylan et al. (JGR 2018), eq. 48. "Model with Order 3 Power

    !        Law". The is denoted as the "M2" model by Liu et al. (JPO 2020)

    !        It is a function of ice thickness and wave period.

    !        ki  = ChfM2*h_ice*freq^3

    !        where ChfM2 is a coefficient of proportionality which formally

    !        includes viscosity, density, and gravity parameters, see

    !        Meylan et al. (JGR 2018) for details.

    !        ChfM2 has units of s3/m2

    !        It is equation 53 in Meylan et al. (2018) and equation 16 in

    !        Liu et al. (2020).

    !        This method is functionally the same as the "M2" model in IC5

    !        in WW3 (IC5 w/IC5VEMOD=3) and is redundantly included here as

    !        IC4M8 because it is in the same "family" as IC4M7 and IC4M9,

    !        being in the form of:

    !        ki=Chf * h_ice^m * freq^n .

    !        Calibrations:

    !        * Liu et al. has ChfM2=eta*(2*pi)^3/(1025*9.81^2)

    !        ** eta=14.0 for "Sikuliaq" case of Liu et al., so ChfM2=0.035

    !        ** eta=3.0 for "SIPEX" case of Liu et al., so ChfM2=0.0075

    !        * Rogers et al. (tech rep. 2021, "RYW2021") :

    !        ** Fit to Rogers et al. (CRST 2021 "RMK2021") ChfM2=0.059 (*SD*)

    !        suggested default is marked with "(*SD*)", for consistency

    !          with SWAN (v41.31AB or later)

    !

    !     9) Rogers et al. (tech. rep. 2021, "RYW2021"): the "monomial power

    !        fit" described in section 2.2.3. It is the general form above,

    !        ki=Chf * h_ice^m * freq^n but is constrained such that m=n/2-1.

    !        This constraint is derived by RYW2021 by invoking the scaling from

    !        Yu et al. (2019), which is based on Reynolds number with ice

    !        thickness as the relevant length scale.

    !        This is also given as equation 2 in Yu et al. (CRST 2022).

    !        Some calibrations are as follows:

    !        * RYW2021, calibration to RMK2021: Chf=2.9 and n=4.5  (*SD*)

    !        * Yu et al. (2022) calibration to RMK2021 : Chf=2.4 and n=4.46

    !          (noting that c_n=0.108 and Chf=c_n*(2*pi/sqrt(g))^n)

    !        * Yu (2022) adjusted the prior calibration to get better fit

    !          to higher frequency lab measurements and got:

    !          Chf=7.89 and n=4.8

    !        suggested default is marked with "(*SD*)", for consistency

    !          with SWAN (v41.31AB or later)

    !

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

    !

    !     For all methods, the user can specify namelist

    !     variables IC4FMIN and IC4KIBK such as:

    !     &SIC4 IC4METHOD = [...], IC4FMIN=0.08, IC4KIBK=1.0e-7, [...]

    !     This accomodates the situation where the empirically-derived

    !     dissipation is uncertain for the lowest frequencies, which can be

    !     the case if estimated dissipation rate is so small that it falls

    !     in the noise level for the estimation method. (This is common,

    !     since some ice types cause only very weak dissipation

    !     to low frequencies.) In the example above, the amplitude

    !     dissipation rate ki is set to some low background level

    !     dissipation IC4KIBK=1.0e-7 1/m when model frequency is less than

    !     0.08 Hz.

    !

    !     More verbose description of implementation of Sice in WW3:

    !      See documentation for IC1

    !

    !     Notes regarding numerics:

    !      See documentation for IC1

    !

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

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

    !     For questions, comments and/or corrections, please refer to:

    !        Method 1 : C. Collins

    !        Method 2 : C. Collins

    !        Method 3 : C. Collins

    !        Method 4 : C. Collins

    !        Method 5 : E. Rogers

    !        Method 6 : E. Rogers

    !        Method 7 : E. Rogers

    !

    !     ALPHA = 2 * WN_I

    !     Though it may seem redundant/unnecessary to have *both* in the

    !       code, we do it this way to make the code easier to read and

    !       relate to other codes and source material, and hopefully avoid

    !       mistakes.

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

    !

    !  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

    USE w3odatmd, ONLY: ndse

    USE w3servmd, ONLY: extcde

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

                        ic4_ki, ic4_fc, ic4_cn, nic4, ic4_fmin,    &

                        ic4_kibk

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

                        mudt, mudv, mudd, 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),   a(nspec), depth

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

    INTEGER, INTENT(IN)     :: ix, iy

    !/

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

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: ient = 0

#endif

#ifdef W3_T0

    INTEGER                 :: ith

    REAL                    :: dout(nk,nth)

#endif

    INTEGER                 :: ikth, ik, ith, ic4method, ifc

    REAL                    :: d1d(nk), eb(nk)

    REAL                    :: icecoef1, icecoef2, icecoef3, &

                               icecoef4, icecoef5, icecoef6, &

                               icecoef7, icecoef8

    REAL                    :: cice1,cice2,cice3,cice4,cice5 ! temporary variables

    REAL                    :: ki1,ki2,ki3,ki4,fc5,fc6,fc7

    REAL                    :: hs, emean, hice

    REAL                    :: chf,mpow,npow

    REAL, ALLOCATABLE       :: wn_i(:)  ! exponential decay rate for amplitude

    REAL, ALLOCATABLE       :: alpha(:) ! exponential decay rate for energy

    REAL, ALLOCATABLE       :: freq(:) ! wave frequency

    REAL, ALLOCATABLE       :: marg1(:), marg2(:) ! Arguments for M2

    REAL, ALLOCATABLE       :: karg1(:), karg2(:), karg3(:) !Arguments for M3

    LOGICAL                 :: nml_input ! if using namelist input for M2


    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SIC4')

#endif

    !

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

    !

    d        = 0.0

    !

    ALLOCATE(wn_i(0:nk+1))

    ALLOCATE(alpha(0:nk+1))

    ALLOCATE(marg1(0:nk+1))

    ALLOCATE(marg2(0:nk+1))

    ALLOCATE(karg1(0:nk+1))

    ALLOCATE(karg2(0:nk+1))

    ALLOCATE(karg3(0:nk+1))

    ALLOCATE(freq(0:nk+1))

    marg1    = 0.0

    marg2    = 0.0

    karg1    = 0.0

    karg2    = 0.0

    karg3    = 0.0

    wn_i     = 0.0

    alpha    = 0.0

    icecoef1 = 0.0

    icecoef2 = 0.0

    icecoef3 = 0.0

    icecoef4 = 0.0

    icecoef5 = 0.0

    icecoef6 = 0.0

    icecoef7 = 0.0

    icecoef8 = 0.0

    hs       = 0.0

    hice     = 0.0

    emean    = 0.0

    freq=sig/tpi

    !

    !     IF (.NOT.INFLAGS2(-7))THEN

    !        WRITE (NDSE,1001) 'ICE PARAMETER 1'

    !        CALL EXTCDE(201)

    !     ENDIF

    !

    !   We cannot remove the other use of INFLAGS below,

    !   because we would get 'array not allocated' error for the methods

    !   that don't use MUDV, etc. and don't have MUDV allocated.


    IF (inflags2(-7)) icecoef1 = icep1(ix,iy) ! a.k.a. IC1

    IF (inflags2(-6)) icecoef2 = icep2(ix,iy) ! etc.

    IF (inflags2(-5)) icecoef3 = icep3(ix,iy)

    IF (inflags2(-4)) icecoef4 = icep4(ix,iy)

    IF (inflags2(-3)) icecoef5 = icep5(ix,iy)


    ! Borrow from Smud (error if BT8 or BT9)

#ifdef W3_BT8

    WRITE (ndse,*) 'DUPLICATE USE OF MUD PARAMETERS'

    CALL extcde(202)

#endif

#ifdef W3_BT9

    WRITE (ndse,*) 'DUPLICATE USE OF MUD PARAMETERS'

    CALL extcde(202)

#endif

    IF (inflags2(-2)) icecoef6 = mudd(ix,iy) ! a.k.a. MDN

    IF (inflags2(-1)) icecoef7 = mudt(ix,iy) ! a.k.a. MTH

    IF (inflags2(0 )) icecoef8 = mudv(ix,iy) ! a.k.a. MVS


    ic4method = ic4pars(1)

    !

#ifdef W3_T38

    !     Write test output ---------------------------------------------- /

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

#endif

    !

    ! 1.  Make calculations ---------------------------------------------- /

    !

    ! 1.a Calculate WN_I


    SELECT CASE (ic4method)


    CASE (1) ! IC4M1 : Exponential fit to Wadhams et al. 1988

      alpha = exp(-icecoef1 * tpi / sig - icecoef2)

      wn_i = 0.5 * alpha


    CASE (2) ! IC4M2 : Polynomial fit, Eq. 3 from Meylan et al. 2014

      !NB: Eq. 3 only includes T^2 and T^4 terms,

      !  which correspond to ICECOEF3, ICECOEF5, so in

      !  regtest: ICECOEF1=ICECOEF2=ICECOEF4=0


      nml_input=.true.

      IF (inflags2(-7).OR.inflags2(-6).OR.inflags2(-5).OR. &

          inflags2(-4).OR.inflags2(-3)) nml_input=.false.


      IF(nml_input)THEN ! get from namelist array


         cice1=ic4_cn(1)

         cice2=ic4_cn(2)

         cice3=ic4_cn(3)

         cice4=ic4_cn(4)

         cice5=ic4_cn(5)


      ELSE ! get from input-field array (ICEP1 etc.)


         cice1=icecoef1

         cice2=icecoef2

         cice3=icecoef3

         cice4=icecoef4

         cice5=icecoef5


      ENDIF


      ! CICE1 is C_{ice,1} in Collins and Rogers (2017), for example.

      marg1 = cice1 + cice2*freq + cice3*freq**2

      marg2 = cice4*freq**3 + cice5*freq**4

      alpha = marg1 + marg2

      wn_i = 0.5 * alpha


    CASE (3) ! IC4M3 : Quadratic fit to Kohout & Meylan'08 in Horvat & Tziperman'15

      hice=icecoef1 ! For this method, ICECOEF1=ice thickness

      karg1 = -0.3203 + 2.058*hice - 0.9375*(tpi/sig)

      karg2 = -0.4269*hice**2 + 0.1566*hice*(tpi/sig)

      karg3 =  0.0006 * (tpi/sig)**2

      alpha  = exp(karg1 + karg2 + karg3)

      wn_i = 0.5 * alpha


    CASE (4) !Eq. 1 from Kohout et al. 2014

      !Calculate HS

      DO ik=1, nk

        eb(ik) = 0.

        DO ith=1, nth

          eb(ik) = eb(ik) + a(ith+(ik-1)*nth)

        END DO

      END DO

      DO ik=1, nk

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

        emean  = emean + eb(ik)

      END DO

      hs = 4.*sqrt( max(0.,emean) )

      ! If Hs < 3 m then do Hs dependent calc, otherwise dH/dx is a constant

      IF (hs <= 3) THEN

        wn_i=icecoef1 ! from: DHDX=ICECOEF1*HS and WN_I=DHDX/HS

      ELSE IF (hs > 3) THEN

        wn_i=icecoef2/hs ! from: DHDX=ICECOEF2 and WN_I=DHDX/HS

      END IF


    CASE (5) ! Simple step function (time- and/or space-varying)

      ! rename variables for clarity

      ki1=icecoef1

      ki2=icecoef2

      ki3=icecoef3

      ki4=icecoef4

      fc5=icecoef5

      fc6=icecoef6

      fc7=icecoef7

      IF((ki1.EQ.0.0).OR.(ki2.EQ.0.0).OR.(ki3.EQ.0.0).OR. &

           (ki4.EQ.0.0).OR.(fc5.EQ.0.0).OR.(fc6.EQ.0.0).OR. &

           (fc7.EQ.0.0))THEN

        WRITE (ndse,1001)'ICE PARAMETERS'

        CALL extcde(201)

      END IF

      DO ik=1, nk

        ! select ki

        IF(freq(ik).LT.fc5)THEN

          wn_i(ik)=ki1

        ELSEIF(freq(ik).LT.fc6)THEN

          wn_i(ik)=ki2

        ELSEIF(freq(ik).LT.fc7)THEN

          wn_i(ik)=ki3

        ELSE

          wn_i(ik)=ki4

        ENDIF

      END DO


    CASE (6) ! Simple step function (from namelist)


      ! error checking: require at least 3 steps

      IF((ic4_ki(1).EQ.0.0).OR.(ic4_ki(2).EQ.0.0).OR. &

           (ic4_ki(3).EQ.0.0).OR.(ic4_fc(1).EQ.0.0).OR. &

           (ic4_fc(2).EQ.0.0) )THEN

        WRITE (ndse,1001)'ICE PARAMETERS'

        CALL extcde(201)

      END IF


      DO ik=1, nk

        ! select ki

        DO ifc=1,nic4

          IF(freq(ik).LT.ic4_fc(ifc))THEN

            wn_i(ik)=ic4_ki(ifc)

            EXIT

          END IF

        END DO

      END DO


    CASE (7) ! Doble et al. (GRL 2015)


      hice=icecoef1 ! For this method, ICECOEF1=ice thickness

      DO ik=1,nk

        alpha(ik)  = 0.2*(freq(ik)**2.13)*hice

      END DO

      wn_i= 0.5 * alpha


    CASE (8) ! Meylan et al. (JGR 2018), Liu et al. (JPO 2020)


      nml_input=.true.

      IF (inflags2(-6)) nml_input=.false.


      IF(nml_input)THEN ! get from namelist array


        chf=ic4_cn(1) ! Denoted "ChfM2" in documentation


      ELSE ! get from input-field array (ICEP1 etc.)


        chf=icecoef2 ! Denoted "ChfM2" in documentation


      ENDIF


      ! Rename variable, for clarity

      hice=icecoef1 ! For this method, ICECOEF1 is ice thickness


      DO ik=1,nk

        wn_i(ik)  = chf*hice*(freq(ik)**3)

      END DO


    CASE (9) ! Rogers et al. (2021) (RYW2021), Yu et al. (JGR 2022)


      nml_input=.true.

      IF (inflags2(-6).OR.inflags2(-5)) nml_input=.false.


      IF(nml_input)THEN ! get from namelist array


        chf=ic4_cn(1) ! Denoted as same in documentation

        npow=ic4_cn(2) ! Denoted "n" in documentation


      ELSE ! get from input-field array (ICEP1 etc.)


        chf=icecoef2 ! Denoted as same in documentation

        npow=icecoef3 ! Denoted "n" in documentation


      ENDIF


      ! Rename variable, for clarity

      hice=icecoef1 ! For this method, ICECOEF1 is ice thickness

      ! Compute

      mpow=0.5*npow-1.0 ! Denoted "m" in documentation

      DO ik=1,nk

        wn_i(ik)  = chf*(hice**mpow)*(freq(ik)**npow)

      END DO


    CASE DEFAULT

      wn_i = icecoef1 !Default to IC1: Uniform in k


    END SELECT


    !

    ! 1.b Calculate DID

    !

    DO ik=1, nk

      !   SBT1 has: D1D(IK) = FACTOR *  MAX(0., (CG(IK)*WN(IK)/SIG(IK)-0.5) )

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

      IF(freq(ik).LT.ic4_fmin)wn_i(ik)=ic4_kibk

      !           write(*,*)freq(ik),wn_i(ik),ICECOEF1,' % :: freq,ki,hice' ! temporary code: do not commit to repo uncommented

      d1d(ik) = -2. * cg(ik) * wn_i(ik)


    END DO


    !

    ! 1.c Fill diagional matrix

    !

    DO ikth=1, nspec

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

    END DO

    !

    !      END IF

    !

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

        '     ',a,' REQUIRED BUT NOT SELECTED'/)

    !

#ifdef W3_T

9000 FORMAT (' TEST W3SIC4 : DEPTH,ICECOEF1  : ',2e10.3)

#endif

    !/

    !/ End of W3SIC4 --------------------------------------------------- /

    !/

  END SUBROUTINE w3sic4

  !/

  !/ End of module W3SIC4MD ------------------------------------------ /

  !/

END MODULE w3sic4md