c*********************************************************************** SUBROUTINE etat_irwin(pp,tt,xchim,deriv,ro,drop,drot,drox,u,dup,dut,dux, 1 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,entropie, 2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) c subroutine private du module mod_etat c interface avec le package libre d'équation d'état free_eos de A. Irwin c téléchargeable sur http://freeeos.sourceforge.net/ c Yveline LEBRETON - novembre 2005 c P.Morel version - novembre 2009 c entrées : c pp : pression c tt : température c xchim : composition chimique c deriv=.FALSE. : evite le calcul de certaines derivées c sorties : c ro : densité et dérivées c u : énergie interne et dérivées c---------------------------------------------------------------- USE mod_donnees, ONLY : aal27, ac12, ac13, afe56, ah, ahe3, ahe4, ah2, 1 amg24, ana23, ane20, ane22, an13, an14, an15, ao16, ao17, ao18, ap31, 2 aradia, asi28, as32, eps_ini, langue, nchim, nom_etat, nom_nuc USE mod_kind IMPLICIT NONE REAL (kind=dp), INTENT(in), DIMENSION(:) :: xchim REAL (kind=dp), INTENT(in) :: pp, tt REAL (kind=dp), INTENT(out) :: ro, drop, drot, drox, u, dup, dut, dux, 1 delta, deltap, deltat, deltax, cp, dcpp, dcpt, dcpx, entropie, 2 gradad, dgradadp, dgradadt, dgradadx, alfa, beta, gamma1 LOGICAL, INTENT(in) :: deriv REAL (kind=dp), PARAMETER :: dx=1.d-3, unpdx=1.d0+dx REAL (kind=dp), SAVE, DIMENSION(20) :: eps REAL (kind=dp), DIMENSION(3) :: degeneracy, 1 pressure, density, energy, enthalpy, entropy REAL (kind=dp), SAVE :: aradias3 REAL(kind=dp) :: cf, cp_out, dstor, eta, fhe2, fhe3, fh2, fl_out, 1 gamma1_out, gamma2, gamma3, gamma_e, gradad_out, h2plusrat, h2rat, 2 lambda, pl_in, pl_out, p_out, qf, qp, rho_out, rl_out, rmue, rtp, 3 tl_in, sf, st, stor, stor0, t_out, xmu1, xmu3 INTEGER, SAVE :: ifion, ifmodified, ifoption, kif_in, neps INTEGER :: iteration_count LOGICAL, SAVE :: init=.TRUE. LOGICAL :: ok NAMELIST/nl_eos/ifoption, ifmodified, ifion c---------------------------------------------------------------------- 2000 FORMAT(8es10.3) c------ INITIALISATIONS : c------ eps (abondances/mole) c eps(1) : H, eps(2): He, eps(3): C, eps(4): N, eps(5): O c eps(6) : Ne, eps(7): Na, eps(8): Mg, eps(9): Al, eps(10): Si c eps(11): P, eps(12): S, eps(13): Cl, eps(14):A, eps(15): Ca c eps(16): Ti, eps(17): Cr, eps(18): Mn, eps(19): Fe, eps(20): Ni c input parameters for free_eos subroutine : ifoption, ifmodified, ifion, c kif_in, eps, neps c recommended choice : ifoption=3, ifmodified=1, ifion=-2 corresponds to EOS1 c EOS1 constrained by OPAL and SCVH fits, ifion=-2 : all 295 ionization states c of the 20 elements treated in detail, more details are given in subroutine c free_eos IF(init)THEN init=.FALSE. ; aradias3=aradia/3.d0 c recherche d'un éventuel fichier opteos dans l'environnement qui permet c d'ajuster les options du calcul, voir le fichier etat_irwin_explik INQUIRE(file='opteos',exist=ok) IF(ok)THEN OPEN(unit=30,form='formatted',status='old',file='opteos') READ(30,nl_eos) ; CLOSE(unit=30) ELSE c les options du calcul sont déduites de nom_etat, du fichier de données c Ex: dans le fichier de données *.don, en codant NOM_ETAT='etat_irwin4' c on se placera dans le cas EOS4. c Extrait du README de la source téléchargeable: c The EOS1, EOS1a, EOS2, and EOS3 option suites respectively take about 80, c 30, 8, and 2 times as long to compute as EOS4 for the locus of points in a c solar model. They represent various compromises between speed and quality c (see paper). For solar vibrational frequencies or generating EOS tables to c be subsequently interpolated I recommend using EOS1. For most stellar work c involving calculations of luminosities and effective temperatures, EOS4 c provides excellent results when called directly from the stellar interiors' c code, but it might be worthwhile to try EOS1 or EOS1a for a few test cases c (especially along the LMS where differences are the maximum) to make sure c the changes are within the observational errors. Also, note that c interpolating EOS tables is a non-trivial task since the tables have to be c quite large to reduce interpolation errors to negligible proportions. IF(INDEX(nom_etat,'1') == 11)THEN ifoption=3 ; ifmodified=1 ; ifion=-2 !EOS1 PREFERRED ELSEIF(INDEX(nom_etat,'2') == 11)THEN ifoption=2 ; ifmodified=1 ; ifion=-1 !EOS2 ELSEIF(INDEX(nom_etat,'3') == 11)THEN ifoption=2 ; ifmodified=1 ; ifion=0 !EOS3 ELSEIF(INDEX(nom_etat,'4') == 11)THEN ifoption=1 ; ifmodified=101 ; ifion=0 !EOS4 + rapide ELSEIF(INDEX(nom_etat,'5') == 11)THEN ifoption=5 ; ifmodified=1 ; ifion=2 !fully ionized EOS5 ELSEIF(INDEX(nom_etat,'6') == 11)THEN ifoption=3 ; ifmodified=1 ; ifion=-1 !EOS1a ELSE ifoption=3 ; ifmodified=1 ; ifion=-2 !EOS1 PREFERRED ENDIF ENDIF WRITE(2,1)nom_etat,ifion,ifmodified,ifoption WRITE(*,1)nom_etat,ifion,ifmodified,ifoption 1 FORMAT(/,'équation d''état de Irwin/ Irwin''s equation of state:',a11,/, 1 'Paramètres: ifion=',i2,', ifmodified=',i3,', ifoption=',i2,/) c autres initialisations kif_in=1 ! ln p and tl are independant variables neps=20 eps(1:2)=eps_ini(1:2) eps(3:5)=eps_ini(6:8) eps(6)=eps_ini(9)+eps_ini(10) !on ajoute F à Ne eps(7:13)=eps_ini(11:17) eps(14)=eps_ini(18)+eps_ini(19) !on ajoute K à Ar eps(15)=eps_ini(20)+eps_ini(21) !on ajoute Sc à Ca eps(16)=eps_ini(22) eps(17)=eps_ini(23)+eps_ini(24) !on ajoute V à Cr eps(18:19)=eps_ini(25:26) eps(20)=eps_ini(27)+eps_ini(28) !on ajoute Co à Ni c write(88,*) 'eps initiaux', eps(1:20) ENDIF c initialisations fin---------------------------- c CALCULATION of current eps------------------------------- SELECT CASE(nom_nuc) CASE ('pp1') eps(1)=xchim(1)/ah !H eps(2)=(1.d0-xchim(1))/ahe4 !He CASE ('pp3') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He CASE ('ppcno9') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He eps(3)=xchim(4)/ac12+xchim(5)/ac13 !C eps(4)=xchim(6)/an14+xchim(7)/an15 !N eps(5)=xchim(8)/ao16+xchim(9)/ao17 !O CASE ('ppcno9Fe') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He eps(3)=xchim(4)/ac12+xchim(5)/ac13 !C eps(4)=xchim(6)/an14+xchim(7)/an15 !N eps(5)=xchim(8)/ao16+xchim(9)/ao17 !O eps(19)=xchim(11)/afe56 !Fe CASE ('ppcno3a9') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He eps(3)=xchim(4)/ac12+xchim(5)/ac13 !C eps(4)=xchim(6)/an14+xchim(7)/an15 !N eps(5)=xchim(8)/ao16+xchim(9)/ao17 !O eps(6)=xchim(11)/ane20 !Ne CASE('ppcno10','ppcno10K') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He eps(3)=xchim(5)/ac12+xchim(6)/ac13 !C eps(4)=xchim(7)/an14+xchim(8)/an15 !N eps(5)=xchim(9)/ao16+xchim(10)/ao17 !O CASE('ppcno10Fe','ppcno10BeBFe') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He eps(3)=xchim(5)/ac12+xchim(6)/ac13 !C eps(4)=xchim(7)/an14+xchim(8)/an15 !N eps(5)=xchim(9)/ao16+xchim(10)/ao17 !O eps(19)=xchim(11)/afe56 !Fe CASE('ppcno11') eps(1)=xchim(1)/ah+xchim(2)/ah2 !H eps(2)=xchim(3)/ahe3+xchim(4)/ahe4 !He eps(3)=xchim(6)/ac12+xchim(7)/ac13 !C eps(4)=xchim(8)/an14+xchim(9)/an15 !N eps(5)=xchim(10)/ao16+xchim(11)/ao17 !O CASE('ppcno12','ppcno12Be','ppcno12Li') eps(1)=xchim(1)/ah+xchim(2)/ah2 !H eps(2)=xchim(3)/ahe3+xchim(4)/ahe4 !He eps(3)=xchim(7)/ac12+xchim(8)/ac13 !C eps(4)=xchim(9)/an14+xchim(10)/an15 !N eps(5)=xchim(11)/ao16+xchim(12)/ao17 !O CASE('ppcno12BeBFe') eps(1)=xchim(1)/ah+xchim(2)/ah2 !H eps(2)=xchim(3)/ahe3+xchim(4)/ahe4 !He eps(3)=xchim(7)/ac12+xchim(8)/ac13 !C eps(4)=xchim(9)/an14+xchim(10)/an15 !N eps(5)=xchim(11)/ao16+xchim(12)/ao17 !O eps(19)=xchim(16)/afe56 !Fe CASE('ppcno3a12Ne') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He eps(3)=xchim(4)/ac12+xchim(5)/ac13 !C eps(4)=xchim(6)/an14+xchim(7)/an15 !N eps(5)=xchim(8)/ao16+xchim(9)/ao17+xchim(12)/ao18 !O eps(6)=xchim(10)/ane20+eps_ini(9) !Ne CASE('nrj31_YL') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He eps(3)=xchim(4)/ac12+xchim(5)/ac13 !C eps(4)=xchim(6)/an14+xchim(7)/an15 !N eps(5)=xchim(8)/ao16+xchim(9)/ao17+xchim(12)/ao18 !O eps(6)=xchim(10)/ane20+eps_ini(9) !Ne eps(7)=xchim(19)/ana23 !Na eps(8)=xchim(15)/amg24 !Mg eps(9)=xchim(20)/aal27 !Al eps(10)=xchim(21)/asi28 !Si eps(11)=xchim(22)/ap31 !P eps(12)=xchim(23)/as32 !S CASE('ppcno3aco') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He eps(3)=xchim(4)/ac12+xchim(5)/ac13 !C eps(4)=xchim(6)/an14+xchim(7)/an15 !N eps(5)=xchim(8)/ao16+xchim(9)/ao17 !O eps(6)=xchim(11)/ane20 !Ne eps(7)=xchim(12)/ana23 !Na eps(8)=xchim(13)/amg24 !Mg eps(9)=xchim(14)/aal27 !Al eps(10)=xchim(15)/asi28 !Si eps(11)=xchim(16)/ap31 !P eps(12)=xchim(17)/as32 !S CASE('ppcno3acos') eps(1)=xchim(1)/ah !H eps(2)=xchim(2)/ahe3+xchim(3)/ahe4 !He eps(3)=xchim(4)/ac12+xchim(5)/ac13 !C eps(4)=xchim(6)/an14+xchim(7)/an15 !N eps(5)=xchim(8)/ao16+xchim(9)/ao17+xchim(12)/ao18 !O eps(6)=xchim(10)/ane20+xchim(17)/ane22 !Ne eps(7)=xchim(16)/ana23 !Na eps(8)=xchim(13)/amg24 !Mg eps(9)=xchim(15)/aal27 !Al eps(10)=xchim(14)/asi28 !Si eps(11)=xchim(18)/ap31 !P eps(12)=xchim(19)/as32 !S CASE ('iben') eps(1)=xchim(1)/ah !H eps(2)=(1.d0-xchim(1))/ahe4 !He CASE DEFAULT SELECT CASE (langue) CASE('english') WRITE(2,1017)nom_nuc ; WRITE(*,1017)nom_nuc 1017 FORMAT('STOP, unknown routine of thermonuclear reactions : ',a,/, 1 'known routines :') CASE DEFAULT WRITE(2,17)nom_nuc ; WRITE(*,17)nom_nuc 17 FORMAT('ARRET, routine de réactions nucléaires inconnue : ',a,/, 1 'routines connues :') END SELECT WRITE(2,18) ; WRITE(*,18) ; STOP 18 FORMAT('pp1, pp3, ppcno9, ppcno10, ppcno10Fe, ppcno10K, ppcno11',/, 1 'ppcno12, ppcno12Be, ppcno12Li, ppcno3a9, ppcno3ac10') END SELECT c write(88,*) 'eps courants', eps(1:20) c APPEL de l'EOS de IRWIN-------------------------------------------------- pl_in=LOG(pp) tl_in=LOG(tt) CALL free_eos(ifoption, ifmodified, ifion, kif_in, eps, neps, 1 pl_in, tl_in, fl_out, t_out, rho_out, rl_out, p_out, pl_out, cf, 2 cp_out, qf, qp, sf, st, gradad_out, 3 rtp, rmue, fh2, fhe2, fhe3, xmu1, xmu3, eta, gamma1_out, gamma2, 4 gamma3, h2rat, h2plusrat, lambda, gamma_e, degeneracy, 5 pressure, density, energy, enthalpy, entropy, iteration_count) c----------------------------------------------------------------------------------- ro=EXP(density(1)) drop=ro/pp*density(2) !ro / P alfa drot=ro/tt*density(3) !- ro / T delta u=energy(1) dup=energy(2)/pp dut=energy(3)/tt gradad=gradad_out cp=cp_out delta=-density(3) alfa=density(2) beta=1.d0-aradias3*tt**4/pp gamma1=gamma1_out entropie=entropy(1) IF(deriv)THEN c dérivées / P stor0=pp ; stor=stor0*unpdx IF(stor < dx) stor=dx dstor=stor-stor0 ; pl_in=LOG(stor) CALL free_eos(ifoption, ifmodified, ifion, kif_in, eps, neps, 1 pl_in, tl_in, fl_out, t_out, rho_out, rl_out, p_out, pl_out, cf, 2 cp_out, qf, qp, sf, st, gradad_out, 3 rtp, rmue, fh2, fhe2, fhe3, xmu1, xmu3, eta, gamma1_out, gamma2, 4 gamma3, h2rat, h2plusrat, lambda, gamma_e, degeneracy, 5 pressure, density, energy, enthalpy, entropy, iteration_count) c drop=(dEXP(density(1))-ro)/dstor c dup=(energy(1)-u)/dstor deltap=(-density(3)-delta)/dstor dgradadp=(gradad_out-gradad)/dstor dcpp=(cp_out-cp)/dstor pl_in=LOG(pp) c dérivées / T stor0=tt ; stor=stor0*unpdx IF(stor < dx)stor=dx dstor=stor-stor0 ; tl_in=LOG(stor) CALL free_eos(ifoption, ifmodified, ifion, kif_in, eps, neps, 1 pl_in, tl_in, fl_out, t_out, rho_out, rl_out, p_out, pl_out, cf, 2 cp_out, qf, qp, sf, st, gradad_out, 3 rtp, rmue, fh2, fhe2, fhe3, xmu1, xmu3, eta, gamma1_out, gamma2, 4 gamma3, h2rat, h2plusrat, lambda, gamma_e, degeneracy, 5 pressure, density, energy, enthalpy, entropy, iteration_count) c drot=(dEXP(density(1))-ro)/dstor c dut=(energy(1)-u)/dstor deltat=(-density(3)-delta)/dstor dgradadt=(gradad_out-gradad)/dstor dcpt=(cp_out-cp)/dstor tl_in=LOG(tt) c dérivées / X stor0=xchim(1) stor=stor0*unpdx IF(stor < dx)stor=dx dstor=(stor-stor0) stor0=eps(1) eps(1)=stor/ah CALL free_eos(ifoption, ifmodified, ifion, kif_in, eps, neps, 1 pl_in, tl_in, fl_out, t_out, rho_out, rl_out, p_out, pl_out, cf, 2 cp_out, qf, qp, sf, st, gradad_out, 3 rtp, rmue, fh2, fhe2, fhe3, xmu1, xmu3, eta, gamma1_out, gamma2, 4 gamma3, h2rat, h2plusrat, lambda, gamma_e, degeneracy, 5 pressure, density, energy, enthalpy, entropy, iteration_count) drox=(EXP(density(1))-ro)/dstor dux=(energy(1)-u)/dstor deltax=(-density(3)-delta)/dstor dgradadx=(gradad_out-gradad)/dstor dcpx=(cp_out-cp)/dstor eps(1)=stor0 RETURN !C bon avec deriv ENDIF !deriv c write(88,*) 'IRWIN' c write(88,*) 'rho, drop, drot,drox',ro,drop,drot,drox c write(88,*) 'u, dup, dut,dux',u,dup,dut,dux c write(88,*) 'delta,deltap,deltat,deltax',delta,deltap,deltat,deltax c write(88,*) 'gradad, dgradadp, dgradadt,dgradadx',gradad, dgradadp, c 1 dgradadt,dgradadx c write(88,*) 'cp,dcpp,cpt,cpx',cp,dcpp,dcpt,dcpx c write(88,*) 'alpha,beta,gamma1',alfa,beta,gamma1 c CALL etat_opalZ(pp,tt,xchim,.true., c 1 ro,drop,drot,drox,u,dup,dut,dux, c 2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, c 3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) c write(88,*) 'OPALZ' c write(88,*) 'rho, drop, drot,drox',ro,drop,drot,drox c write(88,*) 'u, dup, dut,dux',u,dup,dut,dux c write(88,*) 'delta,deltap,deltat,deltax',delta,deltap,deltat,deltax c write(88,*) 'gradad, dgradadp, dgradadt,dgradadx',gradad, dgradadp, c 1 dgradadt,dgradadx c write(88,*) 'cp,dcpp,cpt,cpx',cp,dcpp,dcpt,dcpx c write(88,*) 'alpha,beta,gamma1',alfa,beta,gamma1 c Lxchim(1:nchim)=xchim(1:nchim) c CALL etat_eff(pp,tt,Lxchim, c 1 ro,drop,drot,drox,u,dup,dut,dux, c 2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, c 3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) c write(88,*) 'EFF' c write(88,*) 'rho, drop, drot,drox',ro,drop,drot,drox c write(88,*) 'u, dup, dut,dux',u,dup,dut,dux c write(88,*) 'delta,deltap,deltat,deltax',delta,deltap,deltat,deltax c write(88,*) 'gradad, dgradadp, dgradadt,dgradadx',gradad, dgradadp, c 1 dgradadt,dgradadx c write(88,*) 'cp,dcpp,cpt,cpx',cp,dcpp,dcpt,dcpx c write(88,*) 'alpha,beta,gamma1',alfa,beta,gamma1 c CALL etat_ceff(pp,tt,Lxchim, c 1 ro,drop,drot,drox,u,dup,dut,dux, c 2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, c 3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) c write(88,*) 'CEFF' c write(88,*) 'rho, drop, drot,drox',ro,drop,drot,drox c write(88,*) 'u, dup, dut,dux',u,dup,dut,dux c write(88,*) 'delta,deltap,deltat,deltax',delta,deltap,deltat,deltax c write(88,*) 'gradad, dgradadp, dgradadt,dgradadx',gradad, dgradadp, c 1 dgradadt,dgradadx c write(88,*) 'cp,dcpp,cpt,cpx',cp,dcpp,dcpt,dcpx c write(88,*) 'alpha,beta,gamma1',alfa,beta,gamma1 c STOP RETURN END SUBROUTINE etat_irwin