BuildFourierModes Subroutine

    subroutine BuildFourierModes(Vrms,peak_mode,Nk,Vk)
      !> Builds the Fourier modes, scaled with the provided
      ! energy spectrum.
      ! Note that Since the FFT of a real field satisfies the Hermitian
      ! symmetry, one needs to describe only half of the wavenumber space.
      ! Direction 1 is truncated in half to conform with FFTW.
      implicit none
      real(wp),     intent(in)  :: Vrms                                        !! Wavevector magnitude
      integer,      intent(in)  :: peak_mode                                   !! Peak
      integer,      intent(out) :: Nk(3)                                       !! Size of fourier modes
      complex(wp), allocatable,  &
                    intent(out) :: Vk(:,:,:,:)                                 !! Fourier modes
      ! Work variables
      real(wp)             :: A
      real(wp)             :: kp
      real(wp)             :: Ek
      real(wp)             :: kmag
      real(wp)             :: phi
      real(wp)             :: kappa(3)
      real(wp)             :: rand(3)
      real(wp)             :: dk(3)
      integer              :: mi,mj,mk
      real(wp),  parameter :: Pi = 4.0_wp*atan(1.0_wp)
      complex(wp),          &
                 parameter :: ii = (0.0_wp,1.0_wp)

      ! Infinitesimal wavenumbers
      dk = 2.0_wp*Pi/(block%pmax-block%pmin)

      ! Peak wavenumber
      kp = real(peak_mode,wp)*maxval(dk)

      ! Spectrum prefactor
      A = 16.0_wp/sqrt(0.5_wp*Pi)*Vrms**2/kp**5

      ! Number of fourier modes in each direction
      Nk = [Ng(1)/2+1, Ng(2), Ng(3)]

      allocate(Vk(3,Nk(1),Nk(2),Nk(3)),mold=(0.0_wp,0.0_wp))

      ! Build Fourier modes following Rogallo's method
      do mk=1,Nk(3)

        ! Wavenumber in 3-dir
        kappa(3) = real(mk-1,wp)*dk(3)
        if (mk.gt.Ng(3)/2+1) kappa(3) = -real(Ng(3)-(mk-1),wp)*dk(3)

        do mj=1,Nk(2)

          ! Wavenumber in 2-dir
          kappa(2) = real(mj-1,wp)*dk(2)
          if (mj.gt.Ng(2)/2+1) kappa(2) = -real(Ng(2)-(mj-1),wp)*dk(2)

          do mi=1,Nk(1)

            ! Wavenumber in 1-dir
            kappa(1) = real(mi-1,wp)*dk(1)

            ! Wavevector magnitude
            kmag=norm2(kappa)

            if (kmag.gt.1.0e-10_wp) then
              ! Energy at this shell
              Ek = GetEnergySpectrum(kmag,kp,A)

              ! Generate a random phase
              call random_number(phi)
              phi    = 2.0_wp*Pi*phi

              ! Generate a complex vector with a random amplitude and the previous random phase
              call random_number(rand)
              Vk(:,mi,mj,mk) = (rand - 0.5_wp) * exp(ii * phi)

              ! Make it orthogonal to the wavevector to satisfy continuity
              Vk(:,mi,mj,mk) = Vk(:,mi,mj,mk) - (1.0_wp/kmag**2)*dot_product(kappa,Vk(:,mi,mj,mk))*kappa

              ! Rescale to match target energy spectrum
              Vk(:,mi,mj,mk) = Vk(:,mi,mj,mk) * sqrt(Ek / (4.0_wp * pi * kmag**2))

            else
              Vk(:,mi,mj,mk) = (0.0_wp,0.0_wp)
            end if

          end do
        end do
      end do

      return
    end subroutine BuildFourierModes