SetUpCaseFields Subroutine

    subroutine SetUpCaseFields()
      !> Builds and writes initial flow fields.
      implicit none
      ! Work variables
      type(Eulerian_set)   :: fields
      type(eulerian_obj_r) :: V(3)
      type(eulerian_obj_r) :: P
      character(str64)     :: filename
      real(wp)             :: L(3)
      real(wp)             :: Reb
      real(wp)             :: h
      real(wp)             :: Vbulk
      real(wp)             :: Vfric
      real(wp)             :: Vlaminar
      real(wp)             :: rand
      real(wp)             :: amp
      real(wp)             :: rho
      real(wp)             :: mu
      real(wp)             :: Cf
      real(wp)             :: tau_w
      real(wp)             :: pgrad
      real(wp)             :: ymin,ymax
      integer              :: i,j,k
      real(wp), parameter  :: twoPi=8.0_wp*atan(1.0_wp)
      real(wp), parameter  :: Re_transition = 650.0_wp

      ! Get info from parser
      call parser%Get("Fields IC file",        filename)
      call parser%Get("Domain size",           L       )
      call parser%Get('Bulk Reynolds',         Reb     )
      call parser%Get('Fluctuation amplitude', amp     )
      call parser%Get("Fluid density",         rho     )
      call parser%Get("Fluid viscosity",       mu      )

      ! 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   )


      ! Channel half height and mean velocity
      h     = 0.5_wp*L(2)
      Vbulk = Reb*(mu/rho)/h

      associate (lo => block%lo, hi=> block%hi,          &
          x =>block%x , y =>block%y,  z => block%z, &
          xm=>block%xm, ym=>block%ym, zm=> block%zm)

        V(1) = 0.0_wp
        V(2) = 0.0_wp
        V(3) = 0.0_wp
        P    = 0.0_wp

        ! Base laminar flow
        ymin = block%pmin(2)
        ymax = block%pmax(2)
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              ! Base laminar contribution
              Vlaminar = Vbulk*6.0_wp*(ym(j)-ymin)*(ymax-ym(j))/(ymax-ymin)**2
              ! Add perturbation to make transition faster
              call random_number(rand)
              V(1)%cell(i,j,k) = Vlaminar + Vbulk*amp*(rand-0.5_wp)*cos(8.0_wp*twoPi*zm(k)/L(3))
              call random_number(rand)
              V(3)%cell(i,j,k) = 0.0_wp  +  Vbulk*amp*(rand-0.5_wp)*cos(8.0_wp*twoPi*xm(i)/L(1))
            end do
          end do
        end do
      end associate

      ! Determine friction factor
      if (Reb.le.Re_transition) then
        Cf    = 6.0_wp/Reb
      else
        ! Dean's formula
        Cf    = 0.0614_wp/Reb**0.25_wp
      end if

      ! Determine other flow parameters
      tau_w = 0.5_wp*rho*Vbulk**2*Cf
      pgrad = -tau_w/h
      Vfric = sqrt(tau_w/rho)

      ! Write some info to stdout
      if (parallel%RankIsRoot()) then
        write(*,*) "Bulk Reynolds     = ", Reb
        write(*,*) "Friction Reynolds = ", rho*Vfric*h/mu
        write(*,*) "Bulk velocity     = ", Vbulk
        write(*,*) "Friction velocity = ", Vfric
        write(*,*) "Pressure gradient = ", pgrad
        write(*,*) "min dx(1)/del_v   = ", block%dx(1)*(rho*Vfric)/mu
        write(*,*) "min dx(2)/del_v   = ", block%dx(2)*(rho*Vfric)/mu
        write(*,*) "min dx(3)/del_v   = ", block%dx(3)*(rho*Vfric)/mu
      end if

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

      ! Clear data
      call fields%Finalize()

      return
    end subroutine SetUpCaseFields