Getting distances of points in 2D space in an array in Fortran using the concept of broadcasting (Python)

Now I've got a bit more time I'll turn my comment into an answer.

Fortran is a compiled language. Python is an interpreted language. This leads to many differences, but one of the big ones is that in Fortran loops are fast, as they have been turned into optimized machine code, and more often than not are the preferred way to express what you want to do. However in Python loops are slow as they have to be (re-)interpreted every step of the way - hence the encouragement to use (to my mind) cryptic single liners invoking optimized functions from the python library. This "vectorisation" is rarely if ever needed in Fortran for simple operations - just write it (to my mind) the simple way with loops.

While I am here I will also mention if you do want performance you have your array indices the wrong way around - you want the one that you change most often ( 1:2 in the above ) to be the first index. This is ultimately to do with the way Fortran organizes it arrays in memory, if interested the key term to look into is "column major"

Finally I will mention the Fortran way is usually to use Subroutines rather than Functions. This is mostly stylistic, but does lead to a few simplifications I won't go into here.

Anyway putting the above together I would write it something like the below

ijb@ijb-Latitude-5410:~/work/stack$ cat not_python.f90
Program test

  Implicit None

  Real, Dimension( :, :, : ), Allocatable :: array1
  Real, Dimension( :, :, : ), Allocatable :: dr
  
  Integer :: d, N, t_size
  Integer :: t_step
  Integer :: id, iN, it
  
  Write( *, * ) 'N, t_size ?'
  Read ( *, * )  N, t_size

  d = 2

  Allocate( array1( 1:d, 1:N, 1:t_size ) )

  Do it = 1, Size( array1, Dim = 3 )
     Do iN = 1, Size( array1, Dim = 2 )
        Do id = 1, Size( array1, Dim = 1 )
           array1( id, iN, it ) = id + 10.0 * in + 100.0 * it
        End Do
     End Do
  End Do

  Allocate( dr( 1:d, 1:n, 1:n ) )

  t_step = 1
  Call f( array1( :, :, t_step ), dr )

  Call printit( dr )
  

Contains

  Pure Subroutine f( a, dr )

    Implicit None

    Real, Dimension(    :, : ), Intent( In    ) :: a
    Real, Dimension( :, :, : ), Intent(   Out ) :: dr

    Integer :: N
    Integer :: i, j
    
    N = Size( a, Dim = 2 )

    ! Could take advantage of symmetry here, dr( j, i ) = - dr( i, j )
    ! but probably gaines you little 
    Do j = 1, N
       Do i = 1, N
          dr( :, i, j ) = a( :, i ) - a( :, j )
       End Do
    End Do

  End Subroutine f

  Subroutine printit( dr )

    Use iso_fortran_env, Only : stdout => output_unit

    Implicit None
    
    Real, Dimension( :, :, : ), Intent( In ) :: dr

    Integer :: it, iN, id
    
    Do it = 1, Size( dr, Dim = 3 )
       Do iN = 1, Size( dr, Dim = 2 )
          Do id = 1, Size( dr, Dim = 1 )
             Write( stdout, '( 3( i2, 1x ), 6x, f5.1 )' ) &
                  id, iN, it, dr( id, iN, it )
          End Do
       End Do
    End Do

  End Subroutine printit

End Program test
ijb@ijb-Latitude-5410:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -O -std=f2018 not_python.f90 -o not_python
ijb@ijb-Latitude-5410:~/work/stack$ ./not_python 
 N, t_size ?
3 2
 1  1  1         0.0
 2  1  1         0.0
 1  2  1        10.0
 2  2  1        10.0
 1  3  1        20.0
 2  3  1        20.0
 1  1  2       -10.0
 2  1  2       -10.0
 1  2  2         0.0
 2  2  2         0.0
 1  3  2        10.0
 2  3  2        10.0
 1  1  3       -20.0
 2  1  3       -20.0
 1  2  3       -10.0
 2  2  3       -10.0
 1  3  3         0.0
 2  3  3         0.0