anchor.diff:     REAL(ReKi), PARAMETER        :: SmllNmbr  = 9.999E-4                            ! A small number representing epsilon for taking numerical derivatives.
      REAL(ReKi)                   :: SQRTNStepWave2                                  ! SQRT( NStepWave/2 )
      REAL(ReKi)                   :: SWaveDir                                        ! SIN( WaveDir )
  -   REAL(ReKi), ALLOCATABLE      :: WaveAcc0H (:,:)                                 ! Instantaneous horizontal acceleration of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s^2)
  +   REAL(ReKi), ALLOCATABLE      :: WaveAcc0H (:,:)                                 ! Instantaneous horizontal acceleration of incident waves before applying stretching at the zi-coordinates for points (m/s^2)
      REAL(ReKi)                   :: WaveAcc0HExtrap                                 ! Temporary value extrapolated from the WaveAcc0H  (:,:) array (m/s^2)
      REAL(ReKi)                   :: WaveAcc0HInterp                                 ! Temporary value interpolated from the WaveAcc0H  (:,:) array (m/s^2)
  -   REAL(ReKi), ALLOCATABLE      :: WaveAcc0V (:,:)                                 ! Instantaneous vertical   acceleration of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s^2)
  +   REAL(ReKi), ALLOCATABLE      :: WaveAcc0V (:,:)                                 ! Instantaneous vertical   acceleration of incident waves before applying stretching at the zi-coordinates for points (m/s^2)
      REAL(ReKi)                   :: WaveAcc0VExtrap                                 ! Temporary value extrapolated from the WaveAcc0V  (:,:) array (m/s^2)
      REAL(ReKi)                   :: WaveAcc0VInterp                                 ! Temporary value interpolated from the WaveAcc0V  (:,:) array (m/s^2)
  -   REAL(ReKi), ALLOCATABLE      :: WaveDynP0B(:,:)                                 ! Instantaneous dynamic pressure        of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (N/m^2)
  +   REAL(ReKi), ALLOCATABLE      :: WaveDynP0B(:,:)                                 ! Instantaneous dynamic pressure        of incident waves before applying stretching at the zi-coordinates for points (N/m^2)
      REAL(ReKi)                   :: WaveDynP0BExtrap                                ! Temporary value extrapolated from the WaveDynP0B (:,:) array (N/m^2)
      REAL(ReKi)                   :: WaveDynP0BInterp                                ! Temporary value interpolated from the WaveDynP0B (:,:) array (N/m^2)
      REAL(ReKi)                   :: WaveElev_Max                                    ! Maximum expected value of the instantaneous elevation of incident waves (meters)
      REAL(ReKi)                   :: WaveElev_Min                                    ! Minimum expected value of the instantaneous elevation of incident waves (meters)
      REAL(ReKi), ALLOCATABLE      :: WaveElevxiPrime(:)                              ! Locations along the wave heading direction for points where the incident wave elevations can be output (meters)
  -   REAL(ReKi), ALLOCATABLE      :: WaveKinzi0Prime(:)                              ! zi-coordinates for points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed before applying stretching; these are relative to the mean see level (meters)
  +   REAL(ReKi)                   :: WaveElevxiPrime0
  +   REAL(ReKi), ALLOCATABLE      :: WaveKinzi0Prime(:)                              ! zi-coordinates for points where the incident wave kinematics will be computed before applying stretching; these are relative to the mean see level (meters)
  +   INTEGER   , ALLOCATABLE      :: WaveKinPrimeMap(:)
      REAL(ReKi), ALLOCATABLE      :: WaveKinzi0St   (:)                              ! Array of elevations found by stretching the elevations in the WaveKinzi0Prime(:) array using the instantaneous wave elevation; these are relative to the mean see level (meters)
      REAL(ReKi)                   :: WaveNmbr                                        ! Wavenumber of the current frequency component (1/meter)
      REAL(ReKi)                   :: WaveS1Sdd                                       ! One-sided power spectral density of the wave spectrum per unit time for the current frequency component (m^2/(rad/s))
      REAL(ReKi)                   :: WaveS2Sdd                                       ! Two-sided power spectral density of the wave spectrum per unit time for the current frequency component (m^2/(rad/s))
      REAL(DbKi)                   :: WaveTMax                                        ! Analysis time for incident wave calculations (sec)
  -   REAL(ReKi), ALLOCATABLE      :: WaveVel0H (:,:)                                 ! Instantaneous horizontal   velocity   of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s  )
  -   REAL(ReKi), ALLOCATABLE      :: WaveVel0Hxi    (:,:)                            ! Instantaneous xi-direction velocity   of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s  )
  +   REAL(ReKi), ALLOCATABLE      :: WaveVel0H (:,:)                                 ! Instantaneous horizontal   velocity   of incident waves before applying stretching at the zi-coordinates for points (m/s  )
  +   REAL(ReKi), ALLOCATABLE      :: WaveVel0Hxi    (:,:)                            ! Instantaneous xi-direction velocity   of incident waves before applying stretching at the zi-coordinates for points (m/s  )
      REAL(ReKi)                   :: WaveVel0HxiExtrap                               ! Temporary value extrapolated from the WaveVel0Hxi(:,:) array (m/s  )
      REAL(ReKi)                   :: WaveVel0HxiInterp                               ! Temporary value interpolated from the WaveVel0Hxi(:,:) array (m/s  )
  -   REAL(ReKi), ALLOCATABLE      :: WaveVel0Hyi    (:,:)                            ! Instantaneous yi-direction velocity   of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s  )
  +   REAL(ReKi), ALLOCATABLE      :: WaveVel0Hyi    (:,:)                            ! Instantaneous yi-direction velocity   of incident waves before applying stretching at the zi-coordinates for points (m/s  )
      REAL(ReKi)                   :: WaveVel0HyiExtrap                               ! Temporary value extrapolated from the WaveVel0Hyi(:,:) array (m/s  )
      REAL(ReKi)                   :: WaveVel0HyiInterp                               ! Temporary value interpolated from the WaveVel0Hyi(:,:) array (m/s  )
  -   REAL(ReKi), ALLOCATABLE      :: WaveVel0V (:,:)                                 ! Instantaneous vertical     velocity   of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s  )
  +   REAL(ReKi), ALLOCATABLE      :: WaveVel0V (:,:)                                 ! Instantaneous vertical     velocity   of incident waves before applying stretching at the zi-coordinates for points (m/s  )
      REAL(ReKi)                   :: WaveVel0VExtrap                                 ! Temporary value extrapolated from the WaveVel0V  (:,:) array (m/s  )
      REAL(ReKi)                   :: WaveVel0VInterp                                 ! Temporary value interpolated from the WaveVel0V  (:,:) array (m/s  )
      REAL(ReKi)                   :: zi_Max                                          ! Maximum elevation where the wave kinematics are to be applied using      stretching to the instantaneous free surface (meters)
  
positive.code:    REAL(ReKi)                   :: ziPrime_Max                                     ! Maximum elevation where the wave kinematics are computed before applying stretching to the instantaneous free surface (meters)
     REAL(ReKi)                   :: ziPrime_Min                                     ! Minimum elevation where the wave kinematics are computed before applying stretching to the instantaneous free surface (meters)
  
     INTEGER                      :: GHNStepWave                                     ! Total number of time steps in the GH Bladed wave data files.
     INTEGER                      :: GHNWvDpth                                       ! Number of vertical locations in GH Bladed wave data files.
     INTEGER                      :: I                                               ! Generic index
  
negative.code: !**********************************************************************************************************************************
  ! The Waves and Waves_Types modules make up a template for creating user-defined calculations in the FAST Modularization 
  ! Framework. Wavess_Types will be auto-generated based on a description of the variables for the module.
  !
  ! "Waves" should be replaced with the name of your module. Example: HydroDyn
  ! "Waves" (in Waves_*) should be replaced with the module name or an abbreviation of it. Example: HD
  !..................................................................................................................................
  ! LICENSING
  ! Copyright (C) 2013  National Renewable Energy Laboratory
  !
  !    This file is part of Waves.
  !
  ! Licensed under the Apache License, Version 2.0 (the "License");
  ! you may not use this file except in compliance with the License.
  ! You may obtain a copy of the License at
  !
  !     http://www.apache.org/licenses/LICENSE-2.0
  !
  ! Unless required by applicable law or agreed to in writing, software
  ! distributed under the License is distributed on an "AS IS" BASIS,
  ! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  ! See the License for the specific language governing permissions and
  ! limitations under the License.
  !    
  !**********************************************************************************************************************************
  ! File last committed: $Date: 2013-09-25 10:41:32 -0600 (Wed, 25 Sep 2013) $
  ! (File) Revision #: $Rev: 209 $
  ! URL: $HeadURL: https://windsvn.nrel.gov/HydroDyn/branches/HydroDyn_Modularization/Source/Waves.f90 $
  !**********************************************************************************************************************************
  MODULE Waves
  
     USE Waves_Types   
     USE NWTC_Library
     USE FFT_Module
        
     IMPLICIT NONE
     
     PRIVATE
  
  !   INTEGER(IntKi), PARAMETER            :: DataFormatID = 1   ! Update this value if the data types change (used in Waves_Pack)
     TYPE(ProgDesc), PARAMETER            :: Waves_ProgDesc = ProgDesc( 'Waves', '(v1.00.01, 19-October-2012)', '05-Mar-2013' )
  
     
        ! ..... Public Subroutines ...................................................................................................
     PUBLIC :: WavePkShpDefault                     ! Return the default value of the peak shape parameter of the incident wave spectrum
     PUBLIC :: Waves_Init                           ! Initialization routine
     PUBLIC :: Waves_End                            ! Ending routine (includes clean up)
     
     PUBLIC :: Waves_UpdateStates                   ! Loose coupling routine for solving for constraint states, integrating 
                                                      !   continuous states, and updating discrete states
     PUBLIC :: Waves_CalcOutput                     ! Routine for computing outputs
     
     PUBLIC :: Waves_CalcConstrStateResidual        ! Tight coupling routine for returning the constraint state residual
     PUBLIC :: Waves_CalcContStateDeriv             ! Tight coupling routine for computing derivatives of continuous states
     PUBLIC :: Waves_UpdateDiscState                ! Tight coupling routine for updating discrete states
        
     !PUBLIC :: Waves_JacobianPInput                 ! Routine to compute the Jacobians of the output (Y), continuous- (X), discrete-
     !                                                 !   (Xd), and constraint-state (Z) equations all with respect to the inputs (u)
     !PUBLIC :: Waves_JacobianPContState             ! Routine to compute the Jacobians of the output (Y), continuous- (X), discrete-
     !                                                 !   (Xd), and constraint-state (Z) equations all with respect to the continuous 
     !                                                 !   states (x)
     !PUBLIC :: Waves_JacobianPDiscState             ! Routine to compute the Jacobians of the output (Y), continuous- (X), discrete-
     !                                                 !   (Xd), and constraint-state (Z) equations all with respect to the discrete 
     !                                                 !   states (xd)
     !PUBLIC :: Waves_JacobianPConstrState           ! Routine to compute the Jacobians of the output (Y), continuous- (X), discrete-
     !                                                 !   (Xd), and constraint-state (Z) equations all with respect to the constraint 
     !                                                 !   states (z)
     
     
     PRIVATE:: WheelerStretching                    ! This FUNCTION applies the principle of Wheeler stretching to (1-Forward) find the elevation where the wave kinematics are to be applied using Wheeler stretching or (2-Backword)   
     PRIVATE:: BoxMuller
     PRIVATE:: JONSWAP
     PRIVATE:: WaveNumber
     PRIVATE:: UserWaveSpctrm   
     PRIVATE:: StillWaterWaves_Init
     PRIVATE:: VariousWaves_Init
     PRIVATE:: GH_BladedWaves_Init
    ! PRIVATE:: WhiteNoiseWaves_Init
  
  CONTAINS
  
     FUNCTION WavePkShpDefault ( Hs, Tp )
  
  
        ! This FUNCTION is used to return the default value of the peak shape
        ! parameter of the incident wave spectrum, conditioned on significant
        ! wave height and peak spectral period.
        !
        ! There are several different versions of the JONSWAP spectrum
        ! formula.  This version is based on the one documented in the
        ! IEC61400-3 wind turbine design standard for offshore wind turbines.
  
  
  
     IMPLICIT                        NONE
  
  
        ! Passed Variables:
  
     REAL(ReKi), INTENT(IN )      :: Hs                                              ! Significant wave height (meters)
     REAL(ReKi), INTENT(IN )      :: Tp                                              ! Peak spectral period (sec)
     REAL(ReKi)                   :: WavePkShpDefault                                ! This function = default value of the peak shape parameter of the incident wave spectrum conditioned on significant wave height and peak spectral period (-)
  
  
        ! Local Variables:
  
     REAL(ReKi)                   :: TpOvrSqrtHs                                     ! = Tp/SQRT(Hs) (s/SQRT(m))
  
  
  
        ! Compute the default peak shape parameter of the incident wave spectrum,
        !   conditioned on significant wave height and peak spectral period:
  
     TpOvrSqrtHs = Tp/SQRT(Hs)
  
     IF (     TpOvrSqrtHs <= 3.6 )  THEN
        WavePkShpDefault = 5.0
     ELSEIF ( TpOvrSqrtHs >= 5.0 )  THEN
        WavePkShpDefault = 1.0
     ELSE
        WavePkShpDefault = EXP( 5.75 - 1.15*TpOvrSqrtHs )
     END IF
  
  
  
     RETURN
     END FUNCTION WavePkShpDefault
     
  !=======================================================================
        FUNCTION BoxMuller ( NDAmp, Phase )
  
  
           ! This FUNCTION uses the Box-Muller method to turn two uniformly
           ! distributed randoms into two unit normal randoms, which are
           ! returned as real and imaginary components.
  
  
  
        IMPLICIT                             NONE
  
  
           ! Passed Variables:
  
        COMPLEX(ReKi)                     :: BoxMuller                                  ! This function
  
        REAL(ReKi), INTENT(IN ), OPTIONAL :: Phase                                      ! Optional phase to override random phase (radians)
  
        LOGICAL,    INTENT(IN )           :: NDAmp                                      ! Flag for normally-distributed amplitudes
  
  
           ! Local Variables:
  
        REAL(ReKi)                   :: C1                                              ! Intermediate variable
        REAL(ReKi)                   :: C2                                              ! Intermediate variable
        REAL(ReKi)                   :: U1                                              ! First  uniformly distributed random
        REAL(ReKi)                   :: U2                                              ! Second uniformly distributed random
  
  
  
           ! Compute the two uniformly distributed randoms:
           ! NOTE: The first random, U1, cannot be zero else the LOG() function
           !       below will blow up; there is no restriction on the value of the
           !       second random, U2.
  
        U1 = 0.0
        DO WHILE ( U1 == 0.0 )
           CALL RANDOM_NUMBER(U1)
        END DO
        CALL    RANDOM_NUMBER(U2)
  
  
           ! Compute intermediate variables:
  
        IF ( NDAmp )  THEN            ! Normally-distributed amplitudes
           C1 = SQRT( -2.0*LOG(U1) )
        ELSE                          ! Constant amplitudes (ignore U1); therefore, C1 = SQRT( 2.0 ) = MEAN( SQRT( -2.0*LOG(U1) ) for a uniform distribution of U1 between 0 and 1
           C1 = SQRT(  2.0         )
        END IF
        
        IF ( PRESENT( Phase ) )  THEN ! Specified phase to replace random phase (ignore U2)
           C2 = Phase
        ELSE                          ! Uniformly-distributed phase
           C2 = TwoPi*U2
        END IF
  
  
           ! Compute the unit normal randoms:
  
        BoxMuller = CMPLX( C1*COS(C2), C1*SIN(C2) )
  
  
  
        RETURN
        END FUNCTION BoxMuller
  !=======================================================================
        FUNCTION JONSWAP ( Omega, Hs, Tp, Gamma )
  
  
           ! This FUNCTION computes the JOint North Sea WAve Project
           ! (JONSWAP) representation of the one-sided power spectral density
           ! or wave spectrum given the frequency, Omega, peak shape
           ! parameter, Gamma, significant wave height, Hs, and peak spectral
           ! period, Tp, as inputs.  If the value of Gamma is 1.0, the
           ! Pierson-Moskowitz wave spectrum is returned.
           !
           ! There are several different versions of the JONSWAP spectrum
           ! formula.  This version is based on the one documented in the
           ! IEC61400-3 wind turbine design standard for offshore wind
           ! turbines.
  
  
  
  
        IMPLICIT                        NONE
  
  
           ! Passed Variables:
  
        REAL(ReKi), INTENT(IN )      :: Gamma                                           ! Peak shape parameter (-)
        REAL(ReKi), INTENT(IN )      :: Hs                                              ! Significant wave height (meters)
        REAL(ReKi)                   :: JONSWAP                                         ! This function = JONSWAP wave spectrum, S (m^2/(rad/s))
        REAL(ReKi), INTENT(IN )      :: Omega                                           ! Wave frequency (rad/s)
        REAL(ReKi), INTENT(IN )      :: Tp                                              ! Peak spectral period (sec)
  
  
           ! Local Variables:
  
        REAL(ReKi)                   :: Alpha                                           ! Exponent on Gamma used in the spectral formulation (-)
        REAL(ReKi)                   :: C                                               ! Normalising factor used in the spectral formulation (-)
        REAL(ReKi)                   :: f                                               ! Wave frequency (Hz)
        REAL(ReKi)                   :: fp                                              ! Peak spectral frequency (Hz)
        REAL(ReKi)                   :: fpOvrf4                                         ! (fp/f)^4
        REAL(ReKi)                   :: Sigma                                           ! Scaling factor used in the spectral formulation (-)
  
         REAL(ReKi)                  :: Inv2Pi   =  0.15915494
  
           ! Compute the JONSWAP wave spectrum, unless Omega is zero, in which case,
           !   return zero:
  
        IF ( Omega == 0.0 )  THEN  ! When .TRUE., the formulation below is ill-conditioned; thus, the known value of zero is returned.
  
  
           JONSWAP  = 0.0
  
  
        ELSE                       ! Omega > 0.0; forumulate the JONSWAP spectrum.
  
  
           ! Compute the wave frequency and peak spectral frequency in Hz:
  
           f        = Inv2Pi*Omega
           fp       = 1/Tp
           fpOvrf4  = (fp/f)**4.0
  
  
           ! Compute the normalising factor:
  
           C        = 1.0 - ( 0.287*LOG(GAMMA) )
  
  
           ! Compute Alpha:
  
           IF ( f <= fp )  THEN
              Sigma = 0.07
           ELSE
              Sigma = 0.09
           END IF
  
           Alpha    = EXP( ( -0.5*( ( (f/fp) - 1.0 )/Sigma )**2.0 ) )
  
  
           ! Compute the wave spectrum:
  
           JONSWAP  = Inv2Pi*C*( 0.3125*Hs*Hs*fpOvrf4/f )*EXP( ( -1.25*fpOvrf4 ) )*( GAMMA**Alpha )
  
  
        END IF
  
  
  
        RETURN
        END FUNCTION JONSWAP
        !=======================================================================
  !JASON: MOVE THIS USER-DEFINED ROUTINE (UserWaveSpctrm) TO THE UserSubs.f90 OF HydroDyn WHEN THE PLATFORM LOADING FUNCTIONALITY HAS BEEN DOCUMENTED!!!!!
        SUBROUTINE UserWaveSpctrm ( Omega, WaveDir, DirRoot, WaveS1Sdd )
  
  
           ! This is a dummy routine for holding the place of a user-specified
           ! wave spectrum.  Modify this code to create your own spectrum.
  
  
  
        IMPLICIT                        NONE
  
  
           ! Passed Variables:
  
        REAL(ReKi), INTENT(IN )      :: Omega                                           ! Wave frequency, rad/s.
        REAL(ReKi), INTENT(IN )      :: WaveDir                                         ! Incident wave propagation heading direction, degrees
        REAL(ReKi), INTENT(OUT)      :: WaveS1Sdd                                       ! One-sided power spectral density of the wave spectrum per unit time for the current frequency component and heading direction, m^2/(rad/s).
  
        CHARACTER(1024), INTENT(IN ) :: DirRoot                                         ! The name of the root file including the full path to the current working directory.  This may be useful if you want this routine to write a permanent record of what it does to be stored with the simulation results: the results should be stored in a file whose name (including path) is generated by appending any suitable extension to DirRoot.
  
  
  
        WaveS1Sdd = 0.0
  
  
  
        RETURN
        END SUBROUTINE UserWaveSpctrm
        !=======================================================================
        FUNCTION WaveNumber ( Omega, g, h )
  
  
           ! This FUNCTION solves the finite depth dispersion relationship:
           !
           !                   k*tanh(k*h)=(Omega^2)/g
           !
           ! for k, the wavenumber (WaveNumber) given the frequency, Omega,
           ! gravitational constant, g, and water depth, h, as inputs.  A
           ! high order initial guess is used in conjunction with a quadratic
           ! Newton's method for the solution with seven significant digits
           ! accuracy using only one iteration pass.  The method is due to
           ! Professor J.N. Newman of M.I.T. as found in routine EIGVAL of
           ! the SWIM-MOTION-LINES (SML) software package in source file
           ! Solve.f of the SWIM module.
  
  
  
        IMPLICIT                        NONE
  
  
           ! Passed Variables:
  
        REAL(ReKi), INTENT(IN )      :: g                                               ! Gravitational acceleration (m/s^2)
        REAL(ReKi), INTENT(IN )      :: h                                               ! Water depth (meters)
        REAL(ReKi), INTENT(IN )      :: Omega                                           ! Wave frequency (rad/s)
        REAL(ReKi)                   :: WaveNumber                                      ! This function = wavenumber, k (1/m)
  
  
           ! Local Variables:
  
        REAL(ReKi)                   :: A                                               ! A temporary variable used in the solution.
        REAL(ReKi)                   :: B                                               ! A temporary variable used in the solution.
        REAL(ReKi)                   :: C                                               ! A temporary variable used in the solution.
        REAL(ReKi)                   :: C2                                              ! A temporary variable used in the solution.
        REAL(ReKi)                   :: CC                                              ! A temporary variable used in the solution.
        REAL(ReKi)                   :: E2                                              ! A temporary variable used in the solution.
        REAL(ReKi)                   :: X0                                              ! A temporary variable used in the solution.
  
  
  
           ! Compute the wavenumber, unless Omega is zero, in which case, return
           !   zero:
  
        IF ( Omega == 0.0 )  THEN  ! When .TRUE., the formulation below is ill-conditioned; thus, the known value of zero is returned.
  
  
           WaveNumber = 0.0
  
  
        ELSE                       ! Omega > 0.0; solve for the wavenumber as usual.
  
  
           C  = Omega*Omega*h/g
           CC = C*C
  
  
           ! Find X0:
  
           IF ( C <= 2.0 )  THEN
  
              X0 = SQRT(C)*( 1.0 + C*( 0.169 + (0.031*C) ) )
  
           ELSE
  
              E2 = EXP(-2.0*C)
  
              X0 = C*( 1.0 + ( E2*( 2.0 - (12.0*E2) ) ) )
  
           END IF
  
  
           ! Find the WaveNumber:
  
           IF ( C <= 4.8 )  THEN
  
              C2 = CC - X0*X0
              A  = 1.0/( C - C2 )
              B  = A*( ( 0.5*LOG( ( X0 + C )/( X0 - C ) ) ) - X0 )
  
              WaveNumber = ( X0 - ( B*C2*( 1.0 + (A*B*C*X0) ) ) )/h
  
           ELSE
  
              WaveNumber = X0/h
  
           END IF
  
  
        END IF
  
  
  
        RETURN
        END FUNCTION WaveNumber
        
        !=======================================================================
        FUNCTION COSHNumOvrCOSHDen ( k, h, z )
  
        
           ! This FUNCTION computes the shallow water hyperbolic numerator
           ! over denominator term in the wave kinematics expressions:
           !
           !                    COSH( k*( z + h ) )/COSH( k*h )
           !
           ! given the wave number, k, water depth, h, and elevation z, as
           ! inputs.
  
        IMPLICIT                        NONE
  
  
           ! Passed Variables:
  
        REAL(ReKi)                   :: COSHNumOvrCOSHDen                               ! This function = COSH( k*( z + h ) )/COSH( k*h ) (-)
        REAL(ReKi), INTENT(IN )      :: h                                               ! Water depth ( h      >  0 ) (meters)
        REAL(ReKi), INTENT(IN )      :: k                                               ! Wave number ( k      >= 0 ) (1/m)
        REAL(ReKi), INTENT(IN )      :: z                                               ! Elevation   (-h <= z <= 0 ) (meters)
  
  
  
           ! Compute the hyperbolic numerator over denominator:
  
        IF ( k*h  > 89.4_ReKi )  THEN   ! When .TRUE., the shallow water formulation will trigger a floating point overflow error; however, COSH( k*( z + h ) )/COSH( k*h ) = EXP( k*z ) + EXP( -k*( z + 2*h ) ) for large k*h.  This equals the deep water formulation, EXP( k*z ), except near z = -h, because h > 14.23*wavelength (since k = 2*Pi/wavelength) in this case.
  
           COSHNumOvrCOSHDen = EXP( k*z ) + EXP( -k*( z + 2.0_ReKi*h ) )
  
        ELSE                       ! 0 < k*h <= 89.4; use the shallow water formulation.
  
           COSHNumOvrCOSHDen = COSH( k*( z + h ) )/COSH( k*h )
  
        END IF
  
  
  
        RETURN
        END FUNCTION COSHNumOvrCOSHDen
  !=======================================================================
        FUNCTION COSHNumOvrSINHDen ( k, h, z )
  
  
           ! This FUNCTION computes the shallow water hyperbolic numerator
           ! over denominator term in the wave kinematics expressions:
           !
           !                    COSH( k*( z + h ) )/SINH( k*h )
           !
           ! given the wave number, k, water depth, h, and elevation z, as
           ! inputs.
  
  
  
        IMPLICIT                        NONE
  
  
           ! Passed Variables:
  
        REAL(ReKi)                   :: COSHNumOvrSINHDen                               ! This function = COSH( k*( z + h ) )/SINH( k*h ) (-)
        REAL(ReKi), INTENT(IN )      :: h                                               ! Water depth ( h      >  0 ) (meters)
        REAL(ReKi), INTENT(IN )      :: k                                               ! Wave number ( k      >= 0 ) (1/m)
        REAL(ReKi), INTENT(IN )      :: z                                               ! Elevation   (-h <= z <= 0 ) (meters)
  
  
  
           ! Compute the hyperbolic numerator over denominator:
  
  !bjj: should this be compared with epsilon instead of 0.0 to avoid numerical instability?
        IF (     k   == 0.0  )  THEN  ! When .TRUE., the shallow water formulation is ill-conditioned; thus, HUGE(k) is returned to approximate the known value of infinity.
  
           COSHNumOvrSINHDen = HUGE( k )
  
        ELSEIF ( k*h  > 89.4_ReKi )  THEN  ! When .TRUE., the shallow water formulation will trigger a floating point overflow error; however, COSH( k*( z + h ) )/SINH( k*h ) = EXP( k*z ) + EXP( -k*( z + 2*h ) ) for large k*h.  This equals the deep water formulation, EXP( k*z ), except near z = -h, because h > 14.23*wavelength (since k = 2*Pi/wavelength) in this case.
  
           COSHNumOvrSINHDen = EXP( k*z ) + EXP( -k*( z + 2*h ) )
  
        ELSE                          ! 0 < k*h <= 89.4; use the shallow water formulation.
  
           COSHNumOvrSINHDen = COSH( k*( z + h ) )/SINH( k*h )
  
        END IF
  
  
  
        RETURN
        END FUNCTION COSHNumOvrSINHDen
  !=======================================================================
        FUNCTION COTH ( X )
  
  
           ! This FUNCTION computes the hyperbolic cotangent,
           ! COSH(X)/SINH(X).
  
  
        USE                             Precision
  
  
        IMPLICIT                        NONE
  
  
           ! Passed Variables:
  
        REAL(ReKi)                   :: COTH                                            ! This function = COSH( X )/SINH( X ) (-)
        REAL(ReKi), INTENT(IN )      :: X                                               ! The argument (-)
  
  
  
           ! Compute the hyperbolic cotangent:
  
        IF ( X == 0.0 )  THEN   ! When .TRUE., the formulation below is ill-conditioned; thus, HUGE(X) is returned to approximate the known value of infinity.
  
           COTH = HUGE( X )
  
        ELSE                    ! X /= 0.0; use the numerically-stable computation of COTH(X) by means of TANH(X).
  
           COTH = 1.0/TANH( X ) ! = COSH( X )/SINH( X )
  
        END IF
  
  
  
        RETURN
        END FUNCTION COTH
      
        !=======================================================================
        FUNCTION SINHNumOvrSINHDen ( k, h, z )
  
  
           ! This FUNCTION computes the shallow water hyperbolic numerator
           ! over denominator term in the wave kinematics expressions:
           !
           !                    SINH( k*( z + h ) )/SINH( k*h )
           !
           ! given the wave number, k, water depth, h, and elevation z, as
           ! inputs.
  
  
        IMPLICIT                        NONE
  
  
           ! Passed Variables:
  
        REAL(ReKi)                   :: SINHNumOvrSINHDen                               ! This function = SINH( k*( z + h ) )/SINH( k*h ) (-)
        REAL(ReKi), INTENT(IN )      :: h                                               ! Water depth ( h      >  0 ) (meters)
        REAL(ReKi), INTENT(IN )      :: k                                               ! Wave number ( k      >= 0 ) (1/m)
        REAL(ReKi), INTENT(IN )      :: z                                               ! Elevation   (-h <= z <= 0 ) (meters)
  
  
  
           ! Compute the hyperbolic numerator over denominator:
  
        IF (     k   == 0.0  )  THEN  ! When .TRUE., the shallow water formulation is ill-conditioned; thus, the known value of unity is returned.
  
           SINHNumOvrSINHDen = 1.0
  
        ELSEIF ( k*h >  89.4_ReKi )  THEN  ! When .TRUE., the shallow water formulation will trigger a floating point overflow error; however, SINH( k*( z + h ) )/SINH( k*h ) = EXP( k*z ) - EXP( -k*( z + 2*h ) ) for large k*h.  This equals the deep water formulation, EXP( k*z ), except near z = -h, because h > 14.23*wavelength (since k = 2*Pi/wavelength) in this case.
  
           SINHNumOvrSINHDen = EXP( k*z ) - EXP( -k*( z + 2.0_ReKi*h ) )
  
        ELSE                          ! 0 < k*h <= 89.4; use the shallow water formulation.
  
           SINHNumOvrSINHDen = SINH( k*( z + h ) )/SINH( k*h )
  
        END IF
  
  
  
        RETURN
        END FUNCTION SINHNumOvrSINHDen
        
        
        
  !----------------------------------------------------------------------------------------------------------------------------------
  !SUBROUTINE WhiteNoiseWaves_Init ( InitInp, InitOut, ErrStat, ErrMsg )    
  !! This routine initializes the waves data for WaveMod = 5 , or white noise waves option
  !!----------------------------------------------------------------------------------------------------------------------------------
  ! 
  !
  !   TYPE(Waves_InitInputType),       INTENT(IN   )  :: InitInp     ! Input data for initialization routine    
  !   TYPE(Waves_InitOutputType),      INTENT(INOUT)  :: InitOut     ! Initialization output data      
  !   INTEGER(IntKi),                  INTENT(  OUT)  :: ErrStat     ! Error status of the operation
  !   CHARACTER(*),                    INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
  !
  !   
  !      ! Local Variables
  !      WaveS1Sdd = WaveDT*WaveHs*WaveHs/(8.0*Pi)
  !END SUBROUTINE WhiteNoiseWaves_Init
  
  !----------------------------------------------------------------------------------------------------------------------------------
  SUBROUTINE StillWaterWaves_Init ( InitInp, InitOut, ErrStat, ErrMsg )    
  ! This routine initializes the waves data for WaveMod = 0 , or still water waves option
  !----------------------------------------------------------------------------------------------------------------------------------
   
  
     TYPE(Waves_InitInputType),       INTENT(IN   )  :: InitInp     ! Input data for initialization routine    
     TYPE(Waves_InitOutputType),      INTENT(INOUT)  :: InitOut     ! Initialization output data      
     INTEGER(IntKi),                  INTENT(  OUT)  :: ErrStat     ! Error status of the operation
     CHARACTER(*),                    INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
  
     
        ! Local Variables
     INTEGER                      :: J                                               ! Generic index
     
        ! Initialize ErrStat
           
     ErrStat = ErrID_None         
     ErrMsg  = "" 
     
     
     ! Initialize everything to zero:
  
        InitOut%NStepWave  = 2                ! We must have at least two elements in order to interpolate later on
        InitOut%NStepWave2 = 1
  
        ALLOCATE ( InitOut%WaveTime      (0:InitOut%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= ErrID_None )  THEN
           ErrMsg =' Error allocating memory for the WaveTime array.'
           ErrStat = ErrID_Fatal
           RETURN
        END IF
  
        ALLOCATE ( InitOut%WaveElevC0    (2, 0:InitOut%NStepWave2                  ) , STAT=ErrStat )
        IF ( ErrStat /= ErrID_None )  THEN
           ErrMsg =' Error allocating memory for the WaveElevC0 array.'
           ErrStat = ErrID_Fatal
           RETURN
        END IF   
  
        ALLOCATE ( InitOut%WaveElev      (0:InitOut%NStepWave-1,InitInp%NWaveElev  ) , STAT=ErrStat )
        IF ( ErrStat /= ErrID_None )  THEN
           ErrMsg =' Error allocating memory for the WaveElev array.'
           ErrStat = ErrID_Fatal
           RETURN
        END IF
  
        ALLOCATE ( InitOut%WaveDynP0     (0:InitOut%NStepWave-1,InitInp%NWaveKin0  ) , STAT=ErrStat )
        IF ( ErrStat /= ErrID_None )  THEN
           ErrMsg =' Error allocating memory for the WaveDynP0 array.'
           ErrStat = ErrID_Fatal
           RETURN
        END IF
  
  
        ALLOCATE ( InitOut%WaveVel0      (0:InitOut%NStepWave-1,InitInp%NWaveKin0,3) , STAT=ErrStat )
        IF ( ErrStat /= ErrID_None )  THEN
           ErrMsg =' Error allocating memory for the WaveVel0 array.'
           ErrStat = ErrID_Fatal
           RETURN
        END IF
  
        ALLOCATE ( InitOut%WaveAcc0      (0:InitOut%NStepWave-1,InitInp%NWaveKin0,3) , STAT=ErrStat )
        IF ( ErrStat /= ErrID_None )  THEN
           ErrMsg =' Error allocating memory for the WaveAcc0 array.'
           ErrStat = ErrID_Fatal
           RETURN
        END IF
  
        InitOut%WaveDOmega = 0.0
        InitOut%WaveTime   = (/ 0.0_DbKi, 1.0_DbKi /)   ! We must have at least two different time steps in the interpolation
        InitOut%WaveElevC0(1,0) = 0.0
        InitOut%WaveElevC0(2,0) = 0.0
        InitOut%WaveElev   = 0.0
        InitOut%WaveDynP0  = 0.0
        InitOut%WaveVel0   = 0.0
        InitOut%WaveAcc0   = 0.0
        
        ! Add the current velocities to the wave velocities:
  
        DO J = 1,InitInp%NWaveKin0   ! Loop through all Morison element nodes where the incident wave kinematics will be computed
           
           InitOut%WaveVel0(:,J,1) =  InitInp%CurrVxi(J)  ! xi-direction
           InitOut%WaveVel0(:,J,2) =  InitInp%CurrVyi(J)  ! yi-direction
  
        END DO                ! J - All points where the incident wave kinematics will be computed
  
  END SUBROUTINE StillWaterWaves_Init
  
  
  !----------------------------------------------------------------------------------------------------------------------------------
  SUBROUTINE VariousWaves_Init ( InitInp, InitOut, ErrStat, ErrMsg )
  ! Compute the wave kinematics and related information for  Plane progressive (regular) wave, JONSWAP/Pierson-Moskowitz spectrum  
  ! (irregular) wave, or user-defined spectrum (irregular) wave.
  !----------------------------------------------------------------------------------------------------------------------------------
  
     TYPE(Waves_InitInputType),       INTENT(IN   )  :: InitInp     ! Input data for initialization routine
     TYPE(Waves_InitOutputType),       INTENT(INOUT)  :: InitOut     ! Output data      
     INTEGER(IntKi),                  INTENT(  OUT)  :: ErrStat     ! Error status of the operation
     CHARACTER(*),                    INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
  
     
     
     ! Local Variables
     COMPLEX(ReKi), PARAMETER     :: ImagNmbr = (0.0,1.0)                         ! The imaginary number, SQRT(-1.0)
     COMPLEX(ReKi)                :: ImagOmega                                       ! = ImagNmbr*Omega (rad/s)
     REAL(ReKi), ALLOCATABLE      :: WaveElev0  (:)                            ! Instantaneous elevation of incident waves at the platform reference point (meters)
        
     COMPLEX(ReKi), ALLOCATABLE   :: PWaveAccC0HPz0 (:)                              ! Partial derivative of WaveAccC0H (:) with respect to zi at zi = 0 (1/s^2) 
     COMPLEX(ReKi), ALLOCATABLE   :: PWaveAccC0VPz0 (:)                              ! Partial derivative of WaveAccC0V (:) with respect to zi at zi = 0 (1/s^2) 
     COMPLEX(ReKi), ALLOCATABLE   :: PWaveDynPC0BPz0(:)                              ! Partial derivative of WaveDynPC0B(:) with respect to zi at zi = 0 (N/m  ) 
     COMPLEX(ReKi), ALLOCATABLE   :: PWaveVelC0HPz0 (:)                              ! Partial derivative of WaveVelC0H (:) with respect to zi at zi = 0 (1/s  ) 
     COMPLEX(ReKi), ALLOCATABLE   :: PWaveVelC0VPz0 (:)                              ! Partial derivative of WaveVelC0V (:) with respect to zi at zi = 0 (1/s  ) 
     COMPLEX(ReKi), ALLOCATABLE   :: WaveAccC0H(:,:)                                 ! Discrete Fourier transform of the instantaneous horizontal acceleration of incident waves before applying stretching at the zi-coordinates for points (m/s^2)
     COMPLEX(ReKi), ALLOCATABLE   :: WaveAccC0V(:,:)                                 ! Discrete Fourier transform of the instantaneous vertical   acceleration of incident waves before applying stretching at the zi-coordinates for points (m/s^2)
     COMPLEX(ReKi), ALLOCATABLE   :: WaveDynPC0(:,:)                                 ! Discrete Fourier transform of the instantaneous dynamic pressure        of incident waves before applying stretching at the zi-coordinates for points (N/m^2)
     COMPLEX(ReKi), ALLOCATABLE   :: WaveElevC (:,:)                                 ! Discrete Fourier transform of the instantaneous elevation of incident waves (meters)
     COMPLEX(ReKi), ALLOCATABLE   :: WaveVelC0H(:,:)                                 ! Discrete Fourier transform of the instantaneous horizontal velocity     of incident waves before applying stretching at the zi-coordinates for points (m/s)
     COMPLEX(ReKi), ALLOCATABLE   :: WaveVelC0V(:,:)                                 ! Discrete Fourier transform of the instantaneous vertical   velocity     of incident waves before applying stretching at the zi-coordinates for points (m/s)
     COMPLEX(ReKi)                :: WGNC                                            ! Discrete Fourier transform of the realization of a White Gaussian Noise (WGN) time series process with unit variance for the current frequency component (-)
     REAL(ReKi)                   :: CurrVxi                                         ! xi-component of the current velocity at the instantaneous elevation (m/s)
     REAL(ReKi)                   :: CurrVyi                                         ! yi-component of the current velocity at the instantaneous elevation (m/s)
  