Segmentation fault - invalid memory reference

You haven't posted your code for MY_LIB.f03 so we cannot compile CFDtest.f03 exactly as you have supplied it.

(As an aside, the usual naming convention is that f90 in a .f90 file is not supposed to imply the language version being targeted. Rather, .f90 denotes free format while .f is used for fixed format. By extension, your .f03 files would be better (i.e., more portable if) named as .f90.)

I commented out the USE MY_LIB line and ran your code through nagfor -u -c cfd_test.f90. The output, broken down, is

Extension: cfd_test.f90, line 13: Byte count on numeric data type
           detected at *@8
Extension: cfd_test.f90, line 15: Byte count on numeric data type
           detected at *@8

Byte counts are not portable. The kind value for an 8-byte integer is selected_int_kind(18). (Similarly you might like to use a kind(0.0d0) kind value for your double precision data.)

Error: cfd_test.f90, line 48: Implicit type for I
       detected at 2@I
Error: cfd_test.f90, line 50: Implicit type for J
       detected at 2@J
Error: cfd_test.f90, line 54: Implicit type for K
       detected at 3@K
Error: cfd_test.f90, line 100: Implicit type for N1
       detected at N1@=

You have these implicitly typed, which implies they are 4-byte (default) integers. You should probably declare these explicitly as 8-byte integers (using the 8-byte integer kind value above) if that's what you intend.

Questionable: cfd_test.f90, line 116: Variable COST set but never referenced
Questionable: cfd_test.f90, line 116: Variable N1 set but never referenced
Warning: cfd_test.f90, line 116: Unused local variable QN
Warning: cfd_test.f90, line 116: Unused local variable WKSPCE

You need to decide what you intend to do with these, or whether they are just deletable cruft.

With the implicit integers declared explicitly, there is further output

Warning: cfd_test.f90, line 116: Variable SO referenced but never set

This looks bad.

Obsolescent: cfd_test.f90, line 66: 2 is a shared DO termination label

Your DO loops would probably be better using the modern END DO terminators (not shared!)

Error: cfd_test.f90, line 114: RETURN is only allowed in SUBROUTINEs and FUNCTIONs

This is obviously easy to fix.

For the LAPACK call, one source of explicit interfaces for these routines is the NAG Fortran Library (through the nag_library module). Since your real data is not single precision, you should be using dgesv instead of sgesv. Adding USE nag_library, ONLY: dgesv and switching to call dgesv instead of sgesv, then recompiling as above, reveals

Incorrect data type INTEGER(KIND=4) (expected INTEGER) for argument N (no. 1) of DGESV

so you should indeed be using default (4-byte integers) - at least for the LAPACK build on your system, which will almost certainly be using 4-byte integers. Thus you might want to forget all about kinding your integers and just use the default integer type for all. Correcting this gives

Array supplied for scalar argument LDA (no. 4) of DGESV

so you do need to add this argument. Maybe pass size(DMA2,1)?

With this argument added to the call the code compiles successfully, but without the definitions for your *LEGL functions I couldn't go through any run-time testing.

Here is my modified (and pretty-printed) version of your program

Program cfd_test

! Use my_lib
! Use nag_library, Only: dgesv
  Implicit None
  Integer, Parameter :: wp = kind(0.0D0)
  Real (Kind=wp) :: ef, oper, sigma1, sigma2, tau
  Integer :: i, inf, j, k, n, sum
  Real (Kind=wp) :: dma(0:10, 0:10), dma2(0:10, 0:10), et(0:10), fu(0:10), &
    so(0:10), vn(0:10), wt(0:10)
  Integer :: pivot(10)
  External :: dgesv, dmlegl, welegl, zelegl
  Intrinsic :: kind, size

! SET THE PARAMETERS

  sigma1 = 0._wp

  sigma2 = 0._wp

  tau = 1._wp

  ef = 1._wp

  Do n = 2, 10

!   COMPUATION OF THE NODES, WEIGHTS AND DERIVATIVE MATRIX

    Call zelegl(n, et, vn)

    Call welegl(n, et, vn, wt)

    Call dmlegl(n, 10, et, vn, dma)


!   CONSTRUCTION OF THE MATRIX CORRESPONDING TO THE

!   DIFFERENTIAL OPERATOR

    Do i = 0, n

      Do j = 0, n

        sum = 0._wp

        Do k = 0, n

          sum = sum + dma(i, k)*dma(k, j)

        End Do

        oper = -sum

        If (i==j) oper = -sum + tau

        dma2(i, j) = oper

      End Do

    End Do

!   CHANGE OF THE ENTRIES OF THE MATRIX ACCORDING TO THE

!   BOUNDARY CONDITIONS

    Do j = 0, n

      dma2(0, j) = 0._wp

      dma2(n, j) = 0._wp


    End Do

    dma2(0, 0) = 1._wp

    dma2(n, n) = 1._wp


!   CONSTRUCTION OF THE RIGHT-HAND SIDE VECTOR

    Do i = 1, n - 1

      fu(i) = ef

    End Do

    fu(0) = sigma1

    fu(n) = sigma2

!   SOLUTION OF THE LINEAR SYSTEM

    Call dgesv(n, n, dma2, size(dma2,1), pivot, fu, n, inf)

    Do i = 0, n

      fu(i) = so(i)

    End Do

    Print *, pivot

  End Do

End Program

In general your development experience will be the most pleasant if you use as good a checking compiler as you can get your hands on and if you make sure you ask it to diagnose as much as it can for you.