c***********************************************************************

      SUBROUTINE simq(aa,n,b)

c routine PUBLIC du module mod_numerique

c résolution d'un système linéaire Ax=b, routine de la ssp 360

c entrées:
c     aa(n,n): matrice A

c entrées/sorties:
c     b: second membre/solution

c Auteur: P. Morel, Département J.D. Cassini, O.C.A.,cesam2k

c--------------------------------------------------------------------

      USE mod_kind

      IMPLICIT NONE

      INTEGER, INTENT(in) :: n
      REAL (kind=dp), INTENT(in), DIMENSION(n,n) :: aa
      REAL (kind=dp), INTENT(inout), DIMENSION(n) :: b

      REAL (kind=dp), DIMENSION(n*n) :: a

      REAL (kind=dp) :: tol, biga, sav
      INTEGER :: jj, j, jy, it, i, ij, imax, i1, i2, iqs, ix, ixj, jx, ixjx,
     1 jjx, ny, ia, ib, ic, k

c----------------------------------------------------------------------

c transformation aa(n,n) ==> aa(n*n)
      a=RESHAPE(aa,(/ n*n /))

c forward solution
      tol=0.d0 ; jj=-n
      DO j=1,n
       jy=j+1 ; jj=jj+n+1 ; biga=0.d0 ; it=jj-j
       DO i=j,n

c search for maximum coefficient in column
        ij=it+i
        IF(ABS(biga)-ABS(a(ij)) < 0.d0)THEN
20       biga=a(ij) ; imax=i
        ENDIF
       ENDDO

c test for pivot less than tolerance (singular matrix)
       IF(ABS(biga)-tol <= 0.d0)THEN
        PRINT 36,biga,tol ; STOP
36      FORMAT('dans simq biga=',es10.3,' < tol=',es10.3)
       ENDIF

c interchange rows if necessary
       i1=j+n*(j-2) ; it=imax-j
       DO k=j,n
        i1=i1+n ; i2=i1+it ; sav=a(i1) ; a(i1)=a(i2) ; a(i2)=sav

c divide equation by leading coefficient
        a(i1)=a(i1)/biga
       ENDDO
       sav=b(imax) ; b(imax)=b(j) ; b(j)=sav/biga

c eliminate next variable
       IF(j-n) 55,70,55
55     iqs=n*(j-1)
       DO ix=jy,n
        ixj=iqs+ix ; it=j-ix
        DO jx=jy,n
         ixjx=n*(jx-1)+ix ; jjx=ixjx+it ; a(ixjx)=a(ixjx)-(a(ixj)*a(jjx))
        ENDDO
        b(ix)=b(ix)-(b(j)*a(ixj))
       ENDDO
      ENDDO

c back solution
70    ny=n-1
      it=n*n
      DO j=1,ny
       ia=it-j ; ib=n-j ; ic=n
       DO k=1,j
        b(ib)=b(ib)-a(ia)*b(ic) ; ia=ia-n ; ic=ic-1
       ENDDO
      ENDDO

      RETURN

      END SUBROUTINE simq