importing fortran/c in python without f2py

 To import fortran (probably c too) library without f2py, use following approach. 

(Although it is possible, it seems to be better to f2py if possible.)

(1) modify the code or wrapper as 

subroutine DJ16py_mom(x1mev,x2mev,L1,L2,S,J,res_mev2) &
    bind(C, name="dj16py_mom")

    use iso_c_binding
    implicit none

    real(c_double), intent(in) :: x1mev, x2mev
    integer(c_int), intent(in) :: L1,L2,S,J
    real(c_double), intent(out) :: res_mev2

    ! Original calculation or call of subroutine 
end subroutine

(2) Compile as shared library. 

FC = gfortran
FFLAGS = -O2 -fPIC -ffree-line-length-none
#or FC=ifort FFLAGS="-O2 -fPIC -free"
SRC = DJ16.f90
OBJ = $(SRC:.f90=.o)

all: libDJ16.so

%.o: %.f90
	$(FC) $(FFLAGS) -c $< -o $@

libDJ16.so: $(OBJ)
	$(FC) -shared -o $@ $(OBJ)

clean:
	rm -f *.o *.mod *.so

(3) In python, one calls them 

import ctypes

lib = ctypes.CDLL("./libDJ16.so")

func = lib.dj16py_mom

func.argtypes = (
    [ctypes.POINTER(ctypes.c_double)] * 2
    + [ctypes.POINTER(ctypes.c_int)] * 4
    + [ctypes.POINTER(ctypes.c_double)]
)
func.restype = None

Now one can call "lib.dh16py_mom" 

In case of arrays, 

subroutine calculate_matrix(n, a)
    integer, intent(in) :: n
    real(8), intent(inout) :: a(n,n)
end subroutine

import numpy as np
import ctypes

n = 10
a = np.zeros((n, n), dtype=np.float64, order="F")

func = lib.calculate_matrix_
func.argtypes = [
    ctypes.POINTER(ctypes.c_int),
    ctypes.POINTER(ctypes.c_double)
]
func.restype = None

n_c = ctypes.c_int(n)

func(
    ctypes.byref(n_c),
    a.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
)