anchor.diff: +!**********************************************************************************************************************************
  +! 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
   
  -
  -   ! This MODULE stores variables and routines associated with incident
  -   ! waves and current.
  -
  -   USE                             NWTC_Library
  +   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
   
   
  -      ! Parameter data
  -
  +      ! 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)
  -!   REAL(ReKi),    PARAMETER     :: Inv2Pi   =  0.15915494                       ! 0.5/Pi.
  -   REAL(ReKi),    PARAMETER     :: PiOvr4   =  0.78539816                       ! Pi/4.
  -
  -
  -   TYPE, PUBLIC :: Waves_DataType                                                ! The type is public
  -!      PRIVATE                                                                   ! but the data is private !actually, we're going to make it public for the platform/support structure modules to access this data
  -
  -      COMPLEX(ReKi), ALLOCATABLE   :: WaveElevC0 (:)                            ! Discrete Fourier transform of the instantaneous elevation of incident waves at the platform reference point (meters)
  -      
  -!bjj: not sure what DZNodes is used for in this module:
  -      REAL(ReKi), ALLOCATABLE      :: DZNodes    (:)                            ! Length of variable-length tower or platform elements for the points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (meters)
  -      REAL(ReKi)                   :: Gravity                                   ! Gravitational acceleration (m/s^2)
  -      REAL(ReKi)                   :: RhoXg       = 0.0                         ! = WtrDens*Gravity
  -
  -      REAL(ReKi), ALLOCATABLE      :: WaveAcc0   (:,:,:)                        ! Instantaneous acceleration of incident waves in the xi- (1), yi- (2), and zi- (3) directions, respectively, accounting for stretching, at each of the NWaveKin0 points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (m/s^2)
  -      REAL(ReKi)                   :: WaveDir     = 0.0                         ! Incident wave propagation heading direction (degrees)
  -      REAL(ReKi)                   :: WaveDOmega  = 0.0                         ! Frequency step for incident wave calculations (rad/s)
  -      REAL(ReKi), ALLOCATABLE      :: WaveDynP0 (:,:)                           ! Instantaneous dynamic pressure of incident waves                                                          , accounting for stretching, at each of the NWaveKin0 points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (N/m^2)
  -      REAL(ReKi), ALLOCATABLE      :: WaveElev   (:,:)                          ! Instantaneous elevation of incident waves at each of the NWaveElev points where the incident wave elevations can be output (meters)
  -      REAL(ReKi), ALLOCATABLE      :: WaveElevxi  (:)                           ! xi-coordinates for points where the incident wave elevations can be output (meters)
  -      REAL(ReKi), ALLOCATABLE      :: WaveElevyi  (:)                           ! yi-coordinates for points where the incident wave elevations can be output (meters)
  -
  -      REAL(ReKi), ALLOCATABLE      :: WaveElev0  (:)                            ! Instantaneous elevation of incident waves at the platform reference point (meters)
  -      REAL(ReKi), ALLOCATABLE      :: WaveKinzi0 (:)                            ! zi-coordinates for points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed; these are relative to the mean see level (meters)
  -      REAL(ReKi), ALLOCATABLE      :: WaveTime   (:)                            ! Simulation times at which the instantaneous elevation of, velocity of, acceleration of, and loads associated with the incident waves are determined (sec)
  -      REAL(ReKi), ALLOCATABLE      :: WaveVel0   (:,:,:)                        ! Instantaneous velocity     of incident waves in the xi- (1), yi- (2), and zi- (3) directions, respectively, accounting for stretching, at each of the NWaveKin0 points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (m/s  ) (The values include both the velocity of incident waves and the velocity of current.)
  -      REAL(ReKi)                   :: WtrDens     = 0.0                         ! Water density (kg/m^3)
  -      REAL(ReKi)                   :: WtrDpth     = 0.0                         ! Water depth (meters)
  -
  -      INTEGER                      :: LastIndAcc  = 1                           ! Index into the arrays saved from the last call to WaveAcceleration    as a starting point for next call.
  -      INTEGER                      :: LastIndElev = 1                           ! Index into the arrays saved from the last call to WaveElevation       as a starting point for next call.
  -      INTEGER                      :: LastIndDynP = 1                           ! Index into the arrays saved from the last call to WaveDynamicPressure as a starting point for next call.
  -      INTEGER                      :: LastIndVel  = 1                           ! Index into the arrays saved from the last call to WaveVelocity        as a starting point for next call.
  -      
  -      INTEGER                      :: NStepWave   = 0                           ! Total number of frequency components = total number of time steps in the incident wave (-)
  -      INTEGER                      :: NStepWave2  = 0                           ! NStepWave/2
  -      INTEGER                      :: NWaveElev   = 0                           ! Number of points where the incident wave elevation can be output (-)
  -      INTEGER                      :: NWaveKin0   = 0                           ! Number of points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (-)
  -      INTEGER                      :: WaveMod     = -1                          ! Incident wave kinematics model {0: none=still water, 1: plane progressive (regular), 2: JONSWAP/Pierson-Moskowitz spectrum (irregular), 3: user-defind spectrum from routine UserWaveSpctrm (irregular)}
  -      INTEGER                      :: WaveStMod   = -1                          ! Model for stretching incident wave kinematics to instantaneous free surface {0: none=no stretching, 1: vertical stretching, 2: extrapolation stretching, 3: Wheeler stretching}
  -
  -   END TYPE Waves_DataType
  -
  -   TYPE, PUBLIC :: Current_DataType
  -      REAL(ReKi)                   :: CurrSSV0    = 0.0                         ! Sub-surface current velocity at still water level (m/s)
  -      REAL(ReKi)                   :: CurrSSDir   = 0.0                         ! Sub-surface current heading direction (degrees)
  -      REAL(ReKi)                   :: CurrNSRef   = 0.0                         ! Near-surface current reference depth (meters)
  -      REAL(ReKi)                   :: CurrNSV0    = 0.0                         ! Near-surface current velocity at still water level (m/s)
  -      REAL(ReKi)                   :: CurrNSDir   = 0.0                         ! Near-surface current heading direction (degrees)
  -      REAL(ReKi)                   :: CurrDIV     = 0.0                         ! Depth-independent current velocity (m/s)
  -      REAL(ReKi)                   :: CurrDIDir   = 0.0                         ! Depth-independent current heading direction (degrees)
  -      INTEGER                      :: CurrMod     = 0                           ! Current profile model {0: none=no current, 1: standard, 2: user-defined from routine UserCurrent}
  -   END TYPE Current_DataType
  -
  -   TYPE, PUBLIC :: Waves_InitDataType
  -      REAL(ReKi)                   :: Gravity                                   ! Gravitational acceleration (m/s^2)
  -      REAL(ReKi)                   :: WtrDens                                   ! Water density (kg/m^3)
  -      REAL(ReKi)                   :: WtrDpth                                   ! Water depth (meters)
  -      REAL(DbKi)                   :: WaveTMax    = 0.0                         ! Analysis time for incident wave calculations; the actual analysis time may be larger than this value in order for the maintain an effecient FFT (sec)
  -      REAL(DbKi)                   :: WaveDT      = 0.0                         ! Time step for incident wave calculations (sec)
  -      REAL(ReKi)                   :: WaveHs      = 0.0                         ! Significant wave height of incident waves (meters)
  -      REAL(ReKi)                   :: WaveTp      = 0.0                         ! Peak spectral period of incident waves (sec)
  -      REAL(ReKi)                   :: WavePkShp   = 1.0                         ! Peak shape parameter of incident wave spectrum [1.0 for Pierson-Moskowitz] (-)
  -      REAL(ReKi)                   :: WaveDir     = 0.0                         ! Incident wave propagation heading direction (degrees)
  -
  -      REAL(ReKi), ALLOCATABLE      :: DZNodes     (:)                           ! Length of variable-length tower or platform elements for the points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (meters)
  -      REAL(ReKi), ALLOCATABLE      :: WaveElevxi  (:)                           ! xi-coordinates for points where the incident wave elevations can be output (meters)
  -      REAL(ReKi), ALLOCATABLE      :: WaveElevyi  (:)                           ! yi-coordinates for points where the incident wave elevations can be output (meters)
  -      REAL(ReKi), ALLOCATABLE      :: WaveKinzi0  (:)                           ! zi-coordinates for points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed; these are relative to the mean see level (meters)
  -
  -      REAL(ReKi)                   :: WavePhase   = 0.0                         ! Specified phase for regular waves (radians)
  -      
  -      LOGICAL                      :: WaveNDAmp   = .FALSE.                     ! Flag for normally-distributed amplitudes in incident waves spectrum (flag)      
  -
  -      INTEGER                      :: NWaveElev   = 0                           ! Number of points where the incident wave elevations can be output (-)
  -      INTEGER                      :: NWaveKin0   = 0                           ! Number of points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (-)
  -
  -      INTEGER                      :: WaveSeed (2)= 0                           ! Random seeds of incident waves [-2147483648 to 2147483647]
  -      INTEGER                      :: WaveMod     = 0                           ! Incident wave kinematics model {0: none=still water, 1: plane progressive (regular), 2: JONSWAP/Pierson-Moskowitz spectrum (irregular), 3: user-defind spectrum from routine UserWaveSpctrm (irregular)}
  -      INTEGER                      :: WaveStMod   = 0                           ! Model for stretching incident wave kinematics to instantaneous free surface {0: none=no stretching, 1: vertical stretching, 2: extrapolation stretching, 3: Wheeler stretching}
  -
  -      CHARACTER(1024)              :: GHWvFile    = ""                          ! The root name of GH Bladed files containing wave data.
  -      CHARACTER(1024)              :: DirRoot     = ""                          ! The name of the root file including the full path.  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.
  -   END TYPE Waves_InitDataType
  -
  -
  -
  -CONTAINS
  -!=======================================================================
  -   SUBROUTINE InitWaves ( Waves_InitData, Current_Data, WaveDat, UnSum, ErrStat  )
  -
  -
  -      ! This routine is used to initialize the variables associated with
  -      ! incident waves and current.
  -
  -
  -   USE                             FFT_Module
  -
  -   IMPLICIT                        NONE
  -
  -
  -      ! Passed Variables:
  -
  -   TYPE(Current_DataType),  INTENT(IN) :: Current_Data                             ! the values needed to initialize the current in the waves module
  -   TYPE(Waves_InitDataType),INTENT(IN) :: Waves_InitData                           ! the values needed to initialize the waves module (but don't need to be retained for further reference)
  -
  -   TYPE(Waves_DataType), INTENT(OUT)   :: WaveDat                                  ! The current instance of Waves
  -   INTEGER, INTENT(IN )                :: UnSum                                    ! The unit number for the HydroDyn summary file. This subroutine assumes the file is open.
  -   INTEGER, INTENT(OUT)                :: ErrStat                                  ! a non-zero value indicates an error was encoutered
  -
  -      ! Local Variables:
  -
      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 along a vertical line passing through the platform reference point (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 along a vertical line passing through the platform reference point (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 along a vertical line passing through the platform reference point (N/m^2)
  +   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 along a vertical line passing through the platform reference point (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 along a vertical line passing through the platform reference point (m/s)
  +   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)
  
positive.code:       IF ( nSeeds /= 2 ) THEN
           CALL ProgWarn( ' The random number generator in use differs from the original code provided by NREL. This pRNG uses ' &
                                    //TRIM(Int2LStr(nSeeds))//' seeds instead of the 2 in the HydroDyn input file.')
        END IF
  
        ALLOCATE ( TmpWaveSeeds ( nSeeds ), STAT=ErrStat )
        IF (ErrStat /= 0 ) THEN
           CALL ProgAbort( ' Error allocating space for TmpWaveSeeds array.', TrapErrors = .TRUE. )
           RETURN
        END IF   
  
           ! We'll just populate this with odd seeds = Seed(1) and even seeds = Seed(2)
        DO I = 1,nSeeds,2
           TmpWaveSeeds(I) = Waves_InitData%WaveSeed(1)
        END DO
        DO I = 2,nSeeds,2
           TmpWaveSeeds(I) = Waves_InitData%WaveSeed(2)
        END DO
                       
                    
        CALL RANDOM_SEED ( PUT=TmpWaveSeeds )
        DEALLOCATE(TmpWaveSeeds, STAT=ErrStat)
        IF (ErrStat /= 0 ) THEN
           CALL ProgWarn( ' Error deallocating space for TmpWaveSeeds array.' )
           ErrStat = 0
        END IF            
                    
        WaveDat%NStepWave  = CEILING ( WaveTMax/Waves_InitData%WaveDT )                        ! Set NStepWave to an even integer
        IF ( MOD(WaveDat%NStepWave,2) == 1 )  WaveDat%NStepWave = WaveDat%NStepWave + 1        !   larger or equal to WaveTMax/WaveDT.
        WaveDat%NStepWave2 = MAX( WaveDat%NStepWave/2, 1 )                                     ! Make sure that NStepWave is an even product of small factors (PSF) that is
        WaveDat%NStepWave  = 2*PSF ( WaveDat%NStepWave2, 9 )                                   !   greater or equal to WaveTMax/WaveDT to ensure that the FFT is efficient.
  
        WaveDat%NStepWave2 = WaveDat%NStepWave/2                                               ! Update the value of NStepWave2 based on the value needed for NStepWave.
        WaveTMax           = WaveDat%NStepWave*Waves_InitData%WaveDT                           ! Update the value of WaveTMax   based on the value needed for NStepWave.
        SQRTNStepWave2 = SQRT( REAL( WaveDat%NStepWave2, ReKi ) )                              ! Compute SQRT( NStepWave/2 ).
        WaveDat%WaveDOmega = TwoPi/WaveTMax                                                    ! Compute the frequency step for incident wave calculations.
        I_WaveTp           = NINT ( TwoPi/(WaveDat%WaveDOmega*Waves_InitData%WaveTp) )         ! Compute the index of the frequency component nearest to WaveTp.
        IF ( WaveDat%WaveMod == 2 )  OmegaCutOff = n_Massel*TwoPi/Waves_InitData%WaveTp        ! Compute the cut-off frequency or upper frequency limit of the wave spectrum beyond which the wave spectrum is zeroed.  The TwoPi/WaveTp is the peak spectral frequency in rad/s; the cut-off frequency is a factor of N_Massel above this value based on the suggestion in: Massel, S. R., Ocean Surface Waves: Their Physics and Prediction, Advanced Series on Ocean Engineering - Vol. 11, World Scientific Publishing, Singapore - New Jersey - London - Hong Kong, 1996.
  
        ALLOCATE ( WaveDat%WaveTime  (0:WaveDat%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveTime array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveDat%WaveElevC0(0:WaveDat%NStepWave2                     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveElevC0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveElevC         (0:WaveDat%NStepWave2 ,WaveDat%NWaveElev  ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveElevC array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveDynPC0        (0:WaveDat%NStepWave2 ,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveDynPC0 array.')
           RETURN
        ENDIF
  
        ALLOCATE ( WaveVelC0H        (0:WaveDat%NStepWave2 ,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveVelC0H array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveVelC0V        (0:WaveDat%NStepWave2 ,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveVelC0V array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveAccC0H        (0:WaveDat%NStepWave2 ,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveAccC0H array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveAccC0V        (0:WaveDat%NStepWave2 ,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveAccC0V array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveDynPC0BPz0   (0:WaveDat%NStepWave2                     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveDynPC0BPz0 array.')
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveVelC0HPz0    (0:WaveDat%NStepWave2                     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveVelC0HPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveVelC0VPz0    (0:WaveDat%NStepWave2                     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveVelC0VPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveAccC0HPz0    (0:WaveDat%NStepWave2                     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveAccC0HPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveAccC0VPz0    (0:WaveDat%NStepWave2                     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveAccC0VPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveDat%WaveElev0 (0:WaveDat%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveElev0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveDat%WaveElev  (0:WaveDat%NStepWave-1,WaveDat%NWaveElev  ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveElev array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
        
        ALLOCATE ( WaveDynP0B        (0:WaveDat%NStepWave-1,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveDynP0B array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveVel0H         (0:WaveDat%NStepWave-1,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveVel0H array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveVel0Hxi       (0:WaveDat%NStepWave-1,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveVel0Hxi array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveVel0Hyi       (0:WaveDat%NStepWave-1,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveVel0Hyi array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveVel0V         (0:WaveDat%NStepWave-1,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveVel0V array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveAcc0H         (0:WaveDat%NStepWave-1,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveAcc0H array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveAcc0V         (0:WaveDat%NStepWave-1,NWaveKin0Prime     ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveAcc0V array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
        
        ALLOCATE ( PWaveDynP0BPz0    (0:WaveDat%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveDynP0BPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveVel0HPz0     (0:WaveDat%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveVel0HPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveVel0HxiPz0   (0:WaveDat%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveVel0HxiPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveVel0HyiPz0   (0:WaveDat%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveVel0HyiPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveVel0VPz0     (0:WaveDat%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveVel0VPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveAcc0HPz0     (0:WaveDat%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveAcc0HPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( PWaveAcc0VPz0     (0:WaveDat%NStepWave-1                    ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the PWaveAcc0VPz0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveDat%WaveDynP0 (0:WaveDat%NStepWave-1,WaveDat%NWaveKin0  ) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveDynP0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveDat%WaveVel0  (0:WaveDat%NStepWave-1,WaveDat%NWaveKin0,3) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveVel0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
        ALLOCATE ( WaveDat%WaveAcc0  (0:WaveDat%NStepWave-1,WaveDat%NWaveKin0,3) , STAT=ErrStat )
        IF ( ErrStat /= 0 )  THEN
           CALL ProgAbort(' Error allocating memory for the WaveAcc0 array.', TrapErrors = .TRUE.)
           RETURN
        ENDIF
  
  
  
  !JASON: IMPLEMENT EQUATIONS (2.12 - 2.13) IN MY DISSERTATION SO THAT ONE CAN READ IN EXTERNAL WAVE DATA?<--BETTER YET, IMPLEMENT WaveElevC0 = DFT(WaveElev) WHERE WaveElev CAN BE READ IN AS GH BLADED WAVE DATA.  THAT IS, ADD AN OPTION TO READ IN WAVE DATA FOR FLOATERS!
  
        ! Compute the positive-frequency components (including zero) of the discrete
        !   Fourier transforms of the wave kinematics:
  
        DO I = 0,WaveDat%NStepWave2  ! Loop through the positive frequency components (including zero) of the discrete Fourier transforms
  
  
        ! Compute the frequency of this component and its imaginary value:
  
               Omega = I*       WaveDat%WaveDOmega
           ImagOmega = ImagNmbr*Omega
  
  
  
negative.code:       !   alternative elevations [indicated here by the NWaveKin0Prime-element
        !   array WaveKinzi0Prime(:)]:
  
     
  
  
  
  