c************************************************************************ SUBROUTINE etat_opal01(pp,tp,xchim,deriv,ro,drop,drot,drox, 1 u,dup,dut,dux,delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) c routine public du module mod_etat c interface de l'équation d'état OPAL avec CESAM2k c package OPAL_EOS c appel a MHD en cas de difficultés c Auteur: JP Marquès, Université de Coimbra c adaptation F95: P.Morel, Département J.D. Cassini, O.C.A. c CESAM2k c entrées : c p : pression c t : température c xchim : composition chimique c deriv=.FALSE. : évite le calcul de certaines dérivées c sorties : c ro : densité et derivées c u : énergie interne et dérivées c------------------------------------------------------------------ USE mod_donnees, only : aradia, f_eos, nchim, nom_chemin, z0 USE mod_kind IMPLICIT NONE REAL (kind=dp), INTENT(in), DIMENSION(:) :: xchim REAL (kind=dp), INTENT(in) :: pp, tt LOGICAL, INTENT(in) :: deriv REAL (kind=dp), INTENT(out) :: ro,drop,drot,drox,u,dup,dut,dux, 1 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1 INTEGER, PARAMETER :: mx=5, mv=10, nr=169, nt=191 REAL (kind=dp), PARAMETER, DIMENSION(mx) :: xa=(/ 0.d0, 0.2d0, 1 0.4d0, 0.6d0, 0.8d0 /) REAL (kind=dp), DIMENSION(mv) :: eos REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: Lxchim REAL (kind=dp), PARAMETER :: dx=1.d-3, unpdx=1.d0+dx REAL (kind=dp), SAVE :: aradias3 REAL (kind=dp) :: stor, stor0, dstor, x, ztab, t6, r, p12, p, t REAL (kind=dp) :: esact INTEGER, PARAMETER, DIMENSION(nr) :: nta=(/ (191,i=1,92),190,189, 1 188,187,186,185,184,174,(134,i=101,104),(133,i=105,107), 2 (132,i=108,109),98,92,(85,i=112,113),(85,i=112,113), 3 (77,i=114,115),71,(63,i=117,119),(59,i=120,121),53,51, 4 (46,i=124,125),(44,i=126,134),(38,i=135,137),(33,i=138,143), 5 (29,i=144,159),27,26,25,23,22,20,19,18,17,16 /) INTEGER, PARAMETER, DIMENSION(10) :: index =(/ i,i=1,10 /) INTEGER, PARAMETER :: iorder=9, irad=1 INTEGER :: i LOGICAL, SAVE :: init=.TRUE. LOGICAL :: ok CHARACTER (len=1), SAVE :: blank=' ' c------------------------------------------------------------------ 2000 FORMAT(8es10.3) p=pp ; t=tt IF(init)THEN init=.FALSE. aradias3=aradia/3.d0 WRITE(2,1) ; WRITE(*,1) 1 FORMAT(//,'-------------Equation d''état-----------------',/, 1 'OPAL2001, ou EFF si problème; dérivées numériques',/) ALLOCATE(Lxchim(nchim)) ENDIF c calculate density (ro) for fixed P. t6=t*1.d-6 ; p12=p*1.d-12 ; x=xchim(1) ; ztab=z0 ro=rhoofp(x,ztab,t6,p12,irad,ok) IF(.NOT.ok)THEN CALL etat_eff(p,t,xchim,ro,drop,drot,drox,u,dup,dut,dux, 1 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) RETURN ENDIF CALL esac(x,ztab,t6,ro,iorder,irad) ! calc EOS; use ro from rhoofp drop=r/p/eos(5) !1 / dp/dro drot=-r/t*eos(6)/eos(5) !- dp/dT / dp/dro u=eos(2)*1.d12 duro=-p*(eos(6)-1.d0)/(ro**2) dup=duro*drop !du/dT dro/dp dut=duro*drot+eos(8)*1.d6 !du/dro dro/dT + du/dt delta=eos(6)/eos(5) !Ki T/Ki ro gradad=1.d0/eos(8) ; cp=p/ro/t*eos(6)*(1.d0/eos(9)+delta) alfa=1.d0/eos(5) !1/Ki ro beta=1.d0-aradias3*t**4/p ; gamma1=eos(7) IF(deriv)THEN c dérivées / P stor0=p ; stor=stor0*unpdx ; IF(stor < dx)stor=dx dstor=stor-stor0 ; p=stor ; p12=p*1.d-12 r=rhoofp(x,ztab,t6,p12,irad,ok) ! calculate density (r) for fixed P. IF(.NOT.ok)THEN CALL etat_eff(p,t,xchim,ro,drop,drot,drox,u,dup,dut,dux, 1 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) RETURN ENDIF CALL esac(x,ztab,t6,r,iorder,irad) ! calc EOS; use r from rhoofp drop=(r-ro)/dstor ; dup=(eos(2)*1.d12-u)/dstor deltap=(eos(6)/eos(5)-delta)/dstor dgradadp=(1./eos(8)-gradad)/dstor dcpp=(p/r/t*eos(6)*(1.d0/eos(9)+eos(6)/eos(5))-cp)/dstor p=stor0 ; p12=p*1.d-12 c dérivées / T stor0=t ; stor=stor0*unpdx ; IF(stor < dx)stor=dx dstor=stor-stor0 ; t=stor ; t6=t*1.d-6 r=rhoofp (x,ztab,t6,p12,irad,ok) ! calculate density (r) for fixed P. IF(.NOT.ok)THEN CALL etat_eff(p,t,xchim,ro,drop,drot,drox,u,dup,dut,dux, 1 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) RETURN ENDIF CALL esac(x,ztab,t6,r,iorder,irad) ! calc EOS; use r from rhoofp drot=(r-ro)/dstor ; dut=(eos(2)*1.d12-u)/dstor deltat=(eos(6)/eos(5)-delta)/dstor dgradadt=(1./eos(8)-gradad)/dstor dcpt=(p/r/t*eos(6)*(1.d0/eos(9)+eos(6)/eos(5))-cp)/dstor t=stor0 ; t6=t*1.d-6 c dérivées / X stor0=xchim(1) ; stor=stor0*unpdx ; IF(stor < dx)stor=dx dstor=stor-stor0 ; x=stor r=rhoofp (x,ztab,t6,p12,irad,ok) ! calculate density (r) for fixed P. IF(.NOT.ok)THEN CALL etat_eff(p,t,xchim,ro,drop,drot,drox,u,dup,dut,dux, 1 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) RETURN ENDIF CALL esac(x,ztab,t6,r,iorder,irad) ! calc EOS; use r from rhoofp drox=(r-ro)/dstor ; dux=(eos(2)*1.d12-u)/dstor deltax=(eos(6)/eos(5)-delta)/dstor dgradadx=(1./eos(8)-gradad)/dstor dcpx=(p/r/t*eos(6)*(1.d0/eos(9)+eos(6)/eos(5))-cp)/dstor x=stor0 ENDIF RETURN c END SUBROUTINE etat_opal01 CONTAINS c************************************************************************* SUBROUTINE esac (xh,ztab,t6,r,iorder,irad) c..... The purpose of this SUBROUTINE is to interpolate c the equation of state and its derivatives in X, T6, density c izi=0 recalulate table indices to use; =1 keep previous c xh=hydrogen mass fraction c ztab is metal fraction of the EOSdata tables you are using. c included only for purpose of preventing mismatch c t6=T6=temperature in millions of degrees kelvin c r=rho=Rho=density(g/cm**3) c..... to use esac insert COMMON/eeos/ esact,eos(10) in the CALLing routine. c This COMMON contains the interpolated EOS values for the EOS c c..... eos(i) are obtained from a quadradic interpolation at c fixed T6 at three values of Rho; followed by quadratic c interpolation along T6. Results smoothed by mixing c overlapping quadratics. c definitions: c c T6 is the temperature in units of 10**6 K c c rho is the density in grams/cc c R=Rho/T6**3 c eos(1) is the pressure in megabars (10**12dyne/cm**2) c eos(2) is energy in 10**12 ergs/gm. Zero is zero T6 c eos(3) is the entropy in units of energy/T6 c eos(4) is the specific heat, dE/dT6 at constant V. c eos(5) is dlogP/dlogRho at constant T6. c Cox and Guil1 eq 9.82 c eos(6) is dlogP/dlogT6 at conxtant Rho. c Cox and Guil1 eq 9.81 c eos(7) is gamma1. Eqs. 9.88 Cox and Guili. c eos(8) is gamma2/(gaamma2-1). Eqs. 9.88 Cox and Guili c eos(9) is gamma3-1. Eqs 9.88 Cox and Guili c iorder sets maximum index for eos(i);i.e., iorder=1 c gives just the pressure c c irad IF=0 no radiation correction; IF=1 adds radiation c index(i),i=1,10 sets order in which the equation of state c variables are stored in eos(i). Above order corresponds c to block data statement: c data (index(i),i=1,10)/1,2,3,4,5,6,7,8,9,10/. c If you, for example, only want to RETURN gamma1: set iorder=1 c and set: data (index(i),i=1,10)/8,2,3,4,5,6,7,1,9,10/ c-------------------------------------------------------------------- IMPLICIT NONE CHARACTER, SAVE (len=1) :: blank=' ' INTEGER itime COMMON/lreadco/itime REAL (kind=dp) epl(mx,nt,nr),xx(mx) COMMON/eeeos/epl,xx REAL (kind=dp) q(4),h(4),xxh COMMON/aaeos/q,h,xxh INTEGER m,mf REAL (kind=dp) :: xz(mx,mv,nt,nr),t6list(nr,nt),rho(nr),t6a(nt),esk(nt,nr), 1 esk2(nt,nr),dfsx(mx),dfs(nt),dfsr(nr),xa(mx) COMMON/aeos/xz,t6list,rho,t6a,esk,esk2,dfsx,dfs,dfsr,m,mf,xa INTEGER :: iri(10),index(10),nta(nr) REAL (kind=dp) :: zz(mx) COMMON/beos/iri,index,nta,zz INTEGER l1,l2,l3,l4,k1,k2,k3,k4,ip,iq COMMON/bbeos/l1,l2,l3,l4,k1,k2,k3,k4,ip,iq REAL (kind=dp) :: esact,eos(mv) COMMON/eeos/esact,eos INTEGER iorder,irad,i,j,ilo,ihi,imd,ipu,iqu,mg,mh,mi,mf2,ms, 1 kf2,ir,it,is,iw,iv REAL (kind=dp) slr,quadeos,slt,moles,frac(7),xh,xxi,ztab,t6,r,ri,dixr, 1 p0,z,sum1,sum2,sum23,sum33,tmass,gmass,eground,fracz REAL (kind=dp), SAVE :: aprop/83.14511d0/ save c------------------------------------------------------------------ IF(iorder > 9 )THEN PRINT*,'iorder cannot exceed 9' ENDIF IF(irad /= 0 .AND. irad /= 1)THEN PRINT*,'Irad must be 0 or 1' ; STOP ENDIF xxi=xh ; ri=r ; slt=t6 ; slr=r IF(itime /= 12345678)THEN itime=12345678 DO i=1,10 DO j=1,10 IF (index(i) == j) iri(i)=j ENDDO ENDDO DO i=1,mx xx(i)=xa(i) ENDDO c..... read the data files c CALL readcoeos CALL r_opal_bin z=zz(1) IF(abs(ztab-z) > 0.00001d0)THEN WRITE(*,1)ztab,z ; STOP 1 FORMAT('requested z=',f10.6,' EOSdata is for z=',f10.6 ENDIF IF(z+xh-1.d-6 > 1.d0 ) GOTO 61 ENDIF c..... Determine T6,rho grid points to use in the c interpolation. IF(slt > t6a(1) .OR. slt < t6a(nt)) GOTO 62 IF(slr < rho(1) .OR. slr > rho(nr)) GOTO 62 c.......IF(izi == 0)THEN ! freeze table indices ifnot 0 ilo=2 ; ihi=mx 8 IF(ihi-ilo > 1)THEN imd=(ihi+ilo)/2 IF(xh <= xa(imd)+1.d-7)THEN ihi=imd ELSE ilo=imd ENDIF GOTO 8 ENDIF i=ihi ; mf=i-2 ; mg=i-1 ; mh=i ; mi=i+1; mf2=mi IF(xh < 1.d-6)THEN mf=1 ; mg=1; mh=1 ; mi=2 ; mf2=1 ENDIF IF(xh <= xa(2)+1.d-7 .OR. xh >= xa(mx-2)-1.d-7) mf2=mh c ilo=2 ihi=nr 12 IF(ihi-ilo > 1)THEN imd=(ihi+ilo)/2 IF(slr == rho(imd))THEN ihi=imd GOTO 13 ENDIF IF(slr <= rho(imd))THEN ihi=imd ELSE ilo=imd ENDIF GOTO 12 ENDIF 13 i=ihi l1=i-2 ; l2=i-1 ; l3=i ; l4=l3+1 ; iqu=3 IF(l4 > nr)iqu=2 ilo=nt ; ihi=2 11 IF(ilo-ihi > 1)THEN imd=(ihi+ilo)/2 IF(t6 == t6list(1,imd))THEN ilo=imd GOTO 14 ENDIF IF(t6 <= t6list(1,imd))THEN ihi=imd ELSE ilo=imd ENDIF GOTO 11 ENDIF 14 i=ilo ; k1=i-2 ; k2=i-1 ; k3=i ; k4=k3+1 ; ipu=3 IF(k4 > nt) ipu=2 IF(k3 == 0)THEN WRITE(*,2)ihi,ilo,imd 2 FORMAT('ihi, ilo, imd',3i5) ENDIF c.......ENDIF c check to determine IFinterpolation indices fall within c table boundaries. choose largest allowed size. sum1=0.0d0 ; sum2=0.0d0 ; sum23=0.0d0 ; sum33=0.0d0 DO m=mf,mf2 DO ir=l1,l1+1 DO it=k1,k1+1 sum1=sum1+xz(m,1,it,ir) ENDDO ENDDO DO ir=l1,l1+2 DO it=k1,k1+2 sum2=sum2+xz(m,1,it,ir) ENDDO ENDDO IF(ipu == 3)THEN DO ir=l1,l1+2 DO it=k1,k1+ipu sum23=sum23+xz(m,1,it,ir) ENDDO ENDDO ELSE sum23=2.d+30 ENDIF IF(iqu == 3)THEN DO ir=l1,l1+3 DO it=k1,k1+ipu sum33=sum33+xz(m,1,it,ir) ENDDO ENDDO ELSE sum33=2.d+30 ENDIF ENDDO iq=2 ip=2 IF(sum2 > 1.d+30)THEN IF(sum1 < 1.d+25 )THEN k1=k3-3 ; k2=k1+1 ; k3=k2+1 ; l1=l3-3 ; l2=l1+1 ; l3=l2+1 GOTO 15 ELSE GOTO 65 ENDIF ENDIF IF(sum23 < 1.d+30)ip=3 IF(sum33 < 1.d+30)iq=3 IF(t6 >= t6list(1,2)+1.d-7)ip=2 IF(slr <= rho(2)+1.d-7)iq=2 IF(l3 == nr .OR. k3 == nt)THEN iq=2 ; ip=2 ENDIF 15 CONTINUE DO iv=1,iorder DO m=mf,mf2 is=0 DO ir=l1,l1+iq DO it=k1,k1+ip epl(m,it,ir)=xz(m,iv,it,ir) ; is=1 ENDDO ENDDO ENDDO IF(zz(mg) /= zz(mf) .OR. zz(mh) /= zz(mf))THEN PRINT*,'Z does not match Z in EOSdata files you are using') STOP ENDIF IF(z /= zz(mf)) GOTO 66 is=0 ; iw=1 DO ir=l1,l1+iq DO it=k1,k1+ip IF(mf2 == 1)THEN esk(it,ir)=epl(mf,it,ir) ; GOTO 46 ENDIF esk(it,ir)=quadeos(is,iw,xh,epl(mf,it,ir),epl(mg,it,ir), 1 epl(mh,it,ir),xx(mf),xx(mg),xx(mh)) IF(esk(it,ir) > 1.d+20)THEN WRITE(*,'(" problem it ir,l3,k3,iq,ip=", 6i5)')it,ir,l3,k3,iq,ip WRITE(*,'(3e12.4)') (epl(ms,it,ir),ms=mf,mf+2) ENDIF is=1 46 CONTINUE ENDDO ENDDO IF(mi == mf2)THEN ! interpolate between quadratics is=0 ; iw=1 ; dixr=(xx(mh)-xh)*dfsx(mh) DO ir=l1,l1+iq DO it=k1,k1+ip esk2(it,ir)=quadeos(is,iw,xh,epl(mg,it,ir),epl(mh,it,ir), 1 epl(mi,it,ir),xx(mg),xx(mh),xx(mi)) IF(esk(it,ir) > 1.d+20)THEN WRITE(*,'(" problem it ir,l3,k3,iq,ip=", 6i5)') it,ir,l3,k3,iq,ip WRITE(*,'(3e12.4)') (epl(ms,it,ir),ms=mg,mg+2) ENDIF esk(it,ir)=esk(it,ir)*dixr+esk2(it,ir)*(1.d0-dixr) ; is=1 ENDDO ENDDO ENDIF is=0 c..... completed X interpolation. Now interpolate T6 and rho on a c 4x4 grid. (t6a(i),i=i1,i1+3),rho(j),j=j1,j1+3)).Procedure c mixes overlapping quadratics to obtain smoothed derivatives. CALL t6rinteos(slr,slt) ; eos(iv)=esact ENDDO p0=t6*r eos(iri(1))=eos(iri(1))*p0 ! interpolated in p/po eos(iri(2))=eos(iri(2))*t6 ! interpolated in E/T6 tmass=gmass(xh,z,moles,eground,fracz,frac) IF(irad == 1)THEN CALL radsub (t6,r,moles,tmass) ELSE eos(iri(4))=eos(iri(4))*moles*aprop/tmass ENDIF RETURN 61 WRITE(*,'(" Mass fractions exceed unity (61)")') STOP 62 WRITE(*,'(" T6/LogR outside of table range (62)")') STOP 65 WRITE(*,'("T6/log rho in empty region od table (65)")') WRITE (*,'("xh,t6,r=", 3e12.4)') xh,t6,r STOP 66 WRITE(*,'(" Z does not match Z in EOSdata* files you are using (66)")') WRITE (*,'("mf,zz(mf)=",i5,e12.4)') mf,zz(mf) WRITE (*,'(" iq,ip,k3,l3,xh,t6,r,z= ",4i5,4e12.4)')ip,iq,k3,l3,xh,t6,r,z STOP RETURN END SUBROUTINE esac c*********************************************************************** SUBROUTINE t6rinteos(slr,slt) c The purpose of this SUBROUTINE is to interpolate in T6 and rho IMPLICIT NONE INTEGER iu,is,kx,iw REAL (kind=dp) slr,slt,dix,dix2,esact2,esactq,esactq2,quadeos REAL (kind=dp) epl(mx,nt,nr),xx(mx) COMMON/eeeos/epl,xx REAL (kind=dp) q(4),h(4),xxh COMMON/aaeos/q,h,xxh INTEGER m,mf REAL (kind=dp) xz(mx,mv,nt,nr),t6list(nr,nt),rho(nr),t6a(nt),esk(nt,nr), 1 esk2(nt,nr),dfsx(mx),dfs(nt),dfsr(nr),xa(mx) COMMON/aeos/ xz,t6list,rho,t6a,esk,esk2,dfsx,dfs,dfsr,m,mf,xa INTEGER l1,l2,l3,l4,k1,k2,k3,k4,ip,iq COMMON/bbeos/l1,l2,l3,l4,k1,k2,k3,k4,ip,iq REAL (kind=dp) esact,eos(mv) COMMON/eeos/esact,eos save c----------------------------------------------------------------- iu=0 is=0 DO kx=k1,k1+ip iw=1 iu=iu+1 h(iu)=quadeos(is,iw,slr,esk(kx,l1),esk(kx,l2),esk(kx,l3),rho(l1), 1 rho(l2),rho(l3)) IF(iq. eq. 3)THEN iw=2 q(iu)=quadeos(is,iw,slr,esk(kx,l2),esk(kx,l3),esk(kx,l4),rho(l2), 1 rho(l3),rho(l4)) ENDIF is=1 ENDDO c is=0 iw=1 c..... eos(i) in lower-right 3x3(i=i1,i1+2 j=j1,j1+2) esact=quadeos(is,iw,slt,h(1),h(2),h(3),t6a(k1),t6a(k2),t6a(k3)) IF(iq == 3)THEN c..... eos(i) upper-right 3x3(i=i1+1,i1+3 j=j1,j1+2) esactq=quadeos(is,iw,slt,q(1),q(2),q(3),t6a(k1),t6a(k2),t6a(k3)) ENDIF IF(ip == 3)THEN c..... eos(i) in lower-left 3x3. esact2=quadeos(is,iw,slt,h(2),h(3),h(4),t6a(k2),t6a(k3),t6a(k4)) c..... eos(i) smoothed in left 3x4 dix=(t6a(k3)-slt)*dfs(k3) esact=esact*dix+esact2*(1.-dix) c ENDIF ! moved to loc a IF(iq == 3)THEN c..... eos(i) in upper-right 3x3. esactq2=quadeos(is,iw,slt,q(2),q(3),q(4),t6a(k2),t6a(k3),t6a(k4)) esactq=esactq*dix+esactq2*(1.-dix) ENDIF ENDIF ! loc a c IF(iq == 3)THEN dix2=(rho(l3)-slr)*dfsr(l3) IF(ip == 3)THEN c..... eos(i) smoothed in both log(T6) and log(R) esact=esact*dix2+esactq*(1-dix2) ENDIF ENDIF IF(esact > 1.d+15)THEN WRITE(*,'("Interpolation indices out of range; please report conditions.")') STOP ENDIF RETURN END SUBROUTINE t6rinteos c********************************************************************** SUBROUTINE readcoeos c..... The purpose of this SUBROUTINE is to read the data tables IMPLICIT NONE include 'modele.COMMON' INTEGER :: i,j,k,l,jcs,numtot,iv REAL (kind=dp) :: dum INTEGER :: itime COMMON/lreadco/itime REAL (kind=dp) q(4),h(4),xxh COMMON/aaeos/ q,h,xxh INTEGER m,mf REAL (kind=dp) xz(mx,mv,nt,nr),t6list(nr,nt),rho(nr),t6a(nt),esk(nt,nr), 1 esk2(nt,nr),dfsx(mx),dfs(nt),dfsr(nr),xa(mx) COMMON/aeos/ xz,t6list,rho,t6a,esk,esk2,dfsx,dfs,dfsr,m,mf,xa INTEGER iri(10),index(10),nta(nr) REAL (kind=dp) zz(mx) COMMON/beos/ iri,index,nta,zz REAL (kind=dp) esact,eos(mv) COMMON/eeos/esact,eos REAL (kind=dp) epl(mx,nt,nr),xx(mx) COMMON/eeeos/ epl,xx INTEGER icycuse(mx,nr) REAL (kind=dp) moles(mx),xin(mx),tmass(mx),rhogr(mx,nr),frac(mx,6), 1 alogr(nr,nt) COMMON/eeeeos/moles,xin,tmass,icycuse,rhogr,frac,alogr save external blk_eos IF(itime /= 12345678)THEN DO i=1,mx DO j=1,mv DO k=1,nt DO l=1,nr xz(i,j,k,l)=1.d+35 ENDDO ENDDO ENDDO ENDDO itime=12345678 ENDIF close (2) c..... read tables c.......CALL system (' gunzip EOSdata') open(60, FILE=f_eos(1)) DO 3 m=1,mx read (60,'(3x,f6.4,3x,f12.9,11x,f10.7,17x,f10.7)')xin(m),zz(m),moles(m), 1 tmass(m) read (60,'(21x,e14.7,4x,e14.7,3x,e11.4,3x,e11.4,3x,e11.4,4x,e11.4)') 1 (frac(m,i),i=1,6) read (60,'(a)') blank DO 2 jcs=1,nr read (60,'(2i5,2f12.7,17x,e15.7)')numtot,icycuse(m,jcs),dum,dum, 1 rhogr(m,jcs) IF(numtot /= jcs)THEN WRITE (*,'(" Data file incorrect: numtot,jcs= ",2i5)')numtot,jcs STOP ENDIF read(60,'(a)') blank read(60,'(a)') blank IF(icycuse(m,jcs) < nta(jcs))THEN WRITE (*,'("problem with data files: X=",f6.4," density=",e14.4)') 1 xin(m), rhogr(m,jcs) STOP ENDIF DO i=1,icycuse(m,jcs) IF(i > nta(jcs))THEN read (60,'(a)') blank GOTO 4 ENDIF read (60,'(f9.5,1x,f6.2,3e13.5,6f8.4)')t6list(jcs,i),alogr(jcs,i), 1 (xz(m,index(iv),i,jcs),iv=1,9) 4 CONTINUE ENDDO read(60,'(a)') blank read(60,'(a)') blank read(60,'(a)') blank 2 CONTINUE read(60,'(a)') blank 3 CONTINUE DO i=1,nt IF(t6list(1,i) == 0.0)THEN WRITE(*,'("READCOEOS: Error:",i4,"-th T6 value is zero")') i STOP ENDIF t6a(i)=t6list(1,i) ENDDO DO 12 i=2,nt 12 dfs(i)=1./(t6a(i)-t6a(i-1)) rho(1)=rhogr(1,1) DO 13 i=2,nr rho(i)=rhogr(1,i) 13 dfsr(i)=1./(rho(i)-rho(i-1)) DO i=2,mx dfsx(i)=1./(xx(i)-xx(i-1)) ENDDO c.......CALL system ('gzip EOSdata') CLOSE(60) RETURN END SUBROUTINE readcoeos c c********************************************************************** SUBROUTINE r_opal_bin c adaptation a CESAM de la routine readco du package OPAL_EOS c..... The purpose of this SUBROUTINE is to read the data tables IMPLICIT NONE include 'modele.COMMON' INTEGER mx,mv,nr,nt parameter (mx=5,mv=10,nr=169,nt=191) INTEGER i,j,k,l,ntuse INTEGER itime COMMON/lreadco/itime INTEGER m,mf REAL (kind=dp) xz(mx,mv,nt,nr),t6list(nr,nt),rho(nr),t6a(nt),esk(nt,nr), 1 esk2(nt,nr),dfsx(mx),dfs(nt),dfsr(nr),xa(mx) COMMON/aeos/xz,t6list,rho,t6a,esk,esk2,dfsx,dfs,dfsr,m,mf,xa INTEGER iri(10),index(10),nta(nr) REAL (kind=dp) zz(mx) COMMON/beos/iri,index,nta,zz REAL (kind=dp) epl(mx,nt,nr),xx(mx) COMMON/eeeos/epl,xx INTEGER icycuse(mx,nr) REAL (kind=dp) moles(mx),xin(mx),tmass(mx),rhogr(mx,nr),frac(mx,6), 1 alogr(nr,nt) COMMON/eeeeos/moles,xin,tmass,icycuse,rhogr,frac,alogr save external blk_eos IF(itime /= 12345678)THEN DO i=1,mx DO j=1,mv DO k=1,nt DO l=1,nr xz(i,j,k,l)=1.d+35 ENDDO ENDDO ENDDO ENDDO itime=12345678 ENDIF 2000 FORMAT(1x,1p8d10.3) WRITE(2,1)f_eos(1) WRITE(*,1)f_eos(1) 1 FORMAT(1x,'donnees prises dans le fichier binaire: ',a50) c CALL system('uncompress '//f_eos(1)) c print*,'decompression du fichier' close(unit=60) open(unit=60,file=f_eos(1),form='unFORMATted') print*,'lecture, des donnees EOS opal, et c''est long' read(60)xz,t6list,alogr,rhogr,xin,zz,moles,tmass,frac,icycuse close(unit=60) c print*,'recompression du fichier binaire EOS' c CALL system('compress '//f_eos(1)) c WRITE(*,2000)(zz(i),i=1,mx) c pause'lecture du binaire' DO i=1,nt IF(t6list(1,i) == 0.0)THEN ntuse=i ELSE t6a(i)=t6list(1,i) ENDIF ENDDO DO i=2,nt dfs(i)=1./(t6a(i)-t6a(i-1)) ENDDO rho(1)=rhogr(1,1) DO i=2,nr rho(i)=rhogr(1,i) dfsr(i)=1./(rho(i)-rho(i-1)) endDO DO i=2,mx dfsx(i)=1./(xx(i)-xx(i-1)) ENDDO print*,'fin de lecture des donnees EOS opal en binaire' RETURN END SUBROUTINE r_opal_bin c*********************************************************************** FUNCTION quadeos(ic,i,x,y1,y2,y3,x1,x2,x3) c..... this function performs a quadratic interpolation. IMPLICIT NONE INTEGER ic,i REAL (kind=dp) x,y1,y2,y3,x1,x2,x3,xx(3),yy(3),xx12(30),quadeos, 1 xx13(30),xx23(30),xx1sq(30),xx1pxx2(30),c1,c2,c3 save xx(1)=x1 xx(2)=x2 xx(3)=x3 yy(1)=y1 yy(2)=y2 yy(3)=y3 IF(ic == 0)THEN xx12(i)=1./(xx(1)-xx(2)) xx13(i)=1./(xx(1)-xx(3)) xx23(i)=1./(xx(2)-xx(3)) xx1sq(i)=xx(1)*xx(1) xx1pxx2(i)=xx(1)+xx(2) ENDIF c3=(yy(1)-yy(2))*xx12(i) c3=c3-(yy(2)-yy(3))*xx23(i) c3=c3*xx13(i) c2=(yy(1)-yy(2))*xx12(i)-(xx1pxx2(i))*c3 c1=yy(1)-xx(1)*c2-xx1sq(i)*c3 quadeos=c1+x*(c2+x*c3) RETURN END FUNCTION quadeos c*********************************************************************** SUBROUTINE RADSUB (t6,density,moles,tmass) IMPLICIT NONE REAL (kind=dp) t6,density,moles,tmass,k,molenak,Na,unitf,unitfold,c,sigma, 1 sigmac,sigmacc,aprop REAL (kind=dp) rat,pr,er,sr,revise,fixerror,st,chir,chitt,pt,et, 1 cvtt,gam3pt,gam1t,gam2pt,dedrhoat,gam3pt_norad,gam1t_norad, 2 gam2pt_norad INTEGER mx,mv,nr,nt parameter (mx=5,mv=10,nr=169,nt=191) REAL (kind=dp) esact,eos(mv) COMMON/eeos/esact,eos INTEGER iri(10),index(10),nta(nr) REAL (kind=dp) zz(mx) COMMON/beos/iri,index,nta,zz data Na/6.0221367e+23/,k/1.380658e-16/,unitf/0.9648530/,unitfold/0.965296/, 1 c/2.9979245e+10/,sigma/5.67051e-5/,sigmac/1.8914785e-15/, 2 sigmacc/1.8914785e-3/,aprop/83.14510/ save c Physical constants c Na=6.0221367e+23 c k =1.380658e-16 ! erg/degree K c Na*k=6.0221367E+23*1.380658e-16 erg/degree K=8.314510E+7 erg/degree K c =8.314510e+7*11604.5 erg/eV=0.9648575E+12 erg/eV c Define unitf= Na*k/e+12=0.9648575 c unitf=0.9648530 ! obtained with latest physical constants c unitfold=0.965296 ! old units- still used in the EOS code c In these units energy/density is in units of Mb-CC/gm c Pressure is in units of E+12 bars=Mb c sigma is the Stefan-Boltzmann constant c sigma=5.67051E-5 ! erg /(s*cm**2*K**4) c c=2.99792458E+10 ! cm/sec c rat=sigma/c ! dyne/(cm**2*K**4) c rat=rat*1.d+24 ! Convert degrees K to units 10**6 K (T replaced with T6) rat=sigmacc molenak=moles*aprop ! Mb-cc/unit T6 pr=4./3.*rat*t6**4 ! Mb er=3.*pr/density ! Mb-cc/gm sr=4./3.*er/t6 ! Mb-cc/(gm-unit T6) c----- Calculate EOS without radiation correction pt=eos(iri(1)) et=eos(iri(2)) st=eos(iri(3)) chir=eos(iri(5))*eos(iri(1))/pt chitt=(eos(iri(1))*eos(iri(6)))/pt cvtt=(eos(iri(4))*molenak/tmass) gam3pt_norad=pt*chitt/(cvtt*density*t6) gam1t_norad=chir+chitt*gam3pt_norad gam2pt_norad=gam1t_norad/gam3pt_norad c---- End no radiation calculation c---- Calculate EOS with radiation calculation pt=eos(iri(1))+pr et=eos(iri(2))+er st=eos(iri(3))+sr chir=eos(iri(5))*eos(iri(1))/pt chitt=(eos(iri(1))*eos(iri(6))+4.*pr)/pt c gam1t(jcs,i)=(p(jcs,i)*gam1(jcs,i)+4./3.*pr)/pt(jcs,i) c gam2pt(jcs,i)=(gam2p(jcs,i)*p(jcs,i)+4.*pr)/pt(jcs,i) c gam3pt(jcs,i)=gam1t(jcs,i)/gam2pt(jcs,i) cvtt=(eos(iri(4))*molenak/tmass+4.*er/t6) gam3pt=pt*chitt/(cvtt*density*t6) ! DIRECT gam1t=chir+chitt*gam3pt !eq 16.16 Landau_Lifshitz (Stat. Mech) ! DIRECT gam2pt=gam1t/gam3pt ! DIRECT c----- End Eos calculations with radiation c normalize cvt to 3/2 when gas is ideal,non-degenerate, c fully-ionized, and has no radiation correction c cvt=(eos(5)*molenak/tmass+4.*er/t6)/molenak eos(iri(1))=pt eos(iri(2))=et eos(iri(3))=st eos(iri(4))=cvtt eos(iri(5))=chir eos(iri(6))=chitt c----- Add difference between EOS with and without radiation. cvtt c calculation is not accurate enough to give accurate results using c eq. 16.16 Landau&Lifshitz (SEE line labeled DIRECT) eos(iri(7))=eos(iri(7))+gam1t-gam1t_norad eos(iri(8))=eos(iri(8))+gam2pt-gam2pt_norad eos(iri(9))=eos(iri(9))+gam3pt-gam3pt_norad RETURN END SUBROUTINE radsub c******************************************************************** FUNCTION rhoofp(x,ztab,t6,p,irad,conv) IMPLICIT NONE INTEGER mx,mv,nr,nt parameter (mx=5,mv=10,nr=169,nt=191) INTEGER itime COMMON/lreadco/itime INTEGER m,mf REAL (kind=dp) xz(mx,mv,nt,nr),t6list(nr,nt),rho(nr),t6a(nt),esk(nt,nr), 1 esk2(nt,nr),dfsx(mx),dfs(nt),dfsr(nr),xa(mx) COMMON/aeos/xz,t6list,rho,t6a,esk,esk2,dfsx,dfs,dfsr,m,mf,xa INTEGER iri(10),index(10),nta(nr) REAL (kind=dp) zz(mx) COMMON/beos/iri,index,nta,zz REAL (kind=dp) esact,eos COMMON/eeos/esact,eos(mv) INTEGER i,nra(nt),irad,ilo,ihi,imd,mlo,klo,icount REAL (kind=dp) sigmacc,rat,pr,t6,p,pnr,ztab,x,pmax,pmin,rhog1,rhog2,rhog3,p1, 1 p2,p3,xinit,rinit,tinit data sigmacc/1.8914785e-3/ REAL (kind=dp) rhoofp c-------------------------------------------------------------------- data (nra(i),i=1,nt)/16*169,168,167,166,165,2*164,163,2*162,161,160,2*159, 1 4*143,5*137,6*134,2*125,5*123,2*122,6*121,4*119,8*116,6*115,8*113, 2 7*111,6*110,34*109,107,104,40*100,10*99,98,97,96,95,94,93,92/ logical conv rat=sigmacc pr=0.0 IF(irad == 1) pr=4./3.*rat*t6**4 ! Mb pnr=p-pr IF(itime /= 12345678)THEN xinit=0.5 tinit=1. rinit=1.d-3 CALL esac (xinit,ztab,tinit,rinit,1,0) ENDIF ilo=2 ihi=mx 8 IF(ihi-ilo > 1)THEN imd=(ihi+ilo)/2 IF(x <= xa(imd)+1.d-7)THEN ihi=imd ELSE ilo=imd ENDIF GOTO 8 ENDIF mlo=ilo ilo=nt ihi=2 11 IF(ilo-ihi > 1)THEN imd=(ihi+ilo)/2 IF(t6 == t6list(1,imd))THEN ilo=imd GOTO 14 ENDIF IF(t6 <= t6list(1,imd))THEN ihi=imd ELSE ilo=imd ENDIF GOTO 11 ENDIF 14 klo=ilo pmax=xz(mlo,1,klo,nra(klo))*t6*rho(nra(klo)) pmin=xz(mlo,1,klo,1)*t6*rho(1) IF((pnr > 1.25*pmax) .OR. (pnr < pmin))THEN WRITE (*,'(" The requested pressure-temperature not in table")') c STOP WRITE (*,'("pnr, pmax,pmin=",3e14.4)') pnr,pmax,pmin ENDIF rhog1=rho(nra(klo))*pnr/pmax CALL esac (x,ztab,t6,rhog1,1,0) p1=eos(1) IF(p1 > pnr)THEN p2=p1 rhog2=rhog1 rhog1=0.2*rhog1 IF(rhog1 < 1.d-14) rhog1=1.d-14 CALL esac (x,ztab,t6,rhog1,1,0) p1=eos(1) ELSE rhog2=5.*rhog1 IF(rhog2 > rho(klo)) rhog2=rho(klo) CALL esac (x,ztab,t6,rhog2,1,0) p2=eos(1) ENDIF icount=0 1 CONTINUE icount=icount+1 rhog3=rhog1+(rhog2-rhog1)*(pnr-p1)/(p2-p1) CALL esac (x,ztab,t6,rhog3,1,0) p3=eos(1) IF(abs((p3-pnr)/pnr) < 1.d-5)THEN rhoofp=rhog3 conv=.TRUE. RETURN ENDIF IF(p3 > pnr)THEN rhog2=rhog3 p2=p3 IF(icount < 11) GOTO 1 WRITE(*,1000)t6*1.d6,p*1.d12,x,rhog2 1000 FORMAT(1x,/,'!ATTENTION pas de convergence dans rhoofp, t=',1pe10.3, 1 ' p=',1pe10.3,' X=',1pe10.3,' ro=',1pe10.3) conv=.FALSE. RETURN ELSE rhog1=rhog3 p1=p3 IF(icount < 11) GOTO 1 WRITE(*,1000)t6*1.d6,p*1.d12,x,rhog1 conv=.FALSE. RETURN ENDIF END FUNCTION rhoofp c*********************************************************************** c BLOCKDATA blk_eos c IMPLICIT NONE c INTEGER mx,mv,nr,nt c parameter (mx=5,mv=10,nr=169,nt=191) c INTEGER i c REAL (kind=dp) q(4),h(4),xxh c COMMON/aaeos/q,h,xxh c INTEGER m,mf c REAL (kind=dp) xz(mx,mv,nt,nr),t6list(nr,nt),rho(nr),t6a(nt),esk(nt,nr), c 1 esk2(nt,nr),dfsx(mx),dfs(nt),dfsr(nr),xa(mx) c COMMON/aeos/xz,t6list,rho,t6a,esk,esk2,dfsx,dfs,dfsr,m,mf,xa c INTEGER iri(10),index(10),nta(nr) c REAL (kind=dp) zz(mx) c COMMON/beos/iri,index,nta,zz c data (xa(i),i=1,mx)/0.0,0.2,0.4,0.6,0.8/ c data (index(i),i=1,10)/1,2,3,4,5,6,7,8,9,10/ c data (nta(i),i=1,nr)/92*191,190,189,188,187,186,185,184,174,4*134,3*133, c 1 2*132,98,92,2*85,2*77,71,3*63,2*59,53,51,2*46,9*44,3*38,6*33, c 2 16*29,27,26,25,23,22,20,19,18,17,16/ c END BLOCKDATA blk_eos END SUBROUTINE etat_opal01