SetUpCaseFields Subroutine

    subroutine SetUpCaseFields()
      !> Set the initial flow field
      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)             :: r
      real(wp)             :: a
      real(wp)             :: Vd
      real(wp)             :: Vr, Vt
      real(wp)             :: costheta
      real(wp)             :: sintheta
      real(wp)             :: center(3)
      real(wp)             :: rho
      real(wp)             :: mu
      real(wp)             :: Rey
      integer              :: i,j,k
      real(wp),parameter   :: twoPi=8.0_wp*atan(1.0_wp)
      ! Bessel firt root and constant C
      real(wp),parameter   :: a1 = 3.8317059702075123115_wp
      real(wp)             :: C

      ! Get info from parser
      call parser%Get("Fields IC file",  filename)
      call parser%Get("Domain size",     L       )
      call parser%Get('Dipole radius',   a       )
      call parser%Get('Dipole velocity', Vd      )
      call parser%Get('Dipole center',   center  )
      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   )


      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

        ! Compute the constant
        C = 2.0_wp/J0(a1)

        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)
              ! U - velocity
              r = sqrt((x(i)-center(1))**2+(ym(j)-center(2))**2)
              costheta= (ym(j)-center(2))/(r+epsilon(1.0_wp))
              sintheta= (x (i)-center(1))/(r+epsilon(1.0_wp))

              if ((r/a).le.1.0_wp) then
                r = a1*r/a
                Vr = Vd * ( C*J1(r)/(r+epsilon(1.0_wp)) - 1.0_WP ) * costheta
                Vt = Vd * ( 1.0_WP - C*(J0(r)-J1(r)/(r+epsilon(1.0_wp))) ) * sintheta
              else
                Vr = - Vd/(r/a)**2 * costheta
                Vt = - Vd/(r/a)**2 * sintheta
              end if
              V(1)%cell(i,j,k) = Vr*sintheta+Vt*costheta

              ! V - velocity
              r = sqrt((xm(i)-center(1))**2+(y(j)-center(2))**2)
              costheta= (y (j)-center(2))/(r+epsilon(1.0_wp))
              sintheta= (xm(i)-center(1))/(r+epsilon(1.0_wp))


              if ((r/a).le.1.0_wp) then
                r = a1*r/a
                Vr = Vd * ( C*J1(r)/(r+epsilon(1.0_wp)) - 1.0_WP ) * costheta
                Vt = Vd * ( 1.0_WP - C*(J0(r)-J1(r)/(r+epsilon(1.0_wp))) ) * sintheta
              else
                Vr = - Vd/(r/a)**2 * costheta
                Vt = - Vd/(r/a)**2 * sintheta
              end if
              V(2)%cell(i,j,k) = Vr*costheta-Vt*sintheta
            end do
          end do
        end do

        V(3) = 0.0_wp
        P    = 0.0_wp
      end associate

      ! Write some info to stdout
      if (parallel%RankIsRoot()) then
        ! Reynolds number
        Rey = rho*Vd*a/(mu+epsilon(1.0_wp))
        write(*,*) "Reynolds number = ", Rey
      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