SetUpCaseBlock Subroutine

    subroutine SetUpCaseBlock()
      !> Builds and writes block file.
      implicit none
      ! Work variables
      character(str64) :: filename
      real(wp)         :: L(3)
      integer          :: N(3), ngc, Nb(3)
      integer          :: dir
      integer          :: i,j,k
      integer          :: nwall
      integer          :: hif
      integer          :: lof
      integer          :: Nf

      call parser%Get('Block file',           filename              )
      call parser%Get("Domain size",          L                     )
      call parser%Get("Grid points",          N                     )
      call parser%Get("Ghost cells",          ngc                   )
      call parser%Get("Partition",            Nb                    )
      call parser%Get("Number of wall cells", nwall,    default = 4 )

      ! Initialize the main block
      call block%Initialize(ngc,parallel)

      ! Setup the domain periodicity
      call block%SetPeriodicity([.true.,.false.,.true.])

      ! Create a uniform grid, with extra cells
      ! in wall
      ! ------------------------------------------
      associate (lo=>block%lo, hi=>block%hi)
        ! Setup bounds
        lo = [1,1,1]
        hi = N

        ! Initialize axes
        do dir=1,3
          call block%axis(dir)%Initialize(lo(dir),hi(dir),ngc)
        end do
        call block%SetConveniencePointers()
      end associate

      associate (x=>block%x, y=>block%y, z=>block%z, &
        lo =>block%lo, hi=>block%hi)
        ! Uniform x-axis
        do i=lo(1),hi(1)+1
          x(i)= (i-lo(1))*L(1)/real(hi(1)-lo(1)+1,wp)
        end do

        ! y-axis
        Nf  = N(2)-2*nwall  ! Total number of points within fluid
        lof = lo(2)+nwall   ! Minimum position within fluid
        hif = lof + Nf - 1  ! Maximum position within fluid

        ! Inner channel points
        do j=lof,hif+1
          y(j)= (j-lof)*L(2)/real(Nf,wp) - 0.5_wp*L(2)
          !y(j)= 0.5_wp*L(2)*tanh(r*y(j)/(0.5_wp*L(2)))/tanh(r)
        end do
        ! Bottom wall points
        do j=lo(2),lof-1
          y(j) = y(lof)-(lof-j)*(y(lof+1)-y(lof))
        end do
        ! Top wall points
        do j=hif+2,hi(2)+1
          y(j) = y(hif+1)+(j-(hif+1))*(y(hif+1)-y(hif))
        end do

        ! Store IB positions
        IBy(1) = y(lof)
        IBy(2) = y(hif+1)

        ! Uniform z-axis
        do k=lo(3),hi(3)+1
          z(k)= (k-lo(3))*L(3)/real(hi(3)-lo(3)+1,wp) - 0.5_wp*L(3)
        end do
      end associate

      ! Partition block for parallel initializations
      call block%Partition(Nb)

      ! Write block to disk
      call block%Write(filename)

      return
    end subroutine SetUpCaseBlock