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

	SUBROUTINE bsp2dn(nf,f,x,xt,nx,knotx,mx,lx,y,yt,ny,knoty,my,ly,
	1 init,xx,yy,fxy,dfxydx,dfxydy)

c	interpolation de nf fonctions à 2 variables
c	dans le tableau f(nf,nx,ny) par produit tensoriel
c	de splines polynomiales d'ordres mx{y} > 1, au point
c	(xx,yy) : x(1) <= xx <= x(nx), y(1) <= yy <= y(ny), on aura aussi
c	 xt(lx) <= xx <= xt(lx+1) ; yt(ly) <= yy <= yt(ly+1)
c	en entrée prendre lx{y}=mx{y}
c	au premier appel, init=.FALSE. il y a initialisation :
c	calul de knotx=mx+nx, formation du tableau xt(knotx). idem pour y
c	!! ATTENTION le tableau f(nf,nx,ny) initial est remplacé par le
c	tableau des coefficients des B-splines en xy
c	ce calcul étant fait selon de Boor p.342
c	a une transition près F = tAx^(-1) f  Ay^(-1)

c	version F95 de sbsp2d.f

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

c entrées:
c	init=.true. : calcul des coefficients des splines
c	nf: nombre de fonctions
c	x / y : abscisses orDOnnees
c	nx / ny: nombre  de points
c	mx / my: ordre des splines
c	(xx,yy) : point d'interpolation

c entrées/sorties
c	knotx / knoty: nb de noeuds
c	f(nf,nx,ny): fonctions a interpoler / coefficients des splines
c	xt / yt: points de raccord
c	lx / ly: localisation

c sorties
c	fxy, dfxydx, dfxydy : fonctions, derivees 1-ieres

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

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

	USE mod_kind
		
	IMPLICIT NONE

	REAL (kind=dp), INTENT(in), DIMENSION (:) :: x, y
	REAL (kind=dp), INTENT(in) :: xx, yy
	INTEGER, INTENT(in) :: nx, mx, ny, my, nf	
	logical, INTENT(in) :: init
			
	REAL (kind=dp), INTENT(inout), DIMENSION(:,:,:) :: f			
	REAL (kind=dp), INTENT(inout), DIMENSION (:) :: xt, yt
	INTEGER, INTENT(inout) :: knotx, lx, knoty, ly
				
	REAL (kind=dp), INTENT(out), DIMENSION(:) :: fxy, dfxydx, dfxydy
			
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:,:) :: fv		
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: ax, axm1, ay, aym1, f1		
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: qx, qy, dx, dy
	REAL (kind=dp) :: xx1, yy1

	INTEGER, DIMENSION(:), ALLOCATABLE :: indpc	
	INTEGER :: i, j, k

	LOGICAL :: inversible	

c------------------------------------------------------------------------
	
2000	FORMAT(8es10.3)

	IF(.NOT.init)THEN	!initialisations

c formation de tAy^(-1)
	 ALLOCATE(ay(ny,ny),aym1(ny,ny),qy(my+1),indpc(ny))
	 CALL noein(y,yt,ny,my,knoty) ; IF(no_croiss)RETURN

c ay <- matrice 1-D des interpolations en y
	 ay=0.d0 ; ly=my	 
	 DO j=1,ny
	  CALL linf(y(j),yt,knoty,ly)
	  CALL bval0(y(j),yt,my,ly,qy)
	  DO k=1,my
	   ay(j,ly-my+k)=qy(k)
	  ENDDO
	 ENDDO
	 		 
c Ay^(-1) 
	 aym1=0.d0
	 DO i=1,ny
	  aym1(i,i)=1.d0
	 ENDDO
	 indpc=1	 
	 CALL gauss_band(ay,aym1,indpc,ny,ny,ny,ny,inversible)
	 IF(.NOT.inversible)THEN
	  PRINT*,'dans bsp2dn Ay n''est pas inversible, Arrêt'
	  STOP
	 ENDIF

c nettoyage
	 DEALLOCATE(ay,qy,indpc)
	 
c formation de tAx^(-1)
	 ALLOCATE(ax(nx,nx),axm1(nx,nx),qx(mx+1),indpc(nx)) 
c	 PRINT*,'x' ; WRITE(*,2000)x
	 CALL noein(x,xt,nx,mx,knotx) ; IF(no_croiss)RETURN

c ax <- matrice 1-D des interpolations en x
	 ax=0.d0
	 lx=mx	 
	 DO j=1,nx
	  CALL linf(x(j),xt,knotx,lx)
	  CALL bval0(x(j),xt,mx,lx,qx)
	  DO k=1,mx
	   ax(j,lx-mx+k)=qx(k)
	  ENDDO
	 ENDDO

c tAx^(-1)	 
	 ax=TRANSPOSE(ax)	!tAx	  
	 axm1=0.d0
	 DO i=1,nx
	  axm1(i,i)=1.d0
	 ENDDO
	 indpc=1
	 CALL gauss_band(ax,axm1,indpc,nx,nx,nx,nx,inversible)
	 IF(.NOT.inversible)THEN
	  PRINT*,'dans bsp2dn Ax n''est pas inversible, Arrêt'
	  STOP
	 ENDIF

c nettoyage
	 DEALLOCATE(ax,qx,indpc)

c f=tAx^(-1)*F*Ay^(-1)
	 ALLOCATE(f1(nx,ny),fv(nf,nx,ny)) 
	 fv=f
	 DO  i=1,nf
	  f1(:,:)=fv(i,:,:)
	  f1=MATMUL(axm1,matmul(f1,aym1))	!f=tAx^(-1)*F*Ay^(-1)
	  f(i,:,:)=f1 
	 ENDDO
	 DEALLOCATE(f1,fv,axm1,aym1)	 
	ENDIF

c----------------initialisations fin-----------------------------	
			
c localisation de (xx,yy)
	IF(xx < xt(1) .OR. xx > xt(knotx))THEN
	 WRITE(*,1)xt(1),xx,knotx,xt(knotx)
1	 FORMAT('dans bsp2dn, interpolation 2-D, le point xx sort'
	1 ' de la grille xt(1)=',es15.8,' xx=',es15.8,
	2 ' xt(',i3,')=',es15.8)
	ENDIF	
	IF(yy < yt(1) .OR. yy > yt(knoty))THEN
	 WRITE(*,2)yt(1),yy,knoty,yt(knoty)
2	 FORMAT('dans sbsp2dn, interpolation 2-D, le point yy sort'
	1' de la grille yt(1)=',es15.8,' yy=',es15.8,' yt(',i3,')=',
	2 es15.8)
	ENDIF

c on ne déborde pas!
	xx1=MIN(x(nx),MAX(xx,x(1)))
	CALL linf(xx1,xt,knotx,lx)
	yy1=MIN(y(ny),MAX(yy,y(1)))	 	
	CALL linf(yy1,yt,knoty,ly)

c interpolation
	ALLOCATE(qx(mx+1),qy(my+1),dx(mx),dy(my))
	CALL bval1(xx1,xt,mx,lx,qx,dx)
	CALL bval1(yy1,yt,my,ly,qy,dy)
	
	fxy=0.d0 ; dfxydx=0.d0 ; dfxydy=0.d0
	DO k=1,nf
	 DO i=1,mx
	  DO j=1,my
	   fxy(k)=fxy(k)      +f(k,lx-mx+i,ly-my+j)*qx(i)*qy(j)
	   dfxydx(k)=dfxydx(k)+f(k,lx-mx+i,ly-my+j)*dx(i)*qy(j)
	   dfxydy(k)=dfxydy(k)+f(k,lx-mx+i,ly-my+j)*qx(i)*dy(j)
	  ENDDO
	 ENDDO
	ENDDO
	
c nettoyage	
	DEALLOCATE(qx,qy,dx,dy)

	RETURN
	
	END SUBROUTINE bsp2dn