w3src4md.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
 
MODULE w3src4md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III                SHOM |
  !/                  !            F. Ardhuin             !
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         13-Nov-2013 |
  !/                  +-----------------------------------+
  !/
  !/    30-Aug-2010 : Origination.                        ( version 3.14-Ifremer )
  !/    02-Nov-2010 : Addding fudge factor for low freq.  ( version 4.03 )
  !/    02-Sep-2011 : Clean up and time optimization      ( version 4.04 )
  !/    04-Sep-2011 : Estimation of whitecap stats.       ( version 4.04 )
  !/    13-Nov-2013 : Reduced frequency range with IG     ( version 4.13 )
  !/    01-Mar-2023 : Clean up of SDS4                    ( version 7.14 )
  !/
  !  1. Purpose :
  !
  !     The 'SHOM/Ifremer' source terms based on P.A.E.M. Janssen's wind input
  !     and dissipation functions by Ardhuin et al. (2009,2010)
  !     and Filipot & Ardhuin (2010)
  !     The wind input is converted from the original
  !     WAM codes, courtesy of P.A.E.M. Janssen and J. Bidlot
  !
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SPR4    Subr. Public   Mean parameters from spectrum.
  !      W3SIN4    Subr. Public   WAM4+ input source term.
  !      INSIN4    Subr. Public   Corresponding initialization routine.
  !      TABU_STRESS, TABU_TAUHF, TABU_TAUHF2
  !                Subr. Public   Populate various tables.
  !      CALC_USTAR
  !                Subr. Public   Compute stresses.
  !      W3SDS4    Subr. Public   Dissipation (Ardhuin & al. / Filipot & Ardhuin)
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !  7. Source code :
  !/
  !/ ------------------------------------------------------------------- /
  !/
  PUBLIC
  !/
  !/ Public variables
  !/
  !air kinematic viscosity (used in WAM)
  INTEGER, PARAMETER      :: itaumax=200,jumax=200
  INTEGER, PARAMETER      :: iustar=100,ialpha=200, ilevtail=50
  ! Tables for total stress and H.F. stress as a function of 2 or 3 variables
  REAL, ALLOCATABLE       :: taut(:,:),tauhft(:,:),tauhft2(:,:,:)
  REAL                    :: delust, delalp,deltauw, delu
  ! Table for swell damping
  REAL                    :: deltail
  REAL,    PARAMETER      :: umax    = 50.
  REAL,    PARAMETER      :: tauwmax = 2.2361 !SQRT(5.)
  INTEGER                 :: dikcumul
  !  Size of wave height table for integrating the PDF of wave heights
  INTEGER,    PARAMETER      :: nkhi=100, fac_kd2=1000
  REAL,    PARAMETER      :: fac_kd1=1.01, khsmax=2., khmax=2.
  REAL,    PARAMETER      ::kdmax=200000.
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE w3spr4 (A, CG, WN, EMEAN, FMEAN, FMEAN1, WNMEAN,     &
       AMAX, U, UDIR,                                    &
#ifdef W3_FLX5
       TAUA, TAUADIR, DAIR,                              &
#endif
       USTAR, USDIR,                                     &
       TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS, DLWMEAN)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III                SHOM |
    !/                  !            F. Ardhuin             !
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         22-Feb-2020 |
    !/                  +-----------------------------------+
    !/
    !/    03-Oct-2007 : Origination.                        ( version 3.13 )
    !/    13-Jun-2011 : Adds f_m0,-1 as FMEAN in the outout ( version 4.04 )
    !/    08-Jun-2018 : use STRACE and FLUSH                ( version 6.04 )
    !/    22-Feb-2020 : Merge Romero (2019) and cleanup     ( version 7.06 )
    !/    22-Jun-2021 : Add FLX5 to use stresses with the ST( version 7.14 )
    !/
    !  1. Purpose :
    !
    !     Calculate mean wave parameters for the use in the source term
    !     routines.
    !
    !  2. Method :
    !
    !     See source term routines.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A       R.A.  I   Action density spectrum.
    !       CG      R.A.  I   Group velocities.
    !       WN      R.A.  I   Wavenumbers.
    !       EMEAN   Real  O   Energy
    !       FMEAN1  Real  O   Mean  frequency (fm0,-1) used for reflection
    !       FMEAN   Real  O   Mean  frequency for determination of tail
    !       WNMEAN  Real  O   Mean wavenumber.
    !       AMAX    Real  O   Maximum of action spectrum.
    !       U       Real  I   Wind speed.
    !       UDIR    Real  I   Wind direction.
    !       TAUA    Real  I   Atm. total stress.            ( /!FLX5 )
    !       TAUADIR Real  I   Atm. total stress direction.  ( /!FLX5 )
    !       DAIR    Real  I   Air density.                  ( /!FLX5 )
    !       USTAR   Real I/O  Friction velocity.
    !       USDIR   Real I/O  wind stress direction.
    !       TAUWX-Y Real  I   Components of wave-supported stress.
    !       CD      Real  O   Drag coefficient at wind level ZWND.
    !       Z0      Real  O   Corresponding z0.
    !       CHARN   Real  O   Corresponding Charnock coefficient
    !       LLWS    L.A.  I   Wind sea true/false array for each component
    !       FMEANWS Real  O   Mean frequency of wind sea, used for tail
    !       DLWMEAN Real  O   Mean Long wave direction  (L. Romero 2019)
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3SRCE   Source term integration routine.
    !       W3OUTP   Point output program.
    !       GXEXPO   GrADS point output program.
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !       !/T      Enable test output.
    !       !/FLX5  Direct use of stress from atmoshperic model.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3odatmd, ONLY: iaproc
    USE constants, ONLY: tpiinv, grav, nu_air
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, dth, dden, wwnmeanp, &
         wwnmeanptail, fte, ftf, sstxftf, sstxftwn,&
         sstxftftail, sswellf, esin, ecos, aaircmin, &
         aairgb, aalpha, zzwnd, ssdsc
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
    USE w3odatmd, ONLY: ndst
#endif
    !
#ifdef W3_FLX5
    USE w3flx5md
#endif
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: A(NTH,NK), CG(NK), WN(NK), U, UDIR
#ifdef W3_FLX5
    REAL, INTENT(IN)        :: TAUA, TAUADIR, DAIR
#endif
    REAL, INTENT(IN)        :: TAUWX, TAUWY
    LOGICAL, INTENT(IN)     :: LLWS(NSPEC)
    REAL, INTENT(INOUT)     :: USTAR ,USDIR
    REAL, INTENT(OUT)       :: EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX,  &
         CD, Z0, CHARN, FMEANWS, DLWMEAN
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IS, IK, ITH
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
 
    REAL                    :: TAUW, EBAND, EMEANWS,UNZ,            &
         EB(NK),EB2(NK),ELCS, ELSN, SIGFAC
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SPR4')
#endif
    !
    unz    = max( 0.01 , u )
    ustar  = max( 0.0001 , ustar )
    !
    emean  = 0.
    emeanws= 0.
    fmeanws= 0.
    fmean  = 0.
    fmean1 = 0.
    wnmean = 0.
    amax   = 0.
    dlwmean =0.
    elcs =0.
    elsn =0.
    !
    ! 1.  Integral over directions and maximum --------------------------- *
    !
    DO ik=1, nk
      eb(ik)  = 0.
      eb2(ik) = 0.
      sigfac=sig(ik)**ssdsc(12) * dden(ik) / cg(ik)
      DO ith=1, nth
        is=ith+(ik-1)*nth
        eb(ik) = eb(ik) + a(ith,ik)
        elcs = elcs + a(ith,ik)*ecos(is)*sigfac
        elsn = elsn + a(ith,ik)*esin(is)*sigfac
        IF (llws(is)) eb2(ik) = eb2(ik) + a(ith,ik)
        amax   = max( amax , a(ith,ik) )
      END DO
    END DO
    !
    dlwmean=atan2(elsn,elcs)
    !
    ! 2.  Integrate over directions -------------------------------------- *
    !
    DO ik=1, nk
      eb(ik)   = eb(ik) * dden(ik) / cg(ik)
      eb2(ik)   = eb2(ik) * dden(ik) / cg(ik)
      emean    = emean  + eb(ik)
      fmean    = fmean  + eb(ik) /sig(ik)
      fmean1   = fmean1 + eb(ik) *(sig(ik)**(2.*wwnmeanptail))
      wnmean   = wnmean + eb(ik) *(wn(ik)**wwnmeanp)
      emeanws  = emeanws+ eb2(ik)
      fmeanws  = fmeanws+ eb2(ik)*(sig(ik)**(2.*wwnmeanptail))
    END DO
    !
    ! 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
    fmean1 = fmean1 + eband * sstxftftail
    wnmean = wnmean + eband * sstxftwn
    eband  = eb2(nk) / dden(nk)
    emeanws = emeanws + eband * fte
    fmeanws = fmeanws + eband * sstxftftail
    !
    ! 4.  Final processing
    !
    fmean  = tpiinv * emean / max( 1.e-7 , fmean )
    IF (fmean1.LT.1.e-7) THEN
      fmean1=tpiinv * sig(nk)
    ELSE
      fmean1  = tpiinv *( max( 1.e-7 , fmean1 )                       &
           / max( 1.e-7 , emean ))**(1/(2.*wwnmeanptail))
    ENDIF
    wnmean = ( max( 1.e-7 , wnmean )                              &
         / max( 1.e-7 , emean ) )**(1/wwnmeanp)
    IF (fmeanws.LT.1.e-7.OR.emeanws.LT.1.e-7) THEN
      fmeanws=tpiinv * sig(nk)
    ELSE
      fmeanws  = tpiinv *( max( 1.e-7 , fmeanws )                       &
           / max( 1.e-7 , emeanws ))**(1/(2.*wwnmeanptail))
    END IF
 
    !
    ! 5.  Cd and z0 ----------------------------------------------- *
    !
    tauw = sqrt(tauwx**2+tauwy**2)
    !
#ifdef W3_FLX5
    CALL w3flx5 ( zzwnd, u, udir, taua, tauadir, dair,  &
         ustar, usdir, z0, cd, charn )
#else
    CALL calc_ustar(u,tauw,ustar,z0,charn)
    unz    = max( 0.01 , u )
    cd     = (ustar/unz)**2
    usdir = udir
#endif
    !
    ! 6.  Final test output ---------------------------------------------- *
    !
#ifdef W3_T
    WRITE (ndst,9060) emean, wnmean, tpiinv, ustar, cd, z0
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9060 FORMAT (' TEST W3SPR4 : E,WN MN :',f8.3,f8.4/                   &
         ' TPIINV, USTAR, CD, Z0 :',f8.3,f7.2,1x,2f9.5)
#endif
    !/
    !/ End of W3SPR4 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3spr4
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE w3sin4 (A, CG, K, U, USTAR, DRAT, AS, USDIR, Z0, CD,    &
       TAUWX, TAUWY, TAUWNX, TAUWNY, S, D, LLWS,       &
       IX, IY, BRLAMBDA)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III                SHOM |
    !/                  !            F. Ardhuin             !
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         05-Dec-2013 |
    !/                  +-----------------------------------+
    !/
    !/    09-Oct-2007 : Origination.                        ( version 3.13 )
    !/    24-Jan-2013 : Adding breaking-related input       ( version 4.16 )
    !/    05-Dec-2013 : Cleaning up the ICE input           ( version 4.16 )
    !/
    !  1. Purpose :
    !
    !     Calculate diagonal and input source term for WAM4+ approach.
    !
    !  2. Method :
    !
    !       WAM-4     : Janssen et al.
    !       WAM-"4.5" : gustiness effect (Cavaleri et al. )
    !       SAT       : high-frequency input reduction for balance with
    !                   saturation dissipation (Ardhuin et al., 2008)
    !       SWELL     : negative wind input (Ardhuin et al. 2008)
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A       R.A.  I   Action density spectrum (1-D).
    !       CG      R.A.  I   Group speed                              *)
    !       K       R.A.  I   Wavenumber for entire spectrum.          *)
    !       U       Real  I   WIND SPEED
    !       USTAR   Real  I   Friction velocity.
    !       DRAT    Real  I   Air/water density ratio.
    !       AS      Real  I   Air-sea temperature difference
    !       USDIR   Real  I   wind stress direction
    !       Z0      Real  I   Air-side roughness lengh.
    !       CD      Real  I   Wind drag coefficient.
    !       USDIR   Real  I   Direction of friction velocity
    !       TAUWX-Y Real  I   Components of the wave-supported stress.
    !       TAUWNX  Real  I   Component of the negative wave-supported stress.
    !       TAUWNY  Real  I   Component of the negative wave-supported stress.
    !       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
    !
    !  4. Subroutines used :
    !
    !       STRACE    Subroutine tracing.                 ( !/S switch )
    !       PRT2DS    Print plot of spectrum.             ( !/T0 switch )
    !       OUTMAT    Print out matrix.                   ( !/T1 switch )
    !
    !  5. Called by :
    !
    !       W3SRCE   Source term integration.
    !       W3EXPO   Point output program.
    !       GXEXPO   GrADS point output program.
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable general test output.
    !     !/T0  2-D print plot of source term.
    !     !/T1  Print arrays.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: grav,nu_air,kappa,tpi,fwtable,sizefwtable, &
#ifdef w3_t
         rade,                                      &
#endif
         delab,abmin
    USE w3gdatmd, ONLY: nk, nth, nspec, dden, sig, sig2, th,         &
         esin, ecos, ec2, zzwnd, aalpha, bbeta, zzalp,&
         ttauwshelter, sswellf, dden2, dth, ssinthp,  &
         zz0rat, ssinbr, sintailpar
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
#ifdef W3_T0
    USE w3odatmd, ONLY: ndst
#endif
    USE w3odatmd, ONLY: iaproc
#ifdef W3_T0
    USE w3arrymd, ONLY: prt2ds
#endif
#ifdef W3_T1
    USE w3arrymd, ONLY: outmat
#endif
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: A(NSPEC), BRLAMBDA(NSPEC)
    REAL, INTENT(IN)        :: CG(NK), K(NSPEC),Z0,U, CD
    REAL, INTENT(IN)        :: USTAR, USDIR, AS, DRAT
    REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC), TAUWX, TAUWY, TAUWNX, TAUWNY
    LOGICAL, INTENT(OUT)    :: LLWS(NSPEC)
    INTEGER, INTENT(IN)     :: IX, IY
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IS,IK,ITH
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: FACLN1, FACLN2, LAMBDA
    REAL                    :: COSU, SINU, TAUX, TAUY, USDIRP, USTP
    REAL                    :: TAUPX, TAUPY, UST2, TAUW, TAUWB
    REAL   , PARAMETER      :: EPS1 = 0.00001, eps2 = 0.000001
    REAL                    :: Usigma           !standard deviation of U due to gustiness
    REAL                    :: USTARsigma       !standard deviation of USTAR due to gustiness
    REAL                    :: CM,UCN,ZCN, &
         Z0VISC, Z0NOZ, EB,  &
         EBX, EBY, AORB, AORB1, FW, UORB, TH2, &
         RE, FU, FUD, SWELLCOEFV, SWELLCOEFT
    REAL                   ::  PTURB, PVISC, SMOOTH
    REAL XI,DELI1,DELI2
    REAL XJ,DELJ1,DELJ2
    REAL XK,DELK1,DELK2
    REAL                    :: CONST, CONST0, CONST2, TAU1, TAU1NT, ZINF, TENSK
    REAL X,ZARG,ZLOG,UST
    REAL                    :: COSWIND, XSTRESS, YSTRESS, TAUHF
    REAL TEMP, TEMP2
    INTEGER IND,J,I,ISTAB
    REAL DSTAB(3,NSPEC), DVISC, DTURB
    REAL STRESSSTAB(3,2),STRESSSTABN(3,2)
    !
    INTEGER, PARAMETER      :: JTOT=50
    REAL   , PARAMETER      :: KM=363.,cmm=0.2325  ! K and C at phase speed minimum in rad/m
    REAL                    :: OMEGACC, OMEGA, ZZ0, ZX, ZBETA, USTR, TAUR,  &
         CONST1, LEVTAIL0, X0, Y, DELY, YC, ZMU,      &
         LEVTAIL, CGTAIL, ALPHAM, FM, ALPHAT, FMEAN
 
    REAL, ALLOCATABLE       :: W(:)
#ifdef W3_T0
    REAL                    :: DOUT(NK,NTH)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SIN4')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) bbeta, ustar, usdir*rade
#endif
    !
    ! 1.  Preparations
    !
    !JDM: Initializing values to zero, they shouldn't be used unless
    !set in another place, but seems to solve some bugs with certain
    !compilers.
    dstab =0.
    stressstab =0.
    stressstabn =0.
    !
    ! Coupling coefficient times density ratio DRAT
    !
    const1=bbeta/kappa**2  ! needed for the tail
    const0=const1*drat     ! needed for the resolved spectrum
    !
    ! 1.a  estimation of surface roughness parameters
    !
    z0visc = 0.1*nu_air/max(ustar,0.0001)
    z0noz = max(z0visc,zz0rat*z0)
    facln1 = u / log(zzwnd/z0noz)
    facln2 = log(z0noz)
    !
    ! 1.b  estimation of surface orbital velocity and displacement
    !
    uorb=0.
    aorb=0.
 
    DO ik=1, nk
      eb  = 0.
      ebx = 0.
      eby = 0.
      DO ith=1, nth
        is=ith+(ik-1)*nth
        eb  = eb  + a(is)
      END DO
      !
      !  At this point UORB and AORB are the variances of the orbital velocity and surface elevation
      !
      uorb = uorb + eb *sig(ik)**2 * dden(ik) / cg(ik)
      aorb = aorb + eb             * dden(ik) / cg(ik)  !correct for deep water only
    END DO
    !      FMEAN = SQRT((UORB+1E-6)/(AORB+1E-6))
    uorb  = 2*sqrt(uorb)                  ! significant orbital amplitude
    aorb1 = 2*aorb**(1-0.5*sswellf(6))    ! half the significant wave height ... if SWELLF(6)=1
    re = 4*uorb*aorb1 / nu_air           ! Reynolds number
    !
    ! Defines the swell dissipation based on the "Reynolds number"
    !
    IF (sswellf(4).GT.0) THEN
      IF (sswellf(7).GT.0.) THEN
        smooth = 0.5*tanh((re-sswellf(4))/sswellf(7))
        pturb=(0.5+smooth)
        pvisc=(0.5-smooth)
      ELSE
        IF (re.LE.sswellf(4)) THEN
          pturb =  0.
          pvisc =  1.
        ELSE
          pturb =  1.
          pvisc =  0.
        END IF
      END IF
    ELSE
      pturb=1.
      pvisc=1.
    END IF
 
    !
    IF (sswellf(2).EQ.0) THEN
      fw=max(abs(sswellf(3)),0.)
      fu=0.
      fud=0.
    ELSE
      fu=abs(sswellf(3))
      fud=sswellf(2)
      aorb=2*sqrt(aorb)
      xi=(alog10(max(aorb/z0noz,3.))-abmin)/delab
      ind  = min(sizefwtable-1, int(xi))
      deli1= min(1. ,xi-float(ind))
      deli2= 1. - deli1
      fw =fwtable(ind)*deli2+fwtable(ind+1)*deli1
    END IF
    !
    ! 2.  Diagonal
    !
    ! Here AS is the air-sea temperature difference in degrees. Expression given by
    ! Abdalla & Cavaleri, JGR 2002 for Usigma. For USTARsigma ... I do not see where
    ! I got it from, maybe just made up from drag law ...
    !
#ifdef W3_STAB3
    IF ( isnan(as) ) THEN
      ! AS is typically NaN on land and can propagate into the domain by interpolation
      usigma = 0.
    ELSE
      usigma = max(0.,-0.025*as)
    END IF
    ustarsigma=(1.0+u/(10.+u))*usigma
#endif
#ifdef W3_T
    WRITE (ndst,9003) as, usigma, ustarsigma, u
#endif
    ust=ustar
    istab=3
#ifdef W3_STAB3
    DO istab=1,2
      IF (istab.EQ.1) ust=ustar*(1.-ustarsigma)
      IF (istab.EQ.2) ust=ustar*(1.+ustarsigma)
#endif
      taux = ust**2* cos(usdir)
      tauy = ust**2* sin(usdir)
#ifdef W3_T
      WRITE (ndst,9001) istab, taux, tauy, ust
#endif
      !
      ! Loop over the resolved part of the spectrum
      !
      stressstab(istab,:)=0.
      stressstabn(istab,:)=0.
      !
      DO ik=1, nk
        taupx=taux-abs(ttauwshelter)*stressstab(istab,1)
        taupy=tauy-abs(ttauwshelter)*stressstab(istab,2)
        ! With MIN and MAX the bug should disappear.... but where did it come from?
        ustp=min((taupx**2+taupy**2)**0.25,max(ust,0.3))
        usdirp=atan2(taupy,taupx)
        cosu   = cos(usdirp)
        sinu   = sin(usdirp)
        is=1+(ik-1)*nth
        cm=k(is)/sig2(is) !inverse of phase speed
        ucn=ustp*cm+zzalp  !this is the inverse wave age
        ! the stress is the real stress (N/m^2) divided by
        ! rho_a, and thus comparable to USTAR**2
        ! it is the integral of rho_w g Sin/C /rho_a
        ! (air-> waves momentum flux)
        const2=dden2(is)/cg(ik) &        !Jacobian to get energy in band
             *grav/(sig(ik)/k(is)*drat) ! coefficient to get momentum
        const=sig2(is)*const0
        ! CM parameter is 1 / C_phi
        ! Z0 corresponds to Z0+Z1 of the Janssen eq. 14
        zcn=alog(k(is)*z0)
        !
        ! precomputes swell factors
        !
        swellcoefv=-sswellf(5)*drat*2*k(is)*sqrt(2*nu_air*sig2(is))
        swellcoeft=-drat*sswellf(1)*16*sig2(is)**2/grav
        !
        DO ith=1,nth
          is=ith+(ik-1)*nth
          coswind=(ecos(is)*cosu+esin(is)*sinu)
          IF (coswind.GT.0.01) THEN
            x=coswind*ucn
            ! this ZARG term is the argument of the exponential
            ! in Janssen 1991 eq. 16.
            zarg=kappa/x
            ! ZLOG is ALOG(MU) where MU is defined by Janssen 1991 eq. 15
            ! MU=
            zlog=zcn+zarg
 
            IF (zlog.LT.0.) THEN
              ! The source term Sp is beta * omega * X**2
              ! as given by Janssen 1991 eq. 19
              ! Note that this is slightly diffent from ECWAM code CY45R2 where ZLOG is replaced by ??
              dstab(istab,is) = const*exp(zlog)*zlog**4*ucn*ucn*coswind**ssinthp
 
              ! Below is an example with breaking probability feeding back to the input...
              !DSTAB(ISTAB,IS) = CONST*EXP(ZLOG)*ZLOG**4  &
              !                  *UCN*UCN*COSWIND**SSINTHP *(1+BRLAMBDA(IS)*20*SSINBR)
              llws(is)=.true.
            ELSE
              dstab(istab,is) = 0.
              llws(is)=.false.
            END IF
            !
            !  Added for consistency with ECWAM implsch.F
            !
            IF (28.*cm*ustar*coswind.GE.1) THEN
              llws(is)=.true.
            END IF
          ELSE  ! (COSWIND.LE.0.01)
            dstab(istab,is) = 0.
            llws(is)=.false.
          END IF
          !
          IF ((sswellf(1).NE.0.AND.dstab(istab,is).LT.1e-7*sig2(is)) &
               .OR.sswellf(3).GT.0) THEN
            !
            dvisc=swellcoefv
            dturb=swellcoeft*(fw*uorb+(fu+fud*coswind)*ustp)
            !
            dstab(istab,is) = dstab(istab,is) + pturb*dturb +  pvisc*dvisc
          END IF
          !
          ! Sums up the wave-supported stress
          !
          ! Wave direction is "direction to"
          ! therefore there is a PLUS sign for the stress
          temp2=const2*dstab(istab,is)*a(is)
          IF (dstab(istab,is).LT.0) THEN
            stressstabn(istab,1)=stressstabn(istab,1)+temp2*ecos(is)
            stressstabn(istab,2)=stressstabn(istab,2)+temp2*esin(is)
          ELSE
            stressstab(istab,1)=stressstab(istab,1)+temp2*ecos(is)
            stressstab(istab,2)=stressstab(istab,2)+temp2*esin(is)
          END IF
        END DO
      END DO
      !
      d(:)=dstab(3,:)
      xstress=stressstab(3,1)
      ystress=stressstab(3,2)
      tauwnx =stressstabn(3,1)
      tauwny =stressstabn(3,2)
#ifdef W3_STAB3
    END DO
    d(:)=0.5*(dstab(1,:)+dstab(2,:))
    xstress=0.5*(stressstab(1,1)+stressstab(2,1))
    ystress=0.5*(stressstab(1,2)+stressstab(2,2))
    tauwnx=0.5*(stressstabn(1,1)+stressstabn(2,1))
    tauwny=0.5*(stressstabn(1,2)+stressstabn(2,2))
#endif
#ifdef W3_T
    WRITE (ndst,9002) sum(d), sum(a), xstress, ystress, tauwnx, tauwny
#endif
    s = d * a
    !
    ! ... Test output of arrays
    !
#ifdef W3_T0
    DO ik=1, nk
      DO ith=1, nth
        dout(ik,ith) = d(ith+(ik-1)*nth)
      END DO
    END DO
    CALL prt2ds (ndst, nk, nk, nth, dout, sig(1), '  ', 1.,         &
         0.0, 0.001, 'Diag Sin', ' ', 'NONAME')
#endif
    !
#ifdef W3_T1
    CALL outmat (ndst, d, nth, nth, nk, 'diag Sin')
#endif
    !
    taupx=taux-abs(ttauwshelter)*xstress
    taupy=tauy-abs(ttauwshelter)*ystress
    ustp=(taupx**2+taupy**2)**0.25
    usdirp=atan2(taupy,taupx)
 
    ust=ustp
    !
    ! Computes HF tail
    !
    ! Computes the high-frequency contribution
    ! the difference in spectal density (kx,ky) to (f,theta)
    ! is integrated in this modified CONST0
    const0=dth*sig(nk)**5/((grav**2)*tpi) &
         *tpi*sig(nk) / cg(nk)  !conversion WAM (E(f,theta) to WW3 A(k,theta)
    temp=0.
    DO ith=1,nth
      is=ith+(nk-1)*nth
      coswind=(ecos(is)*cosu+esin(is)*sinu)
      temp=temp+a(is)*(max(coswind,0.))**3
    END DO
    !
    levtail0= const0*temp  ! LEVTAIL is sum over theta of A(k,theta)*cos^3(theta-wind)*DTH*SIG^5/(g^2*2pi)*2*pi*SIG/CG
                           !  which is the same as sum of E(f,theta)*cos^3(theta-wind)*DTH*SIG^5/(g^2*2pi)
                           ! reminder:  sum of E(f,theta)*DTH*SIG^5/(g^2*2pi) is 2*k^3*E(k)
!
! Computation of stress supported by tail: uses table if SINTAILPAR(1)=1 , correspoding to SINTABLE = 1
!
    IF (sintailpar(1).LT.0.5) THEN
      ALLOCATE(w(jtot))
      w(2:jtot-1)=1.
      w(1)=0.5
      w(jtot)=0.5
      x0 = 0.05
      !
      ustr= ust
      zz0=z0
      omegacc  = max(sig(nk),x0*grav/ust)
      yc       = omegacc*sqrt(zz0/grav)
 
      ! DELY     = MAX((1.-YC)/REAL(JTOT),0.)
      ! Changed integration variable from Y to LOG(Y) and to log(K)
      !ZINF      = LOG(YC)
      !DELY     = MAX((1.-ZINF)/REAL(JTOT),0.)
      zinf = log(sig(nk)**2/grav)
      dely = (log(tpi/0.005)-zinf)/real(jtot)
 
      taur=ust**2
      tau1=0.
 
      ! Integration loop over the tail wavenumbers or frequencies ...
      DO j=1,jtot
        !Y        = YC+REAL(J-1)*DELY
        !OMEGA    = Y*SQRT(GRAV/ZZ0)
        !OMEGA    = SQRT(GRAV*Y)
        ! This is the deep water phase speed... No surface tension !!
        !CM       = GRAV/OMEGA
        ! With this form, Y is the wavenumber in the tail;
        y= exp(zinf+real(j-1)*dely)
        tensk    =1+(y/km)**2
        omega    = sqrt(grav*y*tensk)
        cm       = sqrt(grav*tensk/y)
        cgtail   = 0.5*(3*(y/km)**2+1)*sqrt(grav/(y*tensk))
        !this is the inverse wave age, shifted by ZZALP (tuning)
        zx       = ustr/cm +zzalp
        zarg     = min(kappa/zx,20.)
        ! ZMU corresponds to EXP(ZCN)
        zmu      = min(grav*zz0/cm**2*exp(zarg),1.)
        zlog     = min(alog(zmu),0.)
        zbeta        = const1*zmu*zlog**4
        !
        ! Optional addition of capillary wave peak if SINTAIL2=1 
        !
        IF (sintailpar(3).GT.0) THEN
          IF (ustr.LT.cm) THEN
            alpham=max(0.,0.01*(1.+alog(ustr/cm)))
          ELSE
            alpham=0.01*(1+3.*alog(ustr/cm))
          END IF
          fm=exp(-0.25*(y/km-1)**2)
 
          alphat=alpham*(cmm/cm)*fm  ! equivalent to 2*Bh in Elfouhaily et al.
          levtail=levtail0*0.5*(1-tanh((y-20)/5))+sintailpar(3)*0.5*(1+tanh((y-20)/5))*alphat
        ELSE
          levtail=levtail0
        END IF
        ! WRITE(991,*) 'TAIL??',SINTAILPAR(3),LEVTAIL0,LEVTAIL,ALPHAT,Y,Y/KM,OMEGA/(TPI)
 
        !TAU1=TAU1+W(J)*ZBETA*(USTR/UST)**2/Y*DELY              ! integration over LOG(Y)
        tau1=tau1+w(j)*zbeta*ustr**2*levtail*dely*cgtail/cm       ! integration over LOG(K)
 
        ! NB: the factor ABS(TTAUWSHELTER) was forgotten in the TAUHFT2 table
        !TAUR=TAUR-W(J)*ABS(TTAUWSHELTER)*USTR**2*ZBETA*LEVTAIL/Y*DELY
        !TAUR=TAUR-W(J)*USTR**2*ZBETA*LEVTAIL*DELY    ! integration over LOG(Y)
        taur=taur-w(j)*sintailpar(2)*ustr**2*zbeta*levtail*dely*cgtail/cm   ! DK/K*CG/C = D OMEGA / OMEGA
        ustr=sqrt(max(taur,0.))
      END DO
      DEALLOCATE(w)
      tau1nt=tau1
      tauhf = tau1
      !
      ! In this case, uses tables for high frequency contribution to TAUW.
      !
    ELSE
      ! finds the values in the tabulated stress TAUHFT
      xi=ust/delust
      ind  = max(1,min(iustar-1, int(xi)))
      deli1= max(min(1. ,xi-float(ind)),0.)
      deli2= 1. - deli1
      xj=max(0.,(grav*z0/max(ust,0.00001)**2-aalpha) / delalp)
      j    = max(1 ,min(ialpha-1, int(xj)))
      delj1= max(0.,min(1.      , xj-float(j)))
      delj2=1. - delj1
      IF (ttauwshelter.GT.0) THEN
        xk = levtail0/ deltail
        i = min(ilevtail-1, int(xk))
        delk1= min(1. ,xk-float(i))
        delk2=1. - delk1
        tau1 =((tauhft2(ind,j,i)*deli2+tauhft2(ind+1,j,i)*deli1 )*delj2 &
             +(tauhft2(ind,j+1,i)*deli2+tauhft2(ind+1,j+1,i)*deli1)*delj1)*delk2 &
             +((tauhft2(ind,j,i+1)*deli2+tauhft2(ind+1,j,i+1)*deli1 )*delj2 &
             +(tauhft2(ind,j+1,i+1)*deli2+tauhft2(ind+1,j+1,i+1)*deli1)*delj1)*delk1
      ELSE
        tau1 =(tauhft(ind,j)*deli2+tauhft(ind+1,j)*deli1 )*delj2 &
             +(tauhft(ind,j+1)*deli2+tauhft(ind+1,j+1)*deli1)*delj1
      END IF
      !
      tauhf = levtail0*ust**2*tau1
    END IF ! End of test on use of table
 
    tauwx = xstress+tauhf*cos(usdirp)
    tauwy = ystress+tauhf*sin(usdirp)
    !
    ! Reduces tail effect to make sure that wave-supported stress
    ! is less than total stress, this is borrowed from ECWAM Stresso.F
    !
    tauw = sqrt(tauwx**2+tauwy**2)
    ust2   = max(ustar,eps2)**2
    tauwb = min(tauw,max(ust2-eps1,eps2**2))
    IF (tauwb.LT.tauw) THEN
      tauwx=tauwx*tauwb/tauw
      tauwy=tauwy*tauwb/tauw
    END IF
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3SIN4 : COMMON FACT.: ',3e10.3)
9001 FORMAT (' TEST W3SIN4 : ISTAB      :',i2/                    &
             '               TAUX       :',e12.3/                 &
             '               TAUY       :',e12.3/                 &
             '               UST        :',e12.3)
9002 FORMAT (' TEST W3SIN4 : SUM(D)     :',e12.3/                 &
             '               SUM(A)     :',e12.3/                 &
             '               STRESSX    :',e12.3/                 &
             '               STRESSY    :',e12.3/                 &
             '               TAUWNX     :',e12.3/                 &
             '               TAUWNY     :',e12.3)
9003 FORMAT (' TEST W3SIN4 : AS         :',f8.4/                  &
             '               Usigma     :',e12.3/                 &
             '               USTARsigma :',e12.3/                 &
             '               U          :',e12.3)
 
#endif
    !/
    !/ End of W3SIN4 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3sin4
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE insin4(FLTABS)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |                         SHOM      |
    !/                  |            F. Ardhuin             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         30-Aug-2010 |
    !/                  +-----------------------------------+
    !/
    !/    30-Aug-2010 : Origination.                        ( version 3.14-Ifremer )
    !
    !  1. Purpose :
    !
    !     Initialization for source term routine.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     ----------------------------------------------------------------
    !      FLTABS    Logical
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SIN4    Subr. W3SRC3MD Corresponding source term.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: tpiinv, rade, grav
    USE w3odatmd,  ONLY: ndse
    USE w3servmd,  ONLY: extcde
    USE w3dispmd,  ONLY: wavnu2
    USE w3gdatmd,  ONLY: sig, dsip, nk, nth, ttauwshelter,             &
         ssdsdth, ssdscos, th, dth, xfr, ecos, esin,   &
         ssdsc,  ssdsbrf1, ssdsbck, ssdsbint, ssdspbk, &
         ssdsabk, ssdshck, iktab, dcki, satindices,    &
         satweights, cumulw, nkhs, nkd, ndtab, qbi,    &
         sintailpar
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    LOGICAL, INTENT(IN)     :: FLTABS
    !/
    !/ ------------------------------------------------------------------- /
    !/
    INTEGER  SDSNTH, ITH, I_INT, J_INT, IK, IK2, ITH2 , IS, IS2
    INTEGER  IKL, ID, ICON, IKD, IKHS, IKH, TOTO
    REAL     C, C2
    REAL     DIFF1, DIFF2, BINF, BSUP, CGG, PROF
    REAL     KIK, DHS, KD, KHS, KH, XT, GAM, DKH, PR, W, EPS
    REAL     DKD
    REAL, DIMENSION(:,:)   , ALLOCATABLE :: SIGTAB
    REAL, DIMENSION(:,:)   , ALLOCATABLE :: K1, K2
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'INSIN4')
#endif
    !
    ! 1.  Initializations ------------------------------------------------ *
    !
    !
    ! These precomputed tables are written in mod_def.ww3
    !
    IF (sintailpar(1).GT.0.5) THEN
      IF (.NOT. ALLOCATED(taut)) ALLOCATE(taut(0:itaumax,0:jumax))
      IF (.NOT. ALLOCATED(tauhft)) ALLOCATE(tauhft(0:iustar,0:ialpha))
      IF (fltabs) THEN
        CALL tabu_stress
        CALL tabu_tauhf(sig(nk) )      !tabulate high-frequency stress: 2D table
      END IF
      IF (ttauwshelter.GT.0) THEN
        IF (.NOT. ALLOCATED(tauhft2)) ALLOCATE(tauhft2(0:iustar,0:ialpha,0:ilevtail))
        IF (fltabs) CALL tabu_tauhf2(sig(nk) )   !tabulate high-frequency stress: 3D table
      END IF
    END IF
    !
    ! 2.  SPONTANEOUS BREAKING
    ! 2.a Precomputes the indices for integrating the spectrum to get saturation (TEST 4xx )
    !
    IF (ssdsdth.LT.180) THEN
      sdsnth  = min(nint(ssdsdth/(dth*rade)),nth/2-1)
      satindices(:,:)=1
      satweights(:,:)=0.
      DO ith=1,nth
        DO i_int=ith-sdsnth, ith+sdsnth
          j_int=i_int
          IF (i_int.LT.1)  j_int=i_int+nth
          IF (i_int.GT.nth) j_int=i_int-nth
          satindices(i_int-(ith-sdsnth)+1,ith)=j_int
          satweights(i_int-(ith-sdsnth)+1,ith)=          &
               cos(th(ith)-th(j_int))**ssdscos
        END DO
      END DO
    ELSE
      satindices(:,:)=1
      satweights(:,:)=1.
    END IF
    !/ ------------------------------------------------------------------- /
    !
    ! Precomputes QBI and DCKI (TEST 500)
    !
    IF (ssdsbck.GT.0) THEN
      !
      ! Precomputes the indices for integrating the spectrum over frequency bandwidth
      !
      binf=(1-ssdsbint) ! Banner et al 2002: Hp=4*sqrt(int_0.7^1.3fp E df), SSDSBINT=0.3
      bsup=(1+ssdsbint)
      kik=0.
      !
      ! High frequency tail for convolution calculation
      !
      ALLOCATE(k1(nk,ndtab))
      ALLOCATE(k2(nk,ndtab))
      ALLOCATE(sigtab(nk,ndtab))
 
      sigtab=0. !contains frequency for upper windows boundaries
      iktab=0  ! contains indices for upper windows boundaries
 
      DO id=1,ndtab
        toto=0
        prof=real(id)
        DO ikl=1,nk ! last window starts at IK=NK
          CALL wavnu2(sig(ikl), prof, kik, cgg, 1e-7, 15, icon)
          k1(ikl,id)=kik  ! wavenumber lower boundary (is directly related to the frequency indices, IK)
          k2(ikl,id)=((bsup/binf)**2.)*k1(ikl,id)! wavenumber upper boundary
          sigtab(ikl,id)=sqrt(grav*k2(ikl,id)*tanh(k2(ikl,id)*id)) ! corresponding frequency upper boundary
          IF(sigtab(ikl,id) .LE. sig(1)) THEN
            iktab(ikl,id)=1
          END IF
          IF(sigtab(ikl,id) .GT. sig(nk)) THEN
            iktab(ikl,id)=nk+toto       ! in w3sds4 only windows with IKSUP<=NK will be kept
            toto=1
          END IF
          DO ik=1,nk-1
            diff1=0.
            diff2=0.
            IF(sig(ik)<sigtab(ikl,id) .AND. sig(ik+1)>=sigtab(ikl,id)) THEN
              diff1=sigtab(ikl,id)-sig(ik)   ! seeks the indices of the upper boundary
              diff2=sig(ik+1)-sigtab(ikl,id)! the indices of lower boudary = IK
              IF (diff1<diff2) THEN
                iktab(ikl,id)=ik
              ELSE
                iktab(ikl,id)=ik+1
              END IF
            END IF
          END DO
        END DO
      END DO
      !
      ! Tabulates DCKI and QBI
      !
      dhs=khsmax/nkhs ! max value of KHS=KHSMAX
      dkh=khmax/nkhi  ! max value of KH=KHMAX
      dkd=kdmax/nkd
      ALLOCATE(dcki(nkhs,nkd))
      ALLOCATE(qbi(nkhs,nkd))
      dcki=0.
      qbi =0.
      DO ikd=1,nkd
        khs=0.
        kd=(fac_kd1**(ikd-fac_kd2))
        xt=tanh(kd)
        gam=1.0314*(xt**3)-1.9958*(xt**2)+1.5522*xt+0.1885
        gam=gam/2.15
        DO ikhs=1,nkhs  ! max value of KHS=1.
          kh=0.
          khs=khs+dhs
          DO ikh=1,nkhi
            kh=kh+dkh
            pr=(4.*kh/(khs**2.))*exp(-(2*((kh/khs)**2.)))
            !              W=1.5*(((KHS)/(SQRT(2.)*GAM*XT))**2.)*(1-exp(-(((KH)/(GAM*XT))**4.))) !CK2002 parameterization
            w=ssdsabk*(((khs)/(sqrt(2.)*gam*xt))**2.)*(1-exp(-(((kh)/(gam*xt))**ssdspbk)))
            eps=-((((ssdsbck/(xt**ssdshck))*kh)**3.)/4)*sqrt(grav/xt)
            dcki(ikhs, ikd)= dcki(ikhs, ikd)+pr*w*eps*dkh
            qbi(ikhs, ikd) = qbi(ikhs, ikd) +pr*w*    dkh
          END DO
        END DO
      END DO
 
      WHERE ( qbi .GT. 1. )
        qbi = 1.
      END WHERE
 
      DEALLOCATE(k1,k2)
      DEALLOCATE(sigtab)
    ELSE
      iktab(:,:)=1
      dcki(:,:) =0.
      qbi(:,:)  =0.
    END IF
    !
    !/ ------------------------------------------------------------------- /
    !                        CUMULATIVE EFFECT
    !/ ------------------------------------------------------------------- /
    !
    ! Precomputes the weights for the cumulative effect (TEST 441 and 500)
    !
    dikcumul = 0
    IF (ssdsc(3).LT.0.) THEN
      !       DIKCUMUL is the integer difference in frequency bands
      !       between the "large breakers" and short "wiped-out waves"
      dikcumul = nint(ssdsbrf1/(xfr-1.))
      !        WRITE(6,*) 'INSIN4b:',DIKCUMUL
      cumulw(:,:)=0.
      DO ik=1,nk
        c = grav/sig(ik)   ! Valid in deep water only
        !C = SIG(IK)/K(IK) ! Valid in all water depth ???
        DO ith=1,nth
          is=ith+(ik-1)*nth
          DO ik2=1,ik-dikcumul
            c2 = grav/sig(ik2) ! Valid in deep water only
            !C2 = SIG(IK2)/K(IK2) ! Valid in all water depth ???
            DO ith2=1,nth
              is2=ith2+(ik2-1)*nth
              cumulw(is2,is)=sqrt(c**2+c2**2-2*c*c2*ecos(1+abs(ith2-ith))) & ! = deltaC
                   *dsip(ik2)/(0.5*c2) * dth                   ! = dk*dtheta (Valid in deep water only)
            END DO
          END DO
        END DO
      END DO
    ELSE
      cumulw(:,:)=0.
    END IF
    !/
    !/ End of INSIN4 ----------------------------------------------------- /
    !/
  END SUBROUTINE insin4
  ! ----------------------------------------------------------------------
 
  SUBROUTINE tabu_stress
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |            F. Ardhuin             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         17-Oct-2007 |
    !/                  +-----------------------------------+
    !/
    !/    23-Jun-2006 : Origination.                        ( version 3.13 )
    !/     adapted from WAM, original:P.A.E.M. JANSSEN    KNMI AUGUST 1990
    !/     adapted version (subr. STRESS): J. BIDLOT    ECMWF OCTOBER 2004
    !/     Table values were checkes against the original f90 result and found to
    !/     be identical (at least at 0.001 m/s accuracy)
    !/
    !  1. Purpose :
    !     TO GENERATE friction velocity table TAUT(TAUW,U10)=SQRT(TAU).
    !     METHOD.
    !       A STEADY STATE WIND PROFILE IS ASSUMED.
    !       THE WIND STRESS IS COMPUTED USING THE ROUGHNESSLENGTH
    !                  Z1=Z0/SQRT(1-TAUW/TAU)
    !       WHERE Z0 IS THE CHARNOCK RELATION , TAUW IS THE WAVE-
    !       INDUCED STRESS AND TAU IS THE TOTAL STRESS.
    !       WE SEARCH FOR STEADY-STATE SOLUTIONS FOR WHICH TAUW/TAU < 1.
    !       FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990.
    !
    !     Initialization for source term routine.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SIN3    Subr. W3SRC3MD Corresponding source term.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S  Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: kappa, grav
    USE w3gdatmd, ONLY: zzwnd, aalpha, zz0max
    IMPLICIT NONE
    INTEGER, PARAMETER      :: NITER=10
    REAL   , PARAMETER      :: XM=0.50, eps1=0.00001
    !     VARIABLE.   TYPE.     PURPOSE.
    !      *XM*        REAL      POWER OF TAUW/TAU IN ROUGHNESS LENGTH.
    !      *XNU*       REAL      KINEMATIC VISCOSITY OF AIR.
    !      *NITER*     INTEGER   NUMBER OF ITERATIONS TO OBTAIN TOTAL STRESS
    !      *EPS1*      REAL      SMALL NUMBER TO MAKE SURE THAT A SOLUTION
    !                            IS OBTAINED IN ITERATION WITH TAU>TAUW.
    ! ----------------------------------------------------------------------
    INTEGER I,J,ITER
    REAL ZTAUW,UTOP,CDRAG,WCD,USTOLD,TAUOLD
    REAL X,UST,ZZ0,F,DELF,ZZ00
    !
    !
    delu    = umax/float(jumax)
    deltauw = tauwmax/float(itaumax)
    DO i=0,itaumax
      ztauw   = (real(i)*deltauw)**2
      DO j=0,jumax
        utop    = float(j)*delu
        cdrag   = 0.0012875
        wcd     = sqrt(cdrag)
        ustold  = utop*wcd
        tauold  = max(ustold**2, ztauw+eps1)
        DO iter=1,niter
          x   = ztauw/tauold
          ust = sqrt(tauold)
          zz00=aalpha*tauold/grav
          IF (zz0max.NE.0) zz00=min(zz00,zz0max)
          ! Corrects roughness ZZ00 for quasi-linear effect
          zz0 = zz00/(1.-x)**xm
          !ZNU = 0.1*nu_air/UST  ! This was removed by Bidlot in 1996
          !ZZ0 = MAX(ZNU,ZZ0)
          f   = ust-kappa*utop/(alog(zzwnd/zz0))
          delf= 1.-kappa*utop/(alog(zzwnd/zz0))**2*2./ust &
               *(1.-(xm+1)*x)/(1.-x)
          ust = ust-f/delf
          tauold= max(ust**2., ztauw+eps1)
        END DO
        taut(i,j)  = sqrt(tauold)
      END DO
    END DO
    i=itaumax
    j=jumax
    !
    !  Force zero wind to have zero stress (Bidlot 1996)
    !
    DO i=0,itaumax
      taut(i,0)=0.0
    END DO
    RETURN
  END SUBROUTINE tabu_stress
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE tabu_tauhf(SIGMAX)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |            F. Ardhuin             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update 2006/08/14            |
    !/                  +-----------------------------------+
    !/
    !/    27-Feb-2004 : Origination in WW3                  ( version 2.22.SHOM )
    !/     the resulting table was checked to be identical to the original f77 result
    !/    14-Aug-2006 : Modified following Bidlot           ( version 2.22.SHOM )
    !/    18-Aug-2006 : Ported to version 3.09
    !
    !  1. Purpose :
    !
    !     Tabulation of the high-frequency wave-supported stress
    !
    !  2. Method :
    !
    !       SEE REFERENCE FOR WAVE STRESS CALCULATION.
    !       FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990.
    !     See tech. Memo ECMWF 03 december 2003 by Bidlot & Janssen
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       SIGMAX   Real  I   maximum frequency * TPI
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3SIN3   Wind input Source term routine.
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !       !/T      Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: kappa, grav
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    USE w3gdatmd, ONLY: aalpha, bbeta, zzalp, fachfe, zz0max
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, intent(in) :: SIGMAX  !  maximum frequency
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    !       USTARM  R.A.  Maximum friction velocity
    !       ALPHAM  R.A.  Maximum Charnock Coefficient
    !       WLV     R.A.  Water levels.
    !       UA      R.A.  Absolute wind speeds.
    !       UD      R.A.  Absolute wind direction.
    !       U10     R.A.  Wind speed used.
    !       U10D    R.A.  Wind direction used.
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: USTARM, ALPHAM
    REAL                    :: CONST1, OMEGA, OMEGAC
    REAL                    :: UST, ZZ0,OMEGACC, CM
    INTEGER, PARAMETER      :: JTOT=250
    REAL, ALLOCATABLE       :: W(:)
    REAL                    :: ZX,ZARG,ZMU,ZLOG,ZZ00,ZBETA
    REAL                    :: Y,YC,DELY
    INTEGER                 :: J,K,L
    REAL                    :: X0
    !
#ifdef W3_S
    CALL strace (ient, 'TABU_HF')
#endif
    !
    ustarm = 5.
    alpham = 20.*aalpha
    delust = ustarm/real(iustar)
    delalp = alpham/real(ialpha)
    const1 = bbeta/kappa**2
    omegac = sigmax
    !
    tauhft(0:iustar,0:ialpha)=0. !table initialization
    !
    ALLOCATE(w(jtot))
    w(2:jtot-1)=1.
    w(1)=0.5
    w(jtot)=0.5
    x0 = 0.05
    !
    DO l=0,ialpha
      DO k=0,iustar
        ust      = max(real(k)*delust,0.000001)
        zz00       = ust**2*aalpha/grav
        IF (zz0max.NE.0) zz00=min(zz00,zz0max)
        zz0       = zz00*(1+float(l)*delalp/aalpha)
        omegacc  = max(omegac,x0*grav/ust)
        yc       = omegacc*sqrt(zz0/grav)
        dely     = max((1.-yc)/real(jtot),0.)
        ! For a given value of UST and ALPHA,
        ! the wave-supported stress is integrated all the way
        ! to 0.05*g/UST
        DO j=1,jtot
          y        = yc+real(j-1)*dely
          omega    = y*sqrt(grav/zz0)
          ! This is the deep water phase speed
          cm       = grav/omega
          !this is the inverse wave age, shifted by ZZALP (tuning)
          zx       = ust/cm +zzalp
          zarg     = min(kappa/zx,20.)
          zmu      = min(grav*zz0/cm**2*exp(zarg),1.)
          zlog     = min(alog(zmu),0.)
          zbeta        = const1*zmu*zlog**4
          ! Power of Y in denominator should be FACHFE-4  tail applied here
          tauhft(k,l)  = tauhft(k,l)+w(j)*zbeta/y*dely
        END DO
#ifdef W3_T
        WRITE (ndst,9000) l,k,aalpha+float(l)*delalp,ust,tauhft(k,l)
#endif
      END DO
    END DO
    DEALLOCATE(w)
    RETURN
#ifdef W3_T
9000 FORMAT ('TABU_HF, L, K, ALPHA, UST, TAUHFT(K,L) :',(2i4,3f8.3))
#endif
  END SUBROUTINE tabu_tauhf
 
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE tabu_tauhf2(SIGMAX)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |            F. Ardhuin             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update 2006/08/14            |
    !/                  | Last update 2013/01/24            |
    !/                  +-----------------------------------+
    !/
    !/    15-May-2007 : Origination in WW3                  ( version 3.10.SHOM )
    !/    24-Jan-2013 : Allows to read in table             ( version 4.08 )
    !
    !  1. Purpose :
    !
    !     Tabulation of the high-frequency wave-supported stress as a function of
    !     ustar, alpha (modified Charnock), and tail energy level
    !
    !  2. Method :
    !
    !       SEE REFERENCE FOR WAVE STRESS CALCULATION.
    !       FOR QUASILINEAR EFFECT SEE PETER A.E.M. JANSSEN,1990.
    !     See tech. Memo ECMWF 03 december 2003 by Bidlot & Janssen
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       SIGMAX   Real  I   maximum frequency*TPI
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3SIN3   Wind input Source term routine.
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !       !/T      Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: kappa, grav, file_endian
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    USE w3gdatmd, ONLY: aalpha, bbeta, zzalp, fachfe,  &
         ttauwshelter, zz0max
    USE w3odatmd, ONLY: ndse
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, intent(in) :: SIGMAX  !  maximum frequency * TPI
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    !       USTARM  R.A.  Maximum friction velocity
    !       ALPHAM  R.A.  Maximum Charnock Coefficient
    !       WLV     R.A.  Water levels.
    !       UA      R.A.  Absolute wind speeds.
    !       UD      R.A.  Absolute wind direction.
    !       U10     R.A.  Wind speed used.
    !       U10D    R.A.  Wind direction used.
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: USTARM, ALPHAM, LEVTAILM
    REAL                    :: CONST1, OMEGA, OMEGAC, LEVTAIL
    REAL                    :: UST, UST0, ZZ0,OMEGACC, CM
    REAL                    :: TAUW, TAUW0
    INTEGER, PARAMETER      :: JTOT=250
    REAL, ALLOCATABLE       :: W(:)
    REAL                    :: ZX,ZARG,ZMU,ZLOG,ZBETA
    REAL                    :: Y,YC,DELY
    INTEGER                 :: I, J, K, L
    REAL                    :: X0, INSIGMAX, INAALPHA, INBBETA, INZZALP, INKAPPA, INGRAV
    INTEGER                 :: INIUSTAR, INIALPHA, INILEVTAIL, IERR
    CHARACTER(160)          :: FNAMETAB
    LOGICAL                 :: NOFILE
    CHARACTER(LEN=10), PARAMETER :: VERGRD = '2018-06-08'
    CHARACTER(LEN=35), PARAMETER :: IDSTR = 'WAVEWATCH III ST4 TABLE FOR STRESS '
    CHARACTER(LEN=10)       :: VERTST=' '
    CHARACTER(LEN=35)       :: IDTST=' '
    !
#ifdef W3_S
    CALL strace (ient, 'TABU_HF')
#endif
    !
    fnametab='ST4TABUHF2.bin'
    nofile=.true.
    OPEN (993,file=fnametab,form='UNFORMATTED', convert=file_endian,iostat=ierr,status='OLD')
    IF (ierr.EQ.0) THEN
      READ(993,iostat=ierr) idtst, vertst, insigmax, inaalpha, inbbeta, iniustar,  &
           inialpha, inilevtail, inzzalp, inkappa, ingrav
      IF (vertst.EQ.vergrd.AND.idtst.EQ.idstr.AND.ierr.EQ.0             &
           .AND.insigmax.EQ.sigmax.AND.inaalpha.EQ.aalpha.AND.inbbeta.EQ.bbeta) THEN
        IF (iniustar.EQ.iustar.AND.inialpha.EQ.ialpha.AND.inilevtail.EQ.ilevtail.AND. &
             inzzalp.EQ.zzalp.AND.ingrav.EQ.grav.AND.inkappa.EQ.kappa) THEN
          nofile=.false.
        ELSE
          CLOSE(993)
        END IF
      END IF
    END IF
    !
    ustarm = 5.
    alpham = 20.*aalpha
    levtailm = 0.05
    delust  = ustarm/real(iustar)
    delalp  = alpham/real(ialpha)
    deltail = alpham/real(ilevtail)
    const1  = bbeta/kappa**2
    omegac  = sigmax
800 CONTINUE
    IF ( nofile ) THEN
      WRITE(ndse,*) 'Filling 3D look-up table for SIN4. please wait'
      WRITE(ndse,*)  idstr, vergrd, sigmax, aalpha, bbeta, iustar, ialpha,  &
           ilevtail, zzalp, kappa, grav
      !
      tauhft(0:iustar,0:ialpha)=0.  !table initialization
      !
      ALLOCATE(w(jtot))
      w(2:jtot-1)=1.
      w(1)=0.5
      w(jtot)=0.5
      x0 = 0.05
      !
      DO k=0,iustar
        ust0      = max(real(k)*delust,0.000001)
        DO l=0,ialpha
          ust=ust0
          zz0       = ust0**2*(aalpha+float(l)*delalp)/grav
          omegacc  = max(omegac,x0*grav/ust)
          yc       = omegacc*sqrt(zz0/grav)
          dely     = max((1.-yc)/real(jtot),0.)
          ! For a given value of UST and ALPHA,
          ! the wave-supported stress is integrated all the way
          ! to 0.05*g/UST
          DO i=0,ilevtail
            levtail=real(i)*deltail
            tauhft(k,l)=0.
            tauhft2(k,l,i)=0.
            tauw0=ust0**2
            tauw=tauw0
            DO j=1,jtot
              y        = yc+real(j-1)*dely
              omega    = y*sqrt(grav/zz0)
              ! This is the deep water phase speed
              cm       = grav/omega
              !this is the inverse wave age, shifted by ZZALP (tuning)
              zx       = ust0/cm +zzalp
              zarg     = min(kappa/zx,20.)
              zmu      = min(grav*zz0/cm**2*exp(zarg),1.)
              zlog     = min(alog(zmu),0.)
              zbeta        = const1*zmu*zlog**4
              ! Power of Y in denominator should be FACHFE-4
              tauhft(k,l)  = tauhft(k,l)+w(j)*zbeta/y*dely
              zx       = ust/cm +zzalp
              zarg     = min(kappa/zx,20.)
              zmu      = min(grav*zz0/cm**2*exp(zarg),1.)
              zlog     = min(alog(zmu),0.)
              zbeta        = const1*zmu*zlog**4
              ! Power of Y in denominator should be FACHFE-4
              tauhft2(k,l,i)  = tauhft2(k,l,i)+w(j)*zbeta*(ust/ust0)**2/y*dely
              tauw=tauw-w(j)*ust**2*zbeta*levtail/y*dely
              ust=sqrt(max(tauw,0.))
            END DO
#ifdef W3_T
            WRITE (ndst,9000) k,l,i,ust0,aalpha+float(l)*delalp,levtail,tauhft2(k,l,i)
#endif
          END DO
        END DO
      END DO
      DEALLOCATE(w)
      OPEN (993,file=fnametab,form='UNFORMATTED', convert=file_endian,iostat=ierr,status='UNKNOWN')
      WRITE(993) idstr, vergrd, sigmax, aalpha, bbeta, iustar, ialpha, ilevtail, zzalp, kappa, grav
      WRITE(993) tauhft(0:iustar,0:ialpha)
      WRITE(993) tauhft2
      CLOSE(993)
      !DO K=0,IUSTAR
      !  DO L=0,IALPHA
      !    DO I=0,ILEVTAIL
      !      WRITE(995,*) K,L,I,MAX(REAL(K)*DELUST,0.000001),AALPHA+FLOAT(L)*DELALP,REAL(I)*DELTAIL,TAUHFT(K,L),TAUHFT2(K,L,I)
      !      END DO
      !    END DO
      !  END DO
      !
    ELSE
      WRITE(ndse,*) 'Reading 3D look-up table for SIN4 from file.'
      READ(993,err=2000,iostat=ierr ) tauhft(0:iustar,0:ialpha)
      READ(993,err=2000,iostat=ierr ) tauhft2
      CLOSE(993)
    END IF
    !
    GOTO 2001
2000 nofile=.true.
    GOTO 800
2001 CONTINUE
    RETURN
#ifdef W3_T
9000 FORMAT (' TEST TABU_HFT2, K, L, I, UST, ALPHA, LEVTAIL, TAUHFT2(K,L,I) :',(3i4,4f10.5))
#endif
  END SUBROUTINE tabu_tauhf2
 
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE calc_ustar(WINDSPEED,TAUW,USTAR,Z0,CHARN)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |            F. Ardhuin             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update 2006/08/14            |
    !/                  +-----------------------------------+
    !/
    !/    27-Feb-2004 : Origination in WW3                  ( version 2.22-SHOM )
    !/     the resulting table was checked to be identical to the original f77 result
    !/    14-Aug-2006 : Modified following Bidlot           ( version 2.22-SHOM )
    !/    18-Aug-2006 : Ported to version 3.09
    !/    03-Apr-2010 : Adding output of Charnock parameter ( version 3.14-IFREMER )
    !/    03-May-2024 : Optional functional form of         ( version 7.15 )
    !/                  Charnock coefficient and surface drag (UK Met Office).
    !
    !  1. Purpose :
    !
    !     Compute friction velocity based on wind speed U10
    !
    !  2. Method :
    !
    !     Computation of u* based on Quasi-linear theory
    !     uses Charnock relation with modified roughness Z1=Z0/SQRT(1-TAUW/TAU)
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       U10,TAUW,USTAR,Z0
    !     ----------------------------------------------------------------
    !       WINDSPEED Real  I   10-m wind speed ... should be NEUTRAL
    !       TAUW      Real  I   Wave-supported stress
    !       USTAR     Real  O   Friction velocity.
    !       Z0        Real  O   air-side roughness length
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3SIN3   Wind input Source term routine.
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !       !/T      Enable test output.
    !
    ! 10. Source code :
    !-----------------------------------------------------------------------------!
    USE constants, ONLY: grav, kappa, nu_air
    USE w3gdatmd,  ONLY: zzwnd, aalpha, zz0max, sintailpar, capchnk
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
    IMPLICIT NONE
    REAL, intent(in) :: WINDSPEED,TAUW
    REAL, intent(out) :: USTAR, Z0, CHARN
    ! local variables
    REAL             :: SQRTCDM1
    REAL             :: XI,DELI1,DELI2,XJ,delj1,delj2  ! used for table version
    INTEGER          :: IND,J
    REAL             :: TAUW_LOCAL
    REAL             :: TAUOLD,CDRAG,WCD,USTOLD,X,UST,ZZ0,ZNU,ZZ00,F,DELF
    REAL             :: CHATH, XMIN ! used for reduction of high winds
    INTEGER, PARAMETER      :: NITER=10
    REAL   , PARAMETER      :: XM=0.50, eps1=0.00001
    INTEGER                 :: ITER
    !     VARIABLE.   TYPE.     PURPOSE.
    !      *XM*        REAL      POWER OF TAUW/TAU IN ROUGHNESS LENGTH.
    !      *XNU*       REAL      KINEMATIC VISCOSITY OF AIR.
    !      *NITER*     INTEGER   NUMBER OF ITERATIONS TO OBTAIN TOTAL STRESS
    !      *EPS1*      REAL      SMALL NUMBER TO MAKE SURE THAT A SOLUTION
    !                            IS OBTAINED IN ITERATION WITH TAU>TAUW.
 
    chath = aalpha
    !
    IF (sintailpar(1).GT.0.5) THEN
      tauw_local=max(min(tauw,tauwmax),0.)
      xi      = sqrt(tauw_local)/deltauw
      ind     = min( itaumax-1, int(xi)) ! index for stress table
      deli1   = min(1.,xi - real(ind))  !interpolation coefficient for stress table
      deli2   = 1. - deli1
      xj      = windspeed/delu
      j       = min( jumax-1, int(xj) )
      delj1   = min(1.,xj - real(j))
      delj2   = 1. - delj1
      ustar=(taut(ind,j)*deli2+taut(ind+1,j  )*deli1)*delj2 &
           + (taut(ind,j+1)*deli2+taut(ind+1,j+1)*deli1)*delj1
    ELSE
      IF (capchnk(1).EQ.1.) THEN
        ! Computation of sea surface roughness and charnock coefficient based
        ! on Donelan (2018). Determines minimum charnock; reduction for winds
        ! above a particular threshold
        chath  = capchnk(2) + 0.5 * (capchnk(3) - capchnk(2)) * (1 & 
                 - tanh((windspeed-capchnk(4))/capchnk(5)))
        xmin   = 0.15 * (capchnk(3)-chath)
      ELSE
        xmin  = 0.
      END IF
 
      ! This max is for comparison ... to be removed later
      !        TAUW_LOCAL=MAX(MIN(TAUW,TAUWMAX),0.)
      tauw_local=tauw
      cdrag   = 0.0012875
      wcd     = sqrt(cdrag)
      ustold  = windspeed*wcd
      tauold  = max(ustold**2, tauw_local+eps1)
      ! Newton method to solve for ustar in U=ustar*log(Z/Z0)
      DO iter=1,niter
        x   = max(tauw_local/tauold, xmin)
        ust = sqrt(tauold)
        zz00 = chath*tauold/grav
        IF (zz0max.NE.0) zz00=min(zz00,zz0max)
        ! Corrects roughness ZZ00 for quasi-linear effect
        zz0 = zz00/(1.-x)**xm
        znu = 0.11*nu_air/max(ust,1e-6)
        zz0 = sintailpar(5)*znu+zz0
        f   = ust-kappa*windspeed/(alog(zzwnd/zz0))
        delf= 1.-kappa*windspeed/(alog(zzwnd/zz0))**2*2./ust &
             *(1.-(xm+1)*x)/(1.-x)
        ust = ust-f/delf
        tauold= max(ust**2., tauw_local+eps1)
      END DO
      ustar=ust
    END IF
    !
    ! Determines roughness length
    !
    IF (ustar.GT.0.001) THEN
      sqrtcdm1  = min(windspeed/ustar,100.0)
      z0  = zzwnd*exp(-kappa*sqrtcdm1)
      charn = grav*z0/ustar**2
    ELSE
      IF (ustar.GT.0) THEN
        sqrtcdm1  = min(windspeed/ustar,100.0)
        z0  = zzwnd*exp(-kappa*sqrtcdm1)
      ELSE
        z0 = chath*0.001*0.001/grav
      END IF
      charn = chath
    END IF
    IF(capchnk(1) .EQ. 1) THEN
      ! Problem with large values of CHARN for low winds
      charn = min( 0.09 , charn )
      IF(charn.LT.chath) charn = chath
    ENDIF
 
    !  WRITE(6,*) 'CALC_USTAR:',WINDSPEED,TAUW,AALPHA,CHARN,Z0,USTAR
    !
    RETURN
  END SUBROUTINE calc_ustar
  !/ ------------------------------------------------------------------- /
 
  SUBROUTINE w3sds4 (A, K, CG, USTAR, USDIR, DEPTH, DAIR, SRHS,    &
       DDIAG, IX, IY, BRLAMBDA, WHITECAP, DLWMEAN )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  ! F. Ardhuin, F. Leckler, L. Romero !
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Aug-2021 |
    !/                  +-----------------------------------+
    !/
    !/    30-Aug-2010 : Clean up from common ST3-ST4 routine( version 3.14-Ifremer )
    !/    23-Jan-2012 : Add output of lambdas to be used in SIN
    !/    13-Nov-2013 : Reduced frequency range with IG1 switch
    !/    06-Jun-2018 : Add optional DEBUGSRC              ( version 6.04 )
    !/    22-Feb-2020 : Option to use Romero (GRL 2019)    ( version 7.06 )
    !/    13-Aug-2021 : Consider DAIR a variable           ( version 7.14 )
    !/    01-Mar-2023 : Clean up of SDS4                   ( version 7.xx )
    !/
    !  1. Purpose :
    !
    !     Calculate wave dissipation source term and diagonal term of derivative.
    !
    !  2. Method :
    !
    !       This codes does either one or the other of
    !       Ardhuin et al. (JPO 2010)
    !       Filipot & Ardhuin (JGR 2012)
    !       Romero (GRL 2009)
    !       the choice depends on SDSBCHOICE
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       IX, IY    Int   I   Grid Index
    !       A         R.A.  I   Action density spectrum (1-D).
    !       K         R.A.  I   Wavenumber for entire spectrum.          *)
    !       USTAR     Real  I   Friction velocity.
    !       USDIR     Real  I   wind stress direction.
    !       DEPTH     Real  I   Water depth.
    !       DAIR      Real  I   Air density
    !       S         R.A.  O   Source term (1-D version).
    !       D         R.A.  O   Diagonal term of derivative.             *)
    !       BRLAMBDA  R.A.  O   Phillips' Lambdas
    !     ----------------------------------------------------------------
    !                         *) Stored in 1-D array with dimension NTH*NK
    !
    !  4. Subroutines used :
    !
    !       STRACE    Subroutine tracing.                 ( !/S switch )
    !       PRT2DS    Print plot of spectrum.             ( !/T0 switch )
    !       OUTMAT    Print out matrix.                   ( !/T1 switch )
    !
    !  5. Called by :
    !
    !       W3SRCE   Source term integration.
    !       W3EXPO   Point output program.
    !       GXEXPO   GrADS point output program.
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable general test output.
    !     !/T0  2-D print plot of source term.
    !     !/T1  Print arrays.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants,ONLY: grav, dwat, pi, tpi, rade, debug_node
    USE w3gdatmd, ONLY: nspec, nth, nk, ssdsbr, ssdsbt, dden,      &
         ssdsc, ec2, es2, esc,                      &
         sig, ssdsp, ecos, esin, dth, aairgb,       &
         ssdsiso, ssdsdth, ssdsbm, aaircmin,        &
         ssdsbrfdf, ssdsbck, iktab, dcki,           &
         satindices, satweights, cumulw, nkhs, nkd, &
         ndtab, qbi, dsip, ssdsbrf1,xfr
#ifdef W3_IG1
    USE w3gdatmd, ONLY: igpars
#endif
    USE w3odatmd, ONLY: flogrd
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
#ifdef W3_T0
    USE w3odatmd, ONLY: ndst
    USE w3arrymd, ONLY: prt2ds
#endif
#ifdef W3_T1
    USE w3arrymd, ONLY: outmat
#endif
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, OPTIONAL, INTENT(IN) :: IX, IY
    REAL, INTENT(IN)        :: A(NSPEC), K(NK), CG(NK),            &
         DEPTH, DAIR, USTAR, USDIR, DLWMEAN
    REAL, INTENT(OUT)       :: SRHS(NSPEC), DDIAG(NSPEC), BRLAMBDA(NSPEC)
    REAL, INTENT(OUT)       :: WHITECAP(1:4)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IS, IS2, IS0, IKL, IKC, ID, NKL
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER                 :: IK, IK1, ITH, IK2, JTH, ITH2,             &
         IKHS, IKD, SDSNTH, IT, IKM, NKM
    INTEGER                 :: NSMOOTH(NK)
    REAL                    :: C, C2, CUMULWISO, COSWIND, ASUM, SDIAGISO
    REAL                    :: COEF1, COEF2, COEF4(NK),      &
         COEF5(NK)
 
    REAL                    :: FACTURB, FACTURB2, DTURB, DVISC, DIAG2, BREAKFRACTION
    REAL                    :: RENEWALFREQ, EPSR
    REAL                    :: S1(NK), E1(NK)
    INTEGER                 :: NTIMES(NK)
    REAL                    :: GAM, XT
    REAL                    :: DK(NK), HS(NK), KBAR(NK), DCK(NK)
    REAL                    :: EFDF(NK)     ! Energy integrated over a spectral band
    INTEGER                 :: IKSUP(NK)
    REAL                    :: FACSAT, DKHS, FACSTRAINB, FACSTRAINL
    REAL                    :: BTH0(NK)     !saturation spectrum
    REAL                    :: BTH(NSPEC)   !saturation spectrum
    REAL                    :: MSSSUM(NK,5),  FACHF
    REAL                    :: MSSLONG
    REAL                    :: MSSPCS, MSSPC2, MSSPS2, MSSP, MSSD, MSSTH
    REAL                    :: MICHE, X, KLOC
#ifdef W3_T0
    REAL                    :: DOUT(NK,NTH)
#endif
    REAL                    :: QB(NK), S2(NK)
    REAL                    :: TSTR, TMAX, DT, T, MFT, DIRFORCUM
    REAL                    :: PB(NSPEC), PB2(NSPEC), BRM12(NK), BTOVER
    REAL                    :: KO, LMODULATION(NTH)
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SDS4')
#endif
    !
    !
    !----------------------------------------------------------------------
    !
    ! 0.  Pre-Initialization to zero out arrays. All arrays should be reset
    !     within the computation, but these are helping with some bugs
    !     found in certain compilers
    nsmooth=0
    s1=0.; e1=0.
    ntimes=0;iksup=0
    dk=0.; hs=0.; kbar=0.; dck=0.; efdf=0.
    bth0=0.; bth=0.; ddiag=0.; srhs=0.; pb=0.
    msssum(:,:)=0.
#ifdef W3_T0
    dout=0.
#endif
    qb=0.; s2=0.;pb=0.; pb2=0.
    brm12(:)=0.
    !
    ! 1.  Initialization and numerical factors
    !
    facturb=ssdsc(5)*ustar**2/grav*dair/dwat
    dikcumul = nint(ssdsbrf1/(xfr-1.))
    breakfraction=0.
    renewalfreq=0.
    ik1=1
#ifdef W3_IG1
    ik1=nint(igpars(5))+1
#endif
    !
    ! 1.b MSS parameters used for Modulation factors for lambda (Romero )
    !
    IF (ssdsc(8).GT.0.OR.ssdsc(11).GT.0.OR.ssdsc(18).GT.0) THEN
      DO ik=1,nk
        mssp   = 0.
        msspc2 = 0.
        mssps2 = 0.
        msspcs = 0.
        !
        ! Sums the contributions to the directional MSS for all angles
        !
        DO ith=1,nth
          is=ith+(ik-1)*nth
          msslong  = k(ik)**ssdsc(20) * a(is) * dden(ik) / cg(ik) ! contribution to MSS
          msspc2 = msspc2 +msslong*ec2(ith)
          mssps2 = mssps2 +msslong*es2(ith)
          msspcs = msspcs +msslong*esc(ith)
          mssp   = mssp   +msslong
        END DO
        msssum(ik:nk,1) = msssum(ik:nk,1) +mssp
        msssum(ik:nk,3) = msssum(ik:nk,3) +msspc2
        msssum(ik:nk,4) = msssum(ik:nk,4) +mssps2
        msssum(ik:nk,5) = msssum(ik:nk,5) +msspcs
        !
        ! Direction of long wave mss summed up to IK
        !
        mssd=0.5*(atan2(2*msssum(ik,5),msssum(ik,3)-msssum(ik,4)))
        IF (mssd.LT.0) mssd = mssd + pi
        msssum(ik,2)  =  mssd
      END DO
    END IF ! SSDSC(8).GT.0) THEN
    !
    ! 2.   Estimation of spontaneous breaking from local saturation
    !
    !############################################################################################"
    SELECT CASE (nint(ssdsc(1)))
    CASE (1)
      !
      ! 2.a  Case of a direction-dependent breaking term following Ardhuin et al. 2010
      !
      epsr = sqrt(ssdsbr)
      !
      ! 2.a.1 Computes saturation
      !
      bth(:) = 0.
 
      DO  ik=ik1, nk
 
        facsat=sig(ik)*k(ik)**3*dth
        is0=(ik-1)*nth
        bth(is0+1)=0.
        asum = sum(a(is0+1:is0+nth))
        bth0(ik)=asum*facsat
        !
        IF (ssdsdth.GE.180) THEN  ! integrates around full circle
          bth(is0+1:is0+nth)=bth0(ik)
        ELSE
          DO ith=1,nth            ! partial integration
            is=ith+(ik-1)*nth
            bth(is)=dot_product(satweights(:,ith),  a(is0+satindices(:,ith)) ) &
                 *facsat
          END DO
 
          bth0(ik)=maxval(bth(is0+1:is0+nth))
        END IF
        !
      END DO !IK=NK
      !
      !  2.a.2  Computes spontaneous breaking dissipation rate
      !
      DO  ik=ik1, nk
        !
        !  Correction of saturation level for shallow-water kinematics
        !
        IF (ssdsbm(0).EQ.1) THEN
          miche=1.
        ELSE
          x=tanh(min(k(ik)*depth,10.))
          ! Correction of saturation threshold for shallow-water kinematics
          miche=(x*(ssdsbm(1)+x*(ssdsbm(2)+x*(ssdsbm(3)+x*ssdsbm(4)))))**2
        END IF
        coef1=(ssdsbr*miche)
        !
        !  Computes isotropic part
        !
        sdiagiso = ssdsc(2) * sig(ik)*ssdsc(6)*(max(0.,bth0(ik)/coef1-1.))**2
        !
        !  Computes anisotropic part and sums isotropic part
        !
        coef2=ssdsc(2) * sig(ik)*(1-ssdsc(6))/(coef1*coef1)
        ddiag((ik-1)*nth+1:ik*nth) = sdiagiso + &
             coef2*((max(0.,bth((ik-1)*nth+1:ik*nth)-coef1))**ssdsp)
      END DO
      !
      ! Computes Breaking probability
      !
      pb = (max(sqrt(bth)-epsr,0.))**2
      !
      ! Multiplies by 28.16 = 22.0 * 1.6² * 1/2 with
      !  22.0 (Banner & al. 2000, figure 6)
      !  1.6  the coefficient that transforms  SQRT(B) to Banner et al. (2000)'s epsilon
      !  1/2  factor to correct overestimation of Banner et al. (2000)'s breaking probability due to zero-crossing analysis
      !
      pb = pb * 28.16
      ! Compute Lambda = PB* l(k,th)
      ! with l(k,th)=1/(2*pi²)= the breaking crest density
      brlambda = pb / (2.*pi**2.)
      srhs = ddiag * a
 
      !############################################################################################"
    CASE(2)
      !
      ! 2.b             Computes spontaneous breaking for T500 (Filipot et al. JGR 2010)
      !
      e1 = 0.
      hs = 0.
      srhs  = 0.
      ddiag = 0.
      pb2  = 0.
      !
      ! Computes Wavenumber spectrum E1 integrated over direction and computes dk
      !
      DO ik=ik1, nk
        e1(ik)=0.
        DO ith=1,nth
          is=ith+(ik-1)*nth
          e1(ik)=e1(ik)+(a(is)*sig(ik))*dth
        END DO
        dk(ik)=dden(ik)/(dth*sig(ik)*cg(ik))
      END DO
      !
      ! Gets windows indices of IKTAB
      !
      id=min(nint(depth),ndtab)
      IF (id < 1) THEN
        id = 1
      ELSE IF(id > ndtab) THEN
        id = ndtab
      END IF
      !
      ! loop over wave scales
      !
      hs=0.
      efdf=0.
      kbar=0.
      nkl=0. !number of windows
      DO ikl=1,nk
        iksup(ikl)=iktab(ikl,id)
        IF (iksup(ikl) .LE. nk) THEN
          efdf(ikl) = dot_product(e1(ikl:iksup(ikl)-1),dk(ikl:iksup(ikl)-1))
          IF (efdf(ikl) .NE. 0) THEN
            kbar(ikl) = dot_product(k(ikl:iksup(ikl)-1)*e1(ikl:iksup(ikl)-1), &
                 dk(ikl:iksup(ikl)-1)) / efdf(ikl)
          ELSE
            kbar(ikl)=0.
          END IF
          ! estimation of Significant wave height of a given scale
          hs(ikl) = 4*sqrt(efdf(ikl))
          nkl = nkl+1
        END IF
      END DO
      !
      ! Computes Dissipation and breaking probability in each scale
      !
      dck=0.
      qb =0.
      dkhs = khsmax/nkhs
      DO ikl=1, nkl
        IF (hs(ikl) .NE. 0. .AND. kbar(ikl) .NE. 0.)  THEN
          ! gets indices for tabulated dissipation DCKI and breaking probability QBI
          !
          ikd = fac_kd2+anint(log(kbar(ikl)*depth)/log(fac_kd1))
          ikhs= 1+anint(kbar(ikl)*hs(ikl)/dkhs)
          IF (ikd > nkd) THEN    ! Deep water
            ikd = nkd
          ELSE IF (ikd < 1) THEN ! Shallow water
            ikd = 1
          END IF
          IF (ikhs > nkhs) THEN
            ikhs = nkhs
          ELSE IF (ikhs < 1) THEN
            ikhs = 1
          END IF
          xt = tanh(kbar(ikl)*depth)
          !
          ! Gamma corrected for water depth
          !
          gam=1.0314*(xt**3)-1.9958*(xt**2)+1.5522*xt+0.1885
          !
          ! Computes the energy dissipated for the scale IKL
          ! using DCKI which is tabulated in INSIN4
          !
          dck(ikl)=((kbar(ikl)**(-2.5))*(kbar(ikl)/(2*pi)))*dcki(ikhs,ikd)
          !
          ! Get the breaking probability for the scale IKL
          !
          qb(ikl) = qbi(ikhs,ikd) ! QBI is tabulated in INSIN4
        ELSE
          dck(ikl)=0.
          qb(ikl) =0.
        END IF
      END DO
      !
      ! Distributes scale dissipation over the frequency spectrum
      !
      s1 = 0.
      s2 = 0.
      ntimes = 0
      DO ikl=1, nkl
        IF (efdf(ikl) .GT. 0.) THEN
          s1(ikl:iksup(ikl))    = s1(ikl:iksup(ikl)) + &
               dck(ikl)*e1(ikl:iksup(ikl)) / efdf(ikl)
          s2(ikl:iksup(ikl))    = s2(ikl:iksup(ikl)) + &
               qb(ikl) *e1(ikl:iksup(ikl)) / efdf(ikl)
          ntimes(ikl:iksup(ikl)) = ntimes(ikl:iksup(ikl)) + 1
        END IF
      END DO
      !
      ! Finish the average
      !
      WHERE (ntimes .GT. 0)
        s1 = s1 / ntimes
        s2 = s2 / ntimes
      ELSEWHERE
        s1 = 0.
        s2 = 0.
      END WHERE
      ! goes back to action for dissipation source term
      s1(1:nk) = s1(1:nk) / sig(1:nk)
      !
      ! Makes Isotropic distribution
      !
      asum = 0.
      DO ik = 1, nk
        asum = (sum(a(((ik-1)*nth+1):(ik*nth)))*dth)
        IF (asum.GT.1.e-8) THEN
          FORALL (is=1+(ik-1)*nth:ik*nth) ddiag(is)  = s1(ik)/asum
          FORALL (is=1+(ik-1)*nth:ik*nth) pb2(is) = s2(ik)/asum
        ELSE
          FORALL (is=1+(ik-1)*nth:ik*nth) ddiag(is)  = 0.
          FORALL (is=1+(ik-1)*nth:ik*nth) pb2(is) = 0.
        END IF
        IF (pb2(1+(ik-1)*nth).GT.0.001) THEN
          bth0(ik) = 2.*ssdsbr
        ELSE
          bth0(ik) = 0.
        END IF
      END DO
      !
      pb = (1-ssdsc(1))*pb2*a + ssdsc(1)*pb
      ! Compute Lambda = PB* l(k,th)
      ! with l(k,th)=1/(2*pi²)= the breaking crest density
      brlambda = pb / (2.*pi**2.)
      srhs = ddiag * a
      !############################################################################################"
    CASE(3)
      !
      ! 2c Romero (GRL 2019)
      !
      ! directional saturation I
      ! integrate in azimuth
      ko=(grav/(1e-6+ustar**2))/(28./ssdsc(16))**2
      DO ik=1,nk
        is0=(ik-1)*nth
        kloc=k(ik)**(2-ssdsc(20)) ! local wavenumber factor, if mss not used.
        bth(1:nth)=max(a(is0+1:is0+nth)*sig(ik)*k(ik)**3,.00000000000001)
        !
        dirforcum=dlwmean
        IF (ssdsc(11).GT.0) dirforcum=msssum(ik,2)
 
        c=sig(ik)/k(ik)
        bth0(ik)=sum(bth(1:nth)*dth)
        IF (ssdsc(18).GT.0) THEN ! Applies modulation factor on Lambda
          DO ith=1,nth
            facstrainl=1.+ssdsc(18)*((msssum(ik,1)*kloc)**ssdsc(14) *      &   ! Romero
                 (ecos(ith)*cos(dirforcum)+esin(ith)*sin(dirforcum))**2)
            lmodulation(ith)= facstrainl**ssdsc(19)
          END DO
        ELSE
          lmodulation(:)= 1.
        END IF
 
        ! Lambda
        brlambda(is0+1:is0+nth)=ssdsc(9)*exp(-ssdsbr/bth(1:nth))              &
             *( 1.0+ssdsc(13)*max(1.,(k(ik)/ko))**ssdsc(15) ) &
             /(ssdsc(13)+1)*lmodulation(1:nth)
        ! Breaking strength : generalisation of Duncan's b parameter
        btover = sqrt(bth0(ik))-sqrt(ssdsbt)
        brm12(ik)=ssdsc(2)*(max(0.,btover))**(2.5)/sig(ik)  ! not function of direction
        !  For consistency set BRLAMBDA set to zero if b is zero
        brlambda(is0+1:is0+nth)= max(0.,sign(brlambda(is0+1:is0+nth),btover))
        !  Source term / sig2  (action dissipation)
        srhs(is0+1:is0+nth)= brm12(ik)/grav**2*brlambda(is0+1:is0+nth)*c**5
        ! diagonal
        ddiag(is0+1:is0+nth) = srhs(is0+1:is0+nth)*ssdsbr/max(1.e-20,bth(1:nth))/max(1e-20,a(is0+1:is0+nth))  !
      END DO
      !   Breaking probability (Is actually the breaking rate)
      pb = brlambda *c
      !
    END SELECT
    !############################################################################################"
    !
    !
    !/ ------------------------------------------------------------------- /
    !             WAVE-TURBULENCE INTERACTION AND CUMULATIVE EFFECT
    !/ ------------------------------------------------------------------- /
    !
    !
    ! loop over spectrum
    !
    IF ( (ssdsc(3).NE.0.) .OR. (ssdsc(5).NE.0.) .OR. (ssdsc(21).NE.0.) ) THEN
      DO  ik=ik1, nk
        renewalfreq = 0.
        facturb2=-2.*sig(ik)*k(ik)*facturb
        dvisc=-4.*ssdsc(21)*k(ik)*k(ik)
        c = sig(ik)/k(ik) ! phase speed
        !
        IF (ssdsc(3).GT.0 .AND. ik.GT.dikcumul) THEN
          ! this is the cheap isotropic version
          DO ik2=ik1,ik-dikcumul
            c2 = sig(ik2)/k(ik2)
            is2=(ik2-1)*nth
            cumulwiso=abs(c2-c)*dsip(ik2)/(0.5*c2) * dth
            renewalfreq=renewalfreq-cumulwiso*sum(brlambda(is2+1:is2+nth))
          END DO
        END IF
 
        DO ith=1,nth
          is=ith+(ik-1)*nth
          !
          ! Computes cumulative effect from Breaking probability
          !
          IF (ssdsc(3).LT.0 .AND. ik.GT.dikcumul) THEN
            renewalfreq = 0.
            ! this is the expensive and largely useless version
            DO ik2=ik1,ik-dikcumul
              IF (bth0(ik2).GT.ssdsbr) THEN
                is2=(ik2-1)*nth
                renewalfreq=renewalfreq+dot_product(cumulw(is2+1:is2+nth,is),brlambda(is2+1:is2+nth))
              END IF
            END DO
          END IF
          !
          ! Computes wave turbulence interaction
          !
          coswind=(ecos(ith)*cos(usdir)+esin(ith)*sin(usdir))
          dturb=facturb2*max(0.,coswind)  ! Theory -> stress direction
          !
          ! Add effects
          !
          diag2 = (ssdsc(3)*renewalfreq+dturb+dvisc)
          ddiag(is) = ddiag(is) + diag2
          srhs(is)  = srhs(is)  + a(is)* diag2
        END DO
      END DO
    END IF
    !
    !  COMPUTES WHITECAP PARAMETERS
    !
    IF ( .NOT. (flogrd(5,7).OR.flogrd(5,8) ) ) THEN
      RETURN
    END IF
    !
    whitecap(1:4) = 0.
    !
    ! precomputes integration of Lambda over direction
    ! times wavelength times a (a=5 in Reul&Chapron JGR 2003) times dk
    !
    DO ik=1,min(floor(aaircmin),nk)
      c=sig(ik)/k(ik)
      is0=(ik-1)*nth
      coef4(ik) = c*c*sum(brlambda(is0+1:is0+nth))                          &
           *2.*pi/grav*ssdsc(7) * dden(ik)/(sig(ik)*cg(ik))
      coef5(ik) = c**3*sum(brlambda(is0+1:is0+nth)                           &
           *brm12(ik))                                 &
           *aairgb/grav * dden(ik)/(sig(ik)*cg(ik))
      !        COEF4(IK) = SUM(BRLAMBDA((IK-1)*NTH+1:IK*NTH) * DTH) *(2*PI/K(IK)) *  &
      !                    SSDSC(7) * DDEN(IK)/(DTH*SIG(IK)*CG(IK))
      !                   NB: SSDSC(7) is WHITECAPWIDTH
    END DO
    ! Need to extrapolate above NK if necessary ... to be added later.
    DO ik=min(floor(aaircmin),nk),nk
      coef4(ik)=0.
      coef5(ik)=0.
    END DO
 
    !/
    IF ( flogrd(5,7) ) THEN
      !
      ! Computes the Total WhiteCap Coverage (a=5. ; Reul and Chapron, 2003)
      !
      DO ik=ik1,min(floor(aaircmin),nk)
        whitecap(1) = whitecap(1) + coef4(ik) * (1-whitecap(1))
        whitecap(4) = whitecap(4) + coef5(ik)
      END DO
    END IF
    !/
    IF ( flogrd(5,8) ) THEN
      !
      ! Calculates the Mean Foam Thickness for component K(IK) => Fig.3, Reul and Chapron, 2003
      ! ( Copied from ST4 - not yet tested/validated with Romero 2019 (Lambda model)
      !
      DO ik=ik1,nk
        !    Duration of active breaking (TAU*)
        tstr = 0.8 * 2*pi/sig(ik)
        !    Time persistence of foam (a=5.)
        tmax = 5.  * 2*pi/sig(ik)
        dt   = tmax / 50
        mft  = 0.
        DO it = 1, 50
          ! integration over time of foam persistance
          t = float(it) * dt
          ! Eq. 5 and 6 of Reul and Chapron, 2003
          IF ( t .LT. tstr ) THEN
            mft = mft + 0.4 / (k(ik)*tstr) * t * dt
          ELSE
            mft = mft + 0.4 / k(ik) * exp(-1*(t-tstr)/3.8) * dt
          END IF
        END DO
        mft = mft / tmax
        !
        ! Computes foam-layer thickness (Reul and Chapron, 2003)
        !
        whitecap(2) = whitecap(2) + coef4(ik) * mft
      END DO
    END IF
    !
    ! End of output computing
    !
    RETURN
    !
    ! Formats
    !
    !/
    !/ End of W3SDS4 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3sds4
 
 
END MODULE w3src4md