#include "w3macros.h"
!/ ------------------------------------------------------------------- /
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
! ----------------------------------------------------------------
! GRAV Real Global Acc. of gravity (m/s2)
! DWAT Real Global Density of water (kg/m3)
! DAIR Real Global Density of air (kg/m3)
! NU_AIR Real Global Kinematic viscosity of air (m2/s)
! NU_WATER Real Global Kinematic viscosity of water (m2/s)
! SED_SG Real Global Specific gravity of sediments (N.D.)
! KAPPA Real Global von Karman's constant (N.D.)
! PI Real Global pi.
! TPI Real Global 2pi.
! HPI Real Global 0.5pi.
! TPIINV Real Global 1/2pi.
! HPIINV Real Global 2/pi.
! RADE Real Global Conv. factor from radians to degrees.
! DERA Real Global Conv. factor from degrees to radians.
! RADIUS Real Global Radius of the earth. (m)
! TSTOUT Log. Global Flag for generation of test files.
! UNDEF Real Global Value for undefined variable in output
! ----------------------------------------------------------------
!
! 5. Remarks
!
! - 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 !
!
!/ ------------------------------------------------------------------- /
!/
LOGICAL, PARAMETER :: TSTOUT = .FALSE.
!
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)
REAL, PARAMETER :: nu_water = 1.31E-6 !mdo WAS: 3.E-6
REAL, PARAMETER :: sed_sg = 2.65
REAL, PARAMETER :: KAPPA = 0.40 !Von Karman's constant
!
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
!
! 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
!
! 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, 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,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 20 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
20 CONTINUE
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 40 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
40 CONTINUE
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 60 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
60 CONTINUE
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