w3src6md Module Reference

Observation-based wind input and dissipation after Donelan et al (2006), and Babanin et al. More...

Functions/Subroutines
subroutine, public	w3spr6 (A, CG, WN, EMEAN, FMEAN, WNMEAN, AMAX, FP)
	Calculate mean wave parameters. More...
subroutine, public	w3sin6 (A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, TAUWX, TAUWY, TAUNWX, TAUNWY, S, D)
	Observation-based source term for wind input. More...
subroutine, public	w3sds6 (A, CG, WN, S, D)
	Observation-based source term for dissipation. More...
subroutine	tau_wave_atmos (S, CINV, SIG, DSII, TAUNWX, TAUNWY)
	Calculated the stress for the negative part of the input term. More...

Detailed Description

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.

Author

S. Zieger

Q. Liu

Date

26-Jun-2018

Function/Subroutine Documentation

tau_wave_atmos()

subroutine w3src6md::tau_wave_atmos	(	real, dimension(nth,nk), intent(in)	S,
		real, dimension(nk), intent(in)	CINV,
		real, dimension(nk), intent(in)	SIG,
		real, dimension(nk), intent(in)	DSII,
		real, intent(out)	TAUNWX,
		real, intent(out)	TAUNWY
	)

Calculated the stress for the negative part of the input term.

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.

Parameters

[in]	S	Wind-input source term Sin.
[in]	CINV	Inverse phase speed.
[in]	SIG	Relative frequencies.
[in]	DSII	Frequency bandwidths.
[out]	TAUNWX	Stress components (wave->atmos).
[out]	TAUNWY	Stress components (wave->atmos).

Author

S. Zieger

Q. Liu

Date

26-Jun-2018

Definition at line 1067 of file w3src6md.F90.

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

References w3gdatmd::dth, constants::dwat, w3gdatmd::ecos, w3gdatmd::esin, constants::grav, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, w3servmd::strace(), constants::tpi, and w3gdatmd::xfr.

Referenced by w3sin6().

w3sds6()

subroutine, public w3src6md::w3sds6	(	real, dimension(nspec), intent(in)	A,
		real, dimension(nk), intent(in)	CG,
		real, dimension(nk), intent(in)	WN,
		real, dimension(nspec), intent(out)	S,
		real, dimension(nspec), intent(out)	D
	)

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

Parameters

[in]	A	Action density spectrum.
[in]	CG	Group velocities.
[in]	WN	Wavenumbers.
[out]	S	Source term (1-D version).
[out]	D	Diagonal term of derivative.

Author

S. Zieger

Q. Liu

Date

26-Jun-2018

Definition at line 547 of file w3src6md.F90.

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

References constants::dair, w3gdatmd::dden, w3gdatmd::dsii, w3gdatmd::dth, w3gdatmd::ecos, w3gdatmd::esin, w3servmd::extcde(), constants::grav, w3odatmd::iaproc, w3odatmd::naperr, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, w3gdatmd::sds6a1, w3gdatmd::sds6a2, w3gdatmd::sds6et, w3gdatmd::sds6p1, w3gdatmd::sds6p2, w3gdatmd::sig2, w3gdatmd::sin6ws, w3timemd::stme21(), w3servmd::strace(), w3wdatmd::time, constants::tpi, and w3gdatmd::xfr.

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

w3sin6()

subroutine, public w3src6md::w3sin6	(	real, dimension (nspec), intent(in)	A,
		real, dimension(nk), intent(in)	CG,
		real, dimension(nspec), intent(in)	WN2,
		real, intent(in)	UABS,
		real, intent(in)	USTAR,
		real, intent(in)	USDIR,
		real, intent(in)	CD,
		real, intent(in)	DAIR,
		real, intent(out)	TAUWX,
		real, intent(out)	TAUWY,
		real, intent(out)	TAUNWX,
		real, intent(out)	TAUNWY,
		real, dimension(nspec), intent(out)	S,
		real, dimension(nspec), intent(out)	D
	)

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

Parameters

[in]	A	Action density spectrum
[in]	CG	Group velocities.
[in]	WN2	Wavenumbers.
[in]	UABS	Wind speed at 10 m above sea level (U10).
[in]	USTAR	Friction velocity.
[in]	USDIR	Direction of USTAR.
[in]	CD	Drag coefficient.
[in]	DAIR	Air density.
[out]	TAUWX	Component of the wave-supported stress.
[out]	TAUWY	Component of the wave-supported stress.
[out]	TAUWNX	Component of the negative part of the stress.
[out]	TAUWNY	Component of the negative part of the stress.
[out]	S	Source term.
[out]	D	Diagonal term of derivative.

Author

S. Zieger

Q. Liu

Date

13-Aug-2021

Definition at line 292 of file w3src6md.F90.

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

References w3gdatmd::dden2, w3gdatmd::dth, constants::dwat, w3gdatmd::ecos, w3gdatmd::esin, w3servmd::extcde(), constants::grav, w3odatmd::ndse, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, w3gdatmd::sig2, w3gdatmd::sin6a0, w3gdatmd::sin6ws, w3servmd::strace(), tau_wave_atmos(), and constants::tpi.

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

w3spr6()

subroutine, public w3src6md::w3spr6	(	real, dimension(nth,nk), intent(in)	A,
		real, dimension(nk), intent(in)	CG,
		real, dimension(nk), intent(in)	WN,
		real, intent(out)	EMEAN,
		real, intent(out)	FMEAN,
		real, intent(out)	WNMEAN,
		real, intent(out)	AMAX,
		real, intent(out)	FP
	)

Calculate mean wave parameters.

See source term routines.

Parameters

[in,out]	A	Action as a function of direction and wavenumber.
[in,out]	CG	Group velocities.
[in,out]	WN	Wavenumbers.
[in,out]	EMEAN	Mean wave energy.
[in,out]	FMEAN	Mean wave frequency.
[in,out]	WNMEAN	Mean wavenumber.
[in,out]	AMAX	Maximum action density in spectrum.
[in,out]	FP	Peak frequency (rad).

Author

S. Zieger

Date

11-Feb-2011

Definition at line 122 of file w3src6md.F90.

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

References w3gdatmd::dden, w3gdatmd::dsii, w3gdatmd::dth, w3servmd::extcde(), w3gdatmd::fte, w3gdatmd::ftf, w3gdatmd::ftwn, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nth, w3gdatmd::sig, w3servmd::strace(), and constants::tpiinv.

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