WW3/w3src6md_8F90_source.html

#include "w3macros.h"

!/ ------------------------------------------------------------------- /


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

  !/ ------------------------------------------------------------------- /


  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

  !/ ------------------------------------------------------------------- /


  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

  !/ ------------------------------------------------------------------- /


  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

  !/ ------------------------------------------------------------------- /

  !/


  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

  !/ ------------------------------------------------------------------- /

  !/


  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

  !/ ------------------------------------------------------------------- /

  !/


  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

  !/ ------------------------------------------------------------------- /

  !/


  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