fill_corner_region_3d Subroutine

  subroutine fill_corner_region_3d(input, output)
    !! Fill corner 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
    integer :: step, idx_in, idx_out, idx_shift
    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)

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

       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
             step = merge(1, -1, output%indices(1) .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) + step * (i - dest(1,1)) + &
                           (orig(1,2) + step * (j - dest(1,2)) - 1) * input_h + &
                           (orig(1,3) + step * (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
             idx_in = orig(1,1) + &
                  (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)
                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
                      output%val(idx_out, s) = input%val(idx_in, s)
                   end do
                end do
             end do
          end select
       end do
    end do

  end subroutine fill_corner_region_3d