Skip to content

Commit 643cb5a

Browse files
authored
branch=IC4_safety_fix
1 parent 1b77ef2 commit 643cb5a

File tree

3 files changed

+41
-21
lines changed

3 files changed

+41
-21
lines changed

manual/eqs/ICE4.tex

Lines changed: 11 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -5,9 +5,9 @@ \subsubsection{~$S_{ice}$: Empirical/parametric damping by sea ice} \label{sec:I
55
\opthead{IC4}{\ws/NRL}{C. Collins and E. Rogers}
66

77
\noindent
8-
The fourth option ({\code IC4}) for damping of waves by sea ice was introduced by \cite{rep:CR17}. It gives methods to implement one of several simple, empirical/parametric forms for the dissipation of wave energy by sea ice. The motivation for {\code IC4} is to provide a simple, flexible, and efficient source term which reproduces, albeit in a highly parameterized way, some basic physics of wave-ice interaction. The method is set by the integer value (presently 1 to 7) for {\code IC4METHOD} namelist parameter: 1) an exponential fit to the field data of \cite{art:WAD88}, 2) the polynomial fit in \cite{art:MBK14}, 3) a quadratic fit to the calculations of \cite{art:KM08} given in \cite{art:HT15}, 4) Eq. 1 of \cite{art:Ko14}, 5) a simple step function with up to 4 steps (may be nonstationary and non-uniform), and 6) a simple step function with up to 10 steps (must be stationary and uniform), and 7) a formula from \cite{art:Dob15} which uses ice thickness. All but the fourth method of {\code IC4} feature frequency-dependent attenuation. With the fourth method, attenuation varies with waveheight but is uniform in frequency space.
8+
The fourth option ({\code IC4}) for damping of waves by sea ice was introduced (in its earliest form) by \cite{rep:CR17}. It gives methods to implement one of several simple, empirical/parametric forms for the dissipation of wave energy by sea ice. The motivation for {\code IC4} is to provide a simple, flexible, and efficient source term which reproduces, albeit in a highly parameterized way, some basic physics of wave-ice interaction. The method is set by the integer value (presently 1 to 10) for {\code IC4METHOD} namelist parameter. The first six are: 1) an exponential fit to the field data of \cite{art:WAD88}, 2) the polynomial fit in \cite{art:MBK14}, 3) a quadratic fit to the calculations of \cite{art:KM08} given in \cite{art:HT15}, 4) Eq. 1 of \cite{art:Ko14}, 5) a simple step function with up to 4 steps (may be nonstationary and non-uniform), and 6) a simple step function with up to 10 steps (must be stationary and uniform). Methods 7, 8, and 9 use multivariate power fits with dependence on frequency and ice thickness. All but the fourth method of {\code IC4} feature frequency-dependent attenuation. With the fourth method, attenuation varies with waveheight but is uniform in frequency space.
99

10-
In the following discussion we use {\code IC4M1} to denote {\code IC4} method 1, and so forth. {\code IC4} appears in the {\file switch} and namelist {\code IC4METHOD=1} (for example) appears in the file {\file ww3\_grid.inp}. Whereas in {\code IC1}, ${C_{ice,1}}$ is the user-determined attenuation, for {\code IC4M1}, {\code IC4M2}, and {\code IC4M4} ${C_{ice,n}}$ are constants of the equations. For {\code IC4M3}, ${C_{ice,1}}$ is ice thickness. For {\code IC4M5}, ${C_{ice,n}}$ controls the step function. Note that ${C_{ice,n}}$ may be provided by the user as non-stationary and non-uniform using methods analogous to methods used to input water levels.
10+
In the following discussion we use {\code IC4M1} to denote {\code IC4} method 1, and so forth. {\code IC4} appears in the {\file switch} and namelist {\code IC4METHOD=1} (for example) appears in the file {\file ww3\_grid.inp}. Whereas in {\code IC1}, ${C_{ice,1}}$ is the user-determined attenuation, here the meaning of ${C_{ice,n}}$ is context-dependent. With {\code IC4M1}, {\code IC4M2}, and {\code IC4M4} ${C_{ice,n}}$ are constants of the equations. For {\code IC4M3}, {\code IC4M7}, {\code IC4M8}, and {\code IC4M9}, ${C_{ice,1}}$ is ice thickness. For {\code IC4M5}, ${C_{ice,n}}$ controls the step function. For {\code IC4M6}, namelist variables control the step function. Note that, for methods which use it, ${C_{ice,n}}$ may be provided by the user as non-stationary and non-uniform using methods analogous to methods used to input water levels.
1111

1212
{\code IC4M1}: an exponential equation was chosen to fit the data contained in table 2 of \cite{art:WAD88} which results in preferential attenuation of high frequency waves. This parameterizes the well-known low-pass filtering effect of ice. The equation has the following form:
1313
\begin{equation}\label{eq:ice1}
@@ -26,13 +26,15 @@ \subsubsection{~$S_{ice}$: Empirical/parametric damping by sea ice} \label{sec:I
2626

2727
With appropriate coefficients, this polynomial method can be used to reproduce the so-called “roll-over effect” where the attenuation is non-monotonic in frequency space. However, some recent studies do not indicate this effect, e.g. \cite{art:RTS16} and \cite{art:LK17}, and it may just be a spurious artifact in prior observational studies.
2828

29+
Though ${C_{ice,1...5}}$ can be specified to vary in time and space, this feature is rarely used in practice. Most users will prefer to use constant values, and for convenience, these can be set using namelist parameters, where $C_{ice,i}$ is specified as {\code IC4CN(i)} in the {\code SIC4} namelist rather than as an input field using {\file ww3\_prep}, etc.
30+
2931
{\code IC4M3}: \cite{art:HT15} fit a quadratic equation to the attenuation coefficient calculated by \cite{art:KM08} as a function of frequency, $T$, and ice thickness, $h$. Attenuation increases for thicker ice and higher frequencies (lower periods). The number of coefficients of the quadratic equation were prohibitively large to be user-determined, so the equation is hardwired in and the tunable parameter, ${C_{ice,1}}$, is ice thickness $h$. This method is described and applied in \cite{rep:CR17}. For reference, the equation is the following:
3032
\begin{equation}\label{eq:ice3}
3133
{\ln{\alpha(T,h)}} = -0.3203 + 2.058h - 0.9375T - 0.4269h^2 + 0.1566hT + 0.0006T^2
3234
\end{equation}
3335

3436
\noindent
35-
There are two warnings to make about {\code IC4M3}. First, the equation itself was an extrapolation of the original range of $h$ used to calculate the attenuation coefficients in \cite{art:KM08} which was between 0.5 and 3 m, see \cite{art:HT15}. Second, in \cite{art:KM08}, wave attenuation predicted is based on scattering (a conservative process), whereas in {\code IC4M3}, the wave attenuation is treated as dissipation (non-conservative). This is ad hoc and not recommended for general use. Most especially, users should think twice before using {\code IC4M3} in combination with scattering routines {\code IS1} or {\code IS2}, since this is essentially double-counting scattering.
37+
There are four warnings to make about {\code IC4M3}. First, the equation itself was an extrapolation of the original range of $h$ used to calculate the attenuation coefficients in \cite{art:KM08} which was between 0.5 and 3 m, see \cite{art:HT15}. Second, in \cite{art:KM08}, wave attenuation predicted is based on scattering (a conservative process), whereas in {\code IC4M3}, the wave attenuation is treated as dissipation (non-conservative). Third, we recommend against using {\code IC4M3} in combination with scattering routines {\code IS1} or {\code IS2}, since this would include scattering twice. Fourth, the implementation assumes that one floe is encountered per meter. This is a departure from \cite{art:HT15} which allows this length scale to vary. The reader is referred to the inline documentation in {\file w3sic4md.F90} for further detail.
3638

3739
{\code IC4M4}: \cite{art:Ko14} found that attenuation was a function of significant wave height. Attenuation increased linearly with ${H_s}$ until ${H_s} = 3$ m at which point attenuation is capped, thus:
3840
\begin{equation}
@@ -44,11 +46,11 @@ \subsubsection{~$S_{ice}$: Empirical/parametric damping by sea ice} \label{sec:I
4446
\end{equation}
4547
where {$k_i=\frac{\partial H_s}{\partial dx}/H_s$}.
4648

47-
The values given in \cite{art:Ko14} are ${C_{ice,1...2}}=[5.35\times 10^{-6}, 16.05\times 10^{-6}]$. See regression test {\file ww3\_tic1.1/input\_IC4/M4} for examples. This method is described and applied in \cite{rep:CR17}.
49+
The values given in \cite{art:Ko14} are ${C_{ice,1...2}}=[5.35\times 10^{-6}, 16.05\times 10^{-6}]$. See regression test {\file ww3\_tic1.1/input\_IC4\_M4} for examples. This method is described and applied in \cite{rep:CR17}.
4850

49-
{\code IC4M5}: This is a simple step function with up to 4 steps. It is controlled by the optionally nonstationary and non-uniform parameters ${C_{ice,1...7}}$. Parameters ${C_{ice,1...4}}$ control the step levels, which are in terms of dissipation rate, ${k_i}$. Parameters ${C_{ice,5...7}}$ control the step boundaries (given in Hz). See regression test {\file ww3\_tic1.1/input\_IC4/M5} for examples. This method is described in \cite{rep:CR17}.
51+
{\code IC4M5}: This is a simple step function with up to 4 steps. It is controlled by the optionally nonstationary and non-uniform parameters ${C_{ice,1...7}}$. Parameters ${C_{ice,1...4}}$ control the step levels, which are in terms of dissipation rate, ${k_i}$. Parameters ${C_{ice,5...7}}$ control the step boundaries (given in Hz). See regression test {\file ww3\_tic1.1/input\_IC4\_M5} for examples. This method is described in \cite{rep:CR17}.
5052

51-
{\code IC4M6}: This is a simple step function with up to 10 steps. It is controlled by the stationary and uniform namelist parameters {\code IC4KI} and {\code IC4FC}. Array {\code IC4KI} controls the step levels, which are in terms of dissipation rate, ${k_i}$, in radians per meter. Array {\code IC4FC} controls the step boundaries (given in Hz). See regression test {\file ww3\_tic1.1/input\_IC4/M6} for examples.
53+
{\code IC4M6}: This is a simple step function with up to 10 steps. It is controlled by the stationary and uniform namelist parameters {\code IC4KI} and {\code IC4FC}. Array {\code IC4KI} controls the step levels, which are in terms of dissipation rate, ${k_i}$, in radians per meter. Array {\code IC4FC} controls the step boundaries (given in Hz). See regression test {\file ww3\_tic1.1/input\_IC4\_M6} for examples.
5254

5355
{\code IC4M7}: This is a formula for dissipation from \cite{art:Dob15}, developed for a mixture of pancake and frazil ice, using data collected in the Weddell Sea (Antarctica). The formula depends on wave frequency and ice thickness:
5456
\begin{equation}\label{eq:ice7}
@@ -60,12 +62,12 @@ \subsubsection{~$S_{ice}$: Empirical/parametric damping by sea ice} \label{sec:I
6062
\begin{equation}\label{eq:ice8}
6163
{k_i=C_{hf}h^mf^n} \:\:\: .
6264
\end{equation}
63-
The formula is taken from \cite{Meylan2018}, where it is described as a ``Model with Order 3 Power Law''. It is applied by \cite{Liu2020}, where it is referred to as the ``M2'' model. The model specifies $m=1$ and $n=3$, and $C_{hf}$ is a user-specified calibration coefficient. \cite{Liu2020} provide calibration to two field cases and \cite{rep:RYW2021} provides a calibration to a third field case, \cite{art:RMK2021}. The third calibration is set as the default for {\code IC4M8}, $C_{hf}=0.059$, but can be changed in using the namelist parameter (constant and uniform) {\code IC4CN}, or using the spatially and/or temporally variable parameter ${C_{ice,2}}$ . Further details on the calibrations are available in the inline documentation in {\file w3sic4md.F90}. This method is functionally the same as the ``{\code M2}'' model in {\code IC5} (i.e., {\code IC5} with {\code IC5VEMOD=3}) and is redundantly included here as {\code IC4M8} because it is in the same ``family'' as {\code IC4M7} and {\code IC4M9}, being in the form of Eq. (\ref{eq:ice8}).
65+
The formula is taken from \cite{Meylan2018}, where it is described as a ``Model with Order 3 Power Law''. It is applied by \cite{Liu2020}, where it is referred to as the ``M2'' model. The model specifies $m=1$ and $n=3$, and $C_{hf}$ is a user-specified calibration coefficient. \cite{Liu2020} provide calibration to two field cases and \cite{rep:RYW2021} provides a calibration to a third field case, \cite{art:RMK2021}. The third calibration is set as the default for {\code IC4M8}, $C_{hf}=0.059$, but can be changed via the namelist parameter {\code IC4CN} (constant and uniform settings), or by using the spatially and/or temporally variable parameter ${C_{ice,2}}$ . Further details on the calibrations are available in the inline documentation in {\file w3sic4md.F90}. This method is functionally the same as the ``{\code M2}'' model in {\code IC5} (i.e., {\code IC5} with {\code IC5VEMOD=3}) and is redundantly included here as {\code IC4M8} because it is in the same ``family'' as {\code IC4M7} and {\code IC4M9}, being in the form of Eq. (\ref{eq:ice8}).
6466

6567
For an example of setting the namelist parameter, see {\file /regtests/ww3\_tic1.1/input\_IC4\_M8}.
6668

67-
{\code IC4M9}: This formula is taken from the ``monomial power fit'' given in section 2.2.3 of \cite{rep:RYW2021}. Like {\code IC4M7} and {\code IC4M8}, it is a specific case of the general form of Eq. (\ref{eq:ice8}). The specificity is the constraint that $m=n/2-1$. This constraint is derived by \cite{rep:RYW2021} by invoking the scaling from \cite{art:YRW2019}, which is based on Reynolds number with ice thickness as the relevant length scale. This is also given as equation 2 in \cite{art:YRW2022}. The default namelist settings are $C_{hf}=2.9$ and $n=4.5$, from calibration by \cite{rep:RYW2021} to \cite{art:RMK2021}. Further details, including alternative calibrations such as \cite{art:Yu2022}, are available in the inline documentation in {\file w3sic4md.F90}. Constant values can be set using namelist parameters, where $C_{hf}$ and $n$ are {\code IC4CN(1)} and {\code IC4CN(2)}, respectively. Spatially and/or temporally versions of the same can be specified as ${C_{ice,2}}$ and ${C_{ice,3}}$, respectively.
69+
{\code IC4M9}: This formula is taken from the ``monomial power fit'' given in section 2.2.3 of \cite{rep:RYW2021}. Like {\code IC4M7} and {\code IC4M8}, it is a specific case of the general form of Eq. (\ref{eq:ice8}). The specificity is the constraint that $m=n/2-1$. This constraint is derived by \cite{rep:RYW2021} by invoking the scaling from \cite{art:YRW2019}, which is based on Reynolds number with ice thickness as the relevant length scale. This is also given as equation 2 in \cite{art:YRW2022}. The default namelist settings are $C_{hf}=2.9$ and $n=4.5$, from calibration by \cite{rep:RYW2021} and \cite{art:YRW2022} to \cite{art:RMK2021}. Further details, including alternative calibrations such as \cite{art:Yu2022}, are available in the inline documentation in {\file w3sic4md.F90}. Constant values can be set using namelist parameters, where $C_{hf}$ and $n$ are {\code IC4CN(1)} and {\code IC4CN(2)}, respectively. Spatial and/or temporal variability of the same can be specified using ${C_{ice,2}}$ and ${C_{ice,3}}$, respectively.
6870

6971
The namelist default $C_{hf}$ values in {\code IC4M8} and {\code IC4M9} are consistent with those of identical formulae implemented in \cite{man:SWAN4145A}.
7072

71-
73+
{\code IC4M10} is a method for attenuation due to scattering by sea ice floes, depending on ice thickness and floe size, based on \cite{art:MHB21}.

manual/manual.bib

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3226,6 +3226,22 @@ @ARTICLE{art:MBK14
32263226
number = "C14",
32273227
PAGES = "5,046--5,051" }
32283228

3229+
% item art:MHB21
3230+
3231+
@article{art:MHB21,
3232+
author = {Michael H. Meylan and Christopher Horvat and Cecilia M. Bitz and Luke G. Bennetts},
3233+
year = {2021},
3234+
title = {A floe size dependent scattering model in two- and three-dimensions for wave attenuation by ice floes},
3235+
journal=OMOD,
3236+
volume = {161},
3237+
pages = {101779},
3238+
issn = {1463-5003},
3239+
doi = {https://doi.org/10.1016/j.ocemod.2021.101779},
3240+
url = {https://www.sciencedirect.com/science/article/pii/S1463500321000299},
3241+
keywords = {Sea ice, Ocean waves, Scattering},
3242+
publisher={Elsevier}
3243+
}
3244+
32293245
% item art:HT15
32303246
32313247
@ARTICLE{art:HT15,

model/src/w3sic4md.F90

Lines changed: 14 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -140,6 +140,7 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
140140
!/ 11-Jan-2024 : Method 9 added (Rogers et al., 2021)
141141
!/ denoted "RYW2021" (E. Rogers)
142142
!/ 14-Aug-2024 : Method 10 added (Meylan et al. 2021) (E. Thomas)
143+
!/ 15-Aug-2025 : Safety fix for negative hice (E. Rogers)
143144
!/
144145
!/ FIXME : Move field input to W3SRCE and provide
145146
!/ (S.Zieger) input parameter to W3SIC1 to make the subroutine
@@ -370,20 +371,21 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
370371
! 7. Remarks :
371372
!
372373
! If ice parameter 1 is zero, no calculations are made.
373-
! For questions, comments and/or corrections, please refer to:
374-
! Method 1 : C. Collins
375-
! Method 2 : C. Collins
376-
! Method 3 : C. Collins
377-
! Method 4 : C. Collins
378-
! Method 5 : E. Rogers
379-
! Method 6 : E. Rogers
380-
! Method 7 : E. Rogers
374+
! For questions, comments and/or corrections, please contact the
375+
! authors in the updates list, above.
381376
!
382377
! ALPHA = 2 * WN_I
383378
! Though it may seem redundant/unnecessary to have *both* in the
384379
! code, we do it this way to make the code easier to read and
385380
! relate to other codes and source material, and hopefully avoid
386381
! mistakes.
382+
!
383+
! For sub-methods M3, M7, M8, and M9, ICECOEF1 (ICEP1) is used
384+
! to represent ice thickness. When ice thickness is taken from an
385+
! ice model such as CICE, we have encountered cases with spurious,
386+
! small, negative values, which can result in NaNs in WW3. Thus,
387+
! for these sub-methods, we set a lower limit of zero for ICECOEF1.
388+
!
387389
!/ ------------------------------------------------------------------- /
388390
!
389391
! 8. Structure :
@@ -576,7 +578,7 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
576578
WN_I = 0.5 * ALPHA
577579

578580
CASE (3) ! IC4M3 : Quadratic fit to Kohout & Meylan'08 in Horvat & Tziperman'15
579-
HICE=ICECOEF1 ! For this method, ICECOEF1=ice thickness
581+
HICE=MAX(0.0,ICECOEF1) ! For this method, ICECOEF1=ice thickness, which cannot be less than zero
580582
KARG1 = -0.3203 + 2.058*HICE - 0.9375*(TPI/SIG)
581583
KARG2 = -0.4269*HICE**2 + 0.1566*HICE*(TPI/SIG)
582584
KARG3 = 0.0006 * (TPI/SIG)**2
@@ -653,7 +655,7 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
653655

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

656-
HICE=ICECOEF1 ! For this method, ICECOEF1=ice thickness
658+
HICE=MAX(0.0,ICECOEF1) ! For this method, ICECOEF1=ice thickness, which cannot be less than zero
657659
DO IK=1,NK
658660
ALPHA(IK) = 0.2*(FREQ(IK)**2.13)*HICE
659661
END DO
@@ -675,7 +677,7 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
675677
ENDIF
676678

677679
! Rename variable, for clarity
678-
hice=ICECOEF1 ! For this method, ICECOEF1 is ice thickness
680+
hice=MAX(0.0,ICECOEF1) ! For this method, ICECOEF1 is ice thickness, which cannot be less than zero
679681

680682
DO IK=1,NK
681683
WN_I(IK) = Chf*hice*(FREQ(IK)**3)
@@ -699,7 +701,7 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
699701
ENDIF
700702

701703
! Rename variable, for clarity
702-
hice=ICECOEF1 ! For this method, ICECOEF1 is ice thickness
704+
hice=MAX(0.0,ICECOEF1) ! For this method, ICECOEF1 is ice thickness, which cannot be less than zero
703705
! Compute
704706
mpow=0.5*npow-1.0 ! Denoted "m" in documentation
705707
DO IK=1,NK

0 commit comments

Comments
 (0)