anchor.diff:        TYPE(W2_SumData_Type)                              :: SumQTFData           !< Data storage for the full sum QTF method
   
            ! Force arrays
  -      REAL(SiKi)                                         :: MnDriftForce(:)      !< MnDrift force array.   Constant for all time.  First index is force component
  +      REAL(SiKi),    ALLOCATABLE                         :: MnDriftForce(:)      !< MnDrift force array.   Constant for all time.  First index is force component
         REAL(SiKi),    ALLOCATABLE                         :: NewmanAppForce(:,:)  !< NewmanApp force array.  Index 1: Time,    Index 2: force component
         REAL(SiKi),    ALLOCATABLE                         :: DiffQTFForce(:,:)    !< DiffQTF force array.    Index 1: Time,    Index 2: force component
         REAL(SiKi),    ALLOCATABLE                         :: SumQTFForce(:,:)     !< SumQTF force array.     Index 1: Time,    Index 2: force component
  
positive.code: 
  
           ! Now loop through all the dimensions and perform the calculation
        DO I=1,6*p%NBody
  
              ! Set the MnDrift force to 0.0 (Even ones we don't calculate)
           MnDriftForce(I)   = 0.0_SiKi
  
           IF (MnDriftData%DataIs3D) THEN
              TmpFlag = MnDriftData%Data3D%LoadComponents(I)
           ELSE
              TmpFlag = MnDriftData%Data4D%LoadComponents(I)
           END IF
  
              ! Only on the dimensions we requested, and if it is present in the data
           IF ( p%MnDriftDims(I) .AND. TmpFlag ) THEN
  
                 ! Set an initial search index for the 3D and 4D array interpolation
              LastIndex3 = (/0,0,0/)
              LastIndex4 = (/0,0,0,0/)
  
                 ! To make things run slightly quicker, copy the data we will be interpolating over into the temporary arrays
              IF (MnDriftData%DataIs3D) THEN
                 TmpData3D = MnDriftData%Data3D%DataSet(:,:,:,I)
              ELSE
                 TmpData4D = MnDriftData%Data4D%DataSet(:,:,:,:,I)
              END IF
              
  
              DO J=1,InitInp%NStepWave2
  
                    ! First get the wave amplitude -- must be reconstructed from the WaveElevC0 array.  First index is the real (1) or
                    ! imaginary (2) part.  Divide by NStepWave2 to remove the built in normalization in WaveElevC0.
                 aWaveElevC = CMPLX( InitInp%WaveElevC0(1,J), InitInp%WaveElevC0(2,J), SiKi) / InitInp%NStepWave2
  
                    ! Calculate the frequency
                 Omega1 = J * InitInp%WaveDOmega
  
  
                    ! Only get a QTF value if within the range of frequencies we have wave amplitudes for (first order cutoffs).  This
                    ! is done only for efficiency. 
                 IF ( (Omega1 >= InitInp%WvLowCOff) .AND. (Omega1 <= InitInp%WvHiCOff) ) THEN
  
                       ! Now get the QTF value that corresponds to this frequency and wavedirection pair.
                    IF ( MnDriftData%DataIs3D ) THEN
  
                          ! Set the (omega1,beta1,beta2) point we are looking for.
                       Coord3 = (/ REAL(Omega1,SiKi), InitInp%WaveDirArr(J), InitInp%WaveDirArr(J) /)
  
                          ! get the interpolated value for F(omega1,beta1,beta2)
                       CALL WAMIT_Interp3D_Cplx( Coord3, TmpData3D, MnDriftData%Data3D%WvFreq1, &
                                            MnDriftData%Data3D%WvDir1, MnDriftData%Data3D%WvDir2, LastIndex3, QTF_Value, ErrStatTmp, ErrMsgTmp )
  
                    ELSE
  
                          ! Set the (omega1,omega2,beta1,beta2) point we are looking for.
                       Coord4 = (/ REAL(Omega1,SiKi), REAL(Omega1,SiKi), InitInp%WaveDirArr(J), InitInp%WaveDirArr(J) /)
  
                          ! get the interpolated value for F(omega1,omega2,beta1,beta2)
                       CALL WAMIT_Interp4D_Cplx( Coord4, TmpData4D, MnDriftData%Data4D%WvFreq1, MnDriftData%Data4D%WvFreq2, &
                                            MnDriftData%Data4D%WvDir1, MnDriftData%Data4D%WvDir2, LastIndex4, QTF_Value, ErrStatTmp, ErrMsgTmp )
  
  
                    ENDIF !QTF value find
  
  
                 ELSE     ! outside the frequency range
  
                    QTF_Value = CMPLX(0.0,0.0,SiKi)
  
                 ENDIF    ! frequency check
  
  
                    ! Check and make sure nothing bombed in the interpolation that we need to be aware of
                 CALL SetErrStat(ErrStatTmp,ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
                 IF ( ErrStat >= AbortErrLev ) THEN
                    IF (ALLOCATED(TmpData3D))        DEALLOCATE(TmpData3D,STAT=ErrStatTmp)
                    IF (ALLOCATED(TmpData4D))        DEALLOCATE(TmpData4D,STAT=ErrStatTmp)
                    RETURN
                 ENDIF
  
  
                    ! Now we have the value of the QTF.  These values should only be real for the omega1=omega2 case of the mean drift.
                    ! However if the value came from the 4D interpolation routine, it might have some residual complex part to it.  So
                    ! we throw the complex part out.
                 QTF_Value = CMPLX(REAL(QTF_Value,SiKi),0.0,SiKi)
  
  
                    ! Now put it all together... note the frequency stepsize is multiplied after the summation
                 MnDriftForce(I) = MnDriftForce(I) + REAL(QTF_Value * aWaveElevC * CONJG(aWaveElevC)) !bjj: put QTF_Value first so that if it's zero, the rest gets set to zero (to hopefully avoid overflow issues)
  
              ENDDO ! NStepWave2
  
           ENDIF    ! Load component to calculate
  
        ENDDO
  
  
  
  
negative.code:    !!
     !! where \f$ k \f$ indicates the index to the load component,  \f$ {F_{{ex}~k}^{{-}(2)}} \f$ is the resulting time
     !! independent mean drift force, and \f$ a_m \f$ and \f$ a_m^* \f$ are the complex wave amplitude and its complex conjugate for
     !! the \f$ m^{th} \f$ frequency.  _Note the lack of time dependence in this equation:_ the mean drift is the average drift
     !! force over the entire simulation. Note also that \f$ F_k^{-} (\omega_m, \beta_m) \f$ is the dimensionalized real valued
     !! diagonal of the QTF read from the WAMIT file and interpolated for the \f$ m^{th} \f$ wave frequency.  Note that this is
     !! is a numerical integral, so the \f$ \Delta\omega \f$ term is necessary.
     !!
     !! @note The mean drift is a static value representing the mean drift force for the entire simulation.  Therefore, any offset
     !!       in the body location can be ignored since it cancels out (it is in the complex part that is ignored in the summation).
     !!       The orientation of the platform is therefore handled after the summation.
     !!
     !! Since the frequency range of the QTF has not yet been checked, we will do that now.
     !!
  