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

      PROGRAM test_tdetau

c     programme de test des lois T(tau)

c     Auteur: P. Morel, Departement J.D. Cassini, O.C.A., Observatoire de Nice
c     CESAM95

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

      USE mod_donnees, ONLY : nom_chemin, nom_tdetau
      USE mod_quasi_static, ONLY : tdetau

      IMPLICIT NONE

      INTEGER :: k

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

      nom_chemin='C:\SUN_STAR_DATA\' ; k=4

      b1: do
       k=k+1
       SELECT CASE(k)
       CASE(1)
        nom_tdetau='edding'
       CASE(2)
        nom_tdetau='hopf'
       CASE(3)
        nom_tdetau='k5750'
       CASE(4)
        nom_tdetau='k5777'
       CASE(5)
        nom_tdetau='roger00'
       CASE(6)
        nom_tdetau='roger02'
       CASE(7)
        nom_tdetau='roger05'
       CASE(8)
        nom_tdetau='roger10a'
       CASE DEFAULT
       EXIT b1
       END SELECT
       PRINT*,'Loi T(tau): ',nom_tdetau ; PRINT* ; CALL calc_tdetau
      ENDDO b1

      STOP

       CONTAINS

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

      SUBROUTINE calc_tdetau

      USE mod_kind
      USE mod_numerique, ONLY : pause

      IMPLICIT NONE

      REAL (kind=dp), PARAMETER :: dd=1.d-5, unpdd=1.-dd
      REAL (kind=dp) :: teff=5698.d0, grav=2150.d0, tau=2.d0
      REAL (kind=dp) :: t, t0, dtsdtau, dtsdteff, dtsdg, ro_ext0,
     1 ro_ext, dro_grav, dro_teff, f_tau, f_tau0, df_tau, df_tau0, d2f_tau,
     2 stor, stor0, dstor

      INTEGER :: k

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

2000  FORMAT(8es10.3)

      k=1
      b1: DO
       WRITE(*,1)tau,teff,grav
1      FORMAT('tau=',es10.3,', Teff=',es10.3,', grav=',es10.3)
       CALL tdetau(tau,teff,grav,t0,dtsdtau,dtsdteff,dtsdg,
     1  ro_ext0,dro_grav,dro_teff,f_tau0,df_tau0,d2f_tau)
       WRITE(*,2000)tau,teff,grav,t0,f_tau0,df_tau0,ro_ext0
       PRINT*,'derivee/tau'
       stor0=tau ; stor=stor0*unpdd
       IF(stor == 0.d0)stor=dd
       dstor=stor-stor0 ; tau=stor
       CALL tdetau(tau,teff,grav,t,dtsdtau,dtsdteff,dtsdg,
     1  ro_ext,dro_grav,dro_teff,f_tau,df_tau,d2f_tau)
       tau=stor0
       PRINT*,'dtsdtau df_tau d2f_tau'
       WRITE(*,2000)dtsdtau,(t-t0)/dstor,df_tau,(f_tau-f_tau0)/dstor,
     1  d2f_tau,(df_tau-df_tau0)/dstor

       PRINT*,'derivee/teff'
       stor0=teff ; stor=stor0*unpdd
       IF(stor == 0.d0)stor=dd
       dstor=stor-stor0 ; teff=stor
       CALL tdetau(tau,teff,grav,t,dtsdtau,dtsdteff,dtsdg,
     1 ro_ext,dro_grav,dro_teff,f_tau,df_tau,d2f_tau)
       teff=stor0
       PRINT*,'dtsdteff,dro_teff'
       WRITE(*,2000)dtsdteff,(t-t0)/dstor,dro_teff,(ro_ext-ro_ext0)/dstor

       PRINT*,'derivee/grav'
       stor0=grav ; stor=stor0*unpdd
       if(stor == 0.d0)stor=dd
       dstor=stor-stor0 ; grav=stor
       CALL tdetau(tau,teff,grav,t,dtsdtau,dtsdteff,dtsdg,
     1  ro_ext,dro_grav,dro_teff,f_tau,df_tau,d2f_tau)
       grav=stor0
       PRINT*,'dtsdg,dro_grav'
       WRITE(*,2000)dtsdg,(t-t0)/dstor,dro_grav,(ro_ext-ro_ext0)/dstor

       CALL pause('nouveau cas') ; k=k+1

       SELECT CASE(k)
       CASE(1)
        teff=5698.d0 ; grav=36150.d0 ; tau=20.d0 ; k=2
        CYCLE b1
       CASE(2)
        teff=4100.d0 ; grav=100.d0 ; tau=20.d0 ; k=3
        CYCLE b1
       CASE DEFAULT
        EXIT b1
       ENDSELECT
      ENDDO b1
	CALL pause('nouvelle loi T(tau)') ; PRINT*

      RETURN

      END SUBROUTINE calc_tdetau

      END PROGRAM test_tdetau