!> @file
!> @brief Observation-based wind input and dissipation after Donelan et al (2006),
!> and Babanin et al. (2010).
!>
!> @author S. Zieger
!> @author Q. Liu
!> @date 26-Jun-2018
!>
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
!> @brief Observation-based wind input and dissipation after Donelan et al (2006),
!> and Babanin et al. (2010).
!>
!> @details Parameterisation is based on the field
!> data from Lake George, Australia. Initial implementation of input
!> and dissipation is based on work from Tsagareli et al. (2010) and
!> Rogers et al. (2012). Parameterisation extended and account for
!> negative input due to opposing winds (see Donelan et al, 2006) and
!> the vector version of the stress computation.
!>
!> @author S. Zieger
!> @author Q. Liu
!> @date 26-Jun-2018
!>
MODULE W3SRC6MD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP/NOPP |
!/ | S. Zieger |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 26-Jun-2018 |
!/ +-----------------------------------+
!/
!/ 29-May-2009 : Origination (w3srcxmd.ftn) ( version 3.14 )
!/ 10-Feb-2011 : Implementation of source terms ( version 4.04 )
!/ (S. Zieger)
!/ 26-Jun-2017 : Recalibration of ST6 ( verison 6.06 )
!/ (Q. Liu )
!/
!/ 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 :
!
! Observation-based wind input and dissipation after Donelan et al (2006),
! and Babanin et al. (2010). Parameterisation is based on the field
! data from Lake George, Australia. Initial implementation of input
! and dissipation is based on work from Tsagareli et al. (2010) and
! Rogers et al. (2012). Parameterisation extended and account for
! negative input due to opposing winds (see Donelan et al, 2006) and
! the vector version of the stress computation.
!
! References:
! Babanin et al. 2010: JPO 40(4) 667- 683
! Donelan et al. 2006: JPO 36(8) 1672-1689
! Tsagareli et al. 2010: JPO 40(4) 656- 666
! Rogers et al. 2012: JTECH 29(9) 1329-1346
!
! 2. Variables and types :
!
! Not applicable.
!
! 3. Subroutines and functions :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! W3SPR6 Subr. Public Integral parameter calculation following !/ST1.
! W3SIN6 Subr. Public Observation-based wind input.
! W3SDS6 Subr. Public Observation-based dissipation.
!
! IRANGE Func. Private Generate a sequence of integer values.
! LFACTOR Func. Private Calculate reduction factor for Sin.
! TAUWINDS Func. Private Normal stress calculation for Sin.
! ----------------------------------------------------------------
!
! 4. Subroutines and functions used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Remarks :
!
! 6. Switches :
!
! !/S Enable subroutine tracing.
! !/T6 Enable test output for wind input and dissipation subroutines.
!
! 7. Source code :
!/
!/ ------------------------------------------------------------------- /
!/
PUBLIC :: W3SPR6, W3SIN6, W3SDS6
PRIVATE :: LFACTOR, TAUWINDS, IRANGE
CONTAINS
!/ ------------------------------------------------------------------- /
!>
!> @brief Calculate mean wave parameters.
!>
!> @details See source term routines.
!>
!> @param[inout] A Action as a function of direction and wavenumber.
!> @param[inout] CG Group velocities.
!> @param[inout] WN Wavenumbers.
!> @param[inout] EMEAN Mean wave energy.
!> @param[inout] FMEAN Mean wave frequency.
!> @param[inout] WNMEAN Mean wavenumber.
!> @param[inout] AMAX Maximum action density in spectrum.
!> @param[inout] FP Peak frequency (rad).
!>
!> @author S. Zieger
!> @date 11-Feb-2011
!>
SUBROUTINE W3SPR6 (A, CG, WN, EMEAN, FMEAN, WNMEAN, AMAX, FP)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP/NOPP |
!/ | S. Zieger |
!/ | FORTRAN 90 |
!/ | Last update : 11-Feb-2011 |
!/ +-----------------------------------+
!/
!/ 08-Oct-2007 : Origination. ( version 3.13 )
!/ 11-Feb-2011 : Implementation based on W3SPR1 ( version 4.04 )
!/ (S. Zieger)
!/
! 1. Purpose :
! Calculate mean wave parameters.
!
! 2. Method :
! See source term routines.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! A R.A. I Action as a function of direction and wavenumber
! CG R.A. I Group velocities
! WN R.A. I Wavenumbers
! EMEAN REAL O Mean wave energy
! FMEAN REAL O Mean wave frequency
! WNMEAN REAL O Mean wavenumber
! AMAX REAL O Maximum action density in spectrum
! FP REAL O Peak frequency (rad)
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3SRCE Subr. W3SRCEMD Source term integration.
! W3EXPO Subr. N/A Point output post-processor.
! GXEXPO Subr. N/A GrADS point output post-processor.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS, ONLY: TPIINV
USE W3GDATMD, ONLY: NK, NTH, SIG, DTH, DDEN, FTE, FTF, FTWN, DSII
USE W3ODATMD, ONLY: NDST, NDSE
USE W3SERVMD, ONLY: EXTCDE
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(IN) :: A(NTH,NK), CG(NK), WN(NK)
REAL, INTENT(OUT) :: EMEAN, FMEAN, WNMEAN, AMAX, FP
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
INTEGER :: IMAX
REAL :: EB(NK), EBAND
REAL, PARAMETER :: HSMIN = 0.05
REAL :: COEFF(3)
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'W3SPR6')
#endif
!
!
! 1. Integrate over directions -------------------------------------- /
EB = SUM(A,1) * DDEN / CG
AMAX = MAXVAL(A)
!
! 2. Integrate over wavenumbers ------------------------------------- /
EMEAN = SUM(EB)
FMEAN = SUM(EB / SIG(1:NK))
WNMEAN = SUM(EB / SQRT(WN))
!
! 3. Add tail beyond discrete spectrum and get mean pars ------------ /
! ( DTH * SIG absorbed in FTxx )
EBAND = EB(NK) / DDEN(NK)
EMEAN = EMEAN + EBAND * FTE
FMEAN = FMEAN + EBAND * FTF
WNMEAN = WNMEAN + EBAND * FTWN
!
! 4. Final processing
FMEAN = TPIINV * EMEAN / MAX(1.0E-7, FMEAN)
WNMEAN = ( EMEAN / MAX(1.0E-7,WNMEAN) )**2
!
! 5. Determine peak frequency using a weighted integral ------------- /
! Young (1999) p239: integrate f F**4 df / integrate F**4 df ----- /
! TODO: keep in mind that **fp** calculated in this way may not
! work under mixing (wind-sea and swell) sea states (QL)
FP = 0.0
!
IF (4.0*SQRT(EMEAN) .GT. HSMIN) THEN
EB = SUM(A,1) * SIG(1:NK) /CG * DTH
FP = SUM(SIG(1:NK) * EB**4 * DSII) / MAX(1E-10, SUM(EB**4 * DSII))
FP = FP * TPIINV
END IF
!
RETURN
!/
!/ End of W3SPR6 ----------------------------------------------------- /
!/
END SUBROUTINE W3SPR6
!/ ------------------------------------------------------------------- /
!>
!> @brief Observation-based source term for wind input.
!>
!> @details After Donelan, Babanin, Young and Banner (Donelan et al ,2006)
!> following the implementation by Rogers et al. (2012).
!>
!> References:
!> Donelan et al. 2006: JPO 36(8) 1672-1689.
!> Rogers et al. 2012: JTECH 29(9) 1329-1346
!>
!>
!> @param[in] A Action density spectrum
!> @param[in] CG Group velocities.
!> @param[in] WN2 Wavenumbers.
!> @param[in] UABS Wind speed at 10 m above sea level (U10).
!> @param[in] USTAR Friction velocity.
!> @param[in] USDIR Direction of USTAR.
!> @param[in] CD Drag coefficient.
!> @param[in] DAIR Air density.
!> @param[out] TAUWX Component of the wave-supported stress.
!> @param[out] TAUWY Component of the wave-supported stress.
!> @param[out] TAUWNX Component of the negative part of the stress.
!> @param[out] TAUWNY Component of the negative part of the stress.
!> @param[out] S Source term.
!> @param[out] D Diagonal term of derivative.
!>
!> @author S. Zieger
!> @author Q. Liu
!> @date 13-Aug-2021
!>
SUBROUTINE W3SIN6 (A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, &
TAUWX, TAUWY, TAUNWX, TAUNWY, S, D )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP/NOPP |
!/ | S. Zieger |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 13-Aug-2021 |
!/ +-----------------------------------+
!/
!/ 20-Dec-2010 : Origination. ( version 4.04 )
!/ (S. Zieger)
!/
!/ 26-Jun-2018 : UPROXY Update & UABS ( version 6.06 )
!/ (Q. Liu)
!/ 13-Aug-2021 : Consider DAIR a variable ( version x.xx )
!/
! 1. Purpose :
!
! Observation-based source term for wind input after Donelan, Babanin,
! Young and Banner (Donelan et al ,2006) following the implementation
! by Rogers et al. (2012).
!
! References:
! Donelan et al. 2006: JPO 36(8) 1672-1689.
! Rogers et al. 2012: JTECH 29(9) 1329-1346
!
! 2. Method :
!
! Sin = B * E
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! A¹ R.A. I Action density spectrum
! CG R.A. I Group velocities
! WN2¹ R.A. I Wavenumbers
! UABS Real I Wind speed at 10 m above sea level (U10)
! USTAR Real I Friction velocity
! USDIR Real I Direction of USTAR
! CD Real I Drag coefficient
! DAIR Real I Air density
! S¹ R.A. O Source term
! D¹ R.A. O Diagonal term of derivative
! TAUWX-Y Real O Component of the wave-supported stress
! TAUNWX-Y Real O Component of the negative part of the stress
! ¹ Stored as 1-D array with dimension NTH*NK (column by column).
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! LFACTOR Subr. W3SRC6MD
! IRANGE Func. W3SRC6MD
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3SRCE Subr. W3SRCEMD Source term integration.
! W3EXPO Subr. N/A Point output post-processor.
! GXEXPO Subr. N/A GrADS point output post-processor.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See comments in source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
! !/T6 Test and diagnostic output for tail reduction.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS, ONLY: DWAT, TPI, GRAV
USE W3GDATMD, ONLY: NK, NTH, NSPEC, DTH, SIG2, DDEN2
USE W3GDATMD, ONLY: ECOS, ESIN, SIN6A0, SIN6WS
USE W3ODATMD, ONLY: NDSE
USE W3SERVMD, ONLY: EXTCDE
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
REAL, INTENT(IN) :: A (NSPEC), CG(NK), WN2(NSPEC)
REAL, INTENT(IN) :: UABS, USTAR, USDIR, CD, DAIR
REAL, INTENT(OUT) :: TAUWX, TAUWY, TAUNWX, TAUNWY
REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC)
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
INTEGER :: IK, ITH, IKN(NK)
REAL :: COSU, SINU, UPROXY
REAL, DIMENSION(NSPEC) :: CG2, ECOS2, ESIN2, DSII2
REAL, DIMENSION(NK) :: DSII, SIG, WN
REAL :: K(NTH,NK), SDENSIG(NTH,NK) ! 1,2,5)
REAL, DIMENSION(NK) :: ADENSIG, KMAX, ANAR, SQRTBN ! 1,2,3)
REAL, DIMENSION(NSPEC) :: W1, W2, SQRTBN2, CINV2 ! 4,7)
REAL, DIMENSION(NK) :: LFACT, CINV ! 5)
!/ ------------------------------------------------------------------- /
#ifdef W3_S
CALL STRACE (IENT, 'W3SIN6')
#endif
!
!/ 0) --- set up a basic variables ----------------------------------- /
COSU = COS(USDIR)
SINU = SIN(USDIR)
!
TAUNWX = 0.
TAUNWY = 0.
TAUWX = 0.
TAUWY = 0.
!
!/ --- scale friction velocity to wind speed (10m) in
!/ the boundary layer ----------------------------------------- /
!/ Donelan et al. (2006) used U10 or U_{λ/2} in their S_{in}
!/ parameterization. To avoid some disadvantages of using U10 or
!/ U_{λ/2}, Rogers et al. (2012) used the following engineering
!/ conversion:
!/ UPROXY = SIN6WS * UST
!/
!/ SIN6WS = 28.0 following Komen et al. (1984)
!/ SIN6WS = 32.0 suggested by E. Rogers (2014)
!
UPROXY = SIN6WS * USTAR ! Scale wind speed
!
ECOS2 = ECOS(1:NSPEC) ! Only indices from 1 to NSPEC
ESIN2 = ESIN(1:NSPEC) ! are requested.
!
IKN = IRANGE(1,NSPEC,NTH) ! Index vector for elements of 1 ... NK
! ! such that e.g. SIG(1:NK) = SIG2(IKN).
DSII2 = DDEN2 / DTH / SIG2 ! Frequency bandwidths (int.) (rad)
DSII = DSII2(IKN)
SIG = SIG2(IKN)
WN = WN2(IKN)
CINV2 = WN2 / SIG2 ! inverse phase speed
!
DO ITH = 1, NTH
CG2(IKN+(ITH-1)) = CG ! Apply CG to all directions.
END DO
!
!/ 1) --- calculate 1d action density spectrum (A(sigma)) and
!/ zero-out values less than 1.0E-32 to avoid NaNs when
!/ computing directional narrowness in step 4). --------------- /
K = RESHAPE(A,(/ NTH,NK /))
ADENSIG = SUM(K,1) * SIG * DTH ! Integrate over directions.
!
!/ 2) --- calculate normalised directional spectrum K(theta,sigma) --- /
KMAX = MAXVAL(K,1)
DO IK = 1,NK
IF (KMAX(IK).LT.1.0E-34) THEN
K(1:NTH,IK) = 1.
ELSE
K(1:NTH,IK) = K(1:NTH,IK)/KMAX(IK)
END IF
END DO
!
!/ 3) --- calculate normalised spectral saturation BN(IK) ------------ /
ANAR = 1.0/( SUM(K,1) * DTH ) ! directional narrowness
!
SQRTBN = SQRT( ANAR * ADENSIG * WN**3 )
DO ITH = 1, NTH
SQRTBN2(IKN+(ITH-1)) = SQRTBN ! Calculate SQRTBN for
END DO ! the entire spectrum.
!
!/ 4) --- calculate growth rate GAMMA and S for all directions for
!/ following winds (U10/c - 1 is positive; W1) and in 7) for
!/ adverse winds (U10/c -1 is negative, W2). W1 and W2
!/ complement one another. ------------------------------------ /
W1 = MAX( 0., UPROXY * CINV2* ( ECOS2*COSU + ESIN2*SINU ) - 1. )**2
!
D = (DAIR / DWAT) * SIG2 * &
(2.8 - ( 1. + TANH(10.*SQRTBN2*W1 - 11.) )) *SQRTBN2*W1
!
S = D * A
!
!/ 5) --- calculate reduction factor LFACT using non-directional
! spectral density of the wind input ------------------------- /
CINV = CINV2(IKN)
SDENSIG = RESHAPE(S*SIG2/CG2,(/ NTH,NK /))
CALL LFACTOR(SDENSIG, CINV, UABS, USTAR, USDIR, SIG, DSII, &
LFACT, TAUWX, TAUWY )
!
!/ 6) --- apply reduction (LFACT) to the entire spectrum ------------- /
IF (SUM(LFACT) .LT. NK) THEN
DO ITH = 1, NTH
D(IKN+ITH-1) = D(IKN+ITH-1) * LFACT
END DO
S = D * A
END IF
!
!/ 7) --- compute negative wind input for adverse winds. negative
!/ growth is typically smaller by a factor of ~2.5 (=.28/.11)
!/ than those for the favourable winds [Donelan, 2006, Eq. (7)].
!/ the factor is adjustable with NAMELIST parameter in
!/ ww3_grid.inp: '&SIN6 SINA0 = 0.04 /' ----------------------- /
IF (SIN6A0.GT.0.0) THEN
W2 = MIN( 0., UPROXY * CINV2* ( ECOS2*COSU + ESIN2*SINU ) - 1. )**2
D = D - ( (DAIR / DWAT) * SIG2 * SIN6A0 * &
(2.8 - ( 1. + TANH(10.*SQRTBN2*W2 - 11.) )) *SQRTBN2*W2 )
S = D * A
! --- compute negative component of the wave supported stresses
! from negative part of the wind input ---------------------- /
SDENSIG = RESHAPE(S*SIG2/CG2,(/ NTH,NK /))
CALL TAU_WAVE_ATMOS(SDENSIG, CINV, SIG, DSII, TAUNWX, TAUNWY )
END IF
!
!/
!/ End of W3SIN6 ----------------------------------------------------- /
!/
END SUBROUTINE W3SIN6
!/ ------------------------------------------------------------------- /
!>
!> @brief Observation-based source term for dissipation.
!>
!> @details After Babanin et al. (2010) following the implementation by
!> Rogers et al. (2012). The dissipation function Sds accommodates an
!> inherent breaking term T1 and an additional cumulative term T2 at
!> all frequencies above the peak. The forced dissipation term T2 is
!> an integral that grows toward higher frequencies and dominates at
!> smaller scales (Babanin et al. 2010).
!>
!> References:
!> Babanin et al. 2010: JPO 40(4), 667-683
!> Rogers et al. 2012: JTECH 29(9) 1329-1346
!>
!> @param[in] A Action density spectrum.
!> @param[in] CG Group velocities.
!> @param[in] WN Wavenumbers.
!> @param[out] S Source term (1-D version).
!> @param[out] D Diagonal term of derivative.
!>
!> @author S. Zieger
!> @author Q. Liu
!> @date 26-Jun-2018
!>
SUBROUTINE W3SDS6 (A, CG, WN, S, D)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | S. Zieger |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 26-Jun-2018 |
!/ +-----------------------------------+
!/
!/ 23-Jun-2010 : Origination. ( version 4.04 )
!/ (S. Zieger)
!/ 26-Jun-2018 : Revise the width of the last bin ( version 6.06 )
!/ (Q. Liu)
!/
! 1. Purpose :
!
! Observation-based source term for dissipation after Babanin et al.
! (2010) following the implementation by Rogers et al. (2012). The
! dissipation function Sds accommodates an inherent breaking term T1
! and an additional cumulative term T2 at all frequencies above the
! peak. The forced dissipation term T2 is an integral that grows
! toward higher frequencies and dominates at smaller scales
! (Babanin et al. 2010).
!
! References:
! Babanin et al. 2010: JPO 40(4), 667-683
! Rogers et al. 2012: JTECH 29(9) 1329-1346
!
! 2. Method :
!
! Sds = (T1 + T2) * E
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! A¹ R.A. I Action density spectrum
! CG R.A. I Group velocities
! WN R.A. I Wavenumbers
! S¹ R.A. O Source term (1-D version)
! D¹ R.A. O Diagonal term of derivative
! ¹ Stored as 1-D array with dimension NTH*NK (column by column).
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3SRCE Subr. W3SRCEMD Source term integration.
! W3EXPO Subr. N/A Point output post-processor.
! GXEXPO Subr. N/A GrADS point output post-processor.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
! !/T6 Test output for dissipation terms T1 and T2.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS, ONLY: GRAV, TPI
USE W3GDATMD, ONLY: NK, NTH, NSPEC, DDEN, DSII, SIG2, DTH, XFR
USE W3GDATMD, ONLY: SDS6A1, SDS6A2, SDS6P1, SDS6P2, SDS6ET
USE W3ODATMD, ONLY: NDSE
USE W3SERVMD, ONLY: EXTCDE
#ifdef W3_T6
USE W3TIMEMD, ONLY: STME21
USE W3WDATMD, ONLY: TIME
USE W3ODATMD, ONLY: NDST
#endif
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
REAL, INTENT(IN) :: A(NSPEC), CG(NK), WN(NK)
REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC)
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
INTEGER :: IK, ITH, IKN(NK)
REAL :: FREQ(NK) ! frequencies [Hz]
REAL :: DFII(NK) ! frequency bandwiths [Hz]
REAL :: ANAR(NK) ! directional narrowness
REAL :: BNT ! empirical constant for
! wave breaking probability
REAL :: EDENS (NK) ! spectral density E(f)
REAL :: ETDENS(NK) ! threshold spec. density ET(f)
REAL :: EXDENS(NK) ! excess spectral density EX(f)
REAL :: NEXDENS(NK) ! normalised excess spectral density
REAL :: T1(NK) ! inherent breaking term
REAL :: T2(NK) ! forced dissipation term
REAL :: T12(NK) ! = T1+T2 or combined dissipation
REAL :: ADF(NK), XFAC, EDENSMAX ! temp. variables
#ifdef W3_T6
CHARACTER(LEN=23) :: IDTIME
#endif
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'W3SDS6')
#endif
!
!/ 0) --- Initialize essential parameters ---------------------------- /
IKN = IRANGE(1,NSPEC,NTH) ! Index vector for elements of 1,
! ! 2,..., NK such that for example
! ! SIG(1:NK) = SIG2(IKN).
FREQ = SIG2(IKN)/TPI
ANAR = 1.0
BNT = 0.035**2
T1 = 0.0
T2 = 0.0
NEXDENS = 0.0
!
!/ 1) --- Calculate threshold spectral density, spectral density, and
!/ the level of exceedence EXDENS(f) -------------------------- /
ETDENS = ( TPI * BNT ) / ( ANAR * CG * WN**3 )
EDENS = SUM(RESHAPE(A,(/ NTH,NK /)),1) * TPI * SIG2(IKN) * DTH / CG ! E(f)
EXDENS = MAX(0.0,EDENS-ETDENS)
!
!/ --- normalise by a generic spectral density -------------------- /
IF (SDS6ET) THEN ! ww3_grid.inp: &SDS6 SDSET = T or F
NEXDENS = EXDENS / ETDENS ! normalise by threshold spectral density
ELSE ! normalise by spectral density
EDENSMAX = MAXVAL(EDENS)*1E-5
IF (ALL(EDENS .GT. EDENSMAX)) THEN
NEXDENS = EXDENS / EDENS
ELSE
DO IK = 1,NK
IF (EDENS(IK) .GT. EDENSMAX) NEXDENS(IK) = EXDENS(IK) / EDENS(IK)
END DO
END IF
END IF
!
!/ 2) --- Calculate inherent breaking component T1 ------------------- /
T1 = SDS6A1 * ANAR * FREQ * (NEXDENS**SDS6P1)
!
!/ 3) --- Calculate T2, the dissipation of waves induced by
!/ the breaking of longer waves T2 ---------------------------- /
ADF = ANAR * (NEXDENS**SDS6P2)
XFAC = (1.0-1.0/XFR)/(XFR-1.0/XFR)
DO IK = 1,NK
DFII = DSII/TPI
! IF (IK .GT. 1) DFII(IK) = DFII(IK) * XFAC
IF (IK .GT. 1 .AND. IK .LT. NK) DFII(IK) = DFII(IK) * XFAC
T2(IK) = SDS6A2 * SUM( ADF(1:IK)*DFII(1:IK) )
END DO
!
!/ 4) --- Sum up dissipation terms and apply to all directions ------- /
T12 = -1.0 * ( MAX(0.0,T1)+MAX(0.0,T2) )
DO ITH = 1, NTH
D(IKN+ITH-1) = T12
END DO
!
S = D * A
!
!/ 5) --- Diagnostic output (switch !/T6) ---------------------------- /
#ifdef W3_T6
CALL STME21 ( TIME , IDTIME )
WRITE (NDST,270) 'T1*E',IDTIME(1:19),(T1*EDENS)
WRITE (NDST,270) 'T2*E',IDTIME(1:19),(T2*EDENS)
WRITE (NDST,271) SUM(SUM(RESHAPE(S,(/ NTH,NK /)),1)*DDEN/CG)
!
270 FORMAT (' TEST W3SDS6 : ',A,'(',A,')',':',70E11.3)
271 FORMAT (' TEST W3SDS6 : Total SDS =',E13.5)
#endif
!/
!/ End of W3SDS6 ----------------------------------------------------- /
!/
END SUBROUTINE W3SDS6
!/ ------------------------------------------------------------------- /
!/
!>
!> @brief Numerical approximation for the reduction factor.
!>
!> @details Numerical approximation for the reduction factor LFACTOR(f) to
!> reduce energy in the high-frequency part of the resolved part
!> of the spectrum to meet the constraint on total stress (TAU).
!> The constraint is TAU <= TAU_TOT (TAU_TOT = TAU_WAV + TAU_VIS),
!> thus the wind input is reduced to match our constraint.
!>
!> @param[in] S Wind input energy density spectrum.
!> @param[in] CINV Inverse phase speed.
!> @param[in] U10 Wind speed.
!> @param[in] USTAR Friction velocity.
!> @param[in] USDIR Wind direction.
!> @param[in] SIG Relative frequencies (in rad.).
!> @param[in] DSII Frequency bandwidths (in rad.).
!> @param[out] LFACT Factor array.
!> @param[out] TAUWX Component of the wave-supported stress.
!> @param[out] TAUWY Component of the wave-supported stress.
!>
!> @author S. Zieger
!> @author Q. Liu
!> @date 26-Jun-2018
!>
SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
LFACT, TAUWX, TAUWY )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | S. Zieger |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 26-Jun-2018 |
!/ +-----------------------------------+
!/
!/ 15-Feb-2011 : Implemented following Rogers et al. (2012)
!/ (S. Zieger)
!/ 26-Jun-2018 : UPROXY, DSII10Hz Updates ( version 6.06 )
!/ (Q. Liu )
!
! Rogers et al. (2012) JTECH 29(9), 1329-1346
!
! 1. Purpose :
!
! Numerical approximation for the reduction factor LFACTOR(f) to
! reduce energy in the high-frequency part of the resolved part
! of the spectrum to meet the constraint on total stress (TAU).
! The constraint is TAU <= TAU_TOT (TAU_TOT = TAU_WAV + TAU_VIS),
! thus the wind input is reduced to match our constraint.
!
! 2. Method :
!
! 1) If required, extend resolved part of the spectrum to 10Hz using
! an approximation for the spectral slope at the high frequency
! limit: Sin(F) prop. F**(-2) and for E(F) prop. F**(-5).
! 2) Calculate stresses:
! total stress: TAU_TOT = DAIR * USTAR**2
! viscous stress: TAU_VIS = DAIR * Cv * U10**2
! viscous stress (x,y-components):
! TAUV_X = TAU_VIS * COS(USDIR)
! TAUV_Y = TAU_VIS * SIN(USDIR)
! wave supported stress (x,y-components): /10Hz
! TAUW_X,Y = GRAV * DWAT * | [SinX,Y(F)]/C(F) dF
! /
! total stress (input): TAU = SQRT( (TAUW_X + TAUV_X)**2
! + (TAUW_Y + TAUV_Y)**2 )
! 3) If TAU does not meet our constraint reduce the wind input
! using reduction factor:
! LFACT(F) = MIN(1,exp((1-U/C(F))*RTAU))
! Then alter RTAU and repeat 3) until our constraint is matched.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! S R.A. I Wind input energy density spectrum (S_{in}(σ, θ))
! CINV R.A. I Inverse phase speed 1/C(sigma)
! U10 Real I Wind speed (10m)
! USTAR Real I Friction velocity
! USDIR Real I Wind direction
! SIG R.A. I Relative frequencies [in rad.]
! DSII R.A. I Frequency bandwiths [in rad.]
! LFACTOR R.A. O Factor array LFACT(sigma)
! TAUWX-Y Real O Component of the wave-supported stress
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! IRANGE Func. Private Index generator (ie, array addressing)
! TAUWINDS Func. Private Normal stress calculation (TAU_NRM)
! ----------------------------------------------------------------
!
! ----------------------------------------------------------------
!
! 5. Error messages :
!
! A warning is issued to NDST using format 280 if the iteration
! procedure reaches the upper iteration limit (ITERMAX). In this
! case the last approximation for RTAU is used.
!
!/
USE CONSTANTS, ONLY: DAIR, GRAV, TPI
USE W3GDATMD, ONLY: NK, NTH, NSPEC, DTH, XFR, ECOS, ESIN
USE W3GDATMD, ONLY: SIN6WS
USE W3ODATMD, ONLY: NDST, NDSE, IAPROC, NAPERR
USE W3TIMEMD, ONLY: STME21
USE W3WDATMD, ONLY: TIME
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
IMPLICIT NONE
!
!/ ------ I/O parameters --------------------------------------------- /
REAL, INTENT(IN) :: S(NTH,NK) ! wind-input source term Sin
REAL, INTENT(IN) :: CINV(NK) ! inverse phase speed
REAL, INTENT(IN) :: U10 ! wind speed
REAL, INTENT(IN) :: USTAR, USDIR ! friction velocity & direction
REAL, INTENT(IN) :: SIG(NK) ! relative frequencies
REAL, INTENT(IN) :: DSII(NK) ! frequency bandwidths
REAL, INTENT(OUT) :: LFACT(NK) ! correction factor
REAL, INTENT(OUT) :: TAUWX, TAUWY ! normal stress components
!
!/ --- local parameters (in order of appearance) ------------------ /
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
REAL, PARAMETER :: FRQMAX = 10. ! Upper freq. limit to extrapolate to.
INTEGER, PARAMETER:: ITERMAX = 80 ! Maximum number of iterations to
! find numerical solution for LFACT.
INTEGER :: IK, NK10Hz, SIGN_NEW, SIGN_OLD
!
REAL :: ECOS2(NSPEC), ESIN2(NSPEC)
REAL, ALLOCATABLE :: IK10Hz(:), LF10Hz(:), SIG10Hz(:), CINV10Hz(:)
REAL, ALLOCATABLE :: SDENS10Hz(:), SDENSX10Hz(:), SDENSY10Hz(:)
REAL, ALLOCATABLE :: DSII10Hz(:), UCINV10Hz(:)
REAL :: TAU_TOT, TAU, TAU_VIS, TAU_WAV
REAL :: TAUVX, TAUVY, TAUX, TAUY
REAL :: TAU_NND, TAU_INIT(2)
REAL :: UPROXY, RTAU, DRTAU, ERR
LOGICAL :: OVERSHOT
CHARACTER(LEN=23) :: IDTIME
!
!/ ------------------------------------------------------------------- /
#ifdef W3_S
CALL STRACE (IENT, 'LFACTOR')
#endif
!
!/ 0) --- Find the number of frequencies required to extend arrays
!/ up to f=10Hz and allocate arrays --------------------------- /
!/ ALOG is the same as LOG
NK10Hz = CEILING(ALOG(FRQMAX/(SIG(1)/TPI))/ALOG(XFR))+1
NK10Hz = MAX(NK,NK10Hz)
!
ALLOCATE(IK10Hz(NK10Hz))
IK10Hz = REAL( IRANGE(1,NK10Hz,1) )
!
ALLOCATE(SIG10Hz(NK10Hz))
ALLOCATE(CINV10Hz(NK10Hz))
ALLOCATE(DSII10Hz(NK10Hz))
ALLOCATE(LF10Hz(NK10Hz))
ALLOCATE(SDENS10Hz(NK10Hz))
ALLOCATE(SDENSX10Hz(NK10Hz))
ALLOCATE(SDENSY10Hz(NK10Hz))
ALLOCATE(UCINV10Hz(NK10Hz))
!
ECOS2 = ECOS(1:NSPEC)
ESIN2 = ESIN(1:NSPEC)
!
!/ 1) --- Either extrapolate arrays up to 10Hz or use discrete spectral
! grid per se. Limit the constraint to the positive part of the
! wind input only. ---------------------------------------------- /
IF (NK .LT. NK10Hz) THEN
SDENS10Hz(1:NK) = SUM(S,1) * DTH
SDENSX10Hz(1:NK) = SUM(MAX(0.,S)*RESHAPE(ECOS2,(/NTH,NK/)),1) * DTH
SDENSY10Hz(1:NK) = SUM(MAX(0.,S)*RESHAPE(ESIN2,(/NTH,NK/)),1) * DTH
SIG10Hz = SIG(1)*XFR**(IK10Hz-1.0)
CINV10Hz(1:NK) = CINV
CINV10Hz(NK+1:NK10Hz) = SIG10Hz(NK+1:NK10Hz)*0.101978 ! 1/c=σ/g
DSII10Hz = 0.5 * SIG10Hz * (XFR-1.0/XFR)
! The first and last frequency bin:
DSII10Hz(1) = 0.5 * SIG10Hz(1) * (XFR-1.0)
DSII10Hz(NK10Hz) = 0.5 * SIG10Hz(NK10Hz) * (XFR-1.0) / XFR
!
! --- Spectral slope for S_IN(F) is proportional to F**(-2) ------ /
SDENS10Hz(NK+1:NK10Hz) = SDENS10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
SDENSX10Hz(NK+1:NK10Hz) = SDENSX10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
SDENSY10hz(NK+1:NK10Hz) = SDENSY10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
ELSE
SIG10Hz = SIG
CINV10Hz = CINV
DSII10Hz = DSII
SDENS10Hz(1:NK) = SUM(S,1) * DTH
SDENSX10Hz(1:NK) = SUM(MAX(0.,S)*RESHAPE(ECOS2,(/NTH,NK/)),1) * DTH
SDENSY10Hz(1:NK) = SUM(MAX(0.,S)*RESHAPE(ESIN2,(/NTH,NK/)),1) * DTH
END IF
!
!/ 2) --- Stress calculation ----------------------------------------- /
! --- The total stress ------------------------------------------- /
TAU_TOT = USTAR**2 * DAIR
!
! --- The viscous stress and check that it does not exceed
! the total stress. ------------------------------------------ /
TAU_VIS = MAX(0.0, -5.0E-5*U10 + 1.1E-3) * U10**2 * DAIR
! TAU_VIS = MIN(0.9 * TAU_TOT, TAU_VIS)
TAU_VIS = MIN(0.95 * TAU_TOT, TAU_VIS)
!
TAUVX = TAU_VIS * COS(USDIR)
TAUVY = TAU_VIS * SIN(USDIR)
!
! --- The wave supported stress. --------------------------------- /
TAUWX = TAUWINDS(SDENSX10Hz,CINV10Hz,DSII10Hz) ! normal stress (x-component)
TAUWY = TAUWINDS(SDENSY10Hz,CINV10Hz,DSII10Hz) ! normal stress (y-component)
TAU_NND = TAUWINDS(SDENS10Hz, CINV10Hz,DSII10Hz) ! normal stress (non-directional)
TAU_WAV = SQRT(TAUWX**2 + TAUWY**2) ! normal stress (magnitude)
TAU_INIT = (/TAUWX,TAUWY/) ! unadjusted normal stress components
!
TAUX = TAUVX + TAUWX ! total stress (x-component)
TAUY = TAUVY + TAUWY ! total stress (y-component)
TAU = SQRT(TAUX**2 + TAUY**2) ! total stress (magnitude)
ERR = (TAU-TAU_TOT)/TAU_TOT ! initial error
!
!/ 3) --- Find reduced Sin(f) = L(f)*Sin(f) to satisfy our constraint
!/ TAU <= TAU_TOT --------------------------------------------- /
CALL STME21 ( TIME , IDTIME )
LF10Hz = 1.0
IK = 0
!
IF (TAU .GT. TAU_TOT) THEN
OVERSHOT = .FALSE.
RTAU = ERR / 90.
DRTAU = 2.0
SIGN_NEW = INT(SIGN(1.0,ERR))
UPROXY = SIN6WS * USTAR
UCINV10Hz = 1.0 - (UPROXY * CINV10Hz)
!
#ifdef W3_T6
WRITE (NDST,270) IDTIME, U10
WRITE (NDST,271)
#endif
DO IK=1,ITERMAX
LF10Hz = MIN(1.0, EXP(UCINV10Hz * RTAU) )
!
TAU_NND = TAUWINDS(SDENS10Hz *LF10Hz,CINV10Hz,DSII10Hz)
TAUWX = TAUWINDS(SDENSX10Hz*LF10Hz,CINV10Hz,DSII10Hz)
TAUWY = TAUWINDS(SDENSY10Hz*LF10Hz,CINV10Hz,DSII10Hz)
TAU_WAV = SQRT(TAUWX**2 + TAUWY**2)
!
TAUX = TAUVX + TAUWX
TAUY = TAUVY + TAUWY
TAU = SQRT(TAUX**2 + TAUY**2)
ERR = (TAU-TAU_TOT) / TAU_TOT
!
SIGN_OLD = SIGN_NEW
SIGN_NEW = INT(SIGN(1.0, ERR))
#ifdef W3_T6
WRITE (NDST,272) IK, RTAU, DRTAU, TAU, TAU_TOT, ERR, &
TAUWX, TAUWY, TAUVX, TAUVY, TAU_NND
#endif
!
! --- Slow down DRTAU when overshot. -------------------------- /
IF (SIGN_NEW .NE. SIGN_OLD) OVERSHOT = .TRUE.
IF (OVERSHOT) DRTAU = MAX(0.5*(1.0+DRTAU),1.00010)
!
RTAU = RTAU * (DRTAU**SIGN_NEW)
!
IF (ABS(ERR) .LT. 1.54E-4) EXIT
END DO
!
IF (IK .GE. ITERMAX) WRITE (NDST,280) IDTIME(1:19), U10, TAU, &
TAU_TOT, ERR, TAUWX, TAUWY, TAUVX, TAUVY,TAU_NND
END IF
!
LFACT(1:NK) = LF10Hz(1:NK)
!
#ifdef W3_T6
WRITE (NDST,273) 'Sin ', IDTIME(1:19), SDENS10Hz*TPI
WRITE (NDST,273) 'SinR', IDTIME(1:19), SDENS10Hz*LF10Hz*TPI
WRITE (NDST,274) 'Sin ', SUM(SDENS10Hz(1:NK)*DSII)
WRITE (NDST,274) 'SinR ', SUM(SDENS10Hz(1:NK)*LF10Hz(1:NK)*DSII)
WRITE (NDST,274) 'SinR/C', TAUWINDS(SDENS10Hz(1:NK)*LFACT,CINV,DSII)
!
270 FORMAT (' TEST W3SIN6 : LFACTOR SUBROUTINE CALCULATING FOR ', &
A,' U10=',F5.1 )
271 FORMAT (' TEST W3SIN6 : IK RTAU DRTAU TAU TAU_TOT' &
' ERR TAUW_X TAUW_Y TAUV_X TAUV_Y TAU1D' )
272 FORMAT (' TEST W3SIN6 : ',I2,2F9.5,2F8.5,E10.2,4F7.4,F7.3 )
273 FORMAT (' TEST W3SIN6 : ',A,'(',A,'):', 70E11.3 )
#endif
274 FORMAT (' TEST W3SIN6 : Total ',A,' =', E13.5 )
280 FORMAT (' WARNING LFACTOR (TIME,U10,TAU,TAU_TOT,ERR,TAUW_XY,' &
'TAUV_XY,TAU_SCALAR): ',A,F6.1,2F7.4,E10.3,4F7.4,F7.3 )
!
DEALLOCATE(IK10Hz,SIG10Hz,CINV10Hz,DSII10Hz,LF10Hz)
DEALLOCATE(SDENS10Hz,SDENSX10Hz,SDENSY10Hz,UCINV10Hz)
!/
END SUBROUTINE LFACTOR
!/ ------------------------------------------------------------------- /
!/
!>
!> @brief Calculated the stress for the negative part of the input term.
!>
!> @details Calculated the stress for the negative part of the input term,
!> that is the stress from the waves to the atmosphere. Relevant
!> in the case of opposing winds.
!>
!> @param[in] S Wind-input source term Sin.
!> @param[in] CINV Inverse phase speed.
!> @param[in] SIG Relative frequencies.
!> @param[in] DSII Frequency bandwidths.
!> @param[out] TAUNWX Stress components (wave->atmos).
!> @param[out] TAUNWY Stress components (wave->atmos).
!>
!> @author S. Zieger
!> @author Q. Liu
!> @date 26-Jun-2018
!>
SUBROUTINE TAU_WAVE_ATMOS(S, CINV, SIG, DSII, TAUNWX, TAUNWY )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | S. Zieger |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 26-Jun-2018 |
!/ +-----------------------------------+
!/
!/ 24-Oct-2013 : Origination following LFACTOR
!/ (S. Zieger)
!/ 26-Jun-2018 : Updates on DSII10Hz ( version 6.06)
!/ (Q. Liu)
!
! 1. Purpose :
!
! Calculated the stress for the negative part of the input term,
! that is the stress from the waves to the atmosphere. Relevant
! in the case of opposing winds.
!
! 2. Method :
! 1) If required, extend resolved part of the spectrum to 10Hz using
! an approximation for the spectral slope at the high frequency
! limit: Sin(F) prop. F**(-2) and for E(F) prop. F**(-5).
! 2) Calculate stresses:
! stress components (x,y): /10Hz
! TAUNW_X,Y = GRAV * DWAT * | [SinX,Y(F)]/C(F) dF
! /
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! S R.A. I Wind input energy density spectrum
! CINV R.A. I Inverse phase speed 1/C(sigma)
! SIG R.A. I Relative frequencies [in rad.]
! DSII R.A. I Frequency bandwiths [in rad.]
! TAUNWX-Y Real O Component of the negative wave-supported stress
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! IRANGE Func. Private Index generator (ie, array addressing)
! TAUWINDS Func. Private Normal stress calculation (TAU_NRM)
! ----------------------------------------------------------------
!
! 5. Source code :
!
!/
USE CONSTANTS, ONLY: GRAV, TPI
USE W3GDATMD, ONLY: NK, NTH, NSPEC, DTH, XFR, ECOS, ESIN
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
IMPLICIT NONE
!
!/ ------ I/O parameters --------------------------------------------- /
REAL, INTENT(IN) :: S(NTH,NK) ! wind-input source term Sin
REAL, INTENT(IN) :: CINV(NK) ! inverse phase speed
REAL, INTENT(IN) :: SIG(NK) ! relative frequencies
REAL, INTENT(IN) :: DSII(NK) ! frequency bandwidths
REAL, INTENT(OUT) :: TAUNWX, TAUNWY ! stress components (wave->atmos)
!
!/ --- local parameters (in order of appearance) ------------------ /
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
REAL, PARAMETER :: FRQMAX = 10. ! Upper freq. limit to extrapolate to.
INTEGER :: NK10Hz
!
REAL :: ECOS2(NSPEC), ESIN2(NSPEC)
REAL, ALLOCATABLE :: IK10Hz(:), SIG10Hz(:), CINV10Hz(:)
REAL, ALLOCATABLE :: SDENSX10Hz(:), SDENSY10Hz(:)
REAL, ALLOCATABLE :: DSII10Hz(:), UCINV10Hz(:)
!
!/ ------------------------------------------------------------------- /
#ifdef W3_S
CALL STRACE (IENT, 'TAU_WAVE_ATMOS')
#endif
!
!/ 0) --- Find the number of frequencies required to extend arrays
!/ up to f=10Hz and allocate arrays --------------------------- /
NK10Hz = CEILING(ALOG(FRQMAX/(SIG(1)/TPI))/ALOG(XFR))+1
NK10Hz = MAX(NK,NK10Hz)
!
ALLOCATE(IK10Hz(NK10Hz))
IK10Hz = REAL( IRANGE(1,NK10Hz,1) )
!
ALLOCATE(SIG10Hz(NK10Hz))
ALLOCATE(CINV10Hz(NK10Hz))
ALLOCATE(DSII10Hz(NK10Hz))
ALLOCATE(SDENSX10Hz(NK10Hz))
ALLOCATE(SDENSY10Hz(NK10Hz))
ALLOCATE(UCINV10Hz(NK10Hz))
!
ECOS2 = ECOS(1:NSPEC)
ESIN2 = ESIN(1:NSPEC)
!
!/ 1) --- Either extrapolate arrays up to 10Hz or use discrete spectral
! grid per se. Limit the constraint to the positive part of the
! wind input only. ---------------------------------------------- /
IF (NK .LT. NK10Hz) THEN
SDENSX10Hz(1:NK) = SUM(ABS(MIN(0.,S))*RESHAPE(ECOS2,(/NTH,NK/)),1) * DTH
SDENSY10Hz(1:NK) = SUM(ABS(MIN(0.,S))*RESHAPE(ESIN2,(/NTH,NK/)),1) * DTH
SIG10Hz = SIG(1)*XFR**(IK10Hz-1.0)
CINV10Hz(1:NK) = CINV
CINV10Hz(NK+1:NK10Hz) = SIG10Hz(NK+1:NK10Hz)*0.101978
DSII10Hz = 0.5 * SIG10Hz * (XFR-1.0/XFR)
! The first and last frequency bin:
DSII10Hz(1) = 0.5 * SIG10Hz(1) * (XFR-1.0)
DSII10Hz(NK10Hz) = 0.5 * SIG10Hz(NK10Hz) * (XFR-1.0) / XFR
!
! --- Spectral slope for S_IN(F) is proportional to F**(-2) ------ /
SDENSX10Hz(NK+1:NK10Hz) = SDENSX10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
SDENSY10hz(NK+1:NK10Hz) = SDENSY10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
ELSE
SIG10Hz = SIG
CINV10Hz = CINV
DSII10Hz = DSII
SDENSX10Hz(1:NK) = SUM(ABS(MIN(0.,S))*RESHAPE(ECOS2,(/NTH,NK/)),1) * DTH
SDENSY10Hz(1:NK) = SUM(ABS(MIN(0.,S))*RESHAPE(ESIN2,(/NTH,NK/)),1) * DTH
END IF
!
!/ 2) --- Stress calculation ----------------------------------------- /
! --- The wave supported stress (waves to atmosphere) ------------ /
TAUNWX = TAUWINDS(SDENSX10Hz,CINV10Hz,DSII10Hz) ! x-component
TAUNWY = TAUWINDS(SDENSY10Hz,CINV10Hz,DSII10Hz) ! y-component
!/
END SUBROUTINE TAU_WAVE_ATMOS
!/ ------------------------------------------------------------------- /
!/
!>
!> @brief Generate a sequence of linear-spaced integer numbers.
!>
!> @details Used for instance array addressing (indexing).
!>
!> @param X0
!> @param X1
!> @param DX
!> @returns IX
!>
!> @author S. Zieger
!> @date 15-Feb-2011
!>
FUNCTION IRANGE(X0,X1,DX) RESULT(IX)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | S. Zieger |
!/ | FORTRAN 90 |
!/ | Last update : 15-Feb-2011 |
!/ +-----------------------------------+
!/
!/ 15-Feb-2011 : Origination ( version 4.04 )
!/ (S. Zieger)
!/
! 1. Purpose :
! Generate a sequence of linear-spaced integer numbers.
! Used for instance array addressing (indexing).
!
!/
IMPLICIT NONE
INTEGER, INTENT(IN) :: X0, X1, DX
INTEGER, ALLOCATABLE :: IX(:)
INTEGER :: N
INTEGER :: I
!
N = INT(REAL(X1-X0)/REAL(DX))+1
ALLOCATE(IX(N))
DO I = 1, N
IX(I) = X0+ (I-1)*DX
END DO
!/
END FUNCTION IRANGE
!/ ------------------------------------------------------------------- /
!/
!>
!> @brief Wind stress (tau) computation from wind-momentum-input.
!>
!> @verbatim
!> Wind stress (tau) computation from wind-momentum-input
!> function which can be obtained from wind-energy-input (Sin).
!>
!> / FRMAX
!> tau = g * rho_water * | Sin(f)/C(f) df
!> /
!> @endverbatim
!>
!> @param SDENSIG Sin(sigma) in (m2/rad-Hz).
!> @param CINV Inverse phase speed.
!> @param DSII Frequency bandwidths in (rad).
!> @returns TAU_WINDS Wind stress.
!>
!> @author S. Zieger
!> @date 13-Aug-2012
!>
FUNCTION TAUWINDS(SDENSIG,CINV,DSII) RESULT(TAU_WINDS)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | S. Zieger |
!/ | FORTRAN 90 |
!/ | Last update : 13-Aug-2012 |
!/ +-----------------------------------+
!/
!/ 15-Feb-2011 : Origination ( version 4.04 )
!/ (S. Zieger)
!/
! 1. Purpose :
! Wind stress (tau) computation from wind-momentum-input
! function which can be obtained from wind-energy-input (Sin).
!
! / FRMAX
! tau = g * rho_water * | Sin(f)/C(f) df
! /
!/
USE CONSTANTS, ONLY: GRAV, DWAT ! gravity, density of water
IMPLICIT NONE
REAL, INTENT(IN) :: SDENSIG(:) ! Sin(sigma) in [m2/rad-Hz]
REAL, INTENT(IN) :: CINV(:) ! inverse phase speed
REAL, INTENT(IN) :: DSII(:) ! freq. bandwidths in [radians]
REAL :: TAU_WINDS ! wind stress
!
TAU_WINDS = GRAV * DWAT * SUM(SDENSIG*CINV*DSII)
!/
END FUNCTION TAUWINDS
!/ ------------------------------------------------------------------- /
!/
!/ End of module W3SRC6MD -------------------------------------------- /
!/
END MODULE W3SRC6MD