WW3/w3src3md_8F90_source.html

#include "w3macros.h"

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


MODULE w3src3md

  !/

  !/                  +-----------------------------------+

  !/                  | WAVEWATCH III                SHOM |

  !/                  !            F. Ardhuin             !

  !/                  |           H. L. Tolman            |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         02-Sep-2012 |

  !/                  +-----------------------------------+

  !/

  !/    09-Oct-2007 : Origination.                        ( version 3.13 )

  !/    10-Oct-2010 : Adding Janssen-style swell damping  ( version 3.14)

  !/    02-Sep-2012 : Clean up test output T, T0, T1      ( version 4.07 )

  !/

  !  1. Purpose :

  !

  !     The 'WAM4+' source terms based on P.A.E.M. Janssen's work, with

  !     extensions by him and by J.-R. Bidlot. Converted from the original

  !     WAM codes by F. Ardhuin, with further extensions to adapt to a

  !     saturation-based breaking and observation-based swell dissipation

  !

  !

  !  2. Variables and types :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !     ----------------------------------------------------------------

  !

  !  3. Subroutines and functions :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !      W3SPR3    Subr. Public   Mean parameters from spectrum.

  !      W3SIN3    Subr. Public   WAM4+ input source term.

  !      INSIN3    Subr. Public   Corresponding initialization routine.

  !      TABU_STRESS, TABU_TAUHF, TABU_TAUHF2

  !                Subr. Public   Populate various tables.

  !      CALC_USTAR

  !                Subr. Public   Compute stresses.

  !      W3SDS3    Subr. Public   User supplied dissipation.

  !     ----------------------------------------------------------------

  !

  !  4. Subroutines and functions used :

  !

  !      Name      Type  Module   Description

  !     ----------------------------------------------------------------

  !     ----------------------------------------------------------------

  !

  !  5. Remarks :

  !

  !  6. Switches :

  !

  !  7. Source code :

  !/

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

  !/

  USE constants, ONLY: kappa, nu_air


  IMPLICIT NONE


  PUBLIC

  !/

  !/ Public variables

  !/

  !air kinematic viscosity (used in WAM)

  INTEGER, PARAMETER      :: itaumax=200,jumax=200

  INTEGER, PARAMETER      :: iustar=100,ialpha=200, ilevtail=50

  REAL                    :: taut(0:itaumax,0:jumax), deltauw, delu

  ! Table for H.F. stress as a function of 2 variables

  REAL                    :: tauhft(0:iustar,0:ialpha), delust, delalp

  REAL,    PARAMETER      :: umax    = 50.

  REAL,    PARAMETER      :: tauwmax = 2.2361 !SQRT(5.)

  !/

CONTAINS

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


  SUBROUTINE w3spr3 (A, CG, WN, EMEAN, FMEAN, FMEANS, WNMEAN,     &

       AMAX, U, UDIR, USTAR, USDIR, TAUWX, TAUWY, CD, Z0,&

       CHARN, LLWS, FMEANWS)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III                SHOM |

    !/                  !            F. Ardhuin             !

    !/                  |           H. L. Tolman            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         17-Oct-2007 |

    !/                  +-----------------------------------+

    !/

    !/    03-Oct-2007 : Origination.                        ( version 3.13 )

    !/

    !  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

    !       FMEAN   Real  O   Mean  frequency for determination of tail

    !       FMEANS  Real  O   Mean  frequency for dissipation source term

    !       WNMEAN  Real  O   Mean wavenumber.

    !       AMAX    Real  O   Maximum of action spectrum.

    !       U       Real  I   Wind speed.

    !       UDIR    Real  I   Wind direction.

    !       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.

    !       LLWS    L.A.  I   Wind sea true/false array for each component

    !       FMEANWS Real  O   Mean frequency of wind sea, used for tail

    !     ----------------------------------------------------------------

    !

    !  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.

    !

    ! 10. Source code :

    !

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

    USE constants, ONLY: tpiinv

    USE w3gdatmd, ONLY: nk, nth, nspec, sig, dth, dden, wwnmeanp, &

         wwnmeanptail, fte, ftf, sstxftf, sstxftwn,&

         sstxftftail, sswellf

#ifdef W3_T

    USE w3odatmd, ONLY: ndst

#endif

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    REAL, INTENT(IN)        :: A(NTH,NK), CG(NK), WN(NK), U, UDIR

    REAL, INTENT(IN)        :: TAUWX, TAUWY

    LOGICAL, INTENT(IN)     :: LLWS(NSPEC)

    REAL, INTENT(INOUT)     :: USTAR ,USDIR

    REAL, INTENT(OUT)       :: EMEAN, FMEAN, FMEANS, WNMEAN, AMAX,  &

         CD, Z0, CHARN, FMEANWS

    !/

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

    !/ Local parameters

    !/

    INTEGER                 :: IS, IK, ITH

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif


    REAL                    :: TAUW, EBAND, EMEANWS, UNZ,       &

         EB(NK),EB2(NK),ALFA(NK)

    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SPR3')

#endif

    !

    unz    = max( 0.01 , u )

    ustar  = max( 0.0001 , ustar )

    !

    emean  = 0.

    emeanws= 0.

    fmeanws= 0.

    fmean  = 0.

    fmeans = 0.

    wnmean = 0.

    amax   = 0.

    !

    ! 1.  Integral over directions and maximum --------------------------- *

    !

    DO ik=1, nk

      eb(ik)  = 0.

      eb2(ik) = 0.

      DO ith=1, nth

        is=ith+(ik-1)*nth

        eb(ik) = eb(ik) + a(ith,ik)

        IF (llws(is)) eb2(ik) = eb2(ik) + a(ith,ik)

        amax   = max( amax , a(ith,ik) )

      END DO

    END DO

    !

    ! 2.  Integrate over directions -------------------------------------- *

    !

    DO ik=1, nk

      alfa(ik) = 2. * dth * sig(ik) * eb(ik) * wn(ik)**3

      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)**(2.*wwnmeanptail))

      fmeans   = fmeans + eb(ik) *(sig(ik)**(2.*wwnmeanp))

      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 * sstxftftail

    fmeans = fmeans + eband * sstxftf

    wnmean = wnmean + eband * sstxftwn

    eband  = eb2(nk) / dden(nk)

    emeanws = emeanws + eband * fte

    fmeanws = fmeanws + eband * sstxftftail

    !

    ! 4.  Final processing

    !


    IF (fmean.LT.1.e-7) THEN

      fmean=tpiinv * sig(nk)

    ELSE

      fmean  = tpiinv *( max( 1.e-7 , fmean )                       &

           / max( 1.e-7 , emean ))**(1/(2.*wwnmeanptail))

    END IF

    IF (fmeans.LT.1.e-7) THEN

      fmeans=tpiinv * sig(nk)

    ELSE

      fmeans = tpiinv *( max( 1.e-7 , fmeans )                      &

           / max( 1.e-7 , emean ))**(1/(2.*wwnmeanp))

    END IF

    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)


    z0=0.

    CALL calc_ustar(u,tauw,ustar,z0,charn)

    unz    = max( 0.01 , u )

    cd     = (ustar/unz)**2

    usdir = udir

    !

    ! 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 W3SPR3 : E,WN MN :',f8.3,f8.4/                   &

         '         TPIINV, USTAR, CD, Z0 :',f8.3,f7.2,1x,2f9.5)

#endif

    !/

    !/ End of W3SPR3 ----------------------------------------------------- /

    !/

  END SUBROUTINE w3spr3

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


  SUBROUTINE w3sin3 (A, CG, K, U, USTAR, DRAT, AS, USDIR, Z0, CD,    &

       TAUWX, TAUWY, TAUWNX, TAUWNY, ICE, S, D, LLWS, IX, IY)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III                SHOM |

    !/                  !            F. Ardhuin             !

    !/                  |           H. L. Tolman            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         02-Sep-2012 |

    !/                  +-----------------------------------+

    !/

    !/    09-Oct-2007 : Origination.                        ( version 3.13 )

    !/    16-May-2010 : Adding sea ice                      ( version 3.14_Ifremer )

    !/    02-Sep-2012 : Clean up test output T, T0, T1      ( version 4.07 )

    !/

    !  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. )

    !       SWELLF: damping coefficient (=1) for Janssen (2004) theory

    !

    !  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.

    !       ICE     Real  I   Sea ice fraction.                        !/Stefan: ICE is DUMMY argument; remove later.

    !       S       R.A.  O   Source term (1-D version).

    !       D       R.A.  O   Diagonal term of derivative.             *)

    !       LLWS    L.A.  O   Wind sea true/false array for each component

    !     ----------------------------------------------------------------

    !                         *) 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, tpi

#ifdef W3_T

    USE constants, ONLY: rade

#endif

    USE w3gdatmd, ONLY: nk, nth, nspec, xfr, dden, sig, sig2, th,   &

         esin, ecos, ec2, zzwnd, aalpha, bbeta, zzalp,&

         sswellf,                                     &

         dden2, dth, ssinthp,zz0rat

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

#ifdef W3_T

    USE w3odatmd, ONLY: ndst

#endif

#ifdef W3_T0

    USE w3odatmd, ONLY: ndst

#endif

#ifdef W3_T1

    USE w3odatmd, ONLY: ndst

#endif

#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)

    REAL, INTENT(IN)        :: CG(NK), K(NSPEC),Z0,U, CD

    REAL, INTENT(IN)        :: USTAR, USDIR, AS, DRAT, ICE       !/Stefan: ICE is DUMMY argument; remove later.

    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                    :: COSU, SINU, TAUX, TAUY

    REAL                    :: UST2, TAUW, TAUWB

    REAL   , PARAMETER      :: EPS1 = 0.00001, eps2 = 0.000001

#ifdef W3_STAB3

    REAL                    :: Usigma           !standard deviation of U due to gustiness

    REAL                    :: USTARsigma       !standard deviation of USTAR due to gustiness

#endif

    REAL                    :: CM,UCO,UCN,ZCN, &

         Z0VISC

    REAL XI,DELI1,DELI2

    REAL XJ,DELJ1,DELJ2

    REAL                    :: CONST, CONST0, CONST2, CONST3,  TAU1

    REAL                    :: X,ZARG,ZLOG,ZBETA,UST

    REAL COSWIND,XSTRESS,YSTRESS,TAUHF

    REAL TEMP, TEMP2

    INTEGER IND,J,ISTAB

    REAL DSTAB(3,NSPEC)

    REAL STRESSSTAB(3,2),STRESSSTABN(3,2)

#ifdef W3_T0

    REAL                    :: DOUT(NK,NTH)

#endif

    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SIN3')

#endif

    !

#ifdef W3_T

    WRITE (ndst,9000) bbeta, ustar, usdir*rade

#endif

    !

    ! 1.  Preparations

    !

    !      JDM: initializing arrays (shouldn't change answers)

    dstab=0.; stressstab=0. ;stressstabn=0.

    !

    ! 1.a  estimation of surface roughness parameters

    !

    z0visc = 0.1*nu_air/max(ustar,0.0001)

    !

    ! 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

    usigma=max(0.,-0.025*as)

    ustarsigma=(1.0+u/(10.+u))*usigma

#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)

      !

      ! Loop over the resolved part of the spectrum

      !

      stressstab(istab,:)=0.

      stressstabn(istab,:)=0.

      const0=bbeta*drat/(kappa**2)

      ! SSWELLF(1) is IDAMPING in ECMWAM

      const3 = sswellf(1)*2.*kappa*drat

      !

      cosu   = cos(usdir)

      sinu   = sin(usdir)

      DO ik=1, nk

        is=1+(ik-1)*nth

        cm=k(is)/sig2(is) !inverse of phase speed

        ucn=ust*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

        ! this CM parameter is 1 / C_phi

        ! this is the "correct" shallow-water expression

        ! here Z0 corresponds to Z0+Z1 of the Janssen eq. 14

        zcn=alog(k(is)*z0)

        ! commented below is the original WAM version (OK for deep water)

        !       ZCN=ALOG(G*Z0b(I)*CM(I)**2)

        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

            zlog=zcn+zarg

            zbeta = const3*(coswind+kappa/(zcn*ust*cm))*ucn**2

            IF (zlog.LT.0.) THEN

              ! The source term Sp is beta * omega * X**2

              ! as given by Janssen 1991 eq. 19

              dstab(istab,is) = const*exp(zlog)*zlog**4*ucn**2*coswind**ssinthp

              llws(is)=.true.

            ELSE

              dstab(istab,is) = zbeta

              llws(is)=.true.

            END IF

          ELSE

            zbeta = const3*(coswind+kappa/(zcn*ust*cm))*ucn**2

            dstab(istab,is) = zbeta

            llws(is)=.false.

          END IF

          !

          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)

          END IF

          stressstab(istab,1)=stressstab(istab,1)+temp2*ecos(is)

          stressstab(istab,2)=stressstab(istab,2)+temp2*esin(is)

        END DO

      END DO

      !

      d(:)=dstab(3,:)

      xstress=stressstab(3,1)

      ystress=stressstab(3,2)

      tauwnx =stressstabn(3,1)

      tauwny =stressstabn(3,2)

      !             WRITE(995,'(A,11G14.5)') 'NEGSTRESS:    ',TAUWNX,TAUWNY,FW*UORB**3

#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

    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:nk), '  ', 1.,      &

         0.0, 0.001, 'Diag Sin', ' ', 'NONAME')

#endif

    !

#ifdef W3_T1

    CALL outmat (ndst, d, nth, nth, nk, 'diag Sin')

#endif

    !

    ! 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

    !

    ! 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

    tau1 =(tauhft(ind,j)*deli2+tauhft(ind+1,j)*deli1 )*delj2 &

         +(tauhft(ind,j+1)*deli2+tauhft(ind+1,j+1)*deli1)*delj1

    tauhf = const0*temp*ust**2*tau1

    tauwx = xstress+tauhf*cos(usdir)

    tauwy = ystress+tauhf*sin(usdir)

    !

    ! 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 W3SIN3 : COMMON FACT.: ',3e10.3)

#endif

    !/

    !/ End of W3SIN3 ----------------------------------------------------- /

    !/

  END SUBROUTINE w3sin3

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


  SUBROUTINE insin3

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |                         SHOM      |

    !/                  |            F. Ardhuin             |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         23-Jul-2009 |

    !/                  +-----------------------------------+

    !/

    !/    23-Jun-2006 : Origination.                        ( version 3.09 )

    !/    23-Jul-2007 : Cleaning up convolutions            ( version 3.14-SHOM)

    !

    !  1. Purpose :

    !

    !     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: tpiinv

    USE w3gdatmd,  ONLY: sig, nk

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !/

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

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

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'INSIN3')

#endif

    !

    ! 1.  .... ----------------------------------------------------------- *

    !

    CALL tabu_stress

    CALL tabu_tauhf(sig(nk) * tpiinv)   !tabulate high-frequency stress

    !/

    !/ End of INSIN3 ----------------------------------------------------- /

    !/

  END SUBROUTINE insin3

  ! ----------------------------------------------------------------------


  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: 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,ZNU,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(FRMAX)

    !/

    !/                  +-----------------------------------+

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

    !     ----------------------------------------------------------------

    !       FRMAX   Real  I   maximum frequency.

    !     ----------------------------------------------------------------

    !

    !  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, tpi

    USE w3gdatmd, ONLY: aalpha, bbeta, zzalp, xfr, fachfe, zz0max

#ifdef W3_T

    USE w3odatmd, ONLY: ndst

#endif

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

    REAL, intent(in) :: FRMAX  !  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 :

    !

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

    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                 :: I,J,K,L

    REAL                    :: X0

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

    CALL strace (ient, 'TABU_HF')

#endif

    !

    ustarm = 5.

    alpham = 20.*aalpha

    delust = ustarm/real(iustar)

    delalp = alpham/real(ialpha)

    const1 = bbeta/kappa**2

    omegac = tpi*frmax

    !

    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

          tauhft(k,l)  = tauhft(k,l)+w(j)*zbeta/y*dely

        END DO

        !IF (MOD(K,5).EQ.0.AND.MOD(L,5).EQ.0) &

        !WRITE(102,'(2I4,3G16.8)') L,K,UST,AALPHA+FLOAT(L)*DELALP,TAUHFT(K,L)

#ifdef W3_T

        WRITE (ndst,9000) l,k,aalpha+float(l)*delalp,ust,tauhft(k,l)

#endif

      END DO

    END DO

    DEALLOCATE(w)

    !      DO L=0,IALPHA

    !         DO K=0,IUSTAR

    !WRITE(101,'(A,2I4,G16.8)') 'L,K,TAUHFT(K,L):',L,K,TAUHFT(K,L)

    !         END DO

    !        END DO

    !WRITE(101,*) 'TAUHFT:',FRMAX,BBETA,AALPHA,CONST1,OMEGAC,TPI

    !WRITE(101,'(20G16.8)') TAUHFT

    RETURN

#ifdef W3_T

9000 FORMAT (' TABU_HF, L, K :',(2i4,3f8.3)/)

#endif

  END SUBROUTINE tabu_tauhf

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


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


  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 )

    !

    !  1. Purpose :

    !

    !     Compute friction velocity based on wind speed U10

    !

    !  2. Method :

    !

    !     Computation of u* based on Quasi-linear theory

    !

    !  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

    USE w3gdatmd,  ONLY: zzwnd, aalpha

    IMPLICIT NONE

    REAL, intent(in) :: WINDSPEED,TAUW

    REAL, intent(out) :: USTAR, Z0, CHARN

    ! local variables

    REAL SQRTCDM1

    REAL X,XI,DELI1,DELI2,XJ,delj1,delj2

    REAL UST,DELTOLD,TAUW_LOCAL

    INTEGER IND,J

    !

    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

    !

    ! Determines roughness length

    !

    sqrtcdm1  = min(windspeed/ustar,100.0)

    z0  = zzwnd*exp(-kappa*sqrtcdm1)

    IF (ustar.GT.0.001) THEN

      charn = grav*z0/ustar**2

    ELSE

      charn = aalpha

    END IF

    !

    RETURN

  END SUBROUTINE calc_ustar

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

  SUBROUTINE w3sds3 (A, K, CG, EMEAN, FMEAN, WNMEAN, USTAR, USDIR,    &

       DEPTH, S, D, IX, IY)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  !            F. Ardhuin             !

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         23-Jul-2009 |

    !/                  +-----------------------------------+

    !/

    !/    05-Dec-1996 : Final FORTRAN 77                    ( version 1.18 )

    !/    08-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )

    !/    14-Aug-2006 : Generic WAM4+ dissipation term      ( version 2.22SHOM )

    !/    23-Jul-2009 : Addition of Filipot &al convolution ( version 3.14-SHOM )

    !/

    !  1. Purpose :

    !

    !     Calculate whitecapping source term and diagonal term of derivative.

    !

    !  2. Method :

    !

    !       WAM-Cycle 4 and following.

    !       The last update (09-May-2005) follows the redefinition of

    !       the mean wavenumber as in Bidlot et al. (2005).

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       A       R.A.  I   Action density spectrum (1-D).

    !       K       R.A.  I   Wavenumber for entire spectrum.          *)

    !       EMEAN   Real  I   Mean wave energy.

    !       FMEAN   Real  I   Mean wave frequency.

    !       WNMEAN  Real  I   Mean wavenumber.

    !       USTAR   Real  I   Friction velocity.

    !       USDIR   Real  I   wind stress direction.

    !       DEPTH   Real  I   Water depth.

    !       S       R.A.  O   Source term (1-D version).

    !       D       R.A.  O   Diagonal term of derivative.             *)

    !     ----------------------------------------------------------------

    !                         *) 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, tpi

    USE w3gdatmd, ONLY: nspec, nth, nk, ddelta1, ddelta2,   &

#ifdef w3_t0

         sig,                                &

#endif

         ssdsc1

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

#ifdef W3_T

    USE w3odatmd, ONLY: ndst

#endif

#ifdef W3_T1

    USE w3odatmd, ONLY: ndst

#endif

#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), K(NK), CG(NK),                     &

         DEPTH, USTAR, USDIR, EMEAN, FMEAN, WNMEAN

    INTEGER, INTENT(IN)     :: IX, IY

    REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC)

    !/

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

    !/ Local parameters

    !/

    INTEGER                 :: IS, IK, ITH

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    REAL                    :: FACTOR, FACTOR2

    REAL                    :: ALFAMEAN, WNMEAN2

#ifdef W3_T0

    REAL                    :: DOUT(NK,NTH)

#endif

    !/

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

    !/

#ifdef W3_S

    CALL strace (ient, 'W3SDS3')

#endif

    !

    ! 0.  Pre-initialization of arrays, should be set before being used

    !     but this is helping with bit reproducibility

    d=0.


    ! 1.  Common factor

    !

    wnmean2 = max( 1.e-10 , wnmean  )

    alfamean=wnmean**2*emean

    factor = ssdsc1 * tpi*fmean * alfamean**2

    !

#ifdef W3_T

    WRITE (ndst,9000) ssdsc1, fmean, wnmean, emean, factor

#endif

    !

    !----------------------------------------------------------------------

    !

    ! 2.  Source term

    !

    DO  ik=1, nk

      !

      !    Original WAM4/WAM4+ dissipation term

      !

      factor2=factor*(ddelta1*k(ik)/wnmean2 + ddelta2*(k(ik)/wnmean2)**2)

      DO ith=1,nth

        is=ith+(ik-1)*nth

        d(is)=  factor2

      END DO

    END DO

    !

    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:nk), '  ', 1.,      &

         0.0, 0.001, 'Diag Sds', ' ', 'NONAME')

#endif

    !

#ifdef W3_T1

    CALL outmat (ndst, d, nth, nth, nk, 'diag Sds')

#endif

    !

    RETURN

    !

    ! Formats

    !

#ifdef W3_T

9000 FORMAT (' TEST W3SDS3 : COMMON FACT.: ',5e10.3)

#endif

    !/

    !/ End of W3SDS3 ----------------------------------------------------- /

    !/

  END SUBROUTINE w3sds3


END MODULE w3src3md