SUBROUTINE Bspline(X0, lgT0, lgQs0, Ys, dYs)

! perform Bspline interpolation of SAHA-S3 tables for given X0, lgT0, lgQs0.
! Bspline computes all 12 values Ys instead of two (as Bspline2f does).
! Bspline computes also derivatives dYs of some these values with respect to X, lgT, lgQs.

	IMPLICIT NONE

	REAL*8, INTENT(in) :: X0, lgT0, lgQs0
	REAL*8, DIMENSION(mv), INTENT(out) :: Ys
	REAL*8, DIMENSION(mv,3), INTENT(out) :: dYs  ! massiv of derivatives with respect to X, lgT, lgQs

	INTEGER ix,iT,iQs,iv
	INTEGER LCX, LCT,LCQs
	INTEGER LeftX/1/, LeftT/1/, LeftQs/1/
	SAVE LeftX, LeftT, LeftQs
	REAL*8, DIMENSION(Jspline,2) :: dBIATX, dBIATT, dBIATQs
	REAL*8 bttt, bttt_X, bttt_lgT, bttt_lgQs
	INTEGER nderiv, mflag


	nderiv=2
	CALL interv8m ( Xknot, X0, LeftX, mflag )
	CALL BSPLVD8( Xknot, Jspline, X0, LeftX, dBIATX, nderiv)
	LCX=LeftX-Jspline+1

	CALL interv8m ( Tknot, lgT0, LeftT, mflag )
	CALL BSPLVD8( Tknot, Jspline, lgT0, LeftT, dBIATT, nderiv)
	LCT=LeftT-Jspline+1

	CALL interv8m ( Qsknot, lgQs0, LeftQs, mflag )
	CALL BSPLVD8( Qsknot, Jspline, lgQs0, LeftQs, dBIATQs, nderiv)
	LCQs=LeftQs-Jspline+1

! functions Ys
	Ys=0.d0
	DO ix=1,Jspline
	 DO iT=1,Jspline
	  DO iQs=1,Jspline
	   bttt=dBIATX(ix,1)*dBIATT(iT,1)*dBIATQs(iQs,1)
	   DO iv=1,mv
	    Ys(iv)=Ys(iv)+Beos(LCQs+iQs-1,LCT+iT-1,iv,LCX+ix-1)*bttt
	   ENDDO
	  ENDDO
         ENDDO
	ENDDO

! derivatives
	dYs=0.d0
	DO ix=1,Jspline
	 DO iT=1,Jspline
	  DO iQs=1,Jspline
	   bttt_X=dBIATX(ix,2)*dBIATT(iT,1)*dBIATQs(iQs,1)
	   bttt_lgT=dBIATX(ix,1)*dBIATT(iT,2)*dBIATQs(iQs,1)
	   bttt_lgQs=dBIATX(ix,1)*dBIATT(iT,1)*dBIATQs(iQs,2)
	   DO iv=4,8,2
	    dYs(iv,1)=dYs(iv,1)+Beos(LCQs+iQs-1,LCT+iT-1,iv,LCX+ix-1)*bttt_X
	    dYs(iv,2)=dYs(iv,2)+Beos(LCQs+iQs-1,LCT+iT-1,iv,LCX+ix-1)*bttt_lgT
	    dYs(iv,3)=dYs(iv,3)+Beos(LCQs+iQs-1,LCT+iT-1,iv,LCX+ix-1)*bttt_lgQs
	   ENDDO
	  ENDDO
	 ENDDO
	ENDDO

	RETURN

	END SUBROUTINE Bspline