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

	PROGRAM des2k_opa

c dessins d'opacités le long d'un modèle

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

c	var: variables
c	itot: nombre de points
c	nglob: nombre de constantes global
c	nvar: nb. de variables du tableau var sans les éléments chimiques
c	nom_elem: nom des éléments du vecteur de comp. chim. généralisée
c	nchim: nombre d'éléments strictement chimiques

c	glob: variables globales du fichier mon_modele-ad.osc
c		glob(1)=mstar*msol
c		glob(2)=rtot*rsol
c		glob(3)=ltot*lsol
c		glob(4)=z0
c		glob(5)=x0
c		glob(6)=alpha
c		glob(7)=X dans ZC
c		glob(8)=Y dans ZC
c		glob(9)=d2p
c		glob(10)=d2ro
c		glob(11)=age
c		glob(12)=wrot vitesse de rotation globale
c		glob(13)=w_rot initial

c	var: variables locales utilisées

c	nvar=22 pour oscillations adiabatiques
c		var(1,i)=r*rsol
c	 	var(2,i)=log(m/mstar) -1.d38 au centre
c		var(3,i)=t
c		var(4,i)=Ptot
c		var(5,i)=ro
c		var(6,i)=gradient réel d ln T / d ln P
c		var(7,i)=l
c		var(8,i)=kap
c		var(9,i)=énergie thermo+gravifique
c		var(10,i)=grand Gammass
c		var(11,i)=gradient adiabatique
c		var(12,i)=delta
c		var(13,i)=cp
c		var(14,i)=mu elec.
c		var(15,i)=vaissala, 0 au centre
c	 	var(16,i)=vitesse angulaire, radian/sec
c	 	var(17,i)=d ln kappa / d ln T
c	 	var(18,i)=d ln kappa / d ln ro
c	 	var(19,i)=d epsilon(nuc) / d ln T
c	 	var(20,i)=d epsilon(nuc) / d ln ro
c		var(21,i)=Ptot / Pgas
c		var(22,i)=gradient radiatif

c	  	var(22+j,i)=xchim(j)*nucleo(j), j=1,nchim

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

	USE mod_donnees, ONLY : ihe4, nchim, nom_fich2, nom_elem,
	1 nom_opa, nom_output
	USE mod_exploit, ONLY : glob, lit_nl, read_ascii, var
	USE mod_kind
	USE mod_opa, ONLY : opa
	USE mod_variables, ONLY : rstar, wrot
	
	IMPLICIT NONE
	
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: xchim
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: ro, t, ychim
	REAL (kind=dp) :: kap, dkapt, dkapro, dkapx
	
	REAL (kind=sp), ALLOCATABLE, DIMENSION(:) :: absc, kappa, mass, ray	
	REAL (kind=sp) :: m_max, r_max, xmax, xmin, ymax, ymin
	
	INTEGER :: i, ic, itot, j, n_fich, nglob, nopa=1, nvar
	
	LOGICAL :: en_rayon, en_masse, en_q, ok
	
	CHARACTER (len=1) :: oui	
	CHARACTER (len=7) :: nabsc, device
	CHARACTER (len=17) :: htext		
	CHARACTER (len=80), DIMENSION(4) :: abid		
	CHARACTER (len=80) :: chaine, modele
	CHARACTER (len=100) :: titre
		
c-------------------------------------------------------------------

2000	FORMAT(8es10.3)

	PRINT*,'nom, sans extension, du fichier Ex: soleil'
	READ*,nom_fich2
	
	chaine=TRIM(nom_fich2)//'.don'

	INQUIRE(file=TRIM(chaine),exist=ok)
	IF(.NOT.ok)THEN
	 PRINT*,'Arrêt: fichier de données inconnu: ',TRIM(chaine)
	 STOP
	ENDIF
	
c lecture du fichier de données
	CALL lit_nl(wrot)
	
	PRINT*,'pour unfichier créé par:'
	PRINT*,'CESAM2k, Ex: soleil-ad.osc, entrer o'	
	PRINT*,'f037_2k, Ex: soleil_f037-ad.osc, entrer n '
	READ*,oui	
	
	IF(nom_output == 'osc_adia')THEN
	 IF(oui == 'o')THEN
	  chaine=TRIM(nom_fich2)//'-ad.osc'
	 ELSE
	  chaine=TRIM(nom_fich2)//'_037-ad.osc'	
	 ENDIF
	ELSEIF(nom_output == 'osc_nadia')THEN
	 IF(oui == 'o')THEN
	  chaine=TRIM(nom_fich2)//'-nad.osc'
	 ELSE
	  chaine=TRIM(nom_fich2)//'_037-nad.osc'	
	 ENDIF
	ELSEIF(nom_output == 'osc_invers')THEN	
	 IF(oui == 'o')THEN
	  chaine=TRIM(nom_fich2)//'-inv.osc'
	 ELSE
	  chaine=TRIM(nom_fich2)//'_037-inv.osc'	
	 ENDIF
	ELSE
	 chaine=nom_output
	 WRITE(*,4)TRIM(chaine) ; STOP  
4	 FORMAT('Arrêt, pour un dessin nom_output doit être',/,
	1 'osc_adia ou osc_nadia ou osc_invers; nom_output : ',a)
	ENDIF
	
	INQUIRE(file=TRIM(chaine),exist=ok)
	IF(ok)THEN
	 WRITE(*,5)TRIM(chaine)
5	 FORMAT('fichier ASCII utilisé: ',a)	 
	ELSE
	 WRITE(*,3)TRIM(chaine) ; STOP	
3	 FORMAT('Arrêt, fichier ASCII non trouvé: ',a)
	ENDIF
	
	CALL read_ascii(chaine,itot,nglob,nvar,abid)
	rstar=glob(2)
	modele='modèle '//TRIM(chaine)

	ALLOCATE(absc(itot),mass(itot),kappa(itot),ray(itot),ro(itot),
	1 t(itot),xchim(nchim,itot),ychim(nchim))
	DO i=1,itot
	 IF(var(2,i) > -1.d30)THEN
	  mass(i)=exp(var(2,i))
	 ELSE
	  mass(i)=0.
	 ENDIF
	 kappa(i)=var(8,i) ; ray(i)=var(1,i) ; ro(i)=var(5,i) ; t(i)=var(3,i)
	 DO j=1,nchim
	  xchim(j,i)=MAX(var(nvar+j,i),1.d-30)
	 ENDDO	 
	ENDDO
	mass=mass/MAXVAL(mass) ; ray=ray/MAXVAL(ray)

c indice de He4
	DO i=1,nchim
	 IF(nom_elem(i) == 'He4 ')ihe4=i
	ENDDO
	
c écarts d'opacité	
	titre=TRIM(nom_opa)
	WRITE(*,1)
1	FORMAT('entrer le nom de la routine d''opacité, Exemple:',/,
	1 'opa_yveline')
	READ*,nom_opa	
	titre='1 - '//TRIM(nom_opa)//' / '//TRIM(titre)
	DO i=1,itot
	 ychim(:)=xchim(:,i)
	 CALL opa(ychim,t(i),ro(i),kap,dkapt,dkapro,dkapx)
	 kappa(i)=1.-kap/kappa(i)
	ENDDO

	ymax=MAXVAL(kappa) ; ymin=MINVAL(kappa)
		
c type de graphe
	en_rayon=.FALSE. ; en_masse=.FALSE. ; en_q=.FALSE.	
	PRINT*,'dessin en rayon o/n ?'
	READ(*,'(a)')oui
	en_rayon= oui == 'o'	
	IF(.NOT.en_rayon)THEN	
	 PRINT*,'dessin en masse? o/n'
	 READ(*,'(a)')oui
	 en_masse= oui == 'o'
	ELSE
	 en_q=.TRUE.
	ENDIF

c cas en rayon	
	IF(en_rayon)THEN
	 PRINT*,'dessin en rayon'; nabsc='r/Rstar'
	 absc=ray
	 
	 PRINT*,'Zoom? o/n'
	 READ(*,'(a)')oui
	 IF(oui == 'o')THEN
	  PRINT*,'entrer rmin, rmax' ; READ*,xmin,xmax
	  WRITE(*,991)xmin,xmax
991	  FORMAT('dessin entre les rayons',8es10.3)
	 ELSE
	  xmin=-.05 ; xmax=1.05
	 ENDIF
	 
c cas en masse
	ELSEIF(en_masse)THEN
	 PRINT*,'dessin en masse' ; nabsc='m/Mstar'
	 absc=mass
	 PRINT*,'Zoom? o/n' ; READ(*,'(a)')oui
	 IF(oui == 'o')THEN
	  PRINT*,'entrer m min, m max?' ; READ*,xmin,xmax
	  WRITE(*,992)xmin,xmax
992	  FORMAT('dessin entre les masses',8es10.3)
	 ELSE
	  xmin=-.05 ; xmax=1.05
	 ENDIF

	ENDIF
	
	device='/xw'
c	WRITE(6,*)'device ? /xw, /ps, /cps'
c	READ(*,'(a)')device
	CALL pgbegin(0,device,1,1)
	CALL pgscf(2)		!roman font
	CALL pgsch(1.2) ; CALL pgslw(3)

c tracé de d kappa / kappa
c	WRITE(*,2000)kappa(1:itot:50)
c	WRITE(*,2000)absc(1:itot:50)
	CALL pgsls(1) ; ic=1 ; CALL pgsci(ic)
	CALL pgenv(xmin,xmax,ymin,ymax,0,0)
	CALL pglabel(nabsc,titre,modele)
	ic=2 ; CALL pgsci(ic)	
 	CALL pgline(itot,absc,kappa)
	
	CALL pgend
		
	STOP

	END PROGRAM des2k_opa