mode_splines.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     ######spl
MODULE mode_splines
!-------------------------------------------------------------------------------
!!    MODIFICATIONS
!!    -------------
!!
!!      J.Escobar     06/2013  for REAL4/8 add EPSILON management
!!
!-------------------------------------------------------------------------------
USE modi_abor1_sfx
USE modd_splines
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
USE modd_surf_par , ONLY : xsurf_huge , xsurf_tiny
!
INTERFACE sscopy
        MODULE PROCEDURE sscopy_1
        MODULE PROCEDURE sscopy_2
END INTERFACE
INTERFACE mxidml
        MODULE PROCEDURE mxidml
END INTERFACE
INTERFACE mtxaxm
        MODULE PROCEDURE mtxaxm
END INTERFACE
INTERFACE mxmspl
        MODULE PROCEDURE mxmspl_1
        MODULE PROCEDURE mxmspl_2
        MODULE PROCEDURE mxmspl_3
END INTERFACE
INTERFACE smxinv
        MODULE PROCEDURE smxinv
END INTERFACE
!génération de matrices
INTERFACE splie
        MODULE PROCEDURE splie
END INTERFACE
INTERFACE splk
        MODULE PROCEDURE splk_1
        MODULE PROCEDURE splk_2
END INTERFACE
INTERFACE splbfin
        MODULE PROCEDURE splbfin
END INTERFACE
INTERFACE splt
        MODULE PROCEDURE splt_1
        MODULE PROCEDURE splt_2
END INTERFACE
INTERFACE sple
        MODULE PROCEDURE sple
END INTERFACE
INTERFACE splb2e1
        MODULE PROCEDURE splb2e1
END INTERFACE
INTERFACE splb2e
        MODULE PROCEDURE splb2e
END INTERFACE
INTERFACE tred2
        MODULE PROCEDURE tred2
END INTERFACE
INTERFACE tql2_2
        MODULE PROCEDURE tql2_2
END INTERFACE
INTERFACE eisrs1
        MODULE PROCEDURE eisrs1
END INTERFACE
INTERFACE splv
        MODULE PROCEDURE splv
END INTERFACE 
INTERFACE spls2vi
        MODULE PROCEDURE spls2vi
END INTERFACE
INTERFACE splbvm
        MODULE PROCEDURE splbvm
END INTERFACE
INTERFACE spls2v
        MODULE PROCEDURE spls2v
END INTERFACE
INTERFACE splds2v
        MODULE PROCEDURE splds2v
END INTERFACE
INTERFACE splps2v
        MODULE PROCEDURE splps2v
END INTERFACE
INTERFACE splri
        MODULE PROCEDURE splri
END INTERFACE
INTERFACE splrs
        MODULE PROCEDURE splrs
END INTERFACE
INTERFACE spldrs
        MODULE PROCEDURE spldrs
END INTERFACE
INTERFACE splpr
        MODULE PROCEDURE splpr
END INTERFACE
INTERFACE spldv
        MODULE PROCEDURE spldv
END INTERFACE
INTERFACE spld2v
        MODULE PROCEDURE spld2v
END INTERFACE
INTERFACE splpv
        MODULE PROCEDURE splpv
END INTERFACE
INTERFACE splp
        MODULE PROCEDURE splp
END INTERFACE
INTERFACE spltt
        MODULE PROCEDURE spltt
END INTERFACE
INTERFACE splr
        MODULE PROCEDURE splr_1
        MODULE PROCEDURE splr_2
END INTERFACE
INTERFACE splu
        MODULE PROCEDURE splu
END INTERFACE
INTERFACE splw
        MODULE PROCEDURE splw
END INTERFACE
INTERFACE splvpq
        MODULE PROCEDURE splvpq
END INTERFACE
INTERFACE splc
        MODULE PROCEDURE splc 
END INTERFACE
INTERFACE splm
        MODULE PROCEDURE splm
END INTERFACE
INTERFACE sp0nop
        MODULE PROCEDURE sp0nop
END INTERFACE
INTERFACE sp0cvq
        MODULE PROCEDURE sp0cvq
END INTERFACE
INTERFACE splbsd
        MODULE PROCEDURE splbsd
END INTERFACE
INTERFACE splbsel
        MODULE PROCEDURE splbsel
END INTERFACE
INTERFACE splb2c
        MODULE PROCEDURE splb2c
END INTERFACE

 CONTAINS

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!COPIE D'UNE MATRICE DANS UNE AUTRE
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!1. copie de matrices potentiellement de tailles différentes
!SPLR
SUBROUTINE sscopy_1(A,B,IA,IB)
!
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: a
INTEGER, INTENT(IN) :: ia
INTEGER, INTENT(IN) :: ib
REAL, DIMENSION(:), INTENT(OUT) :: b

INTEGER :: i, j, k
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SSCOPY_1',0,zhook_handle)
j = ib
k = ia
DO i=1,SIZE(a)
  b(j) = a(k)
  j = j + ib
  k = k + ia
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SSCOPY_1',1,zhook_handle)

END SUBROUTINE sscopy_1
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!1. copie de matrices potentiellement de tailles différentes
SUBROUTINE sscopy_2(A,B,IA,IB)
!
IMPLICIT NONE

REAL, DIMENSION(:,:), INTENT(IN) :: a
INTEGER, INTENT(IN) :: ia
INTEGER, INTENT(IN) :: ib
REAL, DIMENSION(:,:), INTENT(OUT) :: b

REAL,DIMENSION(SIZE(A,1)*SIZE(A,2)) :: a1
REAL,DIMENSION(SIZE(B,1)*SIZE(B,2)) :: b1

INTEGER :: n
INTEGER :: i, j, k
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SSCOPY_2',0,zhook_handle)
n=SIZE(a,1)*SIZE(a,2)

DO j=1,SIZE(a,2)
  a1((j-1)*SIZE(a,1)+1:j*SIZE(a,1))=a(:,j)
ENDDO

j = ib
k = ia
DO i=1,n-1
  b1(j) = a1(k)
  j = j + ib
  k = k + ia
ENDDO

DO j=1,SIZE(b,2)
  b(:,j)=b1((j-1)*SIZE(b,1)+1:j*SIZE(b,1))
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SSCOPY_2',1,zhook_handle)

END SUBROUTINE sscopy_2
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!OPERATIONS SUR LES MATRICES
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!2. A(I) B(I*J) R(I*J) R(:,K)=B(:,K)/A(:)
SUBROUTINE mxidml(A,B,R)
!SPLTT, SPLR
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: a
REAL, DIMENSION(:,:), INTENT(IN) :: b
REAL, DIMENSION(:,:), INTENT(OUT):: r

INTEGER :: k
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:MXIDML',0,zhook_handle)
IF (SIZE(b,1).EQ.0 .OR. SIZE(b,2).EQ.0 .AND. lhook) &
        CALL dr_hook('MODE_SPLINES:MXIDML',1,zhook_handle)
IF (SIZE(b,1).EQ.0 .OR. SIZE(b,2).EQ.0) RETURN

DO k=1,SIZE(b,2)
  r(:,k)=b(:,k)/a(:)
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:MXIDML',1,zhook_handle)

END SUBROUTINE mxidml
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!3. produit bizarre de matrices
SUBROUTINE mtxaxm(A,B,R)
!SPLR,SPLU,SPLVPQ
IMPLICIT NONE

REAL, DIMENSION(:,:),INTENT(IN) :: a
REAL, DIMENSION(:,:),INTENT(IN) :: b
REAL, DIMENSION(:,:),INTENT(OUT)::r

INTEGER :: i,j,k, l
REAL :: rij
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:MTXAXM',0,zhook_handle)
r(:,:)=0.
DO i=1,SIZE(a,2)
  DO j=1,SIZE(a,1)
    rij=0.
    DO k=1,SIZE(a,1)
        rij=rij+b(j,k)*a(k,i)
    ENDDO
    DO l=1,i
        r(l,i)=r(l,i)+rij*a(j,l)
        r(i,l)=r(l,i)
    ENDDO
  ENDDO
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:MTXAXM',1,zhook_handle)

END SUBROUTINE mtxaxm
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!4. produit simple de matrices (mettre les deux dans la même interface)
SUBROUTINE mxmspl_1(A,B,R)
!SPLE, SPLC, SPLTT, SPLR, SPLW, SPLVPQ
IMPLICIT NONE

REAL, DIMENSION(:,:),INTENT(IN) :: a
REAL, DIMENSION(:,:),INTENT(IN) :: b
REAL, DIMENSION(:,:),INTENT(OUT)::r

INTEGER :: i,j,k
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:MXMSPL_1',0,zhook_handle)
r(:,:)=0.

DO i=1,SIZE(a,1)
  DO j=1,SIZE(b,2)
    DO k=1,SIZE(a,2)
      r(i,j)=r(i,j)+a(i,k)*b(k,j)
    ENDDO
  ENDDO
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:MXMSPL_1',1,zhook_handle)

END SUBROUTINE mxmspl_1
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!4. produit simple de matrices 2
SUBROUTINE mxmspl_2(A,B,R)
!
IMPLICIT NONE

REAL, DIMENSION(:,:),INTENT(IN) :: a
REAL, DIMENSION(:),INTENT(IN) :: b
REAL, DIMENSION(:),INTENT(OUT)::r

INTEGER :: i,j
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:MXMSPL_2',0,zhook_handle)
r(:)=0.

DO i=1,SIZE(a,1)
  DO j=1,SIZE(a,2)
    r(i)=r(i)+a(i,j)*b(j)
  ENDDO
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:MXMSPL_2',1,zhook_handle)

END SUBROUTINE mxmspl_2    
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!4. produit simple de matrices 3
SUBROUTINE mxmspl_3(A,B,RES)
!
IMPLICIT NONE

REAL, DIMENSION(:),INTENT(IN) :: a
REAL, DIMENSION(:),INTENT(IN) :: b
REAL,INTENT(OUT)::res

INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:MXMSPL_3',0,zhook_handle)
res=0.

DO i=1,SIZE(a)
  res=res+a(i)*b(i)
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:MXMSPL_3',1,zhook_handle)

END SUBROUTINE mxmspl_3
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!5. inversion de la matrice A.
SUBROUTINE smxinv(A,IREP)
!SPLTT
IMPLICIT NONE

REAL, DIMENSION(:,:), INTENT(INOUT) :: a
INTEGER, INTENT(OUT) :: irep

REAL, DIMENSION(SIZE(A,1)*SIZE(A,2)) :: a2
INTEGER, DIMENSION(SIZE(A,1)) :: index
REAL, DIMENSION(SIZE(A,1)) :: ri
INTEGER :: n, np, np1, i, j, jj, k, l, kk, kj, jl
INTEGER :: ij0, ji0, ij, ji
REAL :: pivot, elm
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SMXINV',0,zhook_handle)
irep=0
n=SIZE(a,1)
np1=n+1

index(:)=1

DO j=1,n
 !les colonnes sont à la suite
  a2((j-1)*n+1:j*n)=a(:,j)
ENDDO

DO i=1,n
  !find pivot
  pivot=0.
  jj=1
  !dans la diagonale, on retient le plus grand élément
  DO j=1,n
    IF(index(j).NE.0) THEN   
      elm=abs(a2(jj))
      IF (elm.GT.pivot) THEN
        pivot=elm
        k=j !numéro de la ligne / colonne
        !numéro de l'élément dans A2
        kk=jj
      ENDIF
    ENDIF
    jj=jj+np1
  ENDDO
  index(k)=0
  IF (pivot.EQ.0.) THEN 
    irep=-1
    IF (lhook) CALL dr_hook('MODE_SPLINES:SMXINV',1,zhook_handle)
    RETURN 
  ENDIF
  pivot=-a2(kk)
  !elimination
  kj=k
  np=n
  DO j=1,n
    !KJ=K, KJ=K+N, KJ=K+2N ... : avant l'élément pivot, 
    !on se déplace sur les colonnes
    !J=K=3: KJ=K+2N+1, KJ=K+2N+2....!une fois qu'on a passé
    !le pivot, on se déplace sur la ligne, à droite de l'élément pivot
    !si on est sur le premier élément de la colonne pivot
    IF (j==k) THEN
      !KJ réfère à l'élément diagonal du pivot
      a2(kj)=1./pivot
      ri(j)=0.
      np=1
    ELSE
      !si on est sur un autre élément de du pivot
      elm=-a2(kj)
      ri(j)=elm/pivot
      IF (elm.NE.0.) THEN
        jl=j
        !JL prend pour valeur les élements de la colonne 
        !du premier en haut au diagonal -1
        DO l=1,j
          !à chaque fois A2 de cet élément est incrémenté
          a2(jl)=a2(jl)+elm*ri(l)
          jl=jl+n
        ENDDO
      ENDIF
      !élement véritablement en cours
      a2(kj)=ri(j)
    ENDIF
    !KJ=K+JN si on n'a pas encore passé l'élément de la colonne pivot
    !KJ=K+J0N+J si on l'a passé
    kj=kj+np
  ENDDO
ENDDO

!change the sign and provisional fill-up
DO j=1,n
  a(:,j)= a2((j-1)*n+1:j*n)
ENDDO

DO i=1,n
  DO j=1,i
    a(i,j)=-a(i,j)
    a(j,i)=a(i,j)
  ENDDO
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SMXINV',1,zhook_handle)

END SUBROUTINE smxinv

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! CALCULS PLUS SPECIFIQUES
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!6. IE(ND,M): calcul de la matrice IE en fonction de NDEG.
SUBROUTINE splie(NDEG,IE)
!SPLE, SPLR
IMPLICIT NONE

INTEGER, INTENT(IN) :: ndeg
INTEGER, DIMENSION(:,:), INTENT(OUT) :: ie

INTEGER, DIMENSION(SIZE(IE,1)) :: nv
INTEGER :: i, j, k, t, n, n0, nc
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLIE',0,zhook_handle)
nv(:)=1
ie(:,1)=0
n=1
n0=1

!boucle d'itération
DO t=1,ndeg
!boucle sur les indices
  DO i=1,SIZE(ie,1)
    nc=0
    !NV(I) varie d'itération T en itération
    DO k=nv(i),n0
      n=n+1
      nc=nc+1
      !boucle sur les indices
      DO j=1,SIZE(ie,1)
        ie(j,n)=ie(j,k)
        IF (j.EQ.i) ie(j,n)=ie(j,n)+1
      ENDDO
    ENDDO
    nv(i)=n-nc+1
  ENDDO
  n0=n
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLIE',1,zhook_handle)

END SUBROUTINE splie
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!7. calcul de RK en fonction des autres RK(N1,N2)
SUBROUTINE splk_1(NORD,X1,X2,G,RK)
!SPLE,SPLU,
IMPLICIT NONE

INTEGER, INTENT(IN) :: nord
REAL, DIMENSION(:), INTENT(IN) :: x1
REAL, DIMENSION(:,:), INTENT(IN) :: x2
REAL, DIMENSION(:,:), INTENT(IN) :: g
REAL, DIMENSION(:), INTENT(OUT):: rk

INTEGER :: i, id, jd, exp, isigne
INTEGER :: nd, n, di, dj
REAL :: d2
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLK_1',0,zhook_handle)
nd=SIZE(x1)
n=SIZE(x2,2)

!Exposant du noyau
exp=(2*nord-nd)/2.

IF (mod(nd,2).EQ.0) THEN
  !Calcul de K dans le cas d'un espace de dimension paire
  isigne=(-1)**(1+nord+nd/2)
ELSE
  !Calcul de K dans le cas d'un espace de dimension impair 
  isigne=(-1)**(nord+nd/2)
ENDIF

rk(:)=0.

DO i=1,n
  d2=0.
  DO id=1,nd
    DO jd=1,nd
      di=x1(id)-x2(id,i)
      dj=x1(jd)-x2(jd,i)
      d2=d2+g(id,jd)*di*dj
    ENDDO
  ENDDO
  IF (d2.NE.0.) THEN
    rk(i) = isigne*d2**exp
    IF (mod(nd,2).EQ.0) rk(i)=rk(i)*0.5*log(d2)
  ENDIF
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLK_1',1,zhook_handle)

END SUBROUTINE splk_1
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!7. calcul de RK en fonction des autres RK(N1,N2)
SUBROUTINE splk_2(NORD,X1,X2,G,RK)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nord
REAL, DIMENSION(:,:), INTENT(IN) :: x1
REAL, DIMENSION(:,:), INTENT(IN) :: x2
REAL, DIMENSION(:,:), INTENT(IN) :: g
REAL, DIMENSION(:,:), INTENT(OUT):: rk

INTEGER :: i1, i2, id, jd, isigne
INTEGER :: nd, n1, n2
REAL :: exp, d2, di, dj
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLK_2',0,zhook_handle)
nd=SIZE(x1,1)
n1=SIZE(x1,2)
n2=SIZE(x2,2)

!Exposant du noyau
exp=(2*nord-nd)/2.

IF (mod(nd,2).EQ.0) THEN
  !Calcul de K dans le cas d'un espace de dimension paire
  isigne=(-1)**(1+nord+nd/2)
ELSE
  !Calcul de K dans le cas d'un espace de dimension impair 
  isigne=(-1)**(nord+nd/2)
ENDIF

rk(:,:)=0.

DO i1=1,n1
  DO i2=1,n2
    d2=0.
    DO id=1,nd
      DO jd=1,nd
        di=x1(id,i1)-x2(id,i2)
        dj=x1(jd,i1)-x2(jd,i2)
        d2=d2+g(id,jd)*di*dj
      ENDDO
    ENDDO
    IF (d2.NE.0.) THEN
      rk(i1,i2) = isigne*d2**exp
      IF (mod(nd,2).EQ.0) rk(i1,i2)=rk(i1,i2)*0.5*log(d2)
    ENDIF
  ENDDO
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLK_2',1,zhook_handle)

END SUBROUTINE splk_2
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!8. on cherche l'indice I de la 3ème dimension de AI tel que XE(ID) soit compris 
!entre AI(ID,1,I) et AI(ID,2,I+1). Il est compris entre 1 et NSD(ID), 
!par défaut il est égal à NDSD(ID)
SUBROUTINE splbfin(ID,NSD,XE,AI,IDX,IDX1)
!SPLB2E
IMPLICIT NONE

INTEGER, INTENT(IN) :: id
INTEGER, DIMENSION(:), INTENT(IN) :: nsd
REAL, DIMENSION(:), INTENT(IN) :: xe
REAL, DIMENSION(:,:,:), INTENT(IN) :: ai
INTEGER, INTENT(OUT) :: idx
INTEGER, INTENT(OUT) :: idx1

!INTEGER, PARAMETER: NSDMAX=20
INTEGER :: i, i1
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBFIN',0,zhook_handle)
idx1=0
DO i=1,nsd(id)-1
  IF (xe(id).LE.ai(id,2,i)) THEN
    idx=i
    i1=i+1
    IF (xe(id).GT.ai(id,1,i1)) idx1=i1
    IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBFIN',1,zhook_handle)
    RETURN
  ENDIF
ENDDO
idx=nsd(id)
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBFIN',1,zhook_handle)
     
END SUBROUTINE splbfin
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!9. X(ND,N) IE(ND,M) : T(N,M) T=produit(X**IE) quand IE est non nul
SUBROUTINE splt_1(X,IE,T)
!SPLE, SPLTT
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: x
INTEGER, DIMENSION(:,:), INTENT(IN) :: ie
REAL, DIMENSION(:), INTENT(OUT) :: t
!
INTEGER :: j,k
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLT_1',0,zhook_handle)
t(:)=1.

DO j=1,SIZE(ie,2)
    DO k=1,SIZE(x,1)
      IF (ie(k,j).NE.0) t(j)=t(j)*x(k)**ie(k,j)
    ENDDO
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLT_1',1,zhook_handle)

END SUBROUTINE splt_1
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!9. X(ND,N) IE(ND,M) : T(N,M) T=produit(X**IE) quand IE est non nul
SUBROUTINE splt_2(X,IE,T)
!
IMPLICIT NONE

REAL, DIMENSION(:,:), INTENT(IN) :: x
INTEGER, DIMENSION(:,:), INTENT(IN) :: ie
REAL, DIMENSION(:,:), INTENT(OUT) :: t
!
INTEGER :: i,j,k
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLT_2',0,zhook_handle)
t(:,:)=1.

DO j=1,SIZE(ie,2)
  DO i=1,SIZE(x,2)
    DO k=1,SIZE(x,1)
      IF (ie(k,j).NE.0) t(i,j)=t(i,j)*x(k,i)**ie(k,j)
    ENDDO
  ENDDO
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLT_2',1,zhook_handle)

END SUBROUTINE splt_2

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! PARTIE "E"
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!10.
SUBROUTINE sple(NORD,X,G,C,XE,ZE)
!SPLB2E1
IMPLICIT NONE

INTEGER, INTENT(IN) :: nord
REAL, DIMENSION(:,:),INTENT(IN) :: x !ND,N
REAL, DIMENSION(:,:), INTENT(IN) :: g !ND,ND
REAL, DIMENSION(:),INTENT(IN) :: c ! N+M
REAL, DIMENSION(:),INTENT(IN) :: xe !ND
REAL, INTENT(OUT) :: ze

INTEGER, DIMENSION(SIZE(X,1),SIZE(C)-SIZE(X,2)) :: iw !ND,M
REAL, DIMENSION(SIZE(X,2)) :: we1
REAL, DIMENSION(SIZE(IW,2)) :: we2
REAL :: zout1,zout2
INTEGER :: m, n, ne, m1, m2, m3
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLE',0,zhook_handle)
n=SIZE(x,2)
m=SIZE(c)-n

! Calcul de la matrice Kxe
 CALL splk(nord,xe,x,g,we1)
! Calcul de Kxe.c
 CALL mxmspl(we1,c(1:n),zout1)

IF (m.NE.0) THEN
  ! Generation des exposants des monomes
  CALL splie(nord-1,iw) !IW(ND,M)
  ! Calcul de la matrice Txe
  CALL splt(xe,iw,we2)
  ! Calcul Txe.d
  CALL mxmspl(we2,c(n+1:n+m),zout2)
  ! Calcul de Txe.d + Kxe.c
  ze = zout1+zout2
ENDIF
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLE',1,zhook_handle)

END SUBROUTINE sple
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!11.
SUBROUTINE splb2e1(NORD,M,G,C,XE,IDX,IDX1,IDY,IDY1,ZE)
!SPLB2E

IMPLICIT NONE

INTEGER, INTENT(IN) :: nord
INTEGER, INTENT(IN) :: m
REAL, DIMENSION(:,:),INTENT(IN) :: g
REAL, DIMENSION(:,:,:), INTENT(IN) :: c
REAL, DIMENSION(:), INTENT(IN) :: xe
INTEGER, INTENT(IN) :: idx,idx1,idy,idy1
REAL, INTENT(OUT) :: ze

INTEGER :: nn,nm
REAL :: rx,ry,alphax,alphay,sum_ax,sum_ay,sum_a
REAL :: zxy, zx1y, zxy1, zx1y1
REAL(KIND=JPRB) :: zhook_handle


IF (lhook) CALL dr_hook('MODE_SPLINES:SPLB2E1',0,zhook_handle)

rx=0
ry=0
IF (idx1.NE.0) rx=(xe(1)-ai(1,1,idx1))/(ai(1,2,idx)-ai(1,1,idx1))
IF (idy1.NE.0) ry=(xe(2)-ai(2,1,idy1))/(ai(2,2,idy)-ai(2,1,idy1))
alphax=(rx-1.)*(rx-1.)*(2*rx+1.)
alphay=(ry-1.)*(ry-1.)*(2*ry+1.)
sum_ax=0.
sum_ay=0.
sum_a=0.
zxy=0.
zx1y=0.
zxy1=0.
zx1y1=0.

nm=0

IF (ns(idx,idy).GE.m) THEN
  nm=1
  sum_ax=sum_ax+alphax
  sum_ay=sum_ay+alphay
  sum_a=sum_a+alphax*alphay
  nn=ns(idx,idy)
  !N=NN
  CALL sple(nord,xs(:,1:nn,idx,idy),g,c(1:nn+m,idx,idy),xe,zxy)
ENDIF

IF (idx1.NE.0) THEN
  IF (ns(idx1,idy).GE.m) THEN
    nm=1
    sum_ax=sum_ax+1.-alphax
    sum_a=sum_a+(1-alphax)*alphay
    nn=ns(idx1,idy)
    CALL sple(nord,xs(:,1:nn,idx1,idy),g,c(1:nn+m,idx1,idy),xe,zx1y)
  ENDIF
ENDIF

IF (idy1.NE.0) THEN
  IF (ns(idx,idy1).GE.m) THEN
    nm=1
    sum_ay=sum_ay+1.-alphay
    sum_a=sum_a+alphax*(1-alphay)
    nn=ns(idx,idy1)
    CALL sple(nord,xs(:,1:nn,idx,idy1),g,c(1:nn+m,idx,idy1),xe,zxy1)
  ENDIF
ENDIF

IF (idx1.NE.0 .AND. idy1.NE.0) THEN
  IF (ns(idx1,idy1).GE.m) THEN
    nm=1
    sum_a=sum_a+(1-alphax)*(1-alphay)
    nn=ns(idx1,idy1)
    CALL sple(nord,xs(:,1:nn,idx1,idy1),g,c(1:nn+m,idx1,idy1),xe,zx1y1)
  ENDIF
ENDIF

IF (nm==1) THEN
  IF (idx1.NE.0 .AND. idy1.EQ.0) THEN
    ze=(alphax*zxy+(1-alphax)*zx1y)/sum_ax
  ELSEIF (idx1.EQ.0 .AND. idy1.NE.0) THEN
    ze=(alphay*zxy+(1-alphay)*zxy1)/sum_ay
  ELSEIF (idx1.NE.0 .AND. idy1.NE.0) THEN
    ze=(alphax*alphay    *zxy  + (1-alphax)*    alphay*zx1y + &
         alphax*(1-alphay)*zxy1 + (1-alphax)*(1-alphay)*zx1y1) / sum_a 
  ELSE
    ze=zxy
  ENDIF
ENDIF
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLB2E1',1,zhook_handle)

END SUBROUTINE splb2e1
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  APPELEE PAR INTERPOL_SPLINES.F90
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!12. 
SUBROUTINE splb2e(NORD,M,NSDI,NSDJ,G,C,XE,ZE)

IMPLICIT NONE

INTEGER, INTENT(IN) :: nord
INTEGER, INTENT(IN) :: m
INTEGER, INTENT(IN) :: nsdi,nsdj
REAL, DIMENSION(:,:), INTENT(IN) :: g
REAL,DIMENSION(:,:,:), INTENT(IN) :: c
REAL, DIMENSION(:,:), INTENT(IN) :: xe
REAL, DIMENSION(:), INTENT(OUT) :: ze

INTEGER, DIMENSION(2) :: nsd
INTEGER :: j, id, ne
INTEGER :: idx, idx1, idy, idy1
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLB2E',0,zhook_handle)
ne=SIZE(xe,2)

nsd(1)=nsdi
nsd(2)=nsdj

DO j=1,ne
  id=1
  CALL splbfin(id,nsd,xe(:,j),ai,idx,idx1)
  id=2
  CALL splbfin(id,nsd,xe(:,j),ai,idy,idy1)
  CALL splb2e1(nord,m,g,c,xe(:,j),idx,idx1,idy,idy1,ze(j))
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLB2E',1,zhook_handle)

END SUBROUTINE splb2e

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!         PARTIE "C"
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

! Calculs EISRS1
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!13. 
SUBROUTINE tred2(A,D,E,Z)
!EISRS1
IMPLICIT NONE

REAL, DIMENSION(:,:), INTENT(IN) :: a
REAL, DIMENSION(:), INTENT(OUT)   :: d
REAL, DIMENSION(:), INTENT(OUT)  :: e
REAL, DIMENSION(:,:), INTENT(OUT):: z

INTEGER :: n, i, ii, j, k, l, jp1
REAL :: f, g, h, scale, hh
REAL(KIND=JPRB) :: zhook_handle


IF (lhook) CALL dr_hook('MODE_SPLINES:TRED2',0,zhook_handle)
n=SIZE(a,1)

!Z=A sur la partie gauceh jusqu'à la diagonale
DO i=1,n
  DO j=1,i
    z(i,j)=a(i,j)
  ENDDO
ENDDO

IF (n.NE.1) THEN
  DO ii=2,n
    !I varie de N à 2
    i=n+2-ii
    !L varie de N-1 à 1
    l=i-1
    h=0.0
    scale=0.0
    IF (l.GE.2) THEN
      !scale = somme des éléments de la ligne à l'exclusion du diagonal
      DO k=1,l
        scale=scale+abs(z(i,k))
      ENDDO
      IF (scale.NE.0.0) THEN
        !normalisation de Z par la valeur de la ligne
        DO k=1,l
          z(i,k)=z(i,k)/scale
          !H = somme des carrés de Z / (somme des ABS(Z)) au carré
          h=h+z(i,k)**2
        ENDDO
        f=z(i,l) 
        !SQRT(H) signé par lélément diagonal-1 de la ligne
        g=-sign(sqrt(h),f)
        !moyenne des éléments au carrés de la ligne, signée
        e(i)=scale*g
        h=h-f*g
        z(i,l)=f-g
        f=0.0
        DO j=1,l
          z(j,i)=z(i,j)/(scale*h)
          g=0.0
          DO k=1,j
            g=g+z(j,k)*z(i,k)
          ENDDO
          jp1=j+1
          IF (l.GE.jp1) THEN
            DO k=jp1,l
              g=g+z(k,j)*z(i,k)
            ENDDO
          ENDIF
          e(j)=g/h
          f=f+e(j)*z(i,j)
        ENDDO
        hh=f/(h+h)
        DO j=1,l
          f=z(i,j)
          g=e(j)-hh*f
          e(j)=g
          DO k=1,j
            z(j,k)=z(j,k)-f*e(k)-g*z(i,k)
          ENDDO
        ENDDO
        DO k=1,l
          z(i,k)=scale*z(i,k)
        ENDDO
      ELSE
        e(i)=z(i,l)
      ENDIF
    ELSE
      e(i)=z(i,l)
    ENDIF
    d(i)=h
  ENDDO
ENDIF

d(1)=0.0
e(1)=0.0
DO i=1,n
  l=i-1
  IF (d(i).NE.0.0) THEN
    DO j=1,l
      g=0.0
      DO k=1,l
        g=g+z(i,k)*z(k,j)
      ENDDO
      DO k=1,l
        z(k,j)=z(k,j)-g*z(k,i)
      ENDDO
    ENDDO
  ENDIF
  d(i)=z(i,i)
  z(i,i)=1.0
  IF (l.GE.1) THEN
    DO j=1,l
      z(i,j)=0.0
      z(j,i)=0.0
    ENDDO
  ENDIF
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:TRED2',1,zhook_handle)

END SUBROUTINE tred2
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!14. 
SUBROUTINE tql2_2(D,E,Z,IERR)
!EISRS1
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(INOUT) :: d
REAL, DIMENSION(:), INTENT(INOUT) :: e
REAL, DIMENSION(:,:), INTENT(OUT) :: z
INTEGER, INTENT(OUT) :: ierr

INTEGER :: i, ii, j, k, l, m, n, mml
REAL :: b, c, f, g, h, p, r, s
REAL :: machep
REAL(KIND=JPRB) :: zhook_handle


IF (lhook) CALL dr_hook('MODE_SPLINES:TQL2_2',0,zhook_handle)
machep=2.**(-52)
ierr=0
n=SIZE(d)

IF (n.NE.1) THEN
  DO i=2,n
    e(i-1)=e(i)
  ENDDO
  f=0.0
  b=0.0
  e(n)=0.0
  DO l=1,n
    j=0
    h=machep*(abs(d(l))+abs(e(l)))
    IF (b.LT.h) b=h
    DO m=l,n
      IF (abs(e(m)).LE.b) EXIT
    ENDDO
    IF (m.NE.l) THEN
      DO WHILE (abs(e(l)).GT.b) 
        IF (j.EQ.30) THEN
          ierr=l
          IF (lhook) CALL dr_hook('MODE_SPLINES:TQL2_2',1,zhook_handle)
          RETURN
        ENDIF
        j=j+1
        p=(d(l+1)-d(l)) / (2.0 * e(l))
        r=sqrt(p*p+1.0)
        h=d(l)-e(l) / (p+sign(r,p))
        DO i=l,n
          d(i)=d(i)-h
        ENDDO
        f=f+h
        p=d(m)
        c=1.0
        s=0.0
        mml=m-l
        DO ii=1,mml
          i=m-ii
          g=c*e(i)
          h=c*p
          IF (abs(p).GE.abs(e(i))) THEN
            c=e(i)/p
            r=sqrt(c*c+1.0)
            e(i+1)=s*p*r
            s=c/r
            c=1.0/r
          ELSE
            c=p/e(i)
            r=sqrt(c*c+1.0)
            e(i+1)=s*e(i)*r
            s=1.0/r
            c=c*s
          ENDIF
          p=c*d(i)-s*g
          d(i+1)=h+s*(c*g+s*d(i))
          DO k=1,n
            h=z(k,i+1)
            z(k,i+1)=s*z(k,i)+c*h
            z(k,i)=c*z(k,i)-s*h
          ENDDO
        ENDDO
        e(l)=s*p
        d(l)=c*p
      ENDDO
    ENDIF
    d(l)=d(l)+f
  ENDDO
  DO ii=2,n
    i=ii-1
    k=i
    p=d(i)
    DO j=ii,n
      IF (d(j).LT.p) THEN
        k=j
        p=d(j)
      ENDIF
    ENDDO
    IF (k.NE.i) THEN
      d(k)=d(i)
      d(i)=p
      DO j=1,n
        p=z(j,i)
        z(j,i)=z(j,k)
        z(j,k)=p
      ENDDO
    ENDIF
  ENDDO
ENDIF
IF (lhook) CALL dr_hook('MODE_SPLINES:TQL2_2',1,zhook_handle)

END SUBROUTINE tql2_2
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!15. SPLR,SPLU
SUBROUTINE eisrs1(AR,WR,ZR,WORK,IERR)
!
!     all eigenvalues and corresponding eigenvectors of a real
!     symmetric matrix
!
IMPLICIT NONE
!toutes les dimensions sont N: ok
REAL, DIMENSION(:,:), INTENT(IN) :: ar
REAL, DIMENSION(:), INTENT(OUT) :: wr
REAL, DIMENSION(:,:), INTENT(OUT) :: zr
REAL, DIMENSION(:), INTENT(OUT) :: work
INTEGER, INTENT(OUT) :: ierr
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:EISRS1',0,zhook_handle)
 CALL tred2(ar,wr,work,zr)
 CALL tql2_2(wr,work,zr,ierr)
IF (lhook) CALL dr_hook('MODE_SPLINES:EISRS1',1,zhook_handle)

END SUBROUTINE eisrs1

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!        GROUPES DES PROCEDURES V
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!16.
SUBROUTINE splv(B,W,N,P,RES)
!SP0NOP
IMPLICIT NONE

REAL,DIMENSION(:),INTENT(IN) :: b
REAL,DIMENSION(:), INTENT(IN) :: w
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: p
REAL, INTENT(OUT) :: res
!
INTEGER :: a, d, i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLV',0,zhook_handle)
a=0.
res=0.
DO i=1,SIZE(b)
  d=b(i)+n*p
  a=a+1./d
  res=res+w(i)**2 / d**2
ENDDO
res=n*res/a**2
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLV',1,zhook_handle)

END SUBROUTINE splv
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!              GROUPE DES PROCEDURES P
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!PROCEDURES INTERMEDIAIRES POUR SPLPS2V
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!17. S2VI=somme(W(i)**2)/size(W)
SUBROUTINE spls2vi(W,RES)
!
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: w
REAL, INTENT(OUT) :: res

INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLS2VI',0,zhook_handle)
res = 0.
DO i=1,SIZE(w)
  res = res + w(i)**2
ENDDO
res=res/SIZE(w)
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLS2VI',1,zhook_handle)

END SUBROUTINE spls2vi
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!18. calcul de BN en fonction de W, B, N et P
SUBROUTINE splbvm(B,W,N,P,RES)
!SPLS2V, SP0VPQ
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: b, w
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: p
REAL, INTENT(OUT) :: res
!
REAL :: a
INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBVM',0,zhook_handle)
res=0.
a=0.
DO i=1,SIZE(b)
  res=res+(w(i)**2)/(b(i)+n*p)**2
  a=a+1./(b(i)+n*p)
ENDDO
res=res/a
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBVM',1,zhook_handle)

END SUBROUTINE splbvm
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!19. 
SUBROUTINE spls2v(BI,WI,N,P,RES)
!SPLPS2V, SPLP
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: bi !dimension N-M
REAL, DIMENSION(:), INTENT(IN) :: wi !dimension N-M
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: p
REAL, INTENT(OUT) :: res
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLS2V',0,zhook_handle)
 CALL splbvm(bi,wi,n,p,res)
res=n*p*res
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLS2V',1,zhook_handle)

END SUBROUTINE spls2v
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!20.
SUBROUTINE splds2v(B,W,N,P,RES)
!SPLPS2V
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: b
REAL, DIMENSION(:), INTENT(IN) :: w
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: p
REAL, INTENT(OUT) :: res

REAL :: da, db, a1, a2, a3, b2, c1, c2
INTEGER :: j, i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLDS2V',0,zhook_handle)
da=0.
db=0.
DO i=1,SIZE(b)
  a1=1./(b(i)+n*p)
  a2=a1**2
  a3=a1**3
  da=da+a1
  b2=w(j)**2
  DO j=1,SIZE(b)
    c1=1./(b(j)+n*p)
    c2=c1**2
    db=db+b2*(a2*c1+n*p*(-2*a3*c1+a2*c2))
  ENDDO
ENDDO
res=(n/(da**2))*db
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLDS2V',1,zhook_handle)

END SUBROUTINE splds2v
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!21. SPLPS2V procedure 
SUBROUTINE splps2v(B,W,N,P0,S2,P,IREP)
!SPLP
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: b
REAL, DIMENSION(:), INTENT(IN) :: w
INTEGER, INTENT(IN) :: n
REAL, INTENT(INOUT) :: p0
REAL, INTENT(IN) :: s2
REAL, INTENT(OUT) :: p
INTEGER, INTENT(OUT) :: irep

INTEGER, PARAMETER :: eps=1.e-2
REAL :: p1, rinf, rs, drs
INTEGER :: niter
REAL(KIND=JPRB) :: zhook_handle

! Calcul de la valeur du s2 de Wahba pour p=infini
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPS2V',0,zhook_handle)
 CALL spls2vi(w,rinf)
IF (rinf.LE.s2) THEN
  irep=-6
  p0=0.
  WRITE(*,fmt='(A53,G15.5)') &
    "SPLPS2V : S2 > VALEUR DE LA FONCTION POUR P INFINI = ",rinf
  p=(b(SIZE(b))**2)/b(1)
  IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPS2V',1,zhook_handle)
  RETURN
ENDIF

p1=0.
! Iterations de Newton
DO niter=1,nmax
  p0=p1
  CALL spls2v(b,w,n,p0,rs)
  CALL splds2v(b,w,n,p0,drs)

  p1=p0+(s2-rs)/drs

  IF (abs(p1-p0)/p0.LT.eps) THEN
    IF (p1.GE.0.) THEN
      irep=0
      p=p1
      p0=p1
    ELSE
      irep=-3
      p=0.
      p0=0.
      WRITE(*,fmt='(A34,G15.5)') &
        "SPLPS2V : SOLUTION NEGATIVE : P = ",p
    ENDIF
    IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPS2V',1,zhook_handle)
    RETURN
  ENDIF
ENDDO

irep=-5
p=p1
p0=0.
WRITE(*,fmt='(A48,A23,G15.5)') &
  "SPLPS2V : NOMBRE MAXIMAL D ITERATIONS ATTEINT : ",&
  "VALEUR DE P ATTEINTE : ",p1
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPS2V',1,zhook_handle)

END SUBROUTINE splps2v
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!PROCEDURES INTERMEDIAIRES POUR SPLPR
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!22. Comme SPLS2VI mais on ne divise pas par SIZE(W) mais par N
SUBROUTINE splri(W,N,RES)
!SPLPR
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: w
INTEGER, INTENT(IN) :: n
REAL, INTENT(OUT) :: res

INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLRI',0,zhook_handle)
res=0.
DO i=1,SIZE(w)
  res = res + w(i)**2
ENDDO
res=res/n
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLRI',1,zhook_handle)

END SUBROUTINE splri
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!23. RS=SQRT(S2/(SOMME((W(I)/B(I))**2)*N)
SUBROUTINE splpr0(B,W,N,S2,P)
!SPLPR
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: b
REAL, DIMENSION(:), INTENT(IN) :: w
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: s2
REAL, INTENT(OUT) :: p

REAL :: res
INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPR0',0,zhook_handle)
res=0.
DO i=1,SIZE(w)
  res=res+(w(i)/b(i))**2
ENDDO
p=sqrt(s2/(res*n))
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPR0',1,zhook_handle)

END SUBROUTINE splpr0   
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!24. RS=N*SOMME((W(I)/(B(I)+N*P))**2)*P**2
SUBROUTINE splrs(B,W,N,P,RES)
!SPLP, SPLPR, SP0CVQ 
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: b
REAL, DIMENSION(:), INTENT(IN) :: w
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: p
REAL, INTENT(OUT):: res

INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLRS',0,zhook_handle)
res=0.
DO i=1,SIZE(b)
  res=res+(w(i)/(b(i)+n*p))**2
ENDDO
res=n*res*p**2
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLRS',1,zhook_handle)

END SUBROUTINE splrs
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!25. DRS=2*N*P*(SOMME(W(I)**2*B(I)/(B(I)+N*P)**3)
SUBROUTINE spldrs(B,W,N,P,RES)
!SPLPR
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: b
REAL, DIMENSION(:), INTENT(IN) :: w
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: p
REAL, INTENT(OUT) :: res
!
INTEGER :: j
REAL :: d
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLDRS',0,zhook_handle)
res=0.
DO j=1,SIZE(b)
  d=b(j)+n*p
  res=res+w(j)**2*b(j)/(d**3)
ENDDO
res=2.*n*p*res
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLDRS',1,zhook_handle)

END SUBROUTINE spldrs   
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!26. SPLPR procedure: parallèle de splps2v
SUBROUTINE splpr(N,W,B,P0,S2,P,IREP)
!SPLP
IMPLICIT NONE

INTEGER, INTENT(IN) :: n
REAL, DIMENSION(:), INTENT(IN) :: w !N-M
REAL, DIMENSION(:), INTENT(IN) :: b !N-M
REAL, INTENT(INOUT) :: p0
REAL, INTENT(INOUT) :: s2
REAL, INTENT(INOUT) :: p
INTEGER, INTENT(OUT) :: irep

INTEGER, PARAMETER :: eps=1.e-2
REAL :: rs, drs, p1, rinf
INTEGER :: niter
REAL(KIND=JPRB) :: zhook_handle

! Calcul de la valeur de la fonction de Reinsch pour p=infini

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPR',0,zhook_handle)
 CALL splri(w,n,rinf)
IF (rinf.LE.s2) THEN
  irep=-6
  p0=0.
  WRITE(*,fmt='(A51,G15.5)') &
    "SPLPR : S2 > VALEUR DE LA FONCTION POUR P INFINI = ",rinf
  p=(b(SIZE(b))**2)/b(1)
  IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPR',1,zhook_handle)
  RETURN
ENDIF

! Recherche du point de depart de l'iteration
IF (p0.EQ.0.) CALL splpr0(b,w,n,s2,p0)
p1=p0

! Iterations de Newton
DO niter=1,nsdmax
  p0=p1
  CALL splrs(b,w,n,p0,rs)
  CALL spldrs(b,w,n,p0,drs)

  p1=p0+(s2-rs)/drs

  IF (abs(p1-p0)/p0.LT.eps) THEN
    IF (p1.GE.0.) THEN
      irep=0
      p=p1
      p0=p1
    ELSE
      irep=-3
      p=0.
      p0=0.
      WRITE(*,fmt='(A32,G15.5)') &
        "SPLPR : SOLUTION NEGATIVE : P = ",p
    ENDIF
    IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPR',1,zhook_handle)
    RETURN
  ENDIF
ENDDO

irep=-5
p=0.
p0=0.
WRITE(*,fmt='(A45,A23,G15.5)') &
  "SPLPR : NOMBRE MAXIMAL D ITERATIONS ATTEINT: ",&
  "VALEUR DE P ATTEINTE : ",p1
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPR',1,zhook_handle)

END SUBROUTINE splpr
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!PROCEDURES INTERMEDIAIRES POUR SPLPV
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!27.
SUBROUTINE spldv(B,W,N,P,RES)
!SPLPV
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: b
REAL, DIMENSION(:), INTENT(IN) :: w
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: p
REAL, INTENT(OUT) :: res
!
REAL :: s1, s2, sb2, sb3, a1, a2, b2
INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLDV',0,zhook_handle)
s1=0.
s2=0.
sb2=0.
sb3=0.
DO i=1,SIZE(b)
  a1=1./(b(i)+n*p)
  a2=a1**2
  b2=a2*w(i)**2
  s1=s1+a1
  s2=s2+a2
  sb2=sb2+b2
  sb3=sb3+b2*a1
ENDDO
res=s2*sb2-s1*sb3
res=2.*n**2*res/(s1**3)
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLDV',1,zhook_handle)

END SUBROUTINE spldv
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!28.
SUBROUTINE spld2v(B,W,N,P,RES)
!SPLPV
IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN) :: b
REAL, DIMENSION(:), INTENT(IN) :: w
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: p
REAL, INTENT(OUT) :: res

REAL :: s1, s2, s3, sb2, sb3, sb4
REAL :: a1, a2, b2
INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLD2V',0,zhook_handle)
s1=0.
s2=0.
s3=0.
sb2=0.
sb3=0.
sb4=0.
DO i=1,SIZE(b)
  a1=1./(b(i)+n*p)
  a2=a1**2
  b2=a2*w(i)**2
  s1=s1+a1
  s2=s2+a2
  s3=s3+a2*a1
  sb2=sb2+b2
  sb3=sb3+b2*a1
  sb4=sb4+b2*a2
ENDDO
res=3.*s2**2*sb2-4.*s2*s1*sb3+3.*s1**2*sb4-2.*s3*s1*sb2
res=res/(s1**4)
res=res*2.*n**3
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLD2V',1,zhook_handle)

END SUBROUTINE spld2v
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!29.
SUBROUTINE splpv(N,B,W,P0,P,IREP)
!SPLP
IMPLICIT NONE

INTEGER, INTENT(IN) :: n
REAL, DIMENSION(:), INTENT(IN):: b, w !N-M
REAL, INTENT(INOUT) :: p0
REAL, INTENT(OUT) :: p
INTEGER, INTENT(OUT) :: irep

INTEGER, PARAMETER:: eps=xsurf_tiny, epsr=1.e-2
REAL :: dvm, d2vm, p1
INTEGER :: iflag, niter
REAL(KIND=JPRB) :: zhook_handle


! Calcul de Vm'(0) . Si Vm'(0) > 0 alors popt = 0.
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPV',0,zhook_handle)
 CALL spldv(b,w,n,0.,dvm)

IF (dvm.GT.0.) THEN
  irep=0.
  p=0.
  IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPV',1,zhook_handle)
  RETURN
ENDIF

! Recherche du point de depart p0 de l'iteration de Newton
iflag=0

DO WHILE(iflag==0)

  IF (p0.EQ.0.) THEN

    iflag=1
    p0=b(1)/n

    !on augmente P0 et on calcul D2VM; quand D2VM est supérieur à 0, 
    !on passe à la suite
    DO niter=1,nsdmax
      CALL spld2v(b,w,n,p0,d2vm)
      IF (d2vm.GT.0.) EXIT
      p0=p0*10.
    ENDDO

    IF (d2vm.LE.0.) THEN
      irep=-1
      p0=0.
      p=0.
      WRITE(*,fmt='(A41)') &
        "SPLPV : PAS DE MINIMUM: D2VM < 0 PARTOUT "
      IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPV',1,zhook_handle)
      RETURN
    ENDIF

  ENDIF

  p1=p0

  DO niter=1,nsdmax

    p0=p1
    CALL spldv(b,w,n,p0,dvm)
    CALL spld2v(b,w,n,p0,d2vm)

    !Si D2Vm est très proche de 0, on sort et P=P0
    IF (abs(d2vm).LT.eps) THEN
      p0=0.
      IF (iflag.EQ.0) EXIT
      irep=-2
      p=p0
      WRITE(*,fmt='(A28)') "SPLPV : PASSAGE PAR D2VM = 0"
      IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPV',1,zhook_handle)
      RETURN
    ENDIF

    p1=p0-dvm/d2vm

     !si l'écart entre P0 et P1 est suffisamment petit
    IF (abs(p1-p0)/p0.LT.epsr) THEN
       !Si P1 est négatif, on repart dans la boucle
      IF (p1.LT.0.) THEN
        p0=0.
        IF (iflag.EQ.0) EXIT
        irep=-3
        p=0.
        WRITE(*,fmt='(A33,G15.5)') &
          "SPLPV : SOLUTION NEGATIVE : P = ",p1
      ELSE
        !si P1 et positif et D2VM positif, on garde P1
        IF (d2vm.GE.0.) THEN
          irep=0
          p=p1
          p0=p1
        ELSE
          !si P1 est négatif, on repart dans la boucle
          p0=0.
          IF (iflag.EQ.0) EXIT
          irep=-4
          p=0.
          WRITE(*,fmt='(A28,G15.5)') &
            "SPLPV : MAXIMUM DE VM : P = ",p1
        ENDIF
      ENDIF
      IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPV',1,zhook_handle)
      RETURN
    ENDIF
  ENDDO

  p0=0.
  IF (iflag.NE.0) THEN
    irep=-5
    p=p1
    WRITE(*,fmt='(A46,A23,G15.5)') &
      "SPLPV : NOMBRE MAXIMAL D ITERATIONS ATTEINT : ",&
      "VALEUR DE P ATTEINTE : ",p1
    IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPV',1,zhook_handle)
    RETURN
  ENDIF

ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLPV',1,zhook_handle)

END SUBROUTINE splpv
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!30.
SUBROUTINE splp(LISSAGE,N,B,W,P,S2,IREP)
!SP0NOP
IMPLICIT NONE

INTEGER, INTENT(IN) :: lissage
INTEGER, INTENT(IN) :: n
REAL, DIMENSION(:), INTENT(IN) :: b
REAL, DIMENSION(:), INTENT(IN) :: w
REAL, INTENT(INOUT) :: p
REAL, INTENT(INOUT) :: s2
INTEGER, INTENT(OUT) :: irep

REAL :: p0
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLP',0,zhook_handle)
IF (lissage.EQ.1) THEN
!  Recherche de la valeur de P pour S2 connu, par la methode de Wahba
  p0=0.
  CALL splps2v(b,w,n,p0,s2,p,irep)
  IF (irep.NE.0) THEN
    irep=-2
    IF (lhook) CALL dr_hook('MODE_SPLINES:SPLP',1,zhook_handle)
    RETURN
  ENDIF
ELSEIF (lissage.EQ.10) THEN
!  Recherche de la valeur de P pour S2 connu, par la methode de Reinsch
  p0=0.
  CALL splpr(n,w,b,p0,s2,p,irep)
  IF (irep.NE.0) THEN
    irep=-2
    IF (lhook) CALL dr_hook('MODE_SPLINES:SPLP',1,zhook_handle)
    RETURN
  ENDIF
ELSEIF (lissage.EQ.2) THEN
!  Recherche de la valeur de P par la methode de G.C.V
  p0=0.
  CALL splpv(n,b,w,p0,p,irep)
  IF (irep.NE.0) irep=-3
!  Calcul de la variance d'erreur de mesure correspondante
  CALL spls2v(b,w,n,p,s2)
ELSEIF (lissage.EQ.3) THEN
!  Calcul de l'ajustement correspondant a une valeur de P donnee
  CALL splrs(b,w,n,p,s2)
ENDIF
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLP',1,zhook_handle)

END SUBROUTINE splp
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!          GROUPE DES PROCEDURES T et R
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!31. TT procedure
SUBROUTINE spltt(X,DS,IW,TM,TT,IREP)
!SPLR
IMPLICIT NONE

REAL,DIMENSION(:,:),INTENT(IN) :: x !ND,N
REAL,DIMENSION(:),INTENT(IN) :: ds !N
INTEGER,DIMENSION(:,:),INTENT(IN) :: iw !ND,M
REAL,DIMENSION(:,:),INTENT(OUT) :: tm !N,M
REAL,DIMENSION(:,:),INTENT(OUT) :: tt !M,M
INTEGER, INTENT(OUT) :: irep

REAL(KIND=JPRB) :: zhook_handle

! Calcul de la matrice T des monomes
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLTT',0,zhook_handle)
 CALL splt(x,iw,tm)
! Calcul de Tm=(Ds-1)*T
 CALL mxidml(ds,tm,tm)
! Calcul de Tm'*Tm
 CALL mxmspl(transpose(tm),tm,tt)
! Calcul de (Tm'*Tm)-1
 CALL smxinv(tt,irep)
! Test du compte-rendu de smxinv
IF (irep.NE.0) WRITE(*,fmt='(A27)') "SPLTT: MATRICE TT NON INVERSIBLE"
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLTT',1,zhook_handle)

END SUBROUTINE spltt
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!32. r procedure
SUBROUTINE splr_1(NORD,X,DS,IW,TTT,R,IREP)
!SP0NOP, SP0CVQ
IMPLICIT NONE

INTEGER,INTENT(IN) :: nord
REAL,DIMENSION(:,:), INTENT(IN) :: x !ND,N
REAL,DIMENSION(:),INTENT(IN) :: ds !N
INTEGER,DIMENSION(:),INTENT(OUT) :: iw !ND
REAL,DIMENSION(:),INTENT(OUT) :: ttt !N
REAL,DIMENSION(:,:),INTENT(OUT) :: r !N,N
INTEGER, INTENT(OUT)::irep

INTEGER :: n, i, j
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLR_1',0,zhook_handle)
n=SIZE(x,2)

! R = I
DO i=1,n
  DO j=1,n
    IF (i.EQ.j) THEN
      r(i,i)=1.
    ELSE
      r(i,j)=0.
    ENDIF
  ENDDO
ENDDO

irep=0

!  Calcul de (Ds-1)*R
 CALL mxidml(ds,r,r)

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLR_1',1,zhook_handle)

END SUBROUTINE splr_1
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!32.
SUBROUTINE splr_2(NORD,X,DS,IW,TTT,R,IREP)

IMPLICIT NONE

INTEGER,INTENT(IN) :: nord
REAL,DIMENSION(:,:), INTENT(IN) :: x !ND,N
REAL,DIMENSION(:),INTENT(IN) :: ds !N
INTEGER,DIMENSION(:,:),INTENT(OUT) :: iw !ND,M
REAL,DIMENSION(:,:),INTENT(OUT) :: ttt !M,N
REAL,DIMENSION(:,:),INTENT(OUT) :: r !N,N-M
INTEGER, INTENT(OUT)::irep

REAL,DIMENSION(SIZE(X,2),SIZE(X,2)) :: c !N,N
REAL,DIMENSION(SIZE(X,2),SIZE(IW,2))  :: tm !N,M
REAL,DIMENSION(SIZE(IW,2),SIZE(IW,2)) :: tt !M,M
REAL, DIMENSION(SIZE(X,2),SIZE(X,2)) :: r1
REAL,DIMENSION(SIZE(X,2)) :: vp !N
REAL,DIMENSION(SIZE(X,2)) :: work !N
INTEGER :: m, n, j
INTEGER, DIMENSION(1) :: imin
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLR_2',0,zhook_handle)
n=SIZE(x,2)
m=SIZE(iw,2)

!     Recherche des monomes de degre inferieur ou egal a nord-1
 CALL splie(nord-1,iw)
!----------------------------------------------------------
!     Calcul de la matrice (Tm'*Tm)-1
 CALL spltt(x,ds,iw,tm,tt,irep)
!     Calcul de (T'm*Tm)-1 Tm'
 CALL mxmspl(tt,transpose(tm),ttt)
!R1(:,1:M)=TRANSPOSE(TTT)
!------------------------------------------------------------
!     Calcul de Tm ((Tm'*Tm)-1) Tm'
 CALL mtxaxm(transpose(tm),tt,c)
DO j=1,n
  c(j,j)=c(j,j)-1.
ENDDO
 c(:,:)=-c(:,:)
!-----------------------------------------------------------
!     Calcul des vecteurs propres de I - Tm c ((Tm'*Tm)-1) c Tm'
 CALL eisrs1(c,vp,r1,work,irep)

!     Test du signe des valeurs propres de C
imin=minloc(vp(m+1:n))
IF(vp(m+imin(1)).LE.0) THEN
  irep=-1
  WRITE(*,fmt='(A31)') "SPLR: VALEUR PROPRE < 0 DE PROJ"
  IF (lhook) CALL dr_hook('MODE_SPLINES:SPLR_2',1,zhook_handle)
  RETURN
ENDIF
!     Compte-rendu OK
irep=0
!  Calcul de (Ds-1)*R
r=r1(:,m+1:n)
 CALL mxidml(ds,r,r)
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLR_2',1,zhook_handle)

END SUBROUTINE splr_2
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!             GROUPE DES PROCEDURES U
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!33. U procedure
SUBROUTINE splu(NORD,X,G,R,RK,U,DB)
!SP0NOP, SP0CVQ
IMPLICIT NONE

INTEGER, INTENT(IN) :: nord
REAL,DIMENSION(:,:), INTENT(IN) :: x !ND,N
REAL,DIMENSION(:,:), INTENT(IN) :: g ! ND,ND
REAL,DIMENSION(:,:), INTENT(IN) :: r !N,N-M
REAL,DIMENSION(:,:), INTENT(OUT) :: rk !N,N
REAL,DIMENSION(:,:), INTENT(OUT) :: u !N-M,N-M
REAL,DIMENSION(:), INTENT(OUT) :: db !N-M

REAL,DIMENSION(SIZE(R,2),SIZE(R,2)) :: rkr !N-M,N-M
REAL,DIMENSION(SIZE(R,2)) :: work !N-M

INTEGER :: m, n, irep
INTEGER, DIMENSION(1) :: imin
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLU',0,zhook_handle)
n=SIZE(x,2)
m=SIZE(x,2)-SIZE(r,2)

! Calcul de la matrice K des noyaux
 CALL splk(nord,x,x,g,rk)
! Calcul de R'*K*R
 CALL mtxaxm(r,rk,rkr)
! Calcul de la matrice U des vecteurs propres de R'*K*R
 CALL eisrs1(rkr,db,u,work,irep)
! Test du signe des valeurs propres de R'KR
imin=minloc(db)
IF (db(imin(1)).LE.0)  &
   WRITE(*,fmt='(A31)') "SPLU : VALEUR PROPRE < 0 DE RKR"
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLU',1,zhook_handle)

END SUBROUTINE splu     
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!               GROUPE DES PROCEDURES W
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!34. w procedure
SUBROUTINE splw(Z,R,U,W,RZ)
!SP0NOP, SP0CVQ
IMPLICIT NONE

REAL,DIMENSION(:), INTENT(IN) :: z !N
REAL, DIMENSION(:,:),INTENT(IN) :: r !N,N-M
REAL,DIMENSION(:,:),INTENT(IN) :: u !N-M,N-M
REAL, DIMENSION(:),INTENT(OUT) :: w !N-M
REAL,DIMENSION(:),INTENT(OUT) :: rz !N-M

INTEGER :: n, m
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLW',0,zhook_handle)
n=SIZE(z)
m=SIZE(z)-SIZE(r,2)

! Calcul de R'*Z
 CALL mxmspl(transpose(r),z,rz)
! Calcul de U'R'Z
 CALL mxmspl(transpose(u),rz,w)
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLW',1,zhook_handle)

END SUBROUTINE splw
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!              GROUPE DES PROCEDURES Q
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!35. splvpQ procedure
SUBROUTINE splvpq(P,R,RK,U,DB,Q)
!SP0CVQ
IMPLICIT NONE

REAL, INTENT(IN) :: p
REAL, DIMENSION(:,:), INTENT(IN) :: r !N,N-M
REAL, DIMENSION(:,:), INTENT(IN) :: rk ! N,N
REAL, DIMENSION(:,:), INTENT(IN) :: u !N-M,N-M
REAL, DIMENSION(:), INTENT(IN) :: db !N-M
REAL, DIMENSION(:,:), INTENT(OUT) :: q ! N,N

REAL, DIMENSION(SIZE(R,1),SIZE(R,1)):: q1 !N,N
REAL, DIMENSION(SIZE(U,1),SIZE(U,2))::qw1,qw2 !N-M,N-M
INTEGER :: n, j
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLVPQ',0,zhook_handle)
n=SIZE(r,1)

! Calcul de (DB+n*p*I)-1
qw1(:,:)=0.
DO j=1,SIZE(db)
  qw1(j,j)=1./(db(j)+n*p)
ENDDO

! Calcul de U((DB+n*p*I)-1)U'
 CALL mtxaxm(transpose(u),qw1,qw2)
! Calcul de RU((DB+n*p*I)-1))UR'
 CALL mtxaxm(transpose(r),qw2,q)
! Calcul de QK
 CALL mxmspl(q,rk,q1)

! Calcul de I - QK
q1(:,:)=-q1(:,:)
DO j=1,n
  q1(j,j)=q1(j,j)+1
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLVPQ',1,zhook_handle)

END SUBROUTINE splvpq
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!              GROUPE DES PROCEDURES C
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!36. 
SUBROUTINE splc(Q,Z,RK,DS,TTT,RKC,C)
!SP0NOP, SP0CVQ
IMPLICIT NONE

REAL, DIMENSION(:,:), INTENT(IN) :: q !Q(N,N)
REAL, DIMENSION(:), INTENT(IN) :: z !Z(N)
REAL, DIMENSION(:,:), INTENT(IN) :: rk !RK(N,n)  
REAL, DIMENSION(:), INTENT(IN) :: ds !DS(N)
REAL, DIMENSION(:,:), INTENT(IN) :: ttt !TTT(M,N)
REAL, DIMENSION(:), INTENT(OUT)  :: rkc !RKC(N)
REAL, DIMENSION(:), INTENT(OUT) :: c !C(N+M)
REAL(KIND=JPRB) :: zhook_handle

INTEGER :: n, m

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLC',0,zhook_handle)
n=SIZE(z)
m=SIZE(ttt,1)

!     Calcul de c
!     -----------
!     Calcul de c = Qz = R((R'KR+n*p*I)-1)R'z
 CALL mxmspl(q,z,c(1:n))

!     Calcul de d
!     -----------
IF (m.NE.0) THEN
! Calcul de K*c
  CALL mxmspl(rk,c(1:n),rkc)
! Calcul de z - K*c
  rkc(:)=z(:)-rkc(:)   
! Calcul de (Ds-1) * (z-K*c)
  rkc(:)=rkc(:)/ds(:)
! Calcul de d = (Tm'*Tm)-1 * Tm' * (Ds-1) * (z-K*c)
  CALL mxmspl(ttt,rkc,c(n+1:n+m))
ENDIF

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLC',1,zhook_handle)

END SUBROUTINE splc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!                  PROCEDURES GLOBALES
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!37. Calcul de M en fonction de NORD et ND (qu'est-ce que c'est que ce ND?)
SUBROUTINE splm(ND,NORD,M)
!SP0NOP
IMPLICIT NONE

INTEGER, INTENT(IN) :: nd
INTEGER, INTENT(IN) :: nord
INTEGER, INTENT(OUT):: m
REAL(KIND=JPRB) :: zhook_handle


IF (lhook) CALL dr_hook('MODE_SPLINES:SPLM',0,zhook_handle)
IF (nd==1) THEN
  m=nord
ELSEIF (nd==2) THEN
  m=nord*(nord+1)/2.
ELSEIF (nd==3) THEN
  m=nord*(nord+1)*(nord+2)/6.
ELSEIF (nd==4) THEN
  m=nord*(nord+1)*(nord+2)*(nord+3)/24.
ENDIF
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLM',1,zhook_handle)

END SUBROUTINE splm
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!38. 
SUBROUTINE sp0nop(X,G,Z,DS,LISSAGE,LORDRE,IOPT,NORDOPT,M,S2,P,&
         iw,ttt,r,rk,u,db,w,rz,irep) 
!SP0CVQ
IMPLICIT NONE

REAL, DIMENSION(:,:), INTENT(IN) :: x !ND,N
REAL, DIMENSION(:,:), INTENT(IN) :: g ! ND,ND
REAL, DIMENSION(:), INTENT(IN) :: z ! N
REAL, DIMENSION(:), INTENT(IN) :: ds ! N
INTEGER, INTENT(IN) :: lissage
INTEGER, INTENT(IN) :: lordre
INTEGER, INTENT(IN) :: iopt
INTEGER, INTENT(INOUT) :: nordopt
INTEGER, INTENT(INOUT) :: m
REAL, INTENT(INOUT) :: s2
REAL, INTENT(INOUT) :: p
INTEGER, DIMENSION(:,:), INTENT(OUT):: iw !ND,M
REAL, DIMENSION(:,:), INTENT(OUT) :: ttt ! M,N
REAL, DIMENSION(:,:), INTENT(OUT) :: r ! N,N-M
REAL,DIMENSION(:,:), INTENT(OUT) :: rk !N,N
REAL,DIMENSION(:,:), INTENT(OUT) :: u !N-M,N-M
REAL,DIMENSION(:), INTENT(OUT) :: db !N-M
REAL, DIMENSION(:),INTENT(OUT) :: w !N-M
REAL,DIMENSION(:),INTENT(OUT) :: rz !N-M
INTEGER, INTENT (OUT) :: irep


REAL, DIMENSION(NORDMAX) :: vm, s2save, psave, msave
INTEGER :: n, nd
INTEGER :: nordmi, nordma, nordm, nord
REAL :: vmmin
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SP0NOP',0,zhook_handle)
n=SIZE(x,2)
nd=SIZE(x,1)

IF (lordre.EQ.0) THEN
!
!  Recherche de p seulement avec nord fixe
!  ---------------------------------------
  IF (iopt.EQ.0) THEN
!   Calcul de IW, TTT, R
    CALL splr(nordopt,x,ds,iw,ttt,r,irep)  
!   Calcul de RK, U, DB    
    CALL splu(nordopt,x,g,r,rk,u,db)
  ENDIF
   !calcul de W, RZ
  CALL splw(z,r,u,w,rz)
  !calcul de P
  CALL splp(lissage,n,db,w,p,s2,irep)

ELSE

!  Recherche de nordopt et popt
!  ----------------------------

  nordmi=nd/2+1

!  Calcul de l'ordre maximal pour le nombre de points n
  CALL splm(nd,nord,m)
  DO WHILE(n.GT.m) 
    nord=nord+1
    CALL splm(nd,nord,m)
  ENDDO
  nordma=min(nord,nordmax)
  nordm=nordma

!  Calcul de Vm pour les differents ordre
  DO nord=nordmi,nordma
    !calcul de M
    CALL splm(nd,nord,m)
    CALL splr(nord,x,ds,iw,ttt,r,irep)
    CALL splu(nord,x,g,r,rk,u,db)    
    CALL splw(z,r,u,w,rz)
    CALL splp(lissage,n,db,w,p,s2,irep)
    msave(nord)=m
    s2save(nord)=s2
    psave(nord)=p
    IF (irep.EQ.0) THEN
      !calcul de VM(NORD)
      CALL splv(db,w,n,p,vm(nord))
      WRITE(*,fmt='(A15,I5,A5,G15.5,A6,G15.5)')&
        "CNORD : NORD = ",nord," P = ",p," VM = ",vm(nord)
    ELSE
      IF (nord.EQ.nordmi) THEN
        nordopt=nord
        IF (lhook) CALL dr_hook('MODE_SPLINES:SP0NOP',1,zhook_handle)
        RETURN
      ELSE
        nordm=nord-1
        irep=0
        EXIT
      ENDIF
    ENDIF
  ENDDO

!  Recherche de l'ordre optimal
  vmmin=xsurf_huge
  DO nord=nordmi,nordm
    IF(vm(nord).LT.vmmin) THEN
      vmmin=vm(nord)
      nordopt=nord
      m=msave(nord)
      s2=s2save(nord)
      p=psave(nord)
    ENDIF
  ENDDO
  WRITE(*,fmt='(A24,I5)') "CNORD : ORDRE OPTIMAL = ",nordopt
ENDIF
IF (lhook) CALL dr_hook('MODE_SPLINES:SP0NOP',1,zhook_handle)

END SUBROUTINE sp0nop
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!39. 
SUBROUTINE sp0cvq(NORD,M,X,G,Z,DS,S2,P,LISSAGE,LORDRE,IOPT,C,IREP)
!SPLB2C
IMPLICIT NONE

INTEGER, INTENT(INOUT) :: nord
INTEGER, INTENT(INOUT) :: m
REAL,DIMENSION(:,:), INTENT(IN) :: x !ND,N
REAL, DIMENSION(:,:), INTENT(IN) :: g !ND,ND
REAL,DIMENSION(:), INTENT(IN) :: z ! N
REAL, DIMENSION(:), INTENT(IN) :: ds !N
REAL, INTENT(INOUT) :: s2
REAL, INTENT(INOUT) :: p
INTEGER, INTENT(IN) :: lissage
INTEGER, INTENT(IN) :: lordre
INTEGER, INTENT(IN) :: iopt
REAL, DIMENSION(:), INTENT(OUT) :: c !N+M
INTEGER, INTENT(OUT) :: irep

INTEGER, DIMENSION(SIZE(X,1),SIZE(C)-SIZE(X,2)) :: iw
REAL, DIMENSION(SIZE(C)-SIZE(X,2),SIZE(X,2)) :: ttt ! M,N
REAL, DIMENSION(SIZE(X,2),2*SIZE(X,2)-SIZE(C)) :: r ! N,N-M
REAL,DIMENSION(SIZE(X,2),SIZE(X,2)) :: rk !N,N
REAL,DIMENSION(2*SIZE(X,2)-SIZE(C),2*SIZE(X,2)-SIZE(C)) :: u !N-M,N-M
REAL,DIMENSION(2*SIZE(X,2)-SIZE(C)):: db !N-M
REAL, DIMENSION(2*SIZE(X,2)-SIZE(C)) :: w !N-M
REAL,DIMENSION(SIZE(R,2)) :: rz !N-M
REAL, DIMENSION(SIZE(X,2),SIZE(X,2)) :: q ! N,N
REAL :: rn
REAL, DIMENSION(SIZE(X,2)) :: rkc !RKC(N)


INTEGER :: n
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SP0CVQ',0,zhook_handle)
n=SIZE(x,2)

! Test des differents parametres
IF (nord.LT.1) THEN
  WRITE(*,fmt='(A17)') "SP0CVQ : NORD < 1"
ELSEIF (m.LT.0) THEN
  WRITE(*,'(A14)') "SP0CVQ : M < 0"
ELSEIF (n.LT.m) THEN
  WRITE(*,'(A14)') "SP0CVQ : N < M"
ELSEIF (s2.LT.0.) THEN
  WRITE(*,'(A15)') "SP0CVQ : S2 < 0"
ELSEIF (p.LT.0.) THEN
  WRITE(*,'(A14)') "SP0CVQ : P < 0"
ELSEIF ((lissage.LT.0.OR.lissage.GT.3).AND.lissage.NE.10) THEN
  WRITE(*,'(A27)') "SP0CVQ : LISSAGE < 0 OU > 3"
ELSEIF (lordre.LT.0.OR.lordre.GT.1) THEN
  WRITE(*,'(A26)') "SP0CVQ : LORDRE < 0 OU > 1"
ELSEIF (iopt.LT.0.OR.iopt.GT.1) THEN
  WRITE(*,'(A24)') "SP0CVQ : IOPT < 0 OU > 1"
ELSE

!     Determination de nord opt ou de p pot
!   -----------------------------------------
  IF (lordre.EQ.1 .OR. &
       lissage.EQ.1 .OR. lissage.EQ.2 .OR. lissage.EQ.10)  THEN 
       !calcul de P en passant par SPLR, SPLU, SPLW, si LORDRE = 0
       !calcul de NORD, M, S2, P si LORDRE=1
    CALL sp0nop(x,g,z,ds,lissage,lordre,iopt,nord,m,s2,p,&
         iw,ttt,r,rk,u,db,w,rz,irep) 
  ELSEIF (lissage.EQ.0) THEN
    p=0.
  ENDIF
!---------------------------------------------------------------------
!---------------------------------------------------------------------

  IF (lordre.EQ.1 .OR. lissage.EQ.0 .OR. lissage.EQ.3) THEN
    IF (iopt.EQ.0) THEN
!     Calcul de IW, TTT, R
      CALL splr(nord,x,ds,iw,ttt,r,irep)
      IF (irep.NE.0 .AND. lhook) CALL dr_hook('MODE_SPLINES:SP0CVQ',1,zhook_handle)
      IF (irep.NE.0) RETURN
!     Calcul de  RK, U, DB
      CALL splu(nord,x,g,r,rk,u,db)  
    ENDIF
!---------------------------------------------------------------------
!   Calcul de W,RZ
    CALL splw(z,r,u,w,rz)
!   Calcul de S2 pour lissage=3
    IF (lissage.EQ.3) CALL splrs(db,w,n,p,s2)
  ENDIF

!---------------------------------------------------------------------
!---------------------------------------------------------------------
! Calcul de Q
  CALL splvpq(p,r,rk,u,db,q)
! Calcul de RN
  CALL splbvm(db,w,n,p,rn)
! Calcul des coefficients C et RKC
  CALL splc(q,z,rk,ds,ttt,rkc,c)
  IF (lhook) CALL dr_hook('MODE_SPLINES:SP0CVQ',1,zhook_handle)
  RETURN

ENDIF

irep=-1
IF (lhook) CALL dr_hook('MODE_SPLINES:SP0CVQ',1,zhook_handle)

END SUBROUTINE sp0cvq
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!40. AI(ND,ND,NSDMAX): calcul de AI en fonction de XD, NSD, INTER; 
!remarque: on dirait que la deuxième dimension de AI est plutôt 2 que ND
SUBROUTINE splbsd(NSD,INTER,XD,AI)
!SPLB2C
IMPLICIT NONE

INTEGER, DIMENSION(:), INTENT(IN) :: nsd
INTEGER, INTENT(IN) :: inter
REAL, DIMENSION(:,:), INTENT(IN) :: xd
REAL, DIMENSION(:,:,:), INTENT(OUT) :: ai
REAL(KIND=JPRB) :: zhook_handle

REAL, DIMENSION(SIZE(NSD)) :: dxi, dxr
INTEGER :: j, i

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBSD',0,zhook_handle)
DO j=1,SIZE(nsd)
  dxi(j)=(xd(2,j)-xd(1,j))/nsd(j)
  dxr(j)=dxi(j)/(2*inter-2)
ENDDO

DO j=1,SIZE(nsd)
  ai(j,1,1)=xd(1,j)
  DO i=2,nsd(j)
    ai(j,1,i)=xd(1,j)+(i-1)*dxi(j)-dxr(j)
  ENDDO
  ai(j,2,nsd(j))=xd(2,j)
  DO i=1,nsd(j)-1
    ai(j,2,i)=xd(2,j)-(nsd(j)-i)*dxi(j)+dxr(j)
  ENDDO
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBSD',1,zhook_handle)

END SUBROUTINE splbsd
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!41. récupère les X et Z tels que les valeurs soient toutes comprises entre 
!A et B, dans XS et ZS. NS est la taille finale des tableaux. 
SUBROUTINE splbsel(X,Z,AI,BI,NS,XS,ZS,IREP)
!SPLB2C
IMPLICIT NONE

REAL, DIMENSION(:,:), INTENT(IN) :: x !ND,N
REAL, DIMENSION(:), INTENT(IN) :: z !N
REAL, DIMENSION(:), INTENT(IN) :: ai !ND
REAL, DIMENSION(:), INTENT(IN) :: bi !ND
INTEGER, INTENT(OUT) :: ns
REAL, DIMENSION(:,:), INTENT(OUT) :: xs !ND,NMAX
REAL, DIMENSION(:), INTENT(OUT) :: zs !NMAX
INTEGER, INTENT(OUT) :: irep

INTEGER :: i, j, i0
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBSEL',0,zhook_handle)
irep=0
ns=0
DO j=1,SIZE(x,2)
  i0=0
  DO i=1,SIZE(x,1)
   !il faut que tous les X(I,J) soient compris entre AI et BI
    IF (x(i,j).LT.ai(i).OR.x(i,j).GT.bi(i)) THEN 
      i0=1
      EXIT
    ENDIF
  ENDDO
  IF (i0.NE.1) THEN
    ns=ns+1
    IF (ns.GT.SIZE(zs)) THEN
      irep=-1
      IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBSEL',1,zhook_handle)
      RETURN
    ENDIF
    DO i=1,SIZE(x,1)
      xs(i,ns)=x(i,j)
    ENDDO
    zs(ns)=z(j)
  ENDIF
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLBSEL',1,zhook_handle)

END SUBROUTINE splbsel
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!42. appelée par INTERPOL_SPLINES
SUBROUTINE splb2c(NORD,M,X,G,Z,S2,P,LISSAGE,IOPT,NSDI,NSDJ,&
                        inter,xd,c,irep)
!INTERPOL_SPLINES
IMPLICIT NONE

INTEGER, INTENT(INOUT) :: nord 
INTEGER, INTENT(INOUT) :: m
REAL, DIMENSION(:,:), INTENT(IN) :: x !ND,N
REAL, DIMENSION(:,:), INTENT(IN) :: g ! ND,ND
REAL, DIMENSION(:), INTENT(IN) :: z !N
REAL, INTENT(INOUT) :: s2
REAL, INTENT(INOUT) :: p
INTEGER, INTENT(IN) :: lissage
INTEGER, INTENT(IN) :: iopt
INTEGER, INTENT(IN) :: nsdi
INTEGER, INTENT(IN) :: nsdj
INTEGER, INTENT(IN) :: inter
REAL, DIMENSION(:,:), INTENT(IN) :: xd !2,ND
REAL, DIMENSION(:,:,:), INTENT(OUT) :: c !nmax+mmax,nsdi,nsdj
INTEGER, INTENT(OUT) :: irep
!
INTEGER :: nn, nd
INTEGER :: lordre, i, j
INTEGER, DIMENSION(2) :: nsd
REAL, DIMENSION(NMAX) :: ds
REAL, DIMENSION(2,2) :: ci

REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_SPLINES:SPLB2C',0,zhook_handle)
nd=SIZE(x,1)
IF (.NOT.ALLOCATED(ai)) ALLOCATE(ai(nd,2,nsdmax))
IF (.NOT.ALLOCATED(ns)) ALLOCATE(ns(nsdmax,nsdmax))
IF (.NOT.ALLOCATED(zs)) ALLOCATE(zs(nmax,nsdmax,nsdmax))
IF (.NOT.ALLOCATED(xs)) ALLOCATE(xs(nd,nmax,nsdmax,nsdmax))

! Initialisation de Ds
ds(:)=1.

! Generation des sous-domaines: matrice AI de MODD_SPLINES
nsd(1)=nsdi
nsd(2)=nsdj
 CALL splbsd(nsd,inter,xd,ai)

lordre=0

! Calcul des splines par sous-domaines
DO i=1,nsdi
  DO j=1,nsdj
    ci(1,1)=ai(1,1,i)
    ci(1,2)=ai(1,2,i)
    ci(2,1)=ai(2,1,j)
    ci(2,2)=ai(2,2,j)
    !selection des points
    CALL splbsel(x,z,ci(:,1),ci(:,2),ns(i,j),xs(:,:,i,j),zs(:,i,j),irep)
    IF (irep.NE.0) THEN
      WRITE(*,fmt='(A35)') "SPLB2C: NOMBRE DE POINTS TROP GRAND"
      irep=-4
      CALL abor1_sfx('SPLINES: SPLB2C: NOMBRE DE POINTS TROP GRAND')
    ELSEIF (ns(i,j).LE.m) THEN
      WRITE(*,fmt='(A42)') &
        "SPB2C : NB DE PTS <= M : COEF NON CALCULES"
    ENDIF
    !calcul de la spline
    IF (ns(i,j).GT.m) THEN
       nn=ns(i,j) 
      !X(ND,NN) ZS(NN) C(NN+M) IW(ND,M) W(7*NN*NN)
      CALL sp0cvq(nord,m,xs(:,1:nn,i,j),g,zs(1:nn,i,j),ds(1:nn),s2,p,&
        lissage,lordre,iopt,c(1:nn+m,i,j),irep)   
    ENDIF 
  ENDDO
ENDDO
IF (lhook) CALL dr_hook('MODE_SPLINES:SPLB2C',1,zhook_handle)

END SUBROUTINE splb2c


END MODULE mode_splines