fill_corner_region_2d Subroutine

  subroutine fill_corner_region_2d(input, output)
    !! Fill corner region for 2D padding
    implicit none

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

    ! Local variables
    integer :: i, j, m, s, f
    integer :: step, idx_in, idx_out, idx_shift
    integer :: input_h, input_w, output_h, output_w
    integer :: pad_h, pad_w
    integer, dimension(2,2) :: orig, dest

    input_h = input%shape(1)
    input_w = input%shape(2)
    output_h = output%shape(1)
    output_w = output%shape(2)
    pad_h = output%indices(2)
    pad_w = output%indices(3)

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

       do concurrent( s = 1:size(output%val, dim=2), m = 1:output%shape(3) )
          select case(output%indices(1))
          case(3, 4) ! circular or reflection
             step = merge(1, -1, output%indices(1) .eq. 3)
             do j = dest(1,2), dest(2,2)
                do i = dest(1,1), dest(2,1)
                   idx_out = i + (j-1) * output_h + (m - 1) * output_h * output_w
                   idx_in = orig(1,1) + step * (i - dest(1,1)) + &
                        (orig(1,2) + step * (j - dest(1,2)) - 1) * input_h + &
                        (m - 1) * input_h * input_w
                   output%val(idx_out, s) = input%val(idx_in, s)
                end do
             end do
          case(5) ! replication
             idx_in = orig(1,1) + (orig(1,2) - 1) * input_h + &
                  (m - 1) * input_h * input_w
             do j = dest(1,2), dest(2,2)
                do i = dest(1,1), dest(2,1)
                   idx_out = i + (j-1) * output_h + (m - 1) * output_h * output_w
                   output%val(idx_out, s) = input%val(idx_in, s)
                end do
             end do
          end select
       end do
    end do

  end subroutine fill_corner_region_2d