WW3/constants_8F90_source.html

#include "w3macros.h"


!

#ifndef ENDIANNESS

#define ENDIANNESS "native"

#endif

!

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

MODULE constants

  !/

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

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |           H. L. Tolman            |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         05-Jun-2018 |

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

  !/

  !/    11-Nov-1999 : Fortran 90 version.                 ( version 2.00 )

  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )

  !/    25-Jun-2011 : Adding Kelvin functions.            ( version 4.05 )

  !/    03-Sep-2012 : Adding TSTOUT flag.                 ( version 4.10 )

  !/    28-Feb-2013 : Adding cap at 0.5 in FWTABLE        ( version 4.08 )

  !/    20-Jan-2017 : Add parameters for ESMF             ( version 6.02 )

  !/    01-Mar-2018 : Add UNDEF parameter                 ( version 6.02 )

  !/    05-Jun-2018 : Add PDLIB parameters                ( version 6.04 )

  !/

  !/    Copyright 2009-2012 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 :

  !

  !     Define some much-used constants for global use (all defined

  !     as PARAMETER).

  !

  !  2. Variables and types :

  !

  !      Name      Type  Scope    Description

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

  !      UNDEF     Real  Global   Value for undefined variable in output

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

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

  !/

  LOGICAL, PARAMETER :: tstout = .false.

  !        The flag for generating test output files is included here as

  !        it is needed in both ww3_shel and ww3_multi at the same time.

  !        Make sure that this flag is true if you want to write to the

  !        test output file !

  REAL, PARAMETER :: grav = 9.806

  REAL, PARAMETER :: dwat = 1000.

  REAL, PARAMETER :: dair = 1.225

  REAL, PARAMETER :: nu_air = 1.4e-5

  !mdo  *** Changing nu_water to be consistent with DWAT=1000 (assumes 10degC)

  !mdo   WAS: 3.E-6

  REAL, PARAMETER :: nu_water = 1.31e-6

  REAL, PARAMETER :: sed_sg = 2.65

  REAL, PARAMETER :: kappa = 0.40

  !

  REAL, PARAMETER :: pi = 3.141592653589793

  REAL, PARAMETER :: tpi = 2.0 * pi

  REAL, PARAMETER :: hpi = 0.5 * pi

  REAL, PARAMETER :: tpiinv = 1. / tpi

  REAL, PARAMETER :: hpiinv = 1. / hpi

  REAL, PARAMETER :: rade = 180. / pi

  REAL, PARAMETER :: dera = pi / 180.

  !

  REAL, PARAMETER :: radius = 4.e7 * tpiinv

  !

  REAL, PARAMETER :: g2pi3i = 1. / ( grav**2 * tpi**3 )

  REAL, PARAMETER :: g1pi1i = 1. / ( grav * tpi )

  !

  REAL            :: undef = -999.9


  CHARACTER(*), PARAMETER :: file_endian = endianness

  !

  ! Parameters for friction factor table

  !

  INTEGER, PARAMETER :: sizefwtable=300

  REAL            :: fwtable(0:sizefwtable)

  REAL            :: delab

  REAL, PARAMETER :: abmin = -1.

  REAL, PRIVATE, PARAMETER :: abmax = 8.

  INTEGER, PARAMETER :: srce_direct = 0

  INTEGER, PARAMETER :: srce_imp_post = 1

  INTEGER, PARAMETER :: srce_imp_pre = 2

  INTEGER, PARAMETER :: debug_node = 1014

  INTEGER, PARAMETER :: debug_element = 50

  LOGICAL            :: lpdlib = .false.

  LOGICAL            :: lsetup = .false.

  !

  ! Parameters in support of running as ESMF component

  !

  ! --- Flag indicating whether or not the model has been invoked as an

  !     ESMF Component.  This flag is set to true in the WMESMFMD ESMF

  !     module during initialization.

  LOGICAL :: is_esmf_component = .false.

  !

CONTAINS

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

  SUBROUTINE tabu_fw

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |            F. Ardhuin             |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         28-Feb-2013 |

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

    !/

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

    !/    28-Feb-2013 : Caps the friction factor to 0.5     ( version 4.08 )

    !/

    !  1. Purpose :

    !     TO estimate friction coefficients in oscillatory boundary layers

    !     METHOD.

    !      tabulation on Kelvin functions

    !

    !  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

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

    !      WW3_GRID  Prog. WW3_GRID Model grid initialization

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

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S  Enable subroutine tracing.

    !

    ! 10. Source code :

    !

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

    IMPLICIT NONE

    INTEGER, PARAMETER      :: NITER=100

    REAL   , PARAMETER      :: XM=0.50

    REAL   , PARAMETER      :: EPS1=0.00001

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

    INTEGER I,ITER

    REAL KER, KEI

    REAL ABR,ABRLOG,L10,FACT,FSUBW,FSUBWMEMO,dzeta0,dzeta0memo

    !

    delab   = (abmax-abmin)/real(sizefwtable)

    l10=alog(10.)

    DO i=0,sizefwtable

      !

      !  index I in this table corresponds to a normalized roughness z0/ABR = 10^ABMIN+REAL(I)*DELAB

      !

      abrlog=abmin+real(i)*delab

      abr=exp(abrlog*l10)

      fact=1/abr/(21.2*kappa)

      fsubw=0.05

      dzeta0=0.

      DO iter=1,niter

        fsubwmemo=fsubw

        dzeta0memo=dzeta0

        dzeta0=fact*fsubw**(-.5)

        CALL kerkei(2.*sqrt(dzeta0),ker,kei)

        fsubw=.08/(ker**2+kei**2)

        fsubw=.5*(fsubwmemo+fsubw)

        dzeta0=.5*(dzeta0memo+dzeta0)

      END DO

      !

      ! Maximum value of 0.5 for fe is based on field

      ! and lab experiment by Lowe et al. JGR 2005, 2007

      !

      fwtable(i)  = min(fsubw,0.5)

      !           WRITE(994,*) 'Friction factor:',I,ABR,FWTABLE(I)

    END DO

    RETURN

  END SUBROUTINE tabu_fw


  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  SUBROUTINE kzeone(X, Y, RE0, IM0, RE1, IM1)

    !  June 1999 adaptation to CRESTb, all tests on range of (x,y) have been

    !  bypassed, we implicitly expect X to be positive or |x,y| non zero

    !

    ! This subroutine is copyright by ACM

    ! see http://www.acm.org/pubs/copyright_policy/softwareCRnotice.html

    ! ACM declines any responsibility of any kind

    !

    ! THE VARIABLES X AND Y ARE THE REAL AND IMAGINARY PARTS OF

    ! THE ARGUMENT OF THE FIRST TWO MODIFIED BESSEL FUNCTIONS

    ! OF THE SECOND KIND,K0 AND K1.  RE0,IM0,RE1 AND IM1 GIVE

    ! THE REAL AND IMAGINARY PARTS OF EXP(X)*K0 AND EXP(X)*K1,

    ! RESPECTIVELY.  ALTHOUGH THE REAL NOTATION USED IN THIS

    ! SUBROUTINE MAY SEEM INELEGANT WHEN COMPARED WITH THE

    ! COMPLEX NOTATION THAT FORTRAN ALLOWS, THIS VERSION RUNS

    ! ABOUT 30 PERCENT FASTER THAN ONE WRITTEN USING COMPLEX

    ! VARIABLES.

    ! ACM Libraries

    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    IMPLICIT NONE

    DOUBLE PRECISION X, Y, X2, Y2, RE0, IM0, RE1, IM1, &

         R1, R2, T1, T2, P1, P2, RTERM, ITERM, L

    DOUBLE PRECISION , PARAMETER, DIMENSION(8) :: EXSQ = &

         (/ 0.5641003087264d0,0.4120286874989d0,0.1584889157959d0, &

         0.3078003387255d-1,0.2778068842913d-2,0.1000044412325d-3, &

         0.1059115547711d-5,0.1522475804254d-8 /)

    DOUBLE PRECISION , PARAMETER, DIMENSION(8) :: TSQ = &

         (/ 0.0d0,3.19303633920635d-1,1.29075862295915d0, &

         2.95837445869665d0,5.40903159724444d0,8.80407957805676d0, &

         1.34685357432515d1,2.02499163658709d1 /)

    INTEGER N,M,K

    ! THE ARRAYS TSQ AND EXSQ CONTAIN THE SQUARE OF THE

    ! ABSCISSAS AND THE WEIGHT FACTORS USED IN THE GAUSS-

    ! HERMITE QUADRATURE.

    r2 = x*x + y*y

    IF (r2.GE.1.96d2) GO TO 50

    IF (r2.GE.1.849d1) GO TO 30

    ! THIS SECTION CALCULATES THE FUNCTIONS USING THE SERIES

    ! EXPANSIONS

    x2 = x/2.0d0

    y2 = y/2.0d0

    p1 = x2*x2

    p2 = y2*y2

    t1 = -(dlog(p1+p2)/2.0d0+0.5772156649015329d0)

    ! THE CONSTANT IN THE PRECEDING STATEMENT IS EULER*S

    ! CONSTANT

    t2 = -datan2(y,x)

    x2 = p1 - p2

    y2 = x*y2

    rterm = 1.0d0

    iterm = 0.0d0

    re0 = t1

    im0 = t2

    t1 = t1 + 0.5d0

    re1 = t1

    im1 = t2

    p2 = dsqrt(r2)

    l = 2.106d0*p2 + 4.4d0

    IF (p2.LT.8.0d-1) l = 2.129d0*p2 + 4.0d0

    DO n=1,int(l)

      p1 = n

      p2 = n*n

      r1 = rterm

      rterm = (r1*x2-iterm*y2)/p2

      iterm = (r1*y2+iterm*x2)/p2

      t1 = t1 + 0.5d0/p1

      re0 = re0 + t1*rterm - t2*iterm

      im0 = im0 + t1*iterm + t2*rterm

      p1 = p1 + 1.0d0

      t1 = t1 + 0.5d0/p1

      re1 = re1 + (t1*rterm-t2*iterm)/p1

      im1 = im1 + (t1*iterm+t2*rterm)/p1

    END DO

    r1 = x/r2 - 0.5d0*(x*re1-y*im1)

    r2 = -y/r2 - 0.5d0*(x*im1+y*re1)

    p1 = dexp(x)

    re0 = p1*re0

    im0 = p1*im0

    re1 = p1*r1

    im1 = p1*r2

    RETURN

    ! THIS SECTION CALCULATES THE FUNCTIONS USING THE INTEGRAL

    ! REPRESENTATION, EQN 3, EVALUATED WITH 15 POINT GAUSS-

    ! HERMITE QUADRATURE

30  x2 = 2.0d0*x

    y2 = 2.0d0*y

    r1 = y2*y2

    p1 = dsqrt(x2*x2+r1)

    p2 = dsqrt(p1+x2)

    t1 = exsq(1)/(2.0d0*p1)

    re0 = t1*p2

    im0 = t1/p2

    re1 = 0.0d0

    im1 = 0.0d0

    DO n=2,8

      t2 = x2 + tsq(n)

      p1 = dsqrt(t2*t2+r1)

      p2 = dsqrt(p1+t2)

      t1 = exsq(n)/p1

      re0 = re0 + t1*p2

      im0 = im0 + t1/p2

      t1 = exsq(n)*tsq(n)

      re1 = re1 + t1*p2

      im1 = im1 + t1/p2

    END DO

    t2 = -y2*im0

    re1 = re1/r2

    r2 = y2*im1/r2

    rterm = 1.41421356237309d0*dcos(y)

    iterm = -1.41421356237309d0*dsin(y)

    ! THE CONSTANT IN THE PREVIOUS STATEMENTS IS,OF COURSE,

    ! SQRT(2.0).

    im0 = re0*iterm + t2*rterm

    re0 = re0*rterm - t2*iterm

    t1 = re1*rterm - r2*iterm

    t2 = re1*iterm + r2*rterm

    re1 = t1*x + t2*y

    im1 = -t1*y + t2*x

    RETURN

    ! THIS SECTION CALCULATES THE FUNCTIONS USING THE

    ! ASYMPTOTIC EXPANSIONS

50  rterm = 1.0d0

    iterm = 0.0d0

    re0 = 1.0d0

    im0 = 0.0d0

    re1 = 1.0d0

    im1 = 0.0d0

    p1 = 8.0d0*r2

    p2 = dsqrt(r2)

    l = 3.91d0+8.12d1/p2

    r1 = 1.0d0

    r2 = 1.0d0

    m = -8

    k = 3

    DO n=1,int(l)

      m = m + 8

      k = k - m

      r1 = float(k-4)*r1

      r2 = float(k)*r2

      t1 = float(n)*p1

      t2 = rterm

      rterm = (t2*x+iterm*y)/t1

      iterm = (-t2*y+iterm*x)/t1

      re0 = re0 + r1*rterm

      im0 = im0 + r1*iterm

      re1 = re1 + r2*rterm

      im1 = im1 + r2*iterm

    END DO

    t1 = dsqrt(p2+x)

    t2 = -y/t1

    p1 = 8.86226925452758d-1/p2

    ! THIS CONSTANT IS SQRT(PI)/2.0, WITH PI=3.14159...

    rterm = p1*dcos(y)

    iterm = -p1*dsin(y)

    r1 = re0*rterm - im0*iterm

    r2 = re0*iterm + im0*rterm

    re0 = t1*r1 - t2*r2

    im0 = t1*r2 + t2*r1

    r1 = re1*rterm - im1*iterm

    r2 = re1*iterm + im1*rterm

    re1 = t1*r1 - t2*r2

    im1 = t1*r2 + t2*r1

    RETURN

  END SUBROUTINE kzeone


  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  SUBROUTINE kerkei(X,KER,KEI)

    !**********************************************************************

    ! Computes the values of the zeroth order Kelvin function Ker and Kei

    ! These functions are used to determine the friction factor fw as a

    ! function of the bottom roughness length assuming a linear profile

    ! of eddy viscosity (See Grant and Madsen, 1979)

    !**********************************************************************

    IMPLICIT NONE


    DOUBLE PRECISION ZR,ZI,CYR,CYI,CYR1,CYI1

    REAL X,KER,KEI


    zr=x*.50d0*sqrt(2.0d0)

    zi=zr

    CALL kzeone(zr, zi, cyr, cyi,cyr1,cyi1)

    ker=cyr/exp(zr)

    kei=cyi/exp(zr)

  END SUBROUTINE kerkei

  !/

  !/ End of module CONSTANTS ------------------------------------------- /

  !/

END MODULE constants