w3src2md.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
 
!/ ------------------------------------------------------------------- /
MODULE w3src2md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         29-May-2009 |
  !/                  +-----------------------------------+
  !/
  !/    04-Feb-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
  !/    21-Feb-2004 : Multiple model version.             ( version 3.06 )
  !/    03-Jul-2006 : Extract stress computation.         ( version 3.09 )
  !/    13-Apr-2007 : EMEAN in W3SPR2 par list.           ( version 3.11 )
  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
  !/
  !/    Copyright 2009 National Weather Service (NWS),
  !/       National Oceanic and Atmospheric Administration.  All rights
  !/       reserved.  WAVEWATCH III is a trademark of the NWS.
  !/       No unauthorized use without permission.
  !/
  !  1. Purpose :
  !
  !     Tolman and Chalikov (1996) input and dissipation source terms.
  !     Bundled with interpolation tables.
  !
  !  2. Variables and types :
  !
  !      Interpolation tables :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      NRSIGA    I.P.  Public   Array dimension (SIGA).
  !      NRDRAG    I.P.  Public   Array dimension (drag coefficient).
  !      SIGAMX    R.P.  Public   Maximum nondiensional frequency.
  !      DRAGMX    R.P.  Public   Maximum drag coefficient.
  !      DSIGA     Real  Public   Table increment.
  !      DDRAG     Real  Public   Id.
  !      BETATB    R.A.  Public   Interpolation table.
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3SPR2    Subr. Public   Mean parameters from spectrum.
  !      W3SIN2    Subr. Public   Input source term.
  !      W3SDS2    Subr. Public   Dissipation source term.
  !      INPTAB    Subr. Public   Interpolation table for wind-wave
  !                               interaction parameter.
  !      W3BETA    R.F.  INPTAB   Id. function.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.            ( !/S )
  !      PRT2DS    Subr. W3ARRYMD Print plot of spectra.        ( !/T0 )
  !      OUTMAT    Subr. W3WRRYMD Print out 2D matrix.          ( !/T1 )
  !      ...       Data  W3DISPMD Interpolation tables to solve
  !                               dispersion relation.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !       !/S      Enable subroutine tracing.
  !       !/T(n)   Test output, see subroutines.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  PUBLIC
  !/
  !/ Interpolation table
  !/
  INTEGER, PARAMETER, PRIVATE :: NRSIGA =  400
  INTEGER, PARAMETER, PRIVATE :: NRDRAG =   20
  REAL, PARAMETER, PRIVATE    :: SIGAMX =   40.
  REAL, PARAMETER, PRIVATE    :: DRAGMX =    1.e-2
  !
  REAL, PRIVATE           :: DSIGA, DDRAG,                        &
       BETATB(-NRSIGA:NRSIGA+1,NRDRAG+1)
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3spr2 (A, CG, WN, DEPTH, FPI, U, USTAR,             &
       EMEAN, FMEAN, WNMEAN, AMAX, ALFA, FP )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |            D.Chalikov             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Apr-2007 |
    !/                  +-----------------------------------+
    !/
    !/    06-Dec-1996 : Final version 1.18 / FORTRAN 77 version.
    !/    16-Nov-1999 : Add itteration to section 5. for removal of W3APR2.
    !/    04-Feb-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    21-Dec-2004 : Multiple model version.             ( version 3.06 )
    !/    03-Jul-2006 : Extract stress computation.         ( version 3.09 )
    !/    13-Apr-2007 : EMEAN in parameter list.            ( version 3.11 )
    !
    !  1. Purpose :
    !
    !     Calculate mean wave parameters for the use in the source term
    !     routines. (Tolman and Chalikov)
    !
    !  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.
    !       DEPTH   Real  I   Water depth.
    !       FPI     Real  I   Peak input frequency.
    !       U       Real  I   Wind speed.
    !       USTAR   Real  I   Friction velocity.
    !       EMEAN   Real  O   Total energy (variance).
    !       FMEAN   Real  O   Mean frequency.
    !       WNMEAN  Real  O   Mean wavenumber.
    !       AMAX    Real  O   Maximum of action spectrum.
    !       ALFA    R.A.  O   Phillips' constant.
    !       FP      Real  O   Peak frequency.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SRCE    Subr. W3SRCEMD Source term integration.
    !      W3EXPO    Subr.   N/A    Point output post-processor.
    !      GXEXPO    Subr.   N/A    GrADS point output post-processor.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !       !/T      Enable test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nth, dth, sig, dden, fte, ftf, ftwn,    &
         nittin, zwind, cinxsi
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
 
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    USE w3dispmd, ONLY: nar1d, dfac, n1max, ecg1, ewn1, dsie
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: A(NTH,NK), CG(NK), WN(NK), DEPTH,    &
         FPI, U, USTAR
    REAL, INTENT(OUT)       :: EMEAN, FMEAN, WNMEAN, AMAX,          &
         ALFA(NK), FP
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IK, ITH, I1, ITT
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: EBAND, FPISTR, EB(NK), UST
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SPR2')
#endif
    !
    ust    = max( 0.0001 , ustar )
    !
    emean  = 0.
    fmean  = 0.
    wnmean = 0.
    amax   = 0.
    !
    ! 1.  Integral over directions and maximum --------------------------- *
    !
    DO ik=1, nk
      eb(ik) = 0.
      DO ith=1, nth
        eb(ik) = eb(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)
      emean    = emean  + eb(ik)
      fmean    = fmean  + eb(ik) / sig(ik)
      wnmean   = wnmean + eb(ik) / sqrt(wn(ik))
    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
    wnmean = wnmean + eband * ftwn
    !
    fmean  = tpiinv * emean / max( 1.e-7 , fmean )
    wnmean = ( emean / max( 1.e-7 , wnmean ) )**2
    !
    ! 4.  Estimate peak frequency from FPI ------------------------------- *
    !
    fpistr = max( 0.008 , fpi * ust / grav )
    fp     = ( 3.6e-4 + 0.92*fpistr - 6.3e-10/fpistr**3 )/ust*grav
    fp     = fp * tpiinv
    !
    RETURN
    !/
    !/ End of W3SPR2 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3spr2
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3sin2 ( A, CG, K, U, UDIR, CD, Z0, FPI, S, D )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |            D.Chalikov             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         21-Feb-2004 |
    !/                  +-----------------------------------+
    !/
    !/    14-Jan-1997 : Final FORTRAN 77                    ( version 1.18 )
    !/    04-Feb-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    21-Feb-2004 : Multiple model version.             ( version 3.06 )
    !/
    !  1. Purpose :
    !
    !     Calculate input source term.
    !
    !  2. Method :
    !
    !     Tolman and Chalikov (1996), see manual.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A       R.A.  I   Action density spectrum (1-D).
    !       CG      R.A.  I   Group velocities for k-axis of spectrum.
    !       K       R.A.  I   Wavenumber for entire spectrum (1-D).
    !       U       Real  I   Wind speed at reference height.
    !       UDIR    Real  I   Direction of U.
    !       CD      Real  I   Drag coefficient at wind level ZWIND.
    !       Z0      Real  I   Corresponding z0.
    !       FPI     R.A.  O   Input 'peak' frequency.
    !       S       R.A.  O   Source term (1-D version).
    !       D       R.A.  O   Diagonal term of derivative (1-D version).
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !      PRT2DS    Subr. W3ARRYMD Print plot of spectra.
    !      OUTMAT    Subr. W3WRRYMD Print out 2D matrix.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SRCE    Subr. W3SRCEMD Source term integration.
    !      W3EXPO    Subr.   N/A    Point output post-processor.
    !      GXEXPO    Subr.   N/A    GrADS point output post-processor.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !     - Actual height of wind speed does not need to be 10 m, but is
    !       given by ZWIND.
    !     - Abs(cos) > 0.0087 to asure continuity in beta. Corresponds
    !       to shift of up to half a degree.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !       !/S      Enable subroutine tracing.
    !       !/T      Enable general test output.
    !       !/T0     Print arrays.
    !       !/T1     Calculation of diagonal within spectrum
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nth, nspec, xfr, dden, sig, sig2,       &
         esin, ecos, fte, fttr, fpimin, zwind,       &
         facti1, facti2, fswell
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
#ifdef W3_S
    USE w3servmd, ONLY: strace
#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), CG(NK), K(NSPEC), U, UDIR, &
         CD, Z0
    REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC), FPI
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IS, IK, IOMA, ICL, NKFILT, NKFIL2
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T0
    INTEGER       ITH
#endif
    REAL                    :: COSU, SINU, COSFAC, LAMBDA, ULAM,    &
         CLAM, OMA, M0, M1, RD1, RD2, BETA,   &
         FACLN1, FACLN2, USTAR, TRANS, FPISTR,&
         FP1STR, FP1, SIN1A(NK)
    REAL, PARAMETER         :: TRANSF = 0.75
    REAL, PARAMETER         :: PEAKFC = 0.8
#ifdef W3_T0
    REAL                    :: DOUT(NK,NTH)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SIN2')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) dsiga, ddrag, u, udir*rade, cd, z0
#endif
    !
    ! 1.  Preparations
    !
    cosu   = cos(udir)
    sinu   = sin(udir)
    !
    ! 2.  Loop over spectrum
    !
#ifdef W3_T2
    WRITE (ndst,9020)
#endif
    !
    facln1 = u / log(zwind/z0)
    facln2 = log(z0)
    !
    DO is=1, nspec
      cosfac = ecos(is)*cosu + esin(is)*sinu
      cosfac = sign( max( 0.0087 , abs(cosfac) ) , cosfac )
      lambda = tpi / ( k(is) * abs(cosfac) )
      ulam   = facln1 * ( log(lambda) - facln2 )
      clam   = cd * ( u / ulam )**2
      oma    = k(is) * ulam * cosfac / sig2(is)
      ioma   = int( oma/dsiga ) +                                  &
           min( 0 , int( sign( -1.1 , oma ) ) )
      icl    = int( clam/ddrag )
      rd1    = oma/dsiga - real(ioma)
      rd2    = clam/ddrag - real(icl)
      ioma   = max( -nrsiga , min( nrsiga , ioma ) )
      icl    = max( 1 , min( nrdrag , icl ) )
      beta   = (1.-rd1) * (1.-rd2) * betatb( ioma , icl )           &
           +    rd1   * (1.-rd2) * betatb(ioma+1, icl )           &
           + (1.-rd1) *    rd2   * betatb( ioma ,icl+1)           &
           +    rd1   *    rd2   * betatb(ioma+1,icl+1)
      d(is)  = beta * sig2(is)
      s(is)  = a(is) * d(is)
#ifdef W3_T2
      WRITE (ndst,9021) is, cosfac, lambda, ulam, clam*1.e3,    &
           oma, beta*1.e4
#endif
    END DO
    !
    ! 3.  Calculate FPI
    !
    DO ik=1, nk
      sin1a(ik) = 0.
      DO is=(ik-1)*nth+1, ik*nth
        sin1a(ik) = sin1a(ik) + max( 0. , s(is) )
      END DO
    END DO
    !
    m0     = 0.
    m1     = 0.
    DO ik=1, nk
      sin1a(ik) = sin1a(ik) * dden(ik) / ( cg(ik) * sig(ik)**3 )
      m0        = m0 + sin1a(ik)
      m1        = m1 + sin1a(ik)/sig(ik)
    END DO
    !
    sin1a(nk) = sin1a(nk) / dden(nk)
    m0        = m0 + sin1a(nk) * fte
    m1        = m1 + sin1a(nk) * fttr
    IF ( m1 .LT. 1e-20 ) THEN
      fpi    = xfr * sig(nk)
    ELSE
      fpi    = m0 / m1
    END IF
    !
    ! 4.  Filter for swell
    !
    ustar  = u * sqrt(cd)
    fpistr = max( fpimin , fpi * ustar / grav )
    fp1str = 3.6e-4 + 0.92*fpistr - 6.3e-10/fpistr**3
    fp1    = peakfc * fp1str * grav / ustar
    !
    nkfilt = min( nk , int(facti2+facti1*log(fp1)) )
    nkfil2 = min( nk , int(facti2+facti1*log(transf*fp1)) )
    nkfil2 = max( 0 , nkfil2 )
    !
    DO is=1, nkfil2*nth
      d(is)  = max( d(is) , fswell*d(is) )
      s(is)  = a(is) * d(is)
    END DO
    !
    DO ik=nkfil2+1, nkfilt
      trans  = ( sig(ik)/fp1 - transf ) / (1.-transf)
      DO is=(ik-1)*nth+1, ik*nth
        d(is)  = (1.-trans)*max(d(is),fswell*d(is)) + trans*d(is)
        s(is)  = a(is) * d(is)
      END DO
    END DO
    !
    ! ... 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
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3SIN2 : DSIGA,DDRAG,U,UDIR,CD,Z0(IN) : '/    &
         '              ',f8.4,f9.6,f7.2,f6.1,f8.5,f8.5)
#endif
    !
#ifdef W3_T2
9020 FORMAT (' TEST W3SIN2 : IS, COS, LAMBDA, ULAM, CLAM*1E3, ', &
         'OMA, BETA*1E4')
9021 FORMAT (6x,i6,f7.2,1x,f6.1,2(1x,f5.2),2(1x,f6.2))
#endif
    !/
    !/ End of W3SIN2 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3sin2
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3sds2 (A, CG, K, FPI, USTAR, ALFA, S, D)
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         21-Feb-2004 |
    !/                  +-----------------------------------+
    !/
    !/    12-Jun-1996 : Final FORTRAN 77                    ( version 1.18 )
    !/    04-Feb-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    23-Apr-2002 : Erick Rogers' fix                   ( version 2.19 )
    !/    21-Feb-2004 : Multiple model version.             ( version 3.06 )
    !/
    !  1. Purpose :
    !
    !     Calculate whitecapping source term and diagonal term of der.
    !
    !  2. Method :
    !
    !     Tolman and Chalikov (1995).
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A       R.A.   I   Input action density spectrum.
    !       CG      R.A.   I   Group velocity array.
    !       K       R.A.   I   Wavenumber array.
    !       FPI     Real   I   'Peak frequency' of input (rad/s).
    !       USTAR   Real   I   Friction velocity (m/s).
    !       ALFA    R.A.   I   Phillips' constant.
    !       S       R.A.   O   Source term (1-D version).
    !       D       R.A.   O   Diagonal term of derivative (1-D version).
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !      PRT2DS    Subr. W3ARRYMD Print plot of spectra.
    !      OUTMAT    Subr. W3WRRYMD Print out 2D matrix.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SRCE    Subr. W3SRCEMD Source term integration.
    !      W3EXPO    Subr.   N/A    Point output post-processor.
    !      GXEXPO    Subr.   N/A    GrADS point output post-processor.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S     Enable subroutine tracing.
    !     !/T     Enable general test output.
    !     !/T0    Print arrays.
    !     !/T1    Print filter and constituents.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nth, sig, dden, dth, fte, fpimin,       &
         facti1, facti2, xf1, xf2, xfh, sdsaln,      &
         cdsa0, cdsa1, cdsa2, cdsb0, cdsb1, cdsb2,   &
         cdsb3
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
#ifdef W3_S
    USE w3servmd, ONLY: strace
#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(NTH,NK), CG(NK), K(NK), FPI,       &
         USTAR, ALFA(NK)
    REAL, INTENT(OUT)       :: S(NTH,NK), D(NTH,NK)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IK, ITH, IKHW
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: FHW, XHW, FPIT, PHI, AF1, AF2,       &
         AFILT, BFILT, CDIST, FILT, POW,      &
         CDISH, CDISP, HW, EHIGH, EBD(NK)
#ifdef W3_T
    REAL          POWMAX
#endif
#ifdef W3_T0
    REAL             DOUT(NK,NTH)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SDS2')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) fpi, ustar
#endif
    !
    ! 1.  Preparations
    ! 1.a HW
    !
    fhw    = xfh*fpi
    xhw    = facti2 + facti1*log(fhw)
    ikhw   = min( nk , int( xhw + 0.5 ) )
    DO ik=ikhw, nk
      ebd(ik) = 0.
      DO ith=1, nth
        ebd(ik) = ebd(ik) + a(ith,ik)
      END DO
    END DO
    !
    IF ( fhw .LT. sig(nk+1) ) THEN
      xhw    = 1. - mod( xhw + 0.5 , 1. )
      IF ( ikhw .EQ. nk ) xhw = max( 0. , xhw - 0.5 )
      hw     = xhw * ebd(ikhw)*dden(ikhw)/cg(ikhw)
      DO ik=ikhw+1, nk
        hw     = hw + ebd(ik)*dden(ik)/cg(ik)
      END DO
      hw     = 4. * sqrt( hw + ebd(nk)/cg(nk)*fte )
    ELSE
      ehigh  = ebd(nk)/cg(nk) * sig(nk)*dth * (sig(nk)/fhw)**5
      hw     = 4. * sqrt( 0.25 * fhw * ehigh )
    END IF
    !
    ! 1.b PHI
    !
    fpit   = max( fpimin , fpi*tpiinv*ustar/grav )
    phi    = cdsb0 + cdsb1*fpit + cdsb2/fpit**cdsb3
    !
    ! 1.c Set-up filter
    !
    af2    = xf2*fpi
    af1    = xf1*fpi
    bfilt  = 1. / ( af2 - af1 )
    afilt  = - bfilt * af1
    !
    ! 1.d Constants
    !
    cdist = - 2. * ustar * hw * phi
    cdish = g2pi3i * ustar**2
    cdisp = g1pi1i * ustar
    !
    ! 2.  Combined diagonal factor
    !
#ifdef W3_T2
    WRITE (ndst,9020)
#endif
#ifdef W3_T
    powmax = 0.
#endif
    DO ik=1, nk
      filt    = min( 1., max( 0. , afilt + bfilt*sig(ik) ))
      pow     = min( 25. , cdsa1 / ( cdisp*sig(ik) )**cdsa2 )
      IF ( filt .GT. 0. ) THEN
        d(1,ik) = (1.-filt)  * cdist * k(ik)**2                   &
             - filt * cdsa0 * cdish * sig(ik)**3             &
             * (alfa(ik)/sdsaln)**pow
      ELSE
        d(1,ik) = (1.-filt)  * cdist * k(ik)**2
      END IF
#ifdef W3_T
      powmax = max(pow*filt,powmax)
#endif
#ifdef W3_T2
      WRITE (ndst,9021) ik, filt, alfa(ik)/sdsaln,              &
           cdist*phi*k(ik)**2, cdsa0*cdish*sig(ik)**3         &
           * (alfa(ik)/sdsaln)**pow, d(1,ik)
#endif
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9010) af1, af2, afilt, bfilt, powmax
#endif
    !
    ! 3.  2-D diagonal array
    !
    DO ik=1, nk
      DO ith=2, nth
        d(ith,ik) = d(1,ik)
      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)
      END DO
    END DO
    CALL prt2ds (ndst, nk, nk, nth, dout, sig(1:), '  ', 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 W3SDS2 : FPI, USTAR           : ',2f8.3)
9010 FORMAT (' TEST W3SDS2 : AF1-2, A-BFILT, PMAX : ',4f7.3,e10.3)
#endif
#ifdef W3_T2
9020 FORMAT (' TEST W3SDS2 : IK, FILT, ALFA, DDST, DDSH, DDS')
9021 FORMAT ('           ',i6,2f7.3,3e11.3)
#endif
    !/
    !/ End of W3SDS2 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3sds2
  !/ ------------------------------------------------------------------- /
  SUBROUTINE inptab
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         21-Feb-2004 |
    !/                  +-----------------------------------+
    !/
    !/    03-Jun-1996 : Final version 1.18 / FORTRAN 77 version.
    !/    06-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    21-Feb-2004 : Multiple model version.             ( version 3.06 )
    !/
    !  1. Purpose :
    !
    !     Generate an interpolation table for the air-sea interaction
    !     parameter of Chalikov and Belevich (1993).
    !
    !  2. Method :
    !
    !     The size of the table is set in parameter statements, the range
    !     is set by the input parameters of this routine. The first counter
    !     of the table corresponds to the nondimensional frequency
    !
    !                  SIGMA Ul
    !        SIGA  =  ----------  COS ( THETA - THETA     )           (1)
    !                     g                          wind
    !
    !     The second counter of the table represents the drag coefficient.
    !     The maximum values of both parameters are passed to the routine
    !     through the parameter list.
    !
    !  3. Parameters :
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !      W3BETA    Func. Internal Function to calculate the
    !                               interaction parameter.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3IOGR    Subr. W3IOGRMD Model definition IO routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable test output.
    !     !/T0  Print table.
    !     !/T1  Estimate maximum errors.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3odatmd, ONLY: ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ISIGA, IDRAG
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T0
    INTEGER                 :: I1
#endif
#ifdef W3_T1
    INTEGER                 :: IE1
#endif
    REAL                    :: SIGA, DRAG
#ifdef W3_T0
    REAL                    :: BMIN, BMAX
#endif
#ifdef W3_T1
    REAL                    :: ENORM, ERR(NRDRAG)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'INPTAB')
#endif
    !
    ! 1.  Determine range and increments of table ------------------------ *
 
    !
    dsiga  = sigamx / real(nrsiga)
    ddrag  = dragmx / real(nrdrag)
    !
#ifdef W3_T
    WRITE (ndst,9000) sigamx, dsiga, dragmx, ddrag
#endif
    !
    ! 2.  Fill table ----------------------------------------------------- *
    !
    DO isiga=-nrsiga,nrsiga+1
      siga   = real(isiga) * dsiga
      DO idrag=1, nrdrag+1
        drag   = real(idrag) * ddrag
        betatb(isiga,idrag) = w3beta( siga, drag , ndst )
      END DO
    END DO
    !
    ! 3.  Test output ---------------------------------------------------- *
    !
#ifdef W3_T0
    WRITE (ndst,9010)
    i1     = min(35,nrdrag)
    DO isiga=-nrsiga,nrsiga
      siga   = real(isiga) * dsiga
      bmin   = 0.
      bmax   = 0.
      DO idrag=1, nrdrag
        bmin   = min( bmin , betatb(isiga,idrag) )
        bmax   = max( bmax , betatb(isiga,idrag) )
      END DO
      bmax   = max( bmax , -bmin )
      WRITE (ndst,9011) isiga, siga, bmax,                      &
           (nint(betatb(isiga,idrag)/bmax*100.),idrag=1,i1)
      IF (i1.LT.nrdrag) WRITE (ndst,9012)                       &
           (nint(betatb(isiga,idrag)/bmax*100.),idrag=i1+1,nrdrag)
    END DO
#endif
    !
#ifdef W3_T1
    WRITE (ndst,9020)
    ie1    = min(30,nrdrag-1)
    enorm  = 1000. / abs(betatb(0,nrdrag))
    DO isiga=-nrsiga,nrsiga
      siga   = real(isiga) * dsiga
      IF ( abs(siga) .LT. 5.01 ) THEN
        DO idrag=1, nrdrag-1
          drag   = ( real(idrag) + 0.5 ) * ddrag
          err(idrag) =  - w3beta(siga,drag,ndst) + 0.5 *     &
               ( betatb(isiga,idrag) + betatb(isiga,idrag+1) )
        END DO
        WRITE (ndst,9021) isiga, siga,                        &
             (nint(enorm*err(idrag)),idrag=1,ie1)
        IF (ie1.LT.nrdrag-1) WRITE (ndst,9022)                &
             (nint(enorm*err(idrag)),idrag=ie1+1,nrdrag-1)
      ENDIF
    END DO
    !
    WRITE (ndst,9030)
    ie1    = min(30,nrdrag)
    enorm  = 1000. / abs(betatb(0,nrdrag))
    DO isiga=-nrsiga,nrsiga-1
      siga   = ( real(isiga) + 0.5 ) * dsiga
      IF ( abs(siga) .LT. 5.01 ) THEN
        DO idrag=1, nrdrag
          drag   = real(idrag) * ddrag
          err(idrag) =  - w3beta(siga,drag,ndst) + 0.5 *     &
               ( betatb(isiga,idrag) + betatb(isiga+1,idrag) )
        END DO
        WRITE (ndst,9031) isiga, siga,                        &
             (nint(enorm*err(idrag)),idrag=1,ie1)
        IF (ie1.LT.nrdrag) WRITE (ndst,9032)                  &
             (nint(enorm*err(idrag)),idrag=ie1+1,nrdrag)
      ENDIF
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT ( ' TEST INPTAB : SIGAMX, DSIGA : ',f6.2,f8.2/        &
         '               DRAGMX, DDRAG : ',f8.4,f9.5)
#endif
    !
#ifdef W3_T0
9010 FORMAT (/' TEST INPTAB : TABLE, NORMALIZED WITH ',          &
         'BETATB(ISIGA,NRDRAG)'/                            &
         '               ISIGA, SIGA, BETA_MAX, TABLE (x100)')
9011 FORMAT (1x,i4,f7.2,f6.4,1x,35i3)
9012 FORMAT (19x,35i3)
#endif
    !
#ifdef W3_T1
9020 FORMAT (/' TEST INPTAB : ERROR DUE TO DRAG, NORMALIZED ',   &
         'WITH BETATB(ISIGA,NRDRAG)'/                       &
         '               ISIGA, SIGA, TABLE (x1000)')
9021 FORMAT (1x,i4,f7.2,35i3)
9022 FORMAT (12x,35i3)
9030 FORMAT (/' TEST INPTAB : ERROR DUE TO SIGA, NORMALIZED WITH ',  &
         'BETATB(ISIGA,NRDRAG)'/                            &
         '               ISIGA, SIGA, TABLE (x1000)')
9031 FORMAT (1x,i4,f7.2,35i3)
9032 FORMAT (12x,35i3)
#endif
    !/
    !/    Internal function W3BETA
    !/
  CONTAINS
    !/ ------------------------------------------------------------------- /
    REAL FUNCTION W3BETA ( OMA , CL , NDST )
      !/
      !/                  +-----------------------------------+
      !/                  | WAVEWATCH III           NOAA/NCEP |
      !/                  |           H. L. Tolman            |
      !/                  |            D.Chalikov             |
      !/                  |                        FORTRAN 90 |
      !/                  | Last update :         21-Feb-2004 |
      !/                  +-----------------------------------+
      !/
      !/    06-Dec-1996 : Final version 1.18 / FORTRAN 77 version.
      !/    06-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
      !/    21-Feb-2004 : Multiple model version.             ( version 3.06 )
      !/
      !  1. Purpose :
      !
      !     Calculate wind-wave interaction parameter beta.
      !
      !  2. Method :
      !
      !     Chalikov and Belevich (1992), see also manual.
      !
      !  3. Parameters :
      !
      !     Parameter list
      !     ----------------------------------------------------------------
      !       W3BETA  Real  O   Wind-wave interaction parameter multiplied
      !                         by density ratio.
      !       OMA     Real  I   Non-dimensional apparent frequency.
      !
      !                         OMA = OMEGA | U | cos(theta-theta ) / g
      !                                        l                 w
      !
      !       CL      Real  I   Drag coefficient at height l
      !     ----------------------------------------------------------------
      !
      !  4. Subroutines used :
      !
      !  5. Called by :
      !
      !  6. Error messages :
      !
      !  7. Remarks :
      !
      !  8. Structure :
      !
      !     See source code.
      !
      !  9. Switches :
      !
      !     !/S   Enable subroutine tracing.
      !     !/T0  Enable test output.
      !
      ! 10. Source code :
      !
      !/ ------------------------------------------------------------------- /
      IMPLICIT NONE
      !/
      !/ ------------------------------------------------------------------- /
      !/ Parameter list
      !/
      INTEGER, INTENT(IN)     :: ndst
      REAL, INTENT(IN)        :: oma, cl
      !/
      !/ ------------------------------------------------------------------- /
      !/ Local parameters
      !/
#ifdef W3_S
      INTEGER, SAVE           :: ient = 0
#endif
      REAL                    :: om1, om2, a0, a1, a2, a3, a4, a5,    &
           a6, a7, a8, a9, a10
      !/
      !/ ------------------------------------------------------------------- /
      !/
#ifdef W3_S
      CALL strace (ient, 'W3BETA')
#endif
      !
#ifdef W3_T0
      WRITE (ndst,9000) oma, cl
#endif
      !
      ! calculate Omegas
      !
      om1    =  1.075 +  75.*cl
      om2    =  1.2   + 300.*cl
      !
      ! calculate factors a
      !
      a1     =  0.25  + 395.*cl
      a2     =  0.35  + 150.*cl
      a4     =  0.3   + 300.*cl
      a9     =  0.35  + 240.*cl
      a10    = -0.06  + 470.*cl
      !
      a5     =  a4 * om1
      a0     =  0.25 * a5**2 / a4
      a3     = (a0-a2-a1) / (a0+a4+a5)
      a6     =  a0 * (1.-a3)
      a7     = (a9*(om2-1)**2+a10) / (om2-om1)
      a8     =  a7 * om1
      !
#ifdef W3_T0
      WRITE (ndst,9001) om1, om2
      WRITE (ndst,9002) a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10
#endif
      !
      ! calculate beta * 1.e4
      !
      IF  ( oma .LT. -1. ) THEN
        w3beta = -a1 * oma**2 - a2
      ELSE IF (oma .LT. om1/2.) THEN
        w3beta =  a3 * oma * ( a4 * oma - a5 ) - a6
      ELSE IF (oma .LT. om1) THEN
        w3beta =       oma * ( a4 * oma - a5 )
      ELSE IF (oma .LT. om2) THEN
        w3beta = a7 * oma - a8
      ELSE
        w3beta = a9 * (oma-1.)**2 + a10
      END IF
      !
      ! beta * dwat / dair
      !
      w3beta = w3beta * 1.e-4
#ifdef W3_T0
      WRITE (ndst,9003) w3beta
#endif
      !
      RETURN
      !
      ! Formats
      !
#ifdef W3_T0
9000  FORMAT ( ' TEST W3BETA : INPUT : ',2e10.3)
9001  FORMAT ( ' TEST W3BETA : OM1-2 : ',2e10.3)
9002  FORMAT ( ' TEST W3BETA : A0-10 : ',5e10.3/                   &
           '             ',6e10.3)
9003  FORMAT ( ' TEST W3BETA : BETA  : ',e10.3)
#endif
      !/
      !/ End of W3BETA ----------------------------------------------------- /
      !/
    END FUNCTION w3beta
    !/
    !/ End of INPTAB ----------------------------------------------------- /
    !/
  END SUBROUTINE inptab
  !/
  !/ End of module W3SRC2MD -------------------------------------------- /
  !/
END MODULE w3src2md