fill_face_region_3d Subroutine

  subroutine fill_face_region_3d(input, output)
    !! Fill face region for 3D padding
    implicit none

    ! Arguments
    type(array_type), intent(in) :: input
    type(array_type), intent(inout) :: output

    ! Local variables
    integer :: i, j, k, m, s, f, idim
    integer :: step1, step2, step3, idx_in, idx_out
    integer :: input_h, input_w, input_d
    integer :: output_h, output_w, output_d
    integer, dimension(2,3) :: orig, dest

    input_h = input%shape(1)
    input_w = input%shape(2)
    input_d = input%shape(3)
    output_h = output%shape(1)
    output_w = output%shape(2)
    output_d = output%shape(3)

    do f = 1, output%indices(5)
       idim = output%indices(7 + f)
       orig(1:2,1) = output%adj_ja(1,(f-1)*6 + 1:(f-1)*6 + 2)
       orig(1:2,2) = output%adj_ja(1,(f-1)*6 + 3:(f-1)*6 + 4)
       orig(1:2,3) = output%adj_ja(1,(f-1)*6 + 5:(f-1)*6 + 6)
       dest(1:2,1) = output%adj_ja(2,(f-1)*6 + 1:(f-1)*6 + 2)
       dest(1:2,2) = output%adj_ja(2,(f-1)*6 + 3:(f-1)*6 + 4)
       dest(1:2,3) = output%adj_ja(2,(f-1)*6 + 5:(f-1)*6 + 6)

       do concurrent( s = 1:size(output%val, dim=2), m = 1:output%shape(4) )
          select case(output%indices(1))
          case(3, 4) ! circular or reflection
             step1 = merge(-1, 1, output%indices(1) .eq. 4 .and. idim .eq. 1)
             step2 = merge(-1, 1, output%indices(1) .eq. 4 .and. idim .eq. 2)
             step3 = merge(-1, 1, output%indices(1) .eq. 4 .and. idim .eq. 3)
             do k = dest(1,3), dest(2,3)
                do j = dest(1,2), dest(2,2)
                   do i = dest(1,1), dest(2,1)
                      idx_out = i + (j-1) * output_h + &
                           (k-1) * output_h * output_w + &
                           (m - 1) * output_h * output_w * output_d
                      idx_in = orig(1,1) + step1 * (i - dest(1,1)) + &
                           (orig(1,2) + step2 * (j - dest(1,2)) - 1) * &
                           input_h + &
                           (orig(1,3) + step3 * (k - dest(1,3)) - 1) * &
                           input_h * input_w + &
                           (m - 1) * input_h * input_w * input_d
                      output%val(idx_out, s) = input%val(idx_in, s)
                   end do
                end do
             end do
          case(5) ! replication
             select case(idim)
             case(1) ! Face perpendicular to dimension 1
                do k = dest(1,3), dest(2,3)
                   do j = dest(1,2), dest(2,2)
                      idx_in = orig(1,1) + &
                           (j - dest(1,2) + orig(1,2) - 1) * input_h + &
                           (k - dest(1,3) + orig(1,3) - 1) * input_h * input_w + &
                           (m - 1) * input_h * input_w * input_d
                      do i = dest(1,1), dest(2,1)
                         idx_out = i + (j - 1) * output_h + &
                              (k - 1) * output_h * output_w + &
                              (m - 1) * output_h * output_w * output_d
                         output%val(idx_out, s) = input%val(idx_in, s)
                      end do
                   end do
                end do
             case(2) ! Face perpendicular to dimension 2
                do k = dest(1,3), dest(2,3)
                   do i = dest(1,1), dest(2,1)
                      idx_in = i - dest(1,1) + orig(1,1) + &
                           (orig(1,2) - 1) * input_h + &
                           (k - dest(1,3) + orig(1,3) - 1) * input_h * input_w + &
                           (m - 1) * input_h * input_w * input_d
                      do j = dest(1,2), dest(2,2)
                         idx_out = i + (j - 1) * output_h + &
                              (k - 1) * output_h * output_w + &
                              (m - 1) * output_h * output_w * output_d
                         output%val(idx_out, s) = input%val(idx_in, s)
                      end do
                   end do
                end do
             case(3) ! Face perpendicular to dimension 3
                do j = dest(1,2), dest(2,2)
                   do i = dest(1,1), dest(2,1)
                      idx_in = i - dest(1,1) + orig(1,1) + &
                           (j - dest(1,2) + orig(1,2) - 1) * input_h + &
                           (orig(1,3) - 1) * input_h * input_w + &
                           (m - 1) * input_h * input_w * input_d
                      do k = dest(1,3), dest(2,3)
                         idx_out = i + (j - 1) * output_h + &
                              (k - 1) * output_h * output_w + &
                              (m - 1) * output_h * output_w * output_d
                         output%val(idx_out, s) = input%val(idx_in, s)
                      end do
                   end do
                end do
             end select
          end select
       end do
    end do

  end subroutine fill_face_region_3d