fill_edge_region_2d Subroutine

  subroutine fill_edge_region_2d(input, output)
    !! Fill edge 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, idim
    integer :: step1, step2, idx_in, idx_out
    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)

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

       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
             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)
             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) + step1 * (i - dest(1,1)) + &
                        (orig(1,2) + step2 * (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
             select case(idim)
             case(1)
                do j = dest(1,2), dest(2,2)
                   idx_in = orig(1,1) + (j - dest(1,2)) * input_h + &
                        (m - 1) * input_h * input_w
                   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
             case(2)
                idx_in = (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 + i - dest(1,1) + 1, s)
                   end do
                end do
             end select
          end select
       end do
    end do

  end subroutine fill_edge_region_2d