I'll put some comments inline to the code, and in the end add a complete implementation of the program, with various changes to improve the execution.

    program BA
     implicit none
     integer, parameter :: N = 20

For real numbers it is advisable to introduce a real kind with selected_real_kind:

     integer, parameter :: rk = selected_real_kind(15)
     integer :: A(n,n)

Use the real kind introduced above, to declare the reals:

     real(kind=rk) :: r, p, r2(2), pj
     real(kind=rk) :: sumA_q
     integer :: sumA
     integer :: m0, i, j, L0, N1, t, jnod, k, di

Often it is better to use more than single characters for variables (for example something like `iNode` or `iEdge`, this makes it visible what the iterator is actually used for.

     integer :: edge_number, ii, ie, je, degree
    
     A = 0
     m0 = 5 ! initial nodes
     L0 = 4 ! initial edges
     ii = 0
     ! adding m0 nodes with random connections
     rand_cnx: do ! while is actually deprecated, you anyway use 'if exit' in the
                  ! loop, so I'd stick to that one. Using a block label to allow
                  ! the identification of the loop in the exit statement.
       do ie = 1, m0
         call random_number(r)

You could create a bunch of `m0` random numbers all at once here, and use them then subsequently. That might be a little faster, but probably has no really big impact.

         je = r * (m0 - 1) + 1
         if (A(ie,je) == 0 ) then
           A(ie,je) = 1
           A(je,ie) = 1
           ii = ii + 1
           if (ii>=L0) EXIT rand_cnx

No need to compute another random number here, if the limit was reached, leave the loop. Use a block label to directly leave the outer loop.

         end if
       end do
     end do rand_cnx
    
     edge_number = L0
     N1 = m0
     degree = m0

     ! Precompute sum of submatrix up to m0-1
     sumA = sum(A(1:m0-1, 1:m0-1))

     ! adding remaining nodes with pereferential attachment
     do t = 1,N-m0
       di = 0
       N1 = m0 + t

The summation over the subarray does not depend on the random numbers (only entries outside the subarray at N1 are modified, if I am not mistaken) and can be pre-computed. As N1 is linearly increasing, we can update the sum of the submatrix up to N1, by just adding the newly filled column and row from the previous iteration.

       ! Update sum of submatrix with the previously computed column and row.
       sumA = sumA + sum(A(1:N1-1,N1-1)) + sum(A(N1-1,1:N1-2))
       sumA_q = 1.0_rk / real(sumA, kind=rk)

       do
         call random_number(r)     
         jnod = r * (N1-2)+1
         pj = sum(A(jnod,1:N1-1))*sumA_q
         call random_number(r)
         if ( r<pj .and. A(N1,jnod)==0 ) then        
            A(N1,jnod) = 1
            A(jnod,N1) = 1
            edge_number = edge_number + 1
            di = sum(A(N1,1:N1))
            if (di >= degree) EXIT
         end if
       end do
     end do
    
     do i = 1,N
        write(12,"(*(I5))") (A(i,j),j=1,N)
     end do
    end program BA

With a little more memory effort, you could also precompute all `sum(A(jnod,1:N1-1))`, use them for the summation of the complete submatrix and in the computation of `pj`, instead of recalculating the sum everytime the random number picks a specific `jnod`. This would result in turning the `di` into an array of length `N`. Implementing this led be to find that the first additional node has to be connected to all previous nodes, as you set `degree=m0`. In this case no randomness seems to be involved? I catch this special case in the implementation below:


    program BA
      implicit none
      integer, parameter :: N = 20
    
      integer, parameter :: rk = selected_real_kind(15)
      integer :: A(n,n)
    
      real(kind=rk) :: r, pj
      real(kind=rk) :: sumA_q
      integer :: sumA
      integer :: m0, i, j, L0, N1, t, jnod
      integer :: di(N)
    
      integer :: edge_number, ii, ie, je, degree
    
      A = 0
      di = 0
    
      m0 = 5 ! initial nodes
      L0 = 4 ! initial edges
    
      ! adding m0 nodes with random connections
      ii = 0
      rand_cnx: do
    
        do ie = 1, m0
          call random_number(r)
          je = nint(r * (m0 - 1)) + 1
          if (A(ie,je) == 0 ) then
            A(ie,je) = 1
            A(je,ie) = 1
            di(ie) = di(ie) + 1
            di(je) = di(je) + 1
            ii = ii + 1
            if (ii>=L0) EXIT rand_cnx
          end if
        end do
    
      end do rand_cnx
    
      edge_number = L0
      N1 = m0
      degree = m0
    
      ! Precompute sum of submatrix up to m0-1
      sumA = sum(di(:m0-1)) - sum(A(:m0-1,m0))
    
      ! adding remaining nodes with pereferential attachment
      do t = 1,N-m0
        N1 = m0 + t
    
        ! Update sum of submatrix with the previously computed column and row.
        sumA = sumA + sum(A(:N1-2,N1-1)) + di(N1-1)
        sumA_q = 1.0_rk / real(sumA, kind=rk)
    
        need_all: if (N1-1 == degree) then
    
          ! Need all nodes to fulfill degree, no randomness required?
          A(N1,:N1-1) = 1
          A(:N1-1,N1) = 1
          edge_number = edge_number + N1-1
          di(:N1-1) = di(:N1-1) + 1
          di(N1) = di(N1) + N1-1
    
        else need_all
    
          do
            call random_number(r)
            jnod = nint(r * (N1-2))+1
            ! Only consider not yet connected nodes.
            if (A(jnod,N1) == 0) then
              pj = di(jnod)*sumA_q
              call random_number(r)
              if ( r<pj ) then
                 A(N1,jnod) = 1
                 A(jnod,N1) = 1
                 edge_number = edge_number + 1
                 di(N1) = di(N1)+1
                 di(jnod) = di(jnod) + 1
                 if (di(N1) >= degree) EXIT
              end if
            end if
          end do
    
        end if need_all
    
      end do
    
      do i = 1,N
         write(12,"(*(I5))") (A(i,j),j=1,N)
      end do
    
    end program BA

No guarantees, that this is in any way correct. I guess, this could be further increased by just picking the random nodes out of the not yet connected ones. This would imply that you need to keep an index list of all 0 values for the current node, and then use `jnod` to pick out of those. This would be more involved book-keeping but would reduce the number of guesses, you need to do.