SetUpCaseFields Subroutine

    subroutine SetUpCaseFields()
      !> Builds and writes initial fields.
      use iso_fortran_env, only : stdout => output_unit
      implicit none
      ! Work variables
      type(Eulerian_set)   :: fields
      type(eulerian_obj_r) :: V(3)
      type(eulerian_obj_r) :: P
      character(str64)     :: filename
      real(wp)             :: rho
      real(wp)             :: mu
      real(wp)             :: nu
      real(wp)             :: dissipation_rate
      real(wp)             :: Re_lambda
      integer              :: peak_mode
      real(wp)             :: ell
      real(wp)             :: L(3)
      integer              :: Nk(3)
      complex(wp),          &
               allocatable :: Vk(:,:,:,:)
      real(wp)             :: lambda
      real(wp)             :: Vrms
      real(wp)             :: eta_k
      real(wp)             :: dV
      integer              :: i,j,k
      real(wp), parameter  :: C=0.19_wp

      ! Get info from parser
      call parser%Get("Fields IC file",   filename               )
      call parser%Get("Fluid density",    rho                    )
      call parser%Get("Fluid viscosity",  mu                     )
      call parser%Get("Reynolds number",  Re_lambda              )
      call parser%Get("Peak mode",        peak_mode, default = 4 )

      ! Initialize fields container
      call fields%Initialize(block,parallel)

      ! Add fields to container (this will allocate data)
      call fields%Add('V1', 1, V(1))
      call fields%Add('V2', 2, V(2))
      call fields%Add('V3', 3, V(3))
      call fields%Add('P',  0, P   )

      ! Kinematic viscosity
      nu = mu/rho

      ! Domain size
      L = block%pmax - block%pmin

      ! Integral length scale based on Rosales & Meneveau, Physics of Fluids, 2005
      ell = C*minval(L)

      ! Dissipation rate
      dissipation_rate = (nu/15.0_wp)**3*Re_lambda**6/ell**4

      ! Velocity fluctuation rms
      Vrms = (dissipation_rate*ell)**(1.0_wp/3.0_wp)

      call BuildFourierModes(Vrms,peak_mode,Nk,Vk)
      call InvertFourierModes(Nk,Vk,V)

      deallocate(Vk)

      ! Rescale
      Vrms = 0.0_wp
      do k=block%lo(3),block%hi(3)
        do j=block%lo(2),block%hi(2)
          do i=block%lo(1),block%hi(1)
            dV = (block%xm(i+1)-block%xm(i)) &
               * (block%ym(j+1)-block%ym(j)) &
               * (block%zm(k+1)-block%zm(k))
            Vrms = Vrms +(V(1)%cell(i,j,k)**2 &
                        + V(2)%cell(i,j,k)**2 &
                        + V(3)%cell(i,j,k)**2)*dV
          end do
        end do
      end do
      Vrms = sqrt((1.0_wp/3.0_wp)/(L(1)*L(2)*L(3))*Vrms)

      V(1)%cell = V(1)%cell* sqrt(nu*Re_lambda/ sqrt(15.0_wp*nu/dissipation_rate))/Vrms
      V(2)%cell = V(2)%cell* sqrt(nu*Re_lambda/ sqrt(15.0_wp*nu/dissipation_rate))/Vrms
      V(3)%cell = V(3)%cell* sqrt(nu*Re_lambda/ sqrt(15.0_wp*nu/dissipation_rate))/Vrms

      ! Compute urms
      Vrms = 0.0_wp
      do k=block%lo(3),block%hi(3)
        do j=block%lo(2),block%hi(2)
          do i=block%lo(1),block%hi(1)
            dV = (block%xm(i+1)-block%xm(i)) &
               * (block%ym(j+1)-block%ym(j)) &
               * (block%zm(k+1)-block%zm(k))
            Vrms = Vrms +(V(1)%cell(i,j,k)**2 &
                        + V(2)%cell(i,j,k)**2 &
                        + V(3)%cell(i,j,k)**2)*dV
          end do
        end do
      end do
      Vrms = sqrt((1.0_wp/3.0_wp)/(L(1)*L(2)*L(3))*Vrms)

      lambda    = sqrt(15.0_wp*nu/dissipation_rate)*Vrms
      Re_lambda = lambda*Vrms/nu
      eta_k     = (nu**0.75_wp)/(dissipation_rate**0.25_wp)

      write(stdout,*) "Taylor Microscale Reynolds : ", Re_lambda
      write(stdout,*) "Velocity fluctuation rms   : ", Vrms
      write(stdout,*) "Taylor microscale (lambda) : ", lambda
      write(stdout,*) "Kolomogorov length (eta_k) : ", eta_k
      write(stdout,*) "Dissipation rate           : ", dissipation_rate
      write(stdout,*) "Eddy turounover time       : ", Vrms**2/dissipation_rate
      write(stdout,*) "HIT linear forcing         : ", dissipation_rate/(3.0_wp*Vrms**2)
      write(stdout,*) "Resolution (kmax*eta_k)    : ", 4.0_wp*atan(1.0_wp)*minval(Ng/L)*eta_k


      ! Write data to disk
      call fields%SetWriteFileName(filename)
      call fields%Write(0,0.0_wp)

      ! Clear data
      call fields%Finalize()

      return
    end subroutine SetUpCaseFields