anchor.diff:  
         IF ( .NOT. MnDriftData%DataIs3D .AND. MnDriftData%Data4D%WvFreqDiagComplete ) THEN
            TmpFlag = .FALSE.    ! if this goes true, then we need to convert to 3D data
  -         DO I=1,MnDriftData%Data4D%NumBodies
  -            IF ( p%MnDriftDims(I) ) THEN        ! Flag indicating which dimension we are calculating for
  -               IF ( MnDriftData%Data4D%DataIsSparse(I) .AND. MnDriftData%Data4D%LoadComponents(I) )      TmpFlag = .TRUE.
  -            ENDIF
  +         DO IBody=1,MnDriftData%Data4D%NumBodies
  +            DO ThisDim=1,6
  +               Idx = (IBody-1)*6+ThisDim
  +               IF ( p%MnDriftDims(IBody) ) THEN        ! Flag indicating which dimension we are calculating for
  +                  IF ( MnDriftData%Data4D%DataIsSparse(Idx) .AND. MnDriftData%Data4D%LoadComponents(Idx) )      TmpFlag = .TRUE.
  +               ENDIF
  +            ENDDO
            ENDDO
   
               ! If we need to create the 3D data set, then
  
positive.code:          ! Make sure we have a value for the mean drift for each of the dimensions we were requested to calculate. To
           ! find out we will just check the status of the flags for each dimension in the MnDriftDims against the ones
           ! for the DiffQTFDims
        DO I=1,6
           IF ( p%DiffQTFDims(I) .AND. (.NOT. p%MnDriftDims(I)) ) &
              CALL SetErrStat( ErrID_Fatal,' The DiffQTF method requires the use of the MnDrift method for the first term. '// &
                       'Something went wrong and the MnDrift method returned a different number of load components.', &
                       ErrStat,ErrMsg,RoutineName)
  
negative.code:                         RETURN
                       ENDIF
  
                          ! Set the value of the first term in the frequency domain
                       Term1ArrayC(2*J) = aWaveElevC1 * aWaveElevC1 * QTF_Value
  
  
                    ENDIF    ! Check on the limits
                 ENDDO       ! First term calculation
  
  
  
                 !---------------------------------------------------------------------------------
                 ! Calculate the second term.
                 ! In this term, we are are now stepping through the sum frequencies.  The inner
                 ! sum essentially covers all the off diagonal terms (omega_m /= omega_n).  The limits
                 ! on the outer integral that is the FFT run through the full frequency range that
                 ! we are using
  
  
                    ! Set an initial search index for the 4D array interpolation
                 LastIndex4 = (/0,0,0,0/)
  
                    ! Initialize the temporary arrays for each term to zero.
                 Term2ArrayC = CMPLX(0.0_SiKi,0.0_SiKi,SiKi)
  
  
                    ! Check the limits for the high frequency cutoff. If WvHiCOffS is less than the
                    ! maximum frequency possible with the value of WaveDT (omega_max = pi/WaveDT = NStepWave2*WaveDOmega),
                    ! then we are good.  If the WvHiCOff > 1/2 omega_max, then we will be potentially
                    ! throwing away information.  However, remember the following:
                    !     WaveDT   Omega_max    wavelength
                    !      (s)      (rad/s)        (m)
                    !       .25       4 Pi           < 1
                    !      0.5        2 Pi           1
                    !      1.0          Pi          10
                    ! so, we don't need a really small WaveDT
  
           !This section has been removed since it is kind of annoying.
           !      IF ( InitInp%WvHiCOffS > InitInp%NStepWave2*InitInp%WaveDOmega ) THEN
           !         CALL SetErrStat( ErrID_Warn,' The high frequency cutoff for second order wave forces, WvHiCOffS, '// &
           !                  'is larger than the Nyquist frequency for the given time step of WaveDT. The Nyquist frequency '// &
           !                  '(highest frequency) that can be computed is OmegaMax = PI/WaveDT = '// &
           !                  TRIM(Num2LStr(InitInp%NStepWave2*InitInp%WaveDOmega))// &
           !                  ' radians/second.  If you need those frequencies, decrease WaveDT.  For reference, 2*PI '// &
           !                  'radians/second corresponds to a wavelength of ~1 meter.',&
           !                  ErrStat,ErrMsg,RoutineName)
           !      ENDIF
  
  
  
                    ! Outer loop to create the Term2ArrayC. This is stepwise through the sum frequencies.
                 DO J=1,InitInp%NStepWave2
  
                       ! Calculate the frequency  -- This is the sum frequency.
                    OmegaSum = J * InitInp%WaveDOmega
  
  
  
                       ! Set the \f$ H^+ \f$ term to zero before we start
                    TmpHPlusC = CMPLX(0.0_SiKi,0.0_SiKi,SiKi)
  
  
                       ! Only perform calculations if the difference frequency is in the right range
                    IF ( (OmegaSum >= InitInp%WvLowCOffS) .AND. (OmegaSum <= InitInp%WvHiCOffS) ) THEN
  
                       !> Now do the inner sum.  We are going to perform a sum up to the maximum frequency that we
                       !! can support (Nyquist frequency) for the given WaveDOmega and NStepWave2 (WaveOmegaMax =
                       !! NStepWave2 * WaveDOmega = Pi / WaveDT rad/second). Note that this means the largest diagonal
                       !! frequency we use is \f$ \omega = \Delta\omega * \f$ _NStepwave_/4.
                       !! So, we essentially end up running into the sampling limit.  If we want higher frequency
                       !! terms, we need to use a smaller stepsize.
  
                       DO K=0,FLOOR(Real(J-1)/2.0_SiKi)
  
                             ! Calculate the frequency pair
                          Omega1 =    K  * InitInp%WaveDOmega
                          Omega2 = (J-K) * InitInp%WaveDOmega
  
                             ! Find the wave amplitude at frequency omega.  Remove the NStepWave2 normalization built into WaveElevC0 from Waves module
                          aWaveElevC1 = CMPLX( InitInp%WaveElevC0(1,  K), InitInp%WaveElevC0(2,  K), SiKi )    / InitInp%NStepWave2
                          aWaveElevC2 = CMPLX( InitInp%WaveElevC0(1,J-K), InitInp%WaveElevC0(2,J-K), SiKi )    / InitInp%NStepWave2
  
                             ! Set the (omega1,omega2,beta1,beta2) point we are looking for.
                          Coord4 = (/ REAL(Omega1,SiKi), REAL(Omega2,SiKi), InitInp%WaveDirArr(K), InitInp%WaveDirArr(J-K) /)
  
                             ! get the interpolated value for F(omega1,omega2,beta1,beta2)  --> QTF_Value
                          CALL WAMIT_Interp4D_Cplx( Coord4, TmpData4D, SumQTFData%Data4D%WvFreq1, SumQTFData%Data4D%WvFreq2, &
                                               SumQTFData%Data4D%WvDir1, SumQTFData%Data4D%WvDir2, LastIndex4, QTF_Value, ErrStatTmp, ErrMsgTmp )
                          CALL SetErrStat(ErrStatTmp,ErrMsgTmp,ErrStat,ErrMsg,RoutineName)
                          IF (ErrStat >= AbortErrLev ) THEN
                             IF (ALLOCATED(TmpData4D))        DEALLOCATE(TmpData4D,STAT=ErrStatTmp)
                             RETURN
                          ENDIF
  
                             ! Set the value of the first term in the frequency domain.
                          TmpHPlusC = TmpHPlusC + aWaveElevC1 * aWaveElevC2 * QTF_Value
  
                       ENDDO
  
                          ! Save the value from the summation.
                       Term2ArrayC(J) = TmpHPlusC
  
  
                    ENDIF    ! Check on the limits
  
                 ENDDO       ! Second term calculation -- frequency step on the sum frequency
  
  
  
                    ! Now we apply the FFT to the result of the sum.
                 CALL ApplyFFT_cx(  Term1Array(:),  Term1ArrayC(:), FFT_Data, ErrStatTmp )
                 CALL SetErrStat(ErrStatTmp,'Error occured while applying the FFT to the first term of the Sum QTF.', &
                                ErrStat,ErrMsg,RoutineName)
                 IF ( ErrStat >= AbortErrLev ) THEN
                    IF (ALLOCATED(TmpData4D))        DEALLOCATE(TmpData4D,STAT=ErrStatTmp)
                    RETURN
                 END IF
  
                    ! Now we apply the FFT to the result of the sum.
                 CALL ApplyFFT_cx(  Term2Array(:), Term2ArrayC(:), FFT_Data, ErrStatTmp )
                 CALL SetErrStat(ErrStatTmp,'Error occured while applying the FFT to the second term of the Sum QTF.', &
                                ErrStat,ErrMsg,RoutineName)
                 IF ( ErrStat >= AbortErrLev ) THEN
                    IF (ALLOCATED(TmpData4D))        DEALLOCATE(TmpData4D,STAT=ErrStatTmp)
                    RETURN
                 ENDIF
  
                    ! Now we add the two terms together.  The 0.5 multiplier on is because the double sided FFT was used.
                 DO J=0,InitInp%NStepWave-1  !bjj: Term1Array and Term2Array don't set the last element, so we can get over-flow errors here. SumQTFForce(InitInp%NStepWave,Idx) gets overwritten later, so Idx'm setting the array bounds to be -1.
                    SumQTFForce(J,Idx) = 0.5_SiKi*(REAL(Term1Array(J) + 2*Term2Array(J), SiKi))
                 ENDDO
  
                    ! Copy the last first term to the first so that it is cyclic
                 SumQTFForce(InitInp%NStepWave,Idx) = SumQTFForce(0,Idx)
  
              ENDIF    ! Load component to calculate
  
           ENDDO ! ThisDim -- current dimension
        ENDDO    ! IBody -- current WAMIT body
  
  
  
  