c	INCLUDE '../SOURCE/mod_numerique.f'
	
c****************************************************************

	PROGRAM test_lim_ZC

c	test d'identification des ZC

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

	use mod_kind
	USE mod_numerique, ONLY : bsp1dn, bsp_dis
	
	IMPLICIT NONE

	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: f, fd ,ff, ffd
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: dfxdx, fx, x, xt
	REAL (kind=dp), PARAMETER :: eps=1.d-3
	REAL (kind=dp) :: xx
		
	INTEGER, ALLOCATABLE, DIMENSION(:) :: id
	INTEGER, ALLOCATABLE, DIMENSION(:) :: convd, convf
	INTEGER :: i, izc, j, knot, l=1, m, n, nd, nzc
	
	LOGICAL, ALLOCATABLE, DIMENSION(:) :: ld
	
c---------------------------------------------------------------------------

2000	FORMAT(8es10.3)

c	il y a 3 ZC
	nzc=3 ; ALLOCATE(convd(nzc),convf(nzc))
	convd=(/ 1, 15, 25 /)
	convf=(/ 10, 20, 30 /)
	
c test ZC - ZR - ZC - ZR - ZC
	
	n=30 ; nd=2*nzc
	IF(convd(1) == 1)nd=nd-1
	IF(convf(nzc) == n)nd=nd-1	
	PRINT*,'nd=',nd
	
	ALLOCATE(id(0:nd+1),ld(nd))
	id(0)=1 ; id(nd+1)=n ; j=0
	DO izc=1,nzc
	 IF(convd(izc) /= 1)THEN
	  j=j+1 ; id(j)=convd(izc) ; ld(j)=.TRUE.
	 ENDIF
	 IF(convf(izc) /= n)THEN 
	  j=j+1 ; id(j)=convf(izc) ; ld(j)=.FALSE.
	 ENDIF
	ENDDO
	PRINT*,'id(1:nd)',id(1:nd) ; PRINT*,'ld(1:nd)',ld(1:nd)
	ALLOCATE(f(1,n+nd),fd(1,nd))
	f(1,:)=1.d0 ; fd(1,:)=1.d0
	DO izc=1,nzc
	 DO j=convd(izc)+1,convf(izc)
	  f(1,j)=-1.d0
	 ENDDO
	ENDDO
	IF(convd(1) == 1)f(1,1)=-1.d0
	DO j=1,nd
	 IF(ld(j))fd(1,j)=-1.d0
	ENDDO
	WRITE(*,2000)f(1,:) ; WRITE(*,2000)fd(1,:)
	m=2 ; knot=n+nd+m	
	ALLOCATE(dfxdx(1),fx(1),x(n),xt(knot))
	x=(/ (i,i=1,n) /)
	
	j=1
	DO i=1,n
	 WRITE(*,2000)x(i),f(1,i)
	 IF(i == id(j))THEN
	  PRINT*,'limite j=',j,', fd=',fd(1,j)
	  j=MIN(j+1,nd)
	 ENDIF
	ENDDO
	CALL bsp_dis(1,x,f,nd,id,fd,n,m,xt,knot,.FALSE.)
	
	j=1
	DO i=1,n
	 xx=x(i)
	 CALL bsp1dn(1,f,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	 WRITE(*,2000)xx,fx
	 IF(i == id(j))THEN
	  xx=xx+eps
	  CALL bsp1dn(1,f,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	  PRINT*,'limite j=',j,', xx=', xx,', fx=',fx
	  j=MIN(j+1,nd)
	 ENDIF
	ENDDO
	
c	remplissage

	ALLOCATE(ff(1,n+nd),ffd(1,nd))
	j=1
	DO i=1,n
	 xx=x(i)	
	 CALL bsp1dn(1,f,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)	
	 IF(fx(1) > 0.d0)THEN
c	zone radiative
	  CALL z_rad
	  ff(1,i)=fx(1)
c	est-on sur une limite? si oui  on ne peut être que convectif à droite
	  IF(i == id(j))THEN
	   CALL z_conv
	   ffd(1,j)=fx(1)
	   j=MIN(nd,j+1)
	  ENDIF
	 ELSE	 
c	zone convective
	  CALL z_conv
	  ff(1,i)=fx(1)
c	est-on sur une limite? si oui  on ne peut être que radiatif à droite
	  IF(i == id(j))THEN
	   xx=x(i)+eps
	   CALL z_rad
	   ffd(1,j)=fx(1)
	   j=MIN(nd,j+1)
	  ENDIF
	 ENDIF
	ENDDO
	
	CALL bsp_dis(1,x,ff,nd,id,ffd,n,m,xt,knot,.FALSE.)	

	PRINT*,'test ZC - ZR - ZC - ZR - ZC'	
	DO i=1,n
	 xx=x(i)
	 CALL bsp1dn(1,ff,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	 WRITE(*,2000)xx,fx(1),dfxdx(1)
	 IF(i /= n)THEN
	  xx=x(i)+0.1d0
	  CALL bsp1dn(1,ff,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	  WRITE(*,2000)xx,fx(1),dfxdx(1)
	 ENDIF	 
	ENDDO		
	PAUSE'test ZC - ZR - ZC - ZR - ZC'	

c	n=n+5 test ZC - ZR - ZC - ZR- ZC - ZR

	DEALLOCATE(dfxdx,f,fd,ff,ffd,fx,id,ld,x,xt)	 
	n=n+5 ; nd=2*nzc
	IF(convd(1) == 1)nd=nd-1
	IF(convf(nzc) == n)nd=nd-1	
	PRINT*,'nd=',nd
	
	ALLOCATE(id(0:nd+1),ld(nd))
	id(0)=1 ; id(nd+1)=n ; j=0
	DO izc=1,nzc
	 IF(convd(izc) /= 1)THEN
	  j=j+1 ; id(j)=convd(izc) ; ld(j)=.TRUE.
	 ENDIF
	 IF(convf(izc) /= n)THEN 
	  j=j+1 ; id(j)=convf(izc) ; ld(j)=.FALSE.
	 ENDIF
	ENDDO
	PRINT*,'id(1:nd)',id(1:nd) ; PRINT*,'ld(1:nd)',ld(1:nd)
	ALLOCATE(f(1,n+nd),fd(1,nd))
	f(1,:)=1.d0 ; fd(1,:)=1.d0
	DO izc=1,nzc
	 DO j=convd(izc)+1,convf(izc)
	  f(1,j)=-1.d0
	 ENDDO
	ENDDO
	IF(convd(1) == 1)f(1,1)=-1.d0
	DO j=1,nd
	 IF(ld(j))fd(1,j)=-1.d0
	ENDDO
	WRITE(*,2000)f(1,:) ; WRITE(*,2000)fd(1,:)
	m=4 ; knot=n+nd+m	
	ALLOCATE(dfxdx(1),fx(1),x(n),xt(knot))
	x=(/ (i,i=1,n) /)
	
	j=1
	DO i=1,n
	 WRITE(*,2000)x(i),f(1,i)
	 IF(i == id(j))THEN
	  PRINT*,'limite j=',j,', fd=',fd(1,j)
	  j=MIN(j+1,nd)
	 ENDIF
	ENDDO	  
c	PAUSE'avant bsp_dis'
	
	CALL bsp_dis(1,x,f,nd,id,fd,n,m,xt,knot,.FALSE.)
	
	j=1
	DO i=1,n
	 xx=x(i)
	 CALL bsp1dn(1,f,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	 WRITE(*,2000)xx,fx
	 IF(i == id(j))THEN
	  xx=xx+eps
	  CALL bsp1dn(1,f,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	  PRINT*,'limite j=',j,', xx=', xx,', fx=',fx
	  j=MIN(j+1,nd)
	 ENDIF
	ENDDO
	
c	remplissage

	ALLOCATE(ff(1,n+nd),ffd(1,nd))
	j=1
	DO i=1,n	
	 CALL bsp1dn(1,f,x,xt,n,m,knot,.TRUE.,x(i),l,fx,dfxdx)	
	 IF(fx(1) > 0.d0)THEN
c	zone radiative
	  ff(1,i)=x(i)
c	est-on sur une limite? si oui  on ne peut être que convectif à droite
	  IF(i == id(j))THEN
	   ffd(1,j)=-1.d0
	   j=MIN(nd,j+1)
	  ENDIF
	 ELSE	 
c	zone convective
	  ff(1,i)=-1.d0
c	est-on sur une limite? si oui  on ne peut être que radiatif à droite
	  IF(i == id(j))THEN
	   ffd(1,j)=x(i)+eps
	   j=MIN(nd,j+1)
	  ENDIF
	 ENDIF
	ENDDO
	
	CALL bsp_dis(1,x,ff,nd,id,ffd,n,m,xt,knot,.FALSE.)	

	PRINT*,'test ZC - ZR - ZC - ZR - ZC - ZR'	
	DO i=1,n
	 xx=x(i)
	 CALL bsp1dn(1,ff,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	 WRITE(*,2000)xx,fx(1),dfxdx(1)
	 IF(i /= n)THEN
	  xx=x(i)+0.1d0
	  CALL bsp1dn(1,ff,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	  WRITE(*,2000)xx,fx(1),dfxdx(1)
	 ENDIF	 
	ENDDO	
	
	PAUSE'test ZC - ZR - ZC - ZR - ZR'	

c	test ZR - ZC - ZR - ZC - ZR - ZC

	DEALLOCATE(dfxdx,f,fd,ff,ffd,fx,id,ld,x,xt)	 
	nd=2*nzc ; convd=convd+5 ; convf=convf+5
	IF(convd(1) == 1)nd=nd-1
	IF(convf(nzc) == n)nd=nd-1	
	PRINT*,'nd=',nd
	
	ALLOCATE(id(0:nd+1),ld(nd))
	id(0)=1 ; id(nd+1)=n ; j=0
	DO izc=1,nzc
	 IF(convd(izc) /= 1)THEN
	  j=j+1 ; id(j)=convd(izc) ; ld(j)=.TRUE.
	 ENDIF
	 IF(convf(izc) /= n)THEN 
	  j=j+1 ; id(j)=convf(izc) ; ld(j)=.FALSE.
	 ENDIF
	ENDDO
	PRINT*,'id(1:nd)',id(1:nd) ; PRINT*,'ld(1:nd)',ld(1:nd)
	ALLOCATE(f(1,n+nd),fd(1,nd))
	f(1,:)=1.d0 ; fd(1,:)=1.d0
	DO izc=1,nzc
	 DO j=convd(izc)+1,convf(izc)
	  f(1,j)=-1.d0
	 ENDDO
	ENDDO
	IF(convd(1) == 1)f(1,1)=-1.d0
	DO j=1,nd
	 IF(ld(j))fd(1,j)=-1.d0
	ENDDO
	WRITE(*,2000)f(1,:) ; WRITE(*,2000)fd(1,:)
	m=4 ; knot=n+nd+m	
	ALLOCATE(dfxdx(1),fx(1),x(n),xt(knot))
	x=(/ (i,i=1,n) /)
	
	j=1
	DO i=1,n
	 WRITE(*,2000)x(i),f(1,i)
	 IF(i == id(j))THEN
	  PRINT*,'limite j=',j,', fd=',fd(1,j)
	  j=MIN(j+1,nd)
	 ENDIF
	ENDDO	  
c	PAUSE'avant bsp_dis'
	
	CALL bsp_dis(1,x,f,nd,id,fd,n,m,xt,knot,.FALSE.)
	
	j=1
	DO i=1,n
	 xx=x(i)
	 CALL bsp1dn(1,f,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	 WRITE(*,2000)xx,fx
	 IF(i == id(j))THEN
	  xx=xx+eps
	  CALL bsp1dn(1,f,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	  PRINT*,'limite j=',j,', xx=', xx,', fx=',fx
	  j=MIN(j+1,nd)
	 ENDIF
	ENDDO
	
c	remplissage

	ALLOCATE(ff(1,n+nd),ffd(1,nd))
	j=1
	DO i=1,n	
	 CALL bsp1dn(1,f,x,xt,n,m,knot,.TRUE.,x(i),l,fx,dfxdx)	
	 IF(fx(1) > 0.d0)THEN
c	zone radiative
	  ff(1,i)=x(i)
c	est-on sur une limite? si oui  on ne peut être que convectif à droite
	  IF(i == id(j))THEN
	   ffd(1,j)=-1.d0
	   j=MIN(nd,j+1)
	  ENDIF
	 ELSE	 
c	zone convective
	  ff(1,i)=-1.d0
c	est-on sur une limite? si oui  on ne peut être que radiatif à droite
	  IF(i == id(j))THEN
	   ffd(1,j)=x(i)+eps
	   j=MIN(nd,j+1)
	  ENDIF
	 ENDIF
	ENDDO
	
	CALL bsp_dis(1,x,ff,nd,id,ffd,n,m,xt,knot,.FALSE.)	

	PRINT*,'test ZR - ZC - ZR - ZC - ZR - ZC'	
	DO i=1,n
	 xx=x(i)
	 CALL bsp1dn(1,ff,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	 WRITE(*,2000)xx,fx(1),dfxdx(1)
	 IF(i /= n)THEN
	  xx=x(i)+0.1d0
	  CALL bsp1dn(1,ff,x,xt,n,m,knot,.TRUE.,xx,l,fx,dfxdx)
	  WRITE(*,2000)xx,fx(1),dfxdx(1)
	 ENDIF	 
	ENDDO	
	
	PAUSE'test ZR - ZC - ZR - ZC - ZR - ZC'	

	STOP
	
	CONTAINS
	
c******************************************************************
	
	 SUBROUTINE z_conv
	 	 
c	zone convective	
	  fx(1)=-1.d0
	 END SUBROUTINE z_conv

c*****************************************************************
	 	
	 SUBROUTINE z_rad	

c	zone radiative
	  fx(1)=xx	
	 END SUBROUTINE z_rad
	 
	END PROGRAM test_lim_ZC