From 5c70e8926a8339f29312cdfbd83fbef0e2766dc6 Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Wed, 9 Nov 2022 09:05:27 -0500 Subject: [PATCH 01/15] w3snl1md.F90: add doxygen markup --- model/src/w3snl1md.F90 | 47 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/model/src/w3snl1md.F90 b/model/src/w3snl1md.F90 index 20ba8f5fa..6db0be264 100644 --- a/model/src/w3snl1md.F90 +++ b/model/src/w3snl1md.F90 @@ -1,5 +1,27 @@ +!> @file +!> @brief Bundles routines calculate nonlinear wave-wave interactions +!> according to the Discrete Interaction Approximation (DIA). +!> +!> @author H. L. Tolman +!> @date 03-Sep-2012 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Bundles routines calculate nonlinear wave-wave interactions +!> according to the Discrete Interaction Approximation (DIA) of +!> Hasselmann et al. (JPO, 1985). +!> +!> @author H. L. Tolman +!> @date 03-Sep-2012 +!> +!> @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. +!> + MODULE W3SNL1MD !/ !/ +-----------------------------------+ @@ -59,6 +81,21 @@ MODULE W3SNL1MD !/ CONTAINS !/ ------------------------------------------------------------------- / + +!> +!> @brief Calculate nonlinear interactions and the diagonal term of +!> its derivative. +!> +!> @param[in] A Action spectrum A(ISP) as a function of +!> direction (rad) and wavenumber. +!> @param[in] CG Group velocities (dimension NK). +!> @param[in] KDMEAN Mean relative depth. +!> @param[out] S Source term. +!> @param[out] D Diagonal term of derivative. +!> +!> @author H. L. Tolman +!> @date 06-Jun-2018 +!> SUBROUTINE W3SNL1 (A, CG, KDMEAN, S, D) !/ !/ +-----------------------------------+ @@ -418,7 +455,15 @@ SUBROUTINE W3SNL1 (A, CG, KDMEAN, S, D) !/ End of W3SNL1 ----------------------------------------------------- / !/ END SUBROUTINE W3SNL1 - !/ ------------------------------------------------------------------- / +!/ ------------------------------------------------------------------- / +!> +!> @brief Preprocessing for nonlinear interactions (weights). +!> +!> @param[in] IMOD Model number. +!> +!> @author H. L. Tolman +!> @date 24-Dec-2004 +!> SUBROUTINE INSNL1 ( IMOD ) !/ !/ +-----------------------------------+ From 9e0018b69bb45578cb14055247e1a668aa2c0da9 Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Wed, 9 Nov 2022 09:24:18 -0500 Subject: [PATCH 02/15] w3snl2md.F90: add doxygen markup --- model/src/w3snl2md.F90 | 43 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 42 insertions(+), 1 deletion(-) diff --git a/model/src/w3snl2md.F90 b/model/src/w3snl2md.F90 index 946851927..6008f4a8e 100644 --- a/model/src/w3snl2md.F90 +++ b/model/src/w3snl2md.F90 @@ -1,5 +1,25 @@ +!> @file +!> @brief Interface module to exact nonlinear interactions. +!> +!> @author H. L. Tolman +!> @author G. Ph. van Vledder +!> @date 29-May-2009 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Interface module to exact nonlinear interactions. +!> +!> @author H. L. Tolman +!> @author G. Ph. van Vledder +!> @date 29-May-2009 +!> +!> @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. +!> MODULE W3SNL2MD !/ !/ +-----------------------------------+ @@ -57,7 +77,21 @@ MODULE W3SNL2MD PUBLIC !/ CONTAINS - !/ ------------------------------------------------------------------- / +!/ ------------------------------------------------------------------- / +!> +!> @brief Interface to exact interactions. +!> +!> @param[in] A Action spectrum A(ITH,IK) as a function of +!> direction (rad) and wavenumber. +!> @param[in] CG Group velocities (dimension NK). +!> @param[in] DEPTH Water depth in meters. +!> @param[out] S Source term. +!> @param[out] D Diagonal term of derivative. +!> +!> @author H. L. Tolman +!> @author G. Ph. van Vledder +!> @date 24-Dec-2004 +!> SUBROUTINE W3SNL2 ( A, CG, DEPTH, S, D ) !/ !/ +-----------------------------------+ @@ -250,6 +284,13 @@ SUBROUTINE W3SNL2 ( A, CG, DEPTH, S, D ) !/ END SUBROUTINE W3SNL2 !/ ------------------------------------------------------------------- / + !> + !> @brief Preprocessing for nonlinear interactions (Xnl). + !> + !> @author H. L. Tolman + !> @author G. Ph. van Vledder + !> @date 24-Dec-2004 + !> SUBROUTINE INSNL2 !/ !/ +-----------------------------------+ From 7bf04fa80f0b14ed595cd7004e2c8ab566e78c62 Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Wed, 9 Nov 2022 10:32:05 -0500 Subject: [PATCH 03/15] w3snl3md.F90: add doxygen markup --- model/src/w3snl3md.F90 | 115 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) diff --git a/model/src/w3snl3md.F90 b/model/src/w3snl3md.F90 index 69eec9718..df12bec39 100644 --- a/model/src/w3snl3md.F90 +++ b/model/src/w3snl3md.F90 @@ -1,5 +1,26 @@ +!> @file +!> @brief Generalized and optimized multiple DIA implementation. +!> +!> @author H. L. Tolman +!> @date 13-Jul-2012 +!> + + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Generalized and optimized multiple DIA implementation. +!> +!> @details Expressions in terms of original F(f,theta) spectrum. +!> +!> @author H. L. Tolman +!> @date 13-Jul-2012 +!> +!> @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. +!> MODULE W3SNL3MD !/ !/ +-----------------------------------+ @@ -116,6 +137,46 @@ MODULE W3SNL3MD !/ CONTAINS !/ ------------------------------------------------------------------- / + + !> + !> @brief Multiple Discrete Interaction Parameterization for arbitrary + !> depths with generalized quadruplet layout. + !> + !> @details This is a direct implementation of the Discrete Interaction + !> Paramterization (DIA) with multiple representative quadruplets + !> (MDIA) for arbitrary water depths. + !> + !> The outer loop of the code is over quadruplet realizations, + !> which implies two realizations for a conventional quadruplet + !> definitions and four for extended definitions (with rescaling + !> of the contants for consistency). Within this loop the compu- + !> tations are performed in two stages. First, interactions + !> contributions are computed for the entire spectral space, + !> second all contributions are combined into the actual inter- + !> actions and diagonal contributions. + !> + !> Arbitrary depths are addressed by generating a lookup table + !> for the relative depth. These tables are used for each discrete + !> frequency separately. Efficient memory usages requires relative + !> addressing to reduce the size of the lookup tables. To use this + !> the spectral space is expanded to higher and lower frequencies + !> as well as directional space is expanded/volded. This is done + !> for the input (pseudo-) spectrum (action spectrum devided by the + !> wavenumber) to determine spectral densities at the quadruplet + !> components, and the spectral space describing individual contri- + !> butions before they are combined into the actual interactions. + !> + !> @param[in] A Action spectrum A(ITH,IK) as a function of + !> direction (rad) and wavenumber. + !> @param[in] CG Group velocities (dimension NK). + !> @param[in] WN Wavenumbers (dimension NK). + !> @param[in] DEPTH Water depth in meters. + !> @param[out] S Source term. + !> @param[out] D Diagonal term of derivative. + !> + !> @author H. L. Tolman + !> @date 01-Dec-2009 + !> SUBROUTINE W3SNL3 ( A, CG, WN, DEPTH, S, D ) !/ !/ +-----------------------------------+ @@ -566,6 +627,15 @@ SUBROUTINE W3SNL3 ( A, CG, WN, DEPTH, S, D ) !/ CONTAINS !/ ------------------------------------------------------------------- / + + !> + !> @brief Expand spectrum, subroutine used to simplify addressing. + !> + !> @param[out] SPEC Expanded spectrum. + !> + !> @author H. L. Tolman + !> @date 21-Aug-2009 + !> SUBROUTINE EXPAND ( SPEC ) !/ !/ +-----------------------------------+ @@ -627,6 +697,18 @@ SUBROUTINE EXPAND ( SPEC ) !/ END SUBROUTINE EXPAND !/ ------------------------------------------------------------------- / + + !> + !> @brief Expand spectrum to simplify indirect addressing. + !> + !> @details Done 'in place' with temporary array ( ARIN = AROUT ). + !> + !> @param[in] ARIN Input array. + !> @param[out] AROUT Output array. + !> + !> @author H. L. Tolman + !> @date 16-Jul-2008 + !> SUBROUTINE EXPND2 ( ARIN, AROUT ) !/ !/ +-----------------------------------+ @@ -689,6 +771,15 @@ END SUBROUTINE EXPND2 !/ END SUBROUTINE W3SNL3 !/ ------------------------------------------------------------------- / + !> + !> @brief Initialization for generalized multiple DIA routine. + !> + !> @details Fill storage aryays as described in the main subroutine + !> with interpolation, model and distribution data. + !> + !> @author H. L. Tolman + !> @date 13-Nov-2009 + !> SUBROUTINE INSNL3 !/ !/ +-----------------------------------+ @@ -1424,6 +1515,17 @@ SUBROUTINE INSNL3 !/ CONTAINS !/ ------------------------------------------------------------------- / + + !> + !> @brief Calculate minimum allowed lambda for quadruplet configuration. + !> + !> @param MU Quadruplet parameters. + !> @param THETA Theta in degrees. + !> @returns MINLAM Minimum allowed lambda. + !> + !> @author H. L. Tolman + !> @date 28-Jan-2004 + !> REAL FUNCTION MINLAM ( MU, THETA ) !/ !/ +-----------------------------------+ @@ -1479,6 +1581,19 @@ REAL FUNCTION MINLAM ( MU, THETA ) !/ END FUNCTION MINLAM !/ ------------------------------------------------------------------- / + + !> + !> @attention replaced (likely typo) 'minimum' from original header here + !> with 'maximum'. + !> @brief Calculate maximum allowed lambda for quadruplet configuration. + !> + !> @param MU Quadruplet parameters. + !> @param THETA Theta in degrees. + !> @returns MAXLAM Maximum allowed lambda. + !> + !> @author H. L. Tolman + !> @date 28-Jan-2004 + !> REAL FUNCTION MAXLAM ( MU, THETA ) !/ !/ +-----------------------------------+ From c82880ef15ee6c97ae6151f6a2b30edfd3bb09ab Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Sat, 12 Nov 2022 16:38:56 -0500 Subject: [PATCH 04/15] w3snl4md.F90: add doxygen markup --- model/src/w3snl4md.F90 | 365 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 364 insertions(+), 1 deletion(-) diff --git a/model/src/w3snl4md.F90 b/model/src/w3snl4md.F90 index 2c75b259c..aa9788b15 100644 --- a/model/src/w3snl4md.F90 +++ b/model/src/w3snl4md.F90 @@ -1,6 +1,27 @@ +!> @file +!> @brief Generic shallow-water Boltzmann integral (FBI or TSA). +!> +!> @author Bash Toulany +!> @author Michael Casey +!> @author William Perrie +!> @date 12-Apr-2016 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / !/ + +!> +!> @brief Generic shallow-water Boltzmann integral (FBI or TSA). +!> +!> @details Interface module for TSA type nonlinear interactions. +!> Based on Resio and Perrie (2008) and Perrie and Resio (2009). +!> +!> @author Bash Toulany +!> @author Michael Casey +!> @author William Perrie +!> @date 12-Apr-2016 +!> MODULE W3SNL4MD !/ !/ +-----------------------------------+ @@ -184,6 +205,16 @@ MODULE W3SNL4MD !! !! ------------------------------------------------------------------ !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + !> + !> @brief It reads look-up tables (generated by gridsetr) if file exists + !> otherwise it must generate the look-up tables file. + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE INSNL4 !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! ------------------------------------------------------------------ @@ -557,6 +588,25 @@ END SUBROUTINE INSNL4 !! !! ------------------------------------------------------------------ !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + !> + !> @brief Interface module for TSA type nonlinear interactions. + !> + !> @details Based on Resio and Perrie (2008) and Perrie and Resio (2009). + !> + !> @param[inout] A 2D Action Density A(NTH,NK) as function of + !> direction (rad) and wavenumber (theta,k). + !> @param[inout] CG Group velocities. + !> @param[inout] WN Wavenumbers. + !> @param[inout] DEPTH Water depth (m). + !> @param[inout] S Source term. + !> @param[inout] D Diagonal term of derivative. + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE W3SNL4 ( A, CG, WN, DEPTH, S, D ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !! ------------------------------------------------------------------ @@ -1392,6 +1442,41 @@ END SUBROUTINE W3SNL4 !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + !> + !> @brief This routine sets up the geometric part of the Boltzmann integral. + !> + !> @details This routine sets up the geometric part of the Boltzmann integral ! + !> based on a grid of wave frequencies and directions, with wave- ! + !> numbers related to frequency and depth by linear dispersion. It ! + !> is adapted from Don's original code with changes to modify the ! + !> indexing so there are fewer unused elements, and a number of algo- ! + !> rithmic changes that are mathematically equivalent to Don's but ! + !> take advantage of intrinsic functions to form smooth results with ! + !> less reliance on if statements. ! + !> ! + !> It calls locus-solving routines shloxr and shlocr and coupling ! + !> coefficient routine cplshr. If shlocr does not converge, ierr_gr ! + !> will be something other than 0 and the routine will terminate, ! + !> returning ierr_gr to the calling program (see shlocr). ! + !> ! + !> It returns array grad(,,), which is an estimate of the product ! + !> C(k1,k2,k3,k4)*H(k1,k3,k4)*ds/|dW/dn| (where n and the k's are all ! + !> vectors) as given, for example, by Eq.(7) of 'Nonlinear energy ! + !> fluxes and the finite depth equilibrium range in wave spectra,' ! + !> by Resio, Pihl, Tracy and Vincent (2001, JGR, 106(C4), p. 6985), ! + !> as well as arrays for indexing, interpolating and weighting locus- ! + !> based wavenumber vectors within the discrete solution grid. + !> + !> @param[in] dep Depth (m) + !> @param[in] wka1 Wavenumber array at one depth dim=(nrng) + !> @param[in] cgnrng1 Group velocity at nrng at one depth. + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE gridsetr ( dep, wka1, cgnrng1 ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !/ @@ -1861,6 +1946,67 @@ END SUBROUTINE gridsetr !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + !> + !> @brief General locus solution for input vectors (wk1x,wk1y). + !> + !> @verbatim + !> General locus solution for input vectors (wk1x,wk1y) and + !> (wk3x,wk3y) of the same magnitude but NOT in the same direction + !> (or singularness will occur), output vectors (wk2x,wk2y) and + !> (wk4x,wk4y), and element length ds along locus curve: + !> + !> With wavenumber vector n identified by (wknx,wkny), its magnitude + !> given by wkn = sqrt(wknx**2+wkny**2) and its associated radian + !> frequency given by sign = sqrt[g*wkn*tanh(wkn*dep)], where g is + !> gravitational acceleration and dep is water depth, the four-wave + !> resonance condition is satisfied along a locus of pts defined by + !> + !> [1] (wk1x,wk1y) + (wk2x,wk2y) - (wk3x,wk3y) - (wk4x,wk4y) = 0 + !> + !> [2] sig1 + sig2 - sig3 - sig4 = 0 + !> + !> In the case where k1 [= sqrt(wk1x**2+wk1y**2)] is equal to k3 + !> [= sqrt(wk3x**2+wk3y**2)], we have by dispersion, + !> + !> [3] sig1 = sqrt[g*k1*tanh(k1*h)] = sqrt[g*k3*tanh(k3*h)] = sig3 + !> + !> so sig1 - sig3 = 0 and [2] becomes sig2 = sig4, where, again by + !> dispersion, + !> + !> [4] sig2 = sqrt(g*k2*tanh(k2*h)] = sig4 = sqrt(g*k4*tanh(k4*h)] + !> + !> and consequently k2 = k4. This simplifies the locus solution + !> considerably, and it can be shown that the (wk2x,wk2y) locus is + !> along the perpendicular bisector of the (px,py) vector given by + !> + !> [5] (px,py) = (wk3x-wk1x,wk3y-wk1y) + !> + !> and thereby from [1] + !> [6] (wk4x,wk4y) = (wk2x,wk2y) - (px,py) + !> + !> We note that these loci are independent of depth, although depth + !> is used to set the length of the locus line by requiring that its + !> range on either side of the p vector correspond to a wave with a + !> wkx freq four times that of the k1 vector (the locus line is made + !> up of npts segments of length ds; the outer edges of the terminal + !> segments satisfy the length constraint; vectors k2 and k4 extend + !> to segment centers and will sufficiently approximate the length + !> constraint). As compared to srshlocr.f, we can do all + !> calculations here in dimensional space. + !> @endverbatim + !> + !> @param[in] dep + !> @param[in] wk1x + !> @param[in] wk1y + !> @param[in] wk3x + !> @param[in] wk3y + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE shloxr ( dep, wk1x,wk1y, wk3x,wk3y ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !/ @@ -1925,7 +2071,7 @@ SUBROUTINE shloxr ( dep, wk1x,wk1y, wk3x,wk3y ) !! We note that these loci are independent of depth, although depth ! !! is used to set the length of the locus line by requiring that its ! !! range on either side of the p vector correspond to a wave with a ! - !!wkx freq four times that of the k1 vector (the locus line is made ! + !! wkx freq four times that of the k1 vector (the locus line is made ! !! up of npts segments of length ds; the outer edges of the terminal ! !! segments satisfy the length constraint; vectors k2 and k4 extend ! !! to segment centers and will sufficiently approximate the length ! @@ -2095,6 +2241,67 @@ END SUBROUTINE shloxr !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + !> + !> @brief N/A + !> + !> @verbatim + !> With wavenumber vector n identified by (wknx,wkny), its magnitude + !> given by wkn = sqrt(wknx**2+wkny**2) and its associated radian + !> frequency given by sign = sqrt[g*wkn*tanh(wkn*dep)], where g is + !> gravitational acceleration and dep is water depth, the four-wave + !> resonance condition is satisfied along a locus of pts defined by + !> + !> [1] (wk1x,wk1y) + (wk2x,wk2y) - (wk3x,wk3y) - (wk4x,wk4y) = 0 + !> + !> [2] sig1 + sig2 - sig3 - sig4 = 0 + !> + !> Because of the influence of depth, it is convenient to define new + !> vectors (wnx,wny) = (wknx*dep,wkny*dep) with magnitudes wn = + !> sqrt(wnx**2+wny**2) = wkn*dep such that a dimensionless frequency + !> is sign*sqrt(dep/g) = sqrt[wkn*dep*tanh(wkn*dep)] + !> = sqrt[wn*tanh(wn)]. + !> With these definitions and vectors (wk1x,wk1y) and (wk3x,wk3y) + !> given as input, we can write (with some rearrangement + !> of [1] and [2]) + !> + !> [3] w3x - w1x = px = w2x - w4x + !> [4] w3y - w1y = py = w2y - w4y + !> [5] sqrt[w3*tanh(w3)] - sqrt[w1*tanh(w1)] = q + !> = sqrt[w2*tanh(w2)] - sqrt[w4*tanh(w4)] + !> + !> With dimensionless vector (px,py) = (w3x-w1x,w3y-w1y) [magnitude + !> p = sqrt(px**2+py**2), direction atan2(py,px)] and dimensionless + !> frequency difference q = sqrt(w3*tanh(w3)] - sqrt(w1*tanh(w1)] + !> defined by input parameters, we see from [3] and [4] that + !> (w4x,w4y) = (w2x-px,w2y-py) [magnitude w4 = sqrt((w2x-px)**2 + + !> (w2y-py)**2)] and thus from [5] we must basically find elements + !> w2x and w2y that satisfy + !> + !> [6] sqrt[sqrt(w2x**2+w2y**2)*tanh(sqrt(w2x**2+w2y**2))] - + !> sqrt[sqrt((w2x-px)**2+(w2y-py)**2) * + !> tanh(sqrt((w2x-px)**2+(w2y-py)**2] = q + !> + !> The locus curve defined by the set of pts (w2x,w2y) crosses the + !> p-vector axis at two points; one with magnitude w2=rmin*p with + !> 0.5 < rmin < 1.0 and one with magnitude w2=rmax*p with rmax > 1. + !> We first isolate rmin, rmax using various iterative algorithms, + !> and then find locus pts that are not on the p-vector axis with + !> another iterative scheme. At the end, we un-normalize the w2 + !> and w4 vectors to find the wk2 and wk4 vectors. + !> @endverbatim + !> + !> @param[in] dep + !> @param[in] wk1x + !> @param[in] wk1y + !> @param[in] wk3x + !> @param[in] wk3y + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE shlocr ( dep, wk1x,wk1y, wk3x,wk3y ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !/ @@ -2625,6 +2832,45 @@ END SUBROUTINE shlocr !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + !> + !> @brief Calculates four-wave Boltzmann coupling coefficient in shallow + !> water. + !> + !> @details Calculates four-wave Boltzmann coupling coefficient in shallow + !> water given k1,k2,k3 and following at least Hasselmann (1962) + !> and probably Herterich and Hasselmann (1982). Dimensional + !> wavenumbers are (wnx0,wny0), n = 1,3, h = depth, csq = coupling + !> coefficient. This is the same as Don's cplesh, except within + !> the algorithm, wavenumbers are made dimensionless with h and + !> frequencies with sqrt(h/g), g = gravitational acceleration (the + !> idea is to simplify and speed up the calculations while keeping + !> a reasonable machine resolution of the result). At the end, + !> dimensionless csqhat is redimensioned as csq = csqhat/(h**6) + !> so it is returned as a dimensional entity. + !> + !> This calculation can be a touchy bird, so we use double precision + !> for internal calculations, using single precision for input and + !> output. + !> + !> @param[in] w1x0 + !> @param[in] w1y0 + !> @param[in] w2x0 + !> @param[in] w2y0 + !> @param[in] w3x0 + !> @param[in] w3y0 + !> @param[in] h + !> @param[out] csq + !> @param[in] irng + !> @param[in] krng + !> @param[in] kang + !> @param[in] ipt + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE cplshr ( w1x0,w1y0, w2x0,w2y0, w3x0,w3y0, & h, csq, irng,krng, kang,ipt ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -2986,6 +3232,29 @@ END SUBROUTINE cplshr !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + !> + !> @brief Splits the Action Density into two parts. + !> + !> @verbatim + !> (1) large-scale part dens1(nrng,nang) and + !> (2) small-scale part dens2(nrng,nang) + !> dens1 & dens2 in Polar Action Density (k,theta) space Norm. (in k) + !> @endverbatim + !> + !> @param[in] nrmn Number of first frequency bin in [1,nrng-1]. + !> @param[in] nrmx Number of last frequency bin in [2,nrng]. + !> @param[in] npk Number of peak frequency in [2,nrng-1]. + !> @param[in] fpk Peak frequency [Hz] of initial frequency spectrum. + !> @param[in] nbins number of bins (see more below). + !> @param[in] wka Wavenumbers array [1/m]. + !> @param[in] cga Group velocities array [m/s]. + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE optsa2 ( nrmn,nrmx, npk,fpk, nbins, wka, cga ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !/ @@ -3467,6 +3736,24 @@ END SUBROUTINE optsa2 !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !> + !> @brief For a given Action Density array dens1(k,theta), computes the + !> rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to + !> wave-wave interaction. + !> + !> @details For a given Action Density array dens1(k,theta), computes the + !> rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to + !> wave-wave interaction, as well as some ancillary arrays relating to + !> positive and negative fluxes and their integrals. + !> + !> @param[in] pha Pha = k*dk*dtheta. + !> @param[in] ialt Integer switch. + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE snlr_fbi ( pha, ialt ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !/ @@ -4081,6 +4368,25 @@ END SUBROUTINE snlr_fbi !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + !> + !> @brief For a given Action Density array dens1(k,theta), computes the + !> rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to + !> wave-wave interaction. + !> + !> @details For a given Action Density array dens1(k,theta), computes the + !> rate-of-change array sumint(k,theta) = dN(k,theta)/dt owing to + !> wave-wave interaction, as well as some ancillary arrays + !> relating to positive and negative fluxes and their integrals. + !> + !> @param[in] pha Pha = k*dk*dtheta. + !> @param[in] ialt Integer switch. + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE snlr_tsa ( pha, ialt ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !/ @@ -4695,6 +5001,19 @@ END SUBROUTINE snlr_tsa !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !> + !> @brief Interpolate bi-linearly to fill in tsa/fbi & diag/diag2 arrays. + !> + !> @details Optionally smooth the interior and the corners + !> after alternating the irng, iang, krng & kang loops in snlr's. + !> + !> @param[inout] X Array to be interpolated and smoothed. + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> SUBROUTINE interp2 ( X ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !/ @@ -4954,6 +5273,28 @@ END SUBROUTINE interp2 !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !> + !> @brief Calculate the Wavenumber k (rad/m) as function of + !> frequency 'f' (Hz) and depth 'dep' (m). + !> + !> @verbatim + !> Using what looks like a "Pade approximation" of an inversion + !> of the linear wave dispersion relation. + !> + !> sigma^2 = gk*tanh(kd), sigma = 2*pi*f + !> Wavenumber k (rad/m) is returned in "wkfnc" + !> + !> @endverbatim + !> + !> @param f Frequency (Hz). + !> @param dep Depth (m). + !> @returns wkfnc Wavenumber 'k' + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> REAL FUNCTION wkfnc ( f, dep ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !/ @@ -5042,6 +5383,28 @@ END FUNCTION wkfnc !!============================================================================== !! !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !> + !> @brief Calculate the Group velocity Cg (m/s) as function of + !> frequency 'f' (Hz), depth 'dep' (m) and phase speed 'cvel' (m/s). + !> + !> @verbatim + !> This routine uses the identity + !> sinh(2x) = 2*tanh(x)/(1-tanh(x)**2) + !> to avoid extreme sinh(2x) for large x. + !> thus, 2kd/sinh(2kd) = kd(1-tanh(kd)**2)/tanh(kd) + !> Group velocity Cg (m/s) is returned in "cgfnc". + !> @endverbatim + !> + !> @param f Frequency (Hz). + !> @param dep Depth (m). + !> @param cvel Phase speed (m/s). + !> @returns cgfnc Group velocity (m/s). + !> + !> @author Bash Toulany + !> @author Michael Casey + !> @author William Perrie + !> @date 12-Apr-2016 + !> REAL FUNCTION cgfnc ( f, dep, cvel ) !! -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !/ From a23d1c2daef73756534708a95e83148be51731dc Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Tue, 6 Dec 2022 11:07:04 -0500 Subject: [PATCH 05/15] w3snl5md.F90: add doxygen markup --- model/src/w3snl5md.F90 | 84 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) diff --git a/model/src/w3snl5md.F90 b/model/src/w3snl5md.F90 index 3f96c78f1..e48691378 100644 --- a/model/src/w3snl5md.F90 +++ b/model/src/w3snl5md.F90 @@ -1,5 +1,27 @@ +!> @file +!> @brief Interface module for GKE (resonant & quasi-resonant four-wave +!> interactions). +!> +!> @author O. Gramstad +!> @author Q. Liu +!> @date 07-Jun-2021 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Interface module for GKE (resonant & quasi-resonant four-wave +!> interactions). +!> +!> @author O. Gramstad +!> @author Q. Liu +!> @date 07-Jun-2021 +!> +!> @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. +!> MODULE W3SNL5MD !/ !/ +-----------------------------------+ @@ -54,6 +76,25 @@ MODULE W3SNL5MD ! CONTAINS !/ ------------------------------------------------------------------- / + !> + !> @brief Interface to CalcQRSNL subroutine of the GKE module. + !> + !> @param[in] A + !> @param[in] CG + !> @param[in] WN + !> @param[in] FMEAN + !> @param[in] T1ABS + !> @param[in] U10 + !> @param[in] UDIR + !> @param[in] JSEA + !> @param[out] S + !> @param[out] D + !> @param[out] KURT + !> + !> @author O. Gramstad + !> @author Q. Liu + !> @date 24-Apr-2019 + !> SUBROUTINE W3SNL5(A, CG, WN, FMEAN, T1ABS, U10, UDIR, JSEA, & S, D, KURT) !/ @@ -382,6 +423,14 @@ SUBROUTINE W3SNL5(A, CG, WN, FMEAN, T1ABS, U10, UDIR, JSEA, & !/ END SUBROUTINE W3SNL5 !/ ------------------------------------------------------------------- / + + !> + !> @brief Initialization for the GKE module (Prepare wavenumber grid & kernel + !> coefficients). + !> + !> @author Q. Liu + !> @date 27-Feb-2019 + !> SUBROUTINE INSNL5 !/ !/ +-----------------------------------+ @@ -493,6 +542,23 @@ SUBROUTINE INSNL5 !/ END SUBROUTINE INSNL5 !/ ------------------------------------------------------------------- / + !> + !> @brief Estimate the dominant wave breaking probability. + !> + !> @details Based on the empirical parameterization proposed by + !> Babanin et al. (2001). + !> + !> @param A + !> @param CG + !> @param WN + !> @param DPT + !> @param U10 + !> @param UDIR + !> @returns CALC_WBTv2 + !> + !> @author Q. Liu + !> @date 24-Apr-2019 + !> FUNCTION CALC_WBTv2 (A, CG, WN, DPT, U10, UDIR) !/ !/ +-----------------------------------+ @@ -647,6 +713,12 @@ FUNCTION CALC_WBTv2 (A, CG, WN, DPT, U10, UDIR) END FUNCTION CALC_WBTv2 !/ ------------------------------------------------------------------- / + !> + !> @brief Initialization for point output (Snl). + !> + !> @author Q. Liu + !> @date 25-Mar-2019 + !> SUBROUTINE INPOUT !/ !/ +-----------------------------------+ @@ -799,6 +871,18 @@ SUBROUTINE INPOUT ! END SUBROUTINE INPOUT !/ ------------------------------------------------------------------- / + + !> + !> @brief Check if the 2D array `ARR2D` contains NaN. + !> + !> @param NK + !> @param NTH + !> @param ARR2D + !> @returns HasNaN + !> + !> @author Q. Liu + !> @date 25-Apr-2019 + !> FUNCTION HasNaN(NK, NTH, ARR2D) !/ !/ +-----------------------------------+ From 4f4d471f727e6215f20e6b3d1b666942a38b4d5c Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Tue, 6 Dec 2022 11:24:43 -0500 Subject: [PATCH 06/15] w3snlsmd.F90: add doxygen markup --- model/src/w3snlsmd.F90 | 60 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/model/src/w3snlsmd.F90 b/model/src/w3snlsmd.F90 index 576980553..df7801793 100644 --- a/model/src/w3snlsmd.F90 +++ b/model/src/w3snlsmd.F90 @@ -1,5 +1,23 @@ +!> @file +!> @brief Nonlinear interaction based `smoother' for high frequencies. +!> +!> @author H. L. Tolman +!> @date 13-Jul-2012 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Nonlinear interaction based `smoother' for high frequencies. +!> +!> @author H. L. Tolman +!> @date 13-Jul-2012 +!> +!> @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. +!> MODULE W3SNLSMD !/ !/ +-----------------------------------+ @@ -86,6 +104,30 @@ MODULE W3SNLSMD !/ CONTAINS !/ ------------------------------------------------------------------- / + !> + !> @brief High-frequeny filter based on the nonlinear interactions for + !> an uresolved quadruplet. + !> + !> @details Compute interactions for a quadruplet that is not resolved + !> by the discrete spectral rsolution, and then reduces to a simple + !> five-point stencil. Furthermore interactions are filtered by + !> frequency to allow for high-frequency impact only, and the + !> integration schem is embedded, and reduces to a filter technique + !> for large time steps or strong interactions. + !> + !> @param[in] A Action spectrum A(ITH,IK) as a function of + !> direction (rad) and wavenumber. + !> @param[in] CG Group velocities (dimension NK). + !> @param[in] WN Wavenumbers (dimension NK). + !> @param[in] DEPTH Water depth in meters. + !> @param[in] UABS Wind speed (m/s). + !> @param[in] DT Numerical time step (s). + !> @param[out] SNL Nonlinear source term. + !> @param[out] AA Averaged spectrum. + !> + !> @author H. L. Tolman + !> @date 04-Aug-2008 + !> SUBROUTINE W3SNLS ( A, CG, WN, DEPTH, UABS, DT, SNL, AA ) !/ !/ +-----------------------------------+ @@ -543,6 +585,15 @@ SUBROUTINE W3SNLS ( A, CG, WN, DEPTH, UABS, DT, SNL, AA ) !/ CONTAINS !/ ------------------------------------------------------------------- / + !> + !> @brief Expand spectrum to simplify indirect addressing. + !> + !> @param[out] PSPC Expanded spectrum. + !> @param[out] SPEC Expanded spectrum. + !> + !> @author H. L. Tolman + !> @date 23-Jul-2008 + !> SUBROUTINE EXPAND ( PSPC, SPEC ) !/ !/ +-----------------------------------+ @@ -603,6 +654,15 @@ END SUBROUTINE EXPAND !/ END SUBROUTINE W3SNLS !/ ------------------------------------------------------------------- / + !> + !> @brief Initializations for the Snl / filter source term for high + !> frequencies. + !> + !> @details Precompute weight functions and store in array. + !> + !> @author H. L. Tolman + !> @date 04-Aug-2008 + !> SUBROUTINE INSNLS !/ !/ +-----------------------------------+ From 121497463522ce03d844e123ac2636aa14ea5d6d Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Tue, 6 Dec 2022 12:03:43 -0500 Subject: [PATCH 07/15] w3sis2md.F90: add doxygen markup --- model/src/w3sis2md.F90 | 129 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 127 insertions(+), 2 deletions(-) diff --git a/model/src/w3sis2md.F90 b/model/src/w3sis2md.F90 index 0e01d36ca..3eba61f14 100644 --- a/model/src/w3sis2md.F90 +++ b/model/src/w3sis2md.F90 @@ -1,5 +1,35 @@ +!> @file +!> @brief Floe-size dependant scattering of waves in the marginal ice zone. +!> @details Based on tabulated scattering coefficients for a semi-infinite +!> ice sheet. +!> +!> @author F. Ardhuin +!> @author P. Nicot +!> @author C. Sevigny +!> @author G. Boutin +!> @date 21-Jan-2018 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Floe-size dependant scattering of waves in the marginal ice zone. +!> +!> @details based on tabulated scattering coefficients for a semi-infinite +!> ice sheet. See papers by Dumont et al. (JGR 2011) and Williams et al. +!> (OM 2013) combined with flexural dissipation and ice break-up. +!> +!> @author F. Ardhuin +!> @author P. Nicot +!> @author C. Sevigny +!> @author G. Boutin +!> @date 21-Jan-2018 +!> +!> @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. +!> MODULE W3SIS2MD !/ !/ +-----------------------------------+ @@ -61,6 +91,15 @@ MODULE W3SIS2MD !/ CONTAINS !/ ------------------------------------------------------------------- / + !> + !> @brief Fill tables used for scattering. + !> + !> @details Linear interpolation. + !> + !> @author P. Nicot + !> @author F. Ardhuin + !> @date 21-Jan-2018 + !> SUBROUTINE INSIS2 !/ !/ +-----------------------------------+ @@ -580,6 +619,36 @@ SUBROUTINE INSIS2 !/ END SUBROUTINE INSIS2 !/ ------------------------------------------------------------------- / + !> + !> @brief Wave scattering in the MIZ, adapted from Dumont et al. + !> + !> @details This scattering routine allows the estimation of the + !> maximum floe size and an estimate of the creep-induced dissipation. + !> For the scattering, it is based on the normal incidence results of + !> Kohout and Meylan which are provided in a table. + !> + !> @param[in] A Action density spectrum (1-D). + !> @param[in] DEPTH Water depth. + !> @param[in] CICE Sea ice concentration. + !> @param[in] ICEH Ice thickness. + !> @param[inout] ICEF Maximum floe size (updated). + !> @param[in] ICEDMAX Maximum floe size. + !> @param[in] IX Not used. + !> @param[in] IY Not used. + !> @param[out] S Source term (1-D version). + !> @param[out] D Diagonal part of scattering (1-D version). + !> @param[out] DISSIP Diagonal dissipation term (1-D version) + !> @param[in] WN Wave number. + !> @param[in] CG Group speed. + !> @param[in] WN_R Wave number in ice. + !> @param[in] CG_ICE Group speed in ice. + !> @param[out] R Ratio of energy to wave energy without ice. + !> + !> @author P. Nicot + !> @author F. Ardhuin + !> @author G. Boutin + !> @date 04-May-2016 + !> SUBROUTINE W3SIS2 (A, DEPTH, CICE, ICEH, ICEF, ICEDMAX, IX, IY, & S, D, DISSIP, WN, CG, WN_R, CG_ICE, R) !/ @@ -1016,9 +1085,17 @@ SUBROUTINE W3SIS2 (A, DEPTH, CICE, ICEH, ICEF, ICEDMAX, IX, IY, & !/ End of W3SIS2 ----------------------------------------------------- / !/ END SUBROUTINE W3SIS2 - - !/ ------------------------------------------------------------------- / + !> + !> @brief NA. + !> + !> @param[in] ICEH + !> @param[inout] WN_I + !> @param[inout] DAMPING + !> @param[inout] CG_I + !> + !> @author C. Sevigny + !> @date 27-Aug-2015 SUBROUTINE W3RPWNICE(ICEH,WN_I,DAMPING,CG_I) !/ !/ +-----------------------------------+ @@ -1113,6 +1190,18 @@ END SUBROUTINE W3RPWNICE !/ ------------------------------------------------------------------- / ! SUBROUTINE FINDROOTS_NR(FUNCD,X0,VISC_RP,WN_ROOT) +!> +!> @brief NA. +!> +!> @param[in] GUESS +!> @param[in] VISC_RP +!> @param[in] ICEH +!> @param[in] SIGMA +!> @param[inout] X +!> +!> @author C. Sevigny +!> @date 27-Aug-2015 +!> SUBROUTINE FINDROOTS_NR(GUESS,VISC_RP,X,ICEH,SIGMA) !/ !/ +-----------------------------------+ @@ -1214,6 +1303,18 @@ END SUBROUTINE FINDROOTS_NR !/ ------------------------------------------------------------------- / + !> + !> @brief Computes the mean flow size from the max floe size. + !> + !> @param ICEDMIN + !> @param ICEDMAX + !> @param FRAGILITY + !> @returns W3FSD_DAVE + !> + !> @author C. Sevigny + !> @author F. Ardhuin + !> @date 06-Nov-2015 + !> FUNCTION W3FSD_DAVE(ICEDMIN,ICEDMAX,FRAGILITY) !/ !/ +-----------------------------------+ @@ -1305,6 +1406,18 @@ FUNCTION W3FSD_DAVE(ICEDMIN,ICEDMAX,FRAGILITY) END FUNCTION W3FSD_DAVE ! !/ ------------------------------------------------------------------- / + !> + !> @brief NA. + !> + !> @param WN_GUESS + !> @param VISC_RP + !> @param ICEH + !> @param SIGMA + !> @returns FVAL + !> + !> @author C. Sevigny + !> @date 27-Aug-2015 + !> FUNCTION FUNCD_FVAL(WN_GUESS, VISC_RP,ICEH,SIGMA) RESULT(FVAL) !/ !/ +-----------------------------------+ @@ -1396,6 +1509,18 @@ END FUNCTION FUNCD_FVAL !/ ------------------------------------------------------------------- / ! FUNCTION FUNCD(WN_GUESS, VISC_RP,ICEH) + !> + !> @brief NA. + !> + !> @param WN_GUESS + !> @param VISC_RP + !> @param ICEH + !> @param SIGMA + !> @returns FDERIV + !> + !> @author C. Sevigny + !> @date 27-Aug-2015 + !> FUNCTION FUNCD_FDERIV(WN_GUESS, VISC_RP,ICEH,SIGMA) RESULT(FDERIV) !/ !/ +-----------------------------------+ From 8bc8dfb05ac1669d813a0737e69d7d3f7ec892e8 Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Tue, 6 Dec 2022 19:59:37 -0500 Subject: [PATCH 08/15] w3sis1md.F90: add doxygen markup --- model/src/w3sis1md.F90 | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/model/src/w3sis1md.F90 b/model/src/w3sis1md.F90 index 54e1b87c0..4c7cff92b 100644 --- a/model/src/w3sis1md.F90 +++ b/model/src/w3sis1md.F90 @@ -1,5 +1,24 @@ +!> @file +!> @brief Diffusion source term. +!> +!> @author S. Zieger +!> @date 20-Dec-2013 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / + +!> +!> @brief Diffusion source term. +!> +!> @author S. Zieger +!> @date 20-Dec-2013 +!> +!> @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. +!> MODULE W3SIS1MD !/ !/ +-----------------------------------+ @@ -42,6 +61,16 @@ MODULE W3SIS1MD !/ CONTAINS !/ ------------------------------------------------------------------- / + !> + !> @brief Spectral reflection due to ice. + !> + !> @param[in] A Action density spectrum (1-D). + !> @param[in] ICE Sea ice concentration. + !> @param[out] S Source term (1-D version). + !> + !> @author S. Zieger + !> @date + !> SUBROUTINE W3SIS1 (A, ICE, S) !/ !/ +-----------------------------------+ From 1663a882a12374ba446cec1f4019cf6653abdedc Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Tue, 6 Dec 2022 22:36:09 -0500 Subject: [PATCH 09/15] w3sic5md.F90: add doxygen markup --- model/src/w3sic5md.F90 | 179 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 179 insertions(+) diff --git a/model/src/w3sic5md.F90 b/model/src/w3sic5md.F90 index 2d7ad953f..851bfcfd9 100644 --- a/model/src/w3sic5md.F90 +++ b/model/src/w3sic5md.F90 @@ -1,5 +1,27 @@ +!> @file +!> @brief Calculate ice source term S_{ice} according to different ice +!> models. +!> +!> @author Q. Liu +!> @author E. Rogers +!> @date 19-May-2021 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Calculate ice source term S_{ice} according to different ice +!> models: +!> +!> @author Q. Liu +!> @author E. Rogers +!> @date 19-May-2021 +!> +!> @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. +!> MODULE W3SIC5MD !/ !/ +-----------------------------------+ @@ -124,6 +146,25 @@ MODULE W3SIC5MD CONTAINS !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Calculate ice source term S_{ice} according to 3 different sea ice + !> models. + !> + !> @details (Mosig et al. 2015, Meylan et al. 2018, Liu et al. 2020). + !> + !> @param[in] A Action density spectrum (1-D). + !> @param[in] DEPTH Local water depth. + !> @param[in] CG Group velocities. + !> @param[in] WN Wavenumbers. + !> @param[in] IX Grid index. + !> @param[in] IY Grid index. + !> @param[out] S Source term (1-D version). + !> @param[out] D Diagonal term of derivative (1-D version). + !> + !> @author Q. Liu + !> @author E. Rogers + !> @date 19-May-2021 + !> SUBROUTINE W3SIC5 (A, DEPTH, CG, WN, IX, IY, S, D) !/ !/ +-----------------------------------+ @@ -430,6 +471,25 @@ SUBROUTINE W3SIC5 (A, DEPTH, CG, WN, IX, IY, S, D) END SUBROUTINE W3SIC5 !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Calculation of complex wavenumber arrays for ice-coupled waves. + !> + !> @details Using the Fox-Squire dispersion relations to get (kr, ki) and + !> then get cg by cg = dσ / dk (here dk uses kr). + !> + !> @param[inout] WN_R The real part of the wave number. + !> @param[inout] WN_I The imaginary part of the wave number. + !> @param[inout] CG Group velocity (m s^{-1}). + !> @param[inout] HICE Thickness of ice (m). + !> @param[inout] IVISC Viscosity parameter of ice (m^2 s^{-1}). + !> @param[inout] RHOI The density of ice (kg m^{-3}). + !> @param[inout] ISMODG Effecitive shear modulus G of ice (Pa). + !> @param[inout] HWAT Water depth. + !> + !> @author Q. Liu + !> @author E. Rogers + !> @date 25-Apr-2017 + !> SUBROUTINE W3IC5WNCG(WN_R, WN_I, CG, HICE, IVISC, RHOI, ISMODG, & HWAT) !/ @@ -573,6 +633,22 @@ SUBROUTINE W3IC5WNCG(WN_R, WN_I, CG, HICE, IVISC, RHOI, ISMODG, & END SUBROUTINE W3IC5WNCG !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Calculate the complex wavenumber for waves in ice according to + !> three different sea ice models. + !> + !> @param[in] HICE Thickness of ice (m). + !> @param[in] IVISC Viscosity parameter of ice (m^2 s^{-1}). + !> @param[in] RHOI The density of ice (kg m^{-3}). + !> @param[in] ISMODG Effecitive shear modulus G of ice (Pa). + !> @param[in] HWAT Water depth. + !> @param[in] WT Wave period (s; 1/freq). + !> @param[out] WNR The real part of the wave number. + !> @param[out] WNI The imaginary part of the wave number. + !> + !> @author Q. Liu + !> @date 19-May-2021 + !> SUBROUTINE FSDISP(HICE, IVISC, RHOI, ISMODG, HWAT, WT, WNR, WNI) !/ !/ +-----------------------------------+ @@ -867,6 +943,16 @@ SUBROUTINE FSDISP(HICE, IVISC, RHOI, ISMODG, HWAT, WT, WNR, WNI) END SUBROUTINE FSDISP !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Reducing the sensitivity of eigenvalues to rounding errors during + !> the execution of some algorithms. + !> + !> @param[in] NMAT The size of one dimension of MATRIX. + !> @param[inout] MATRIX A matrix with the shape (NMAT, NMAT). + !> + !> @author Q. Liu + !> @date 15-Mar-2016 + !> SUBROUTINE BALANCING_MATRIX(NMAT, MATRIX) !/ !/ +-----------------------------------+ @@ -1023,6 +1109,17 @@ SUBROUTINE BALANCING_MATRIX(NMAT, MATRIX) END SUBROUTINE BALANCING_MATRIX !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Calculate the eigenvalues of a general matrix. + !> + !> @param[in] NMAT The size of one dimension of HMAT. + !> @param[inout] HMAT The Hessenberg-type matrix (NMAT, NMAT). + !> @param[out] EIGR The real part of the N eigenvalues. + !> @param[out] EIGI The imaginary part of the N eigenvalues. + !> + !> @author Q. Liu + !> @date 17-Mar-2016 + !> SUBROUTINE EIG_HQR (NMAT, HMAT, EIGR, EIGI) !/ !/ +-----------------------------------+ @@ -1323,6 +1420,18 @@ SUBROUTINE EIG_HQR (NMAT, HMAT, EIGR, EIGI) END SUBROUTINE EIG_HQR !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Find the roots of arbitrary polynomials through finding the + !> eigenvalues of companion matrix. + !> + !> @param[in] NPC The # of the Polynomial coefficients. + !> @param[in] PCVEC The 1d vector for the Polynomial coefficients. + !> @param[out] RTRL The real part of all of the roots shape: [NPC-1]. + !> @param[out] RTIM The real part of all of the roots shape: [NPC-1]. + !> + !> @author Q. Liu + !> @date 16-Mar-2016 + !> SUBROUTINE POLYROOTS(NPC, PCVEC, RTRL, RTIM) !/ !/ +-----------------------------------+ @@ -1477,6 +1586,19 @@ SUBROUTINE POLYROOTS(NPC, PCVEC, RTRL, RTIM) END SUBROUTINE POLYROOTS !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Calculate the corrected term in the Newton-Raphson root-finding + !> method (Must use double precision). + !> + !> @param K Complex wave number. + !> @param C1 C1 in FSDISP. + !> @param C2 C2 in FSDISP. + !> @param H Water depth. + !> @returns NR_CORR Newton-Raphson corrected term (DK). + !> + !> @author Q. Liu + !> @date 19-May-2021 + !> FUNCTION NR_CORR(K, C1, C2, H) !/ !/ +-----------------------------------+ @@ -1628,6 +1750,18 @@ FUNCTION NR_CORR(K, C1, C2, H) END FUNCTION NR_CORR !/ ------------------------------------------------------------------- / !/ + !> + !> @brief The iterative procedure of the Newton-Raphson method. + !> + !> @param C1 C1 in FS dipsersion relations. + !> @param C2 C2 in FS dipsersion relations. + !> @param H Water depth. + !> @param GUESS The first guess obtained from POLYROOTS. + !> @returns NR_ROOT The calculated complex wave number. + !> + !> @author Q. Liu + !> @date 19-May-2021 + !> FUNCTION NR_ROOT(C1, C2, H, GUESS) !/ !/ +-----------------------------------+ @@ -1790,6 +1924,19 @@ FUNCTION NR_ROOT(C1, C2, H, GUESS) END FUNCTION NR_ROOT !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Custom written sinh function for complex inputs. + !> + !> @details For a number of compilers, the built-in function sinh, + !> cosh and tanh do not support the complex inputs. So here I write an + !> external one. + !> + !> @param X A double-precision complex variable. + !> @returns CMPLX_SINH Complex sinh(x). + !> + !> @author Q. Liu + !> @date 24-Mar-2016 + !> FUNCTION CMPLX_SINH(X) !/ !/ +-----------------------------------+ @@ -1887,6 +2034,19 @@ FUNCTION CMPLX_SINH(X) END FUNCTION CMPLX_SINH !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Hand written cosh function for complex inputs. + !> + !> @details For a number of compilers, the built-in function sinh, + !> cosh and tanh do not support the complex inputs. So here I write an + !> external one. + !> + !> @param X A double-precision complex variable. + !> @returns CMPLX_COSH Complex cosh(x). + !> + !> @author Q. Liu + !> @date 24-Mar-2016 + !> FUNCTION CMPLX_COSH(X) !/ !/ +-----------------------------------+ @@ -1984,6 +2144,19 @@ FUNCTION CMPLX_COSH(X) END FUNCTION CMPLX_COSH !/ ------------------------------------------------------------------- / !/ + !> + !> @brief An alternate version of tanh to avoid overflow error. + !> + !> @details We may encounter overflow error for the above tanh + !> function as kh becomes huge. This is another version of tanh + !> function. + !> + !> @param X A double-precision complex variable. + !> @returns CMPLX_TANH2 Complex tanh(x). + !> + !> @author Q. Liu + !> @date 24-Mar-2016 + !> FUNCTION CMPLX_TANH2(X) !/ !/ +-----------------------------------+ @@ -2084,6 +2257,12 @@ END FUNCTION CMPLX_TANH2 !/ !/ ------------------------------------------------------------------- / !/ + !> + !> @brief Initialize the random seed based on the system's time. + !> + !> @author Q. Liu + !> @date 24-Mar-2016 + !> SUBROUTINE INIT_RANDOM_SEED() !/ !/ +-----------------------------------+ From a0890b868091cbd7ff75eae60528d066fb7d90a0 Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Wed, 7 Dec 2022 11:36:03 -0500 Subject: [PATCH 10/15] w3sic4md.F90: add doxygen markup --- model/src/w3sic4md.F90 | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/model/src/w3sic4md.F90 b/model/src/w3sic4md.F90 index 9c7cdd27b..7b1c9c67a 100644 --- a/model/src/w3sic4md.F90 +++ b/model/src/w3sic4md.F90 @@ -1,5 +1,29 @@ +!> @file +!> @brief Calculate ice source term S_{ice} according to simple methods. +!> +!> @author C. Collins +!> @author E. Rogers +!> @date 21-Jan-2015 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Calculate ice source term S_{ice} according to simple methods. +!> +!> @details 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 +!> @author 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. +!> MODULE W3SIC4MD !/ !/ +-----------------------------------+ @@ -70,6 +94,21 @@ MODULE W3SIC4MD !/ CONTAINS !/ ------------------------------------------------------------------- / + !> + !> @brief S_{ice} source term using 5 parameters read from input files. + !> + !> @param[in] A Action density spectrum (1-D). + !> @param[in] DEPTH Local water depth. + !> @param[in] CG Group velocities. + !> @param[in] IX Grid index. + !> @param[in] IY Grid index. + !> @param[out] S Source term (1-D version). + !> @param[out] D Diagonal term of derivative (1-D version). + !> + !> @author C. Collins + !> @author E. Rogers + !> @date 24-Feb-2017 + !> SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D) !/ !/ +-----------------------------------+ From 29690403702eb21df95a966d24c5ab9c6f2455f9 Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Wed, 7 Dec 2022 14:19:36 -0500 Subject: [PATCH 11/15] w3sic2md.F90: add doxygen markup --- model/src/w3sic2md.F90 | 48 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/model/src/w3sic2md.F90 b/model/src/w3sic2md.F90 index d15ae7abf..b3b050496 100644 --- a/model/src/w3sic2md.F90 +++ b/model/src/w3sic2md.F90 @@ -1,5 +1,29 @@ +!> @file +!> @brief Calculate ice dissipation source term S_{ice}. +!> +!> @author E. Rogers +!> @author S. Zieger +!> @author F. Ardhuin +!> @author G. Boutin +!> @date 05-Jan-2018 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Calculate ice dissipation source term S_{ice}. +!> +!> @author E. Rogers +!> @author S. Zieger +!> @author F. Ardhuin +!> @author G. Boutin +!> @date 05-Jan-2018 +!> +!> @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. +!> MODULE W3SIC2MD !/ !/ +-----------------------------------+ @@ -69,6 +93,30 @@ MODULE W3SIC2MD !/ CONTAINS !/ ------------------------------------------------------------------- / + !> + !> @brief S_{ice} source term using 5 parameters read from input files. + !> + !> @param[in] A Action density spectrum (1-D). + !> @param[in] DEPTH Local water depth. + !> @param[in] ICEH Ice thickness. + !> @param[in] ICEF Ice Floe diameter. + !> @param[in] CG Group velocities. + !> @param[in] WN Wavenumbers. + !> @param[in] IX Grid index. + !> @param[in] IY Grid index. + !> @param[out] S Source term (1-D version). + !> @param[out] D Diagonal term of derivative (1-D version). + !> @param[in] WN_R Wavenumbers in ice. + !> @param[in] CG_ICE Group velocities in ice. + !> @param[in] ALPHA Exponential decay rate of energy. + !> @param[in] R Ratio of energy to wave energy without ice. + !> + !> @author E. Rogers + !> @author S. Zieger + !> @author F. Ardhuin + !> @author G. Boutin + !> @date 04-Jan-2018 + !> SUBROUTINE W3SIC2 (A, DEPTH, ICEH, ICEF, CG, WN, IX, IY, S, D, WN_R, & CG_ICE, ALPHA, R) !/ From 7f8a574ce90c3060bcecf383aa4f24ce3641f44a Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Wed, 7 Dec 2022 14:27:17 -0500 Subject: [PATCH 12/15] w3sic1md.F90: add doxygen markup --- model/src/w3sic1md.F90 | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/model/src/w3sic1md.F90 b/model/src/w3sic1md.F90 index 16f7b1693..acd749696 100644 --- a/model/src/w3sic1md.F90 +++ b/model/src/w3sic1md.F90 @@ -1,5 +1,25 @@ +!> @file +!> @brief Calculate ice source term S_{ice} according to simple methods. +!> +!> @author E. Rogers +!> @author S. Zieger +!> @date 11-Oct-2013 +!> + #include "w3macros.h" !/ ------------------------------------------------------------------- / +!> +!> @brief Calculate ice source term S_{ice} according to simple methods. +!> +!> @author E. Rogers +!> @author S. Zieger +!> @date 11-Oct-2013 +!> +!> @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. +!> MODULE W3SIC1MD !/ !/ +-----------------------------------+ @@ -55,6 +75,21 @@ MODULE W3SIC1MD !/ CONTAINS !/ ------------------------------------------------------------------- / + !> + !> @brief S_{ice} source term using 5 parameters read from input files. + !> + !> @param[in] A Action density spectrum (1-D). + !> @param[in] DEPTH Local water depth. + !> @param[in] CG Group velocities. + !> @param[in] IX Grid index. + !> @param[in] IY Grid index. + !> @param[out] S Source term (1-D version). + !> @param[out] D Diagonal term of derivative (1-D version). + !> + !> @author E. Rogers + !> @author S. Zieger + !> @date 11-Oct-2013 + !> SUBROUTINE W3SIC1 (A, DEPTH, CG, IX, IY, S, D) !/ !/ +-----------------------------------+ From ca1a841f4c1f00c8e77aca20c07dc1f75303dcc0 Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Wed, 7 Dec 2022 14:35:26 -0500 Subject: [PATCH 13/15] w3snl1md.F90: grammer fix --- model/src/w3snl1md.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/src/w3snl1md.F90 b/model/src/w3snl1md.F90 index 6db0be264..e21349ede 100644 --- a/model/src/w3snl1md.F90 +++ b/model/src/w3snl1md.F90 @@ -9,7 +9,7 @@ #include "w3macros.h" !/ ------------------------------------------------------------------- / !> -!> @brief Bundles routines calculate nonlinear wave-wave interactions +!> @brief Bundles routines to calculate nonlinear wave-wave interactions !> according to the Discrete Interaction Approximation (DIA) of !> Hasselmann et al. (JPO, 1985). !> From b20757256f434003eb474c7e2c95f74dec4a2ff8 Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Wed, 7 Dec 2022 14:40:59 -0500 Subject: [PATCH 14/15] w3snl3md.F90: spelling --- model/src/w3snl3md.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/src/w3snl3md.F90 b/model/src/w3snl3md.F90 index df12bec39..9869a60f6 100644 --- a/model/src/w3snl3md.F90 +++ b/model/src/w3snl3md.F90 @@ -774,7 +774,7 @@ END SUBROUTINE W3SNL3 !> !> @brief Initialization for generalized multiple DIA routine. !> - !> @details Fill storage aryays as described in the main subroutine + !> @details Fill storage arrays as described in the main subroutine !> with interpolation, model and distribution data. !> !> @author H. L. Tolman From e74f356975de9742699896e2ce391d6eff339bab Mon Sep 17 00:00:00 2001 From: Matt Masarik Date: Wed, 14 Dec 2022 11:01:13 -0500 Subject: [PATCH 15/15] w3snl4md.F90: remove stray '\!' characters --- model/src/w3snl4md.F90 | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/model/src/w3snl4md.F90 b/model/src/w3snl4md.F90 index aa9788b15..012e65ba4 100644 --- a/model/src/w3snl4md.F90 +++ b/model/src/w3snl4md.F90 @@ -1446,26 +1446,26 @@ END SUBROUTINE W3SNL4 !> !> @brief This routine sets up the geometric part of the Boltzmann integral. !> - !> @details This routine sets up the geometric part of the Boltzmann integral ! - !> based on a grid of wave frequencies and directions, with wave- ! - !> numbers related to frequency and depth by linear dispersion. It ! - !> is adapted from Don's original code with changes to modify the ! - !> indexing so there are fewer unused elements, and a number of algo- ! - !> rithmic changes that are mathematically equivalent to Don's but ! - !> take advantage of intrinsic functions to form smooth results with ! - !> less reliance on if statements. ! - !> ! - !> It calls locus-solving routines shloxr and shlocr and coupling ! - !> coefficient routine cplshr. If shlocr does not converge, ierr_gr ! - !> will be something other than 0 and the routine will terminate, ! - !> returning ierr_gr to the calling program (see shlocr). ! - !> ! - !> It returns array grad(,,), which is an estimate of the product ! - !> C(k1,k2,k3,k4)*H(k1,k3,k4)*ds/|dW/dn| (where n and the k's are all ! - !> vectors) as given, for example, by Eq.(7) of 'Nonlinear energy ! - !> fluxes and the finite depth equilibrium range in wave spectra,' ! - !> by Resio, Pihl, Tracy and Vincent (2001, JGR, 106(C4), p. 6985), ! - !> as well as arrays for indexing, interpolating and weighting locus- ! + !> @details This routine sets up the geometric part of the Boltzmann integral + !> based on a grid of wave frequencies and directions, with wave- + !> numbers related to frequency and depth by linear dispersion. It + !> is adapted from Don's original code with changes to modify the + !> indexing so there are fewer unused elements, and a number of algo- + !> rithmic changes that are mathematically equivalent to Don's but + !> take advantage of intrinsic functions to form smooth results with + !> less reliance on if statements. + !> + !> It calls locus-solving routines shloxr and shlocr and coupling + !> coefficient routine cplshr. If shlocr does not converge, ierr_gr + !> will be something other than 0 and the routine will terminate, + !> returning ierr_gr to the calling program (see shlocr). + !> + !> It returns array grad(,,), which is an estimate of the product + !> C(k1,k2,k3,k4)*H(k1,k3,k4)*ds/|dW/dn| (where n and the k's are all + !> vectors) as given, for example, by Eq.(7) of 'Nonlinear energy + !> fluxes and the finite depth equilibrium range in wave spectra,' + !> by Resio, Pihl, Tracy and Vincent (2001, JGR, 106(C4), p. 6985), + !> as well as arrays for indexing, interpolating and weighting locus- !> based wavenumber vectors within the discrete solution grid. !> !> @param[in] dep Depth (m)