INCLUDE '../SOURCE_T/mod_donnees.f'
	INCLUDE '../SOURCE_T/mod_conv.f'
	
c**********************************************************************

	PROGRAM test_convection

c	verification des dérivées de la routine de convection

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

c	CESAM2k

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

	USE mod_donnees, ONLY : alpha, ini_ctes, lit_nl, nom_fich2,
	1 print_ctes, rsol, w_rot
	USE mod_kind
	USE mod_conv, ONLY : conv
	USE mod_variables, ONLY : id_conv, if_conv, rstar, r_zc_conv	     

	IMPLICIT NONE

	REAL (kind=dp), DIMENSION(9) ::  dg, dy, y
	REAL (kind=dp), PARAMETER :: dd=1.d-5, unpdd=1.d0+dd
	REAL (kind=dp) :: krad, gravite, delta, cp, ro, hp, taur, gradrad,
	1 gradad, grad, grad0, dgradkra, dgradgra, dgradel, dgradcp,
	2 dgradro, dgradtaur, dgradhp, dgradgrad, dgradgad,
	3 r, stor, stor0, dstor,
	4 gam, gam0, dgamkra, dgamgra, dgamdel, dgamcp, dgamro,
	5 dgamhp, dgamtaur, dgamgrad, dgamgad

	INTEGER :: cas, i

	CHARACTER (len=7), DIMENSION(9) :: deriv=(/ 'krad   ','gravite',
	1 'delta  ','cp     ','ro     ','hp     ','taur   ','gradrad',
	2 'gradad ' /)

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

2000	FORMAT(8es10.3)
2001	FORMAT(10es9.2)

c	lecture de test.don

	nom_fich2='test' ; CALL lit_nl(w_rot) ; CALL ini_ctes

	DO cas=0,6
	 SELECT CASE(cas)
	 CASE(0)
	  id_conv=0 ; if_conv=3
	  r_zc_conv(0)=-1.259d-1 ; r_zc_conv(1)=1.259d-1
	  r_zc_conv(2)=0.6232d0 ; r_zc_conv(3)=0.860789d0
	  krad=5.059d12 ; gravite=3.698d4 ; delta=1.d0 ; cp=1.365d9
	  ro=3.081d-7 ; hp=2.327d7 ; taur=6.489D+02 ; gradrad=2.181d1
	  gradad=1.080d-1 ; r=0.860789d0
	 CASE(1)
	  id_conv=0 ; if_conv=5
	  r_zc_conv(id_conv)=0.0d0 ; r_zc_conv(1)=0.2d0
	  r_zc_conv(2)=0.5d0 ; r_zc_conv(3)=0.8d0
	  r_zc_conv(4)=0.9d0 ; r_zc_conv(if_conv)=1.d0
	  krad=1.153D+13 ; gravite=3.784D+04 ; delta=1.174D+00 ; cp=3.870D+08
	  ro=4.802D-07 ; hp=1.617D+07 ; taur=6.489D+02 ; gradrad=8.399D+00
	  gradad=1.996D-01 ; r=0.99d0
	 CASE(2)
	  id_conv=0 ; if_conv=5
	  r_zc_conv(id_conv)=0.0d0 ; r_zc_conv(1)=0.2d0
	  r_zc_conv(2)=0.5d0 ; r_zc_conv(3)=0.8d0
	  r_zc_conv(4)=0.9d0 ; r_zc_conv(if_conv)=1.d0
	  krad=1.460D+14 ; gravite=3.513D+04  ; delta=1.000D+00 ; cp=1.63D+08
	  ro=4.091D-07 ; hp=1.220D+07 ; taur=1.352D+01 ; gradrad=0.7233
	  gradad=3.983D-01 ; r=0.55d0
	 CASE(3)
	  id_conv=0 ; if_conv=5
	  r_zc_conv(id_conv)=0.0d0 ; r_zc_conv(1)=0.2d0
	  r_zc_conv(2)=0.5d0 ; r_zc_conv(3)=0.8d0
	  r_zc_conv(4)=0.9d0 ; r_zc_conv(if_conv)=1.d0
	  krad=6.548D+15 ; gravite=2.255D+04 ; delta=9.652D-01 ; cp=3.311D+08
	  ro=7.757D+01 ; hp=8.188D+10 ; taur=9.309D+12 ; gradrad=7.657D-01
	  gradad=3.986D-01 ; r=0.01d0
	 CASE(4)
	  id_conv=0 ; if_conv=3
	  r_zc_conv(0)=-1.259d-1 ; r_zc_conv(1)=1.259d-1
	  r_zc_conv(2)=0.6232d0 ; r_zc_conv(3)=0.860789d0
	  krad=6.740d15 ; gravite=3.880d4 ; delta=0.d0 ; cp=3.295d8
	  ro=7.778d1 ; hp=4.771d10 ; taur=0.d0 ; gradrad=8.309d-1
	  gradad=3.985d-1 ; r=3.542d-3
	 CASE(5)    !r en dehors de la ZC 1
	  id_conv=0 ; if_conv=3
	  r_zc_conv(0)=-1.259d-1 ; r_zc_conv(1)=1.259d-1
	  r_zc_conv(2)=0.6232d0 ; r_zc_conv(3)=0.860789d0
	  krad=6.740d15 ; gravite=3.880d4 ; delta=0.d0 ; cp=3.295d8
	  ro=7.778d1 ; hp=4.771d10 ; taur=0.d0 ; gradrad=8.309d-1
	  gradad=3.985d-1 ; r=1.3d-1
	 CASE(6)    !dans l'atmosphère
	  id_conv=0 ; if_conv=3 ; r=1.614293547199790
	  r_zc_conv(0)=-1.939E-01 ; r_zc_conv(1)=1.939E-01
	  r_zc_conv(2)=1.465560530330532	  
	  r_zc_conv(3)=1.614314981893674
	  krad=8.768E+14 ; gravite= 1.535E+04 ; delta=1.d0 ; cp=1.546E+08
	  ro=1.503E-07 ; hp=2.540E+07 ; taur=1.d2 ; gradrad=4.085E-01
	  gradad=3.983E-01 ; rstar=1.7
	 END SELECT
         y(1)=krad ; y(2)=gravite ; y(3)=delta ; y(4)=cp ; y(5)=ro
         y(6)=hp ; y(7)=taur ; y(8)=gradrad ; y(9)=gradad
	 PRINT*,'R, Y, cas:',cas ; WRITE(*,2000)r ; WRITE(*,2001)y
	 CALL conv(r,krad,gravite,delta,cp,ro,hp,taur,gradrad,gradad,.TRUE.,
	1  grad0,dgradkra,dgradgra,dgradel,dgradcp,dgradro,
	2  dgradhp,dgradtaur,dgradgrad,dgradgad,
	3  gam0,dgamkra,dgamgra,dgamdel,dgamcp,dgamro,
	4  dgamhp,dgamtaur,dgamgrad,dgamgad)

	 PRINT*,'grad0',grad0 ; PRINT*,'gam0',gam0

c	PRINT*,'dgradkra,dgradgra,dgradel,dgradcp,dgradro,dgradtaur,dgradhp,dgradgrad,dgradgad'
c	WRITE(*,2001)dgradkra,dgradgra,dgradel,dgradcp,dgradro,
c	1   dgradhp,dgradtaur,dgradgrad,dgradgad
c	PRINT*,'dgamkra,dgamgra,dgamdel,dgamcp,dgamro,dgamhp,dgamtaur,dgamgrad,dgamgad' 
c	WRITE(*,2001)dgamkra,dgamgra,dgamdel,dgamcp,dgamro,
c	1   dgamhp,dgamtaur,dgamgrad,dgamgad

	 dy(1)=dgradkra ; dy(2)=dgradgra ; dy(3)=dgradel ; dy(4)=dgradcp
	 dy(5)=dgradro ; dy(6)=dgradhp ; dy(7)=dgradtaur ; dy(8)=dgradgrad
	 dy(9)=dgradgad

	 dg(1)=dgamkra ; dg(2)=dgamgra ; dg(3)=dgamdel ; dg(4)=dgamcp
	 dg(5)=dgamro ; dg(6)=dgamhp ; dg(7)=dgamtaur ; dg(8)=dgamgrad
	 dg(9)=dgamgad

	 DO i=1,9
	  stor0=y(i) ; stor=stor0*unpdd
	  IF(stor == 0.d0)stor=dd
	  dstor=stor-stor0 ; y(i)=stor
	  krad=y(1) ; gravite=y(2) ; delta=y(3) ; cp=y(4) ; ro=y(5)
	  hp=y(6) ; taur=y(7) ; gradrad=y(8) ; gradad=y(9)
	  CALL conv(r,krad,gravite,delta,cp,ro,hp,taur,gradrad,gradad,.TRUE.,
	1   grad,dgradkra,dgradgra,dgradel,dgradcp,dgradro,
	2   dgradhp,dgradtaur,dgradgrad,dgradgad,
	3   gam,dgamkra,dgamgra,dgamdel,dgamcp,dgamro,
	4   dgamhp,dgamtaur,dgamgrad,dgamgad)
	  y(i)=stor0
	  PRINT*,'derivee ',i,deriv(i)
	  WRITE(*,2001)(grad-grad0)/dstor,dy(i),(gam-gam0)/dstor,dg(i)
	  WRITE(*,'(a)')'dérivée/'//deriv(i) 
	  PAUSE'suivant'
	 ENDDO
	ENDDO

	STOP

	END PROGRAM test_convection