w3sic4md Module Reference

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

Functions/Subroutines
subroutine, public	w3sic4 (A, DEPTH, CG, IX, IY, S, D)
	S_{ice} source term using 5 parameters read from input files. More...

Detailed Description

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

Author

C. Collins

E. Rogers

Date

21-Jan-2015

Copyright

Copyright 2009-2022 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.

Function/Subroutine Documentation

w3sic4()

subroutine, public w3sic4md::w3sic4	(	real, dimension(nspec), intent(in)	A,
		real, intent(in)	DEPTH,
		real, dimension(nk), intent(in)	CG,
		integer, intent(in)	IX,
		integer, intent(in)	IY,
		real, dimension(nspec), intent(out)	S,
		real, dimension(nspec), intent(out)	D
	)

S_{ice} source term using 5 parameters read from input files.

Parameters

[in]	A	Action density spectrum (1-D).
[in]	DEPTH	Local water depth.
[in]	CG	Group velocities.
[in]	IX	Grid index.
[in]	IY	Grid index.
[out]	S	Source term (1-D version).
[out]	D	Diagonal term of derivative (1-D version).

Author

C. Collins

E. Rogers

Date

24-Feb-2017

Definition at line 121 of file w3sic4md.F90.

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

References w3gdatmd::dden, w3servmd::extcde(), w3gdatmd::ic4_cn, w3gdatmd::ic4_fc, w3gdatmd::ic4_fmin, w3gdatmd::ic4_ki, w3gdatmd::ic4_kibk, w3gdatmd::ic4pars, w3idatmd::inflags2, w3gdatmd::mapwn, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nic4, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, w3arrymd::outmat(), w3arrymd::prt2ds(), w3gdatmd::sig, w3servmd::strace(), and constants::tpi.

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