#include "w3macros.h"
!/ ------------------------------------------------------------------- /
      MODULE CONSTANTS
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         20-Jan-2017 |
!/                  +-----------------------------------+
!/
!/    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 )
!/
!/    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.
!     ----------------------------------------------------------------
!
!  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.1415927
      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 )
!
! Parameters for friction factor table 
!
      INTEGER, PARAMETER       :: SIZEFWTABLE=300  
      REAL                     :: FWTABLE(0:SIZEFWTABLE)
      REAL                     :: DELAB
      REAL,    PARAMETER       :: ABMIN = -1. 
      REAL, PRIVATE, PARAMETER  :: ABMAX = 8. 
!
! 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
   INTEGER NZ,IERR
   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
