c**************************************************************** SUBROUTINE kappa_cond(xh,t,ro,kappa,dkapdt,dkapdr,dkapdx) c routine private du module mod_opacite c ajout de l'opacite conductive c adaptation de la partie concernée du code de Genève procurée par c Y. Lebreton et selon ses instructions c expression of Iben(1975 AP.j. 196,525 Appendix A) and numerical c evaluations of the derivatives with respect to density and c temperature assuming dlro=0.001 and dlt=0.0001. Derivatives c are calculated assuming vmye=2./(1.+x)=constante=nb. nucleons par c electron libre c Auteur: P. Morel, Département J.D. Cassini, O.C.A. c CESAM2k c entrées: c xh(1)=X : comp. chim. en fraction de masse c t : température K c ro : densité cgs c entrées/sorties : c kapa : opacité gr / cm2) c dkapdt : kappa / d t c dkapdr : kappa / d densité c dkapdx : kappa / d xchim(1) c-------------------------------------------------------------------------- USE mod_donnees, ONLY : ihe4, langue, nchim, nom_elem, pi, z0 USE mod_kind IMPLICIT NONE REAL (kind=dp), INTENT(in), DIMENSION(:) :: xh REAL (kind=dp), INTENT(in) :: t, ro REAL (kind=dp), INTENT(inout) :: kappa, dkapdt, dkapdr, dkapdx REAL (kind=dp), PARAMETER, DIMENSION(8) :: 1 zr=(/ 1.d0, 2.d0, 6.d0, 7.d0, 8.d0, 10.d0, 10.d0, 12.d0/), 2 a=(/1.d0, 4.d0, 12.d0, 14.d0, 16.d0, 20.d0, 22.d0, 24.d0/) REAL (kind=dp), DIMENSION(8) :: xx REAL (kind=dp), PARAMETER :: dx=0.00001d0, unpdx=1.d0+dx REAL (kind=dp) :: vmye, vmye_log, za, za_log, zb, xz2, xhh, 1 rol6, tl6, dellg, del, eta0, eta2, t1, ro1, 2 a1, a2, b1, b2, rnedne, flg, blamr2, alfa, zc, f, thxlg, 3 thx, penktl, ef, efm, gam, glg, vkcc, thy, thc, thylg, thclg, 4 vcond, vkchr, kappa_c, dkapdr_c, dkapdt_c, dkapdx_c, 5 store, store0, dstore INTEGER, SAVE :: j12=0, j14=0, j16=0, j20=0, j22=0, j24=0 INTEGER :: j, i, k LOGICAL, SAVE :: init=.TRUE. c---------------------------------------------------------------------- 2000 FORMAT(8es10.3) IF(nchim <= 3)RETURN !pas d'opacités conductives IF(init)THEN !indices des elements c WRITE(*,2000)xh(1:nchim),t,ro ; WRITE(*,*)nom_elem(1:nchim) init=.FALSE. DO j=3,nchim IF(nom_elem(j) == 'C12 ')THEN j12=j ELSEIF(nom_elem(j) == 'N14 ')THEN j14=j ELSEIF(nom_elem(j) == 'O16 ')THEN j16=j ELSEIF(nom_elem(j) == 'Ne20')THEN j20=j ELSEIF(nom_elem(j) == 'Ne22')THEN j22=j ELSEIF(nom_elem(j) == 'Mg24')THEN j24=j ENDIF ENDDO SELECT CASE(langue) CASE('english') WRITE(*,1001) ; WRITE(2,1001) 1001 FORMAT(/,'Conductives Opacities (Iben Apj 196,515,1975)') CASE DEFAULT WRITE(*,1) ; WRITE(2,1) 1 FORMAT(/,'Opacités conductives (Iben Apj 196,515,1975)') END SELECT ENDIF c test whether electron conduction is likely to be important c calculation of conductive opacity according to the analytic IF(ro <= 1.d-4 .OR. t <= 2.d4)RETURN c expression of Iben(1975 AP.j. 196,525 Appendix A) and numerical c evaluations of the derivatives with respect to density and c temperature and chemical composition X. ro1=ro ; t1=t ; xhh=xh(1) DO k=1,4 IF(k == 2)THEN store0=t1 store=store0*unpdx dstore=store-store0 t1=store ELSEIF(k == 3)THEN store0=ro1 store=store0*unpdx dstore=store-store0 ro1=store ELSEIF(k == 4)THEN store0=xh(1) store=max(1.d-5,store0*unpdx) dstore=store-store0 xhh=store ENDIF c valeurs de T6 et ro6 en Log10 tl6=log10(t1)-6.d0 rol6=log10(ro1)-6.d0 c abondances en fraction de masse de H, He4, C12, c N14, O16, Ne20, Ne22 et Mg24 xx=0.d0 xx(1)=xhh xx(2)=1.d0-(xhh+SUM(xh(ihe4+1:nchim))) IF(j12 >= 3)xx(3)=xh(j12) IF(j14 >= 3)xx(4)=xh(j14) IF(j16 >= 3)xx(5)=xh(j16) IF(j20 >= 3)xx(6)=xh(j20) IF(j22 >= 3)xx(7)=xh(j22) IF(j24 >= 3)xx(8)=xh(j24) c l'opacité conductive vmye=2.d0/(1.d0+xx(1)) ; vmye_log=log10(vmye) ; za=0.d0 ; zb=0.d0 DO i=1,8 xz2=xx(i)*zr(i)**2 ; za=za+xz2/(a(i)**(1.d0/3.d0)) zb=zb+xz2/a(i) ENDDO za_log=log10(za) IF(rol6 <= 0.3d0) THEN !Hubbart pour ro<2.d6 dellg=rol6+6.d0-1.5d0*tl6-vmye_log ; del=10.d0**(dellg) !A1 eta0=10.d0**(-0.52255d0+2.d0*dellg/3.d0) ; eta2=eta0**2 !A6 c logique de A2 a A5 a1=-3.29243d0+log10(del*(1.d0+0.02804d0*del)) !A3 b1=-4.80946d0+log10(del*del*(1.d0+9.376d0/eta2)) !A4 IF(dellg <= 0.645d0) THEN flg=-3.2862d0+log10(del*(1.d0+0.024417d0*del)) !A2 ELSEIF(dellg <= 2.d0)THEN flg=a1 !A3 ELSEIF(dellg <= 2.5d0) THEN flg=2.d0*a1*(2.5d0-dellg)+2.d0*b1*(dellg-2.d0) !A5 ELSE flg=b1 !A4 ENDIF c logique de A8 a A10 a2=log10(1.d0+0.021876d0*del) !A8 b2=log10(0.4d0*eta0*(1.d0+4.1124d0/eta2)) !A9 IF(dellg <= 1.5d0)THEN penktl=a2 ELSEIF(dellg <= 2.d0) THEN penktl=2.d0*a2*(2.d0-dellg)+2.d0*b2*(dellg-1.5d0) !A10 ELSE flg=b1 !A9 penktl=2.d0*a2*(2.d0-dellg)+2.d0*b2*(dellg-1.5d0) !à vérifier ENDIF c logique de A11, A12 IF(del <= 40.d0)THEN rnedne=1.0d0-0.01d0*del*(2.8966d0-0.034838d0*del) !A11 ELSE rnedne=(1.5d0/eta0)*(1.d0-0.8225d0/eta2) !A12 ENDIF c logique de A13 a A21 blamr2=9.24735d-3*10.d0**(dellg-0.5d0*tl6-penktl) !A7 1 *(vmye*zb+rnedne) !lamba ou lambda_barre alfa = log10(blamr2) ! corrigé IF(alfa <= -3.d0)THEN !pour H thxlg=1.048d0-0.124d0*alfa !A13 ELSEIF(alfa <= -1.d0)THEN thxlg=0.13d0-alfa*(0.745d0+0.105d0*alfa) !A14 ELSE thxlg=0.185d0-0.585d0*alfa !A15 corrigée ENDIF IF(alfa <= -3.d0)THEN !pour He4 thylg=0.937d0-0.111d0*alfa !A16 ELSEIF(alfa <= 0.d0)THEN thylg=0.24d0-alfa*(0.55d0+0.0689d0*alfa) !A17 ELSE thylg=0.24d0-0.6d0*alfa !A18 ENDIF IF(alfa <= -2.5d0)THEN !pour C thclg=1.27d0-0.1d0*alfa !A19 ELSEIF(alfa <= .5)THEN thclg=0.727d0-alfa*(0.511d0+0.0778d0*alfa) !A20 ELSE thclg=0.843d0-0.785d0*alfa !A21 ENDIF thx=10.d0**(thxlg) ; thy=10.d0**(thylg) ; thc=10.d0**(thclg) c zc is an effective abundance of C,N,O,Ne weighted by the c squares of the charges and normalised with respect to C since c electron conduction has been calculated by Hubbard+Lampe only c for pure H,pure He and pure C zc=(3.d0*xx(3)+3.5d0*xx(4)+4.d0*xx(5)+5.d0*xx(6)+4.54d0*xx(7)+ 1 6.d0*xx(8))/3.d0 !A23 vkchr=(xx(1)*thx+xx(2)*thy+zc*thc)*10.d0**(-tl6-flg) !A22 IF(rol6 <= 0.d0)THEN !Hubbart en dessous de ro=1.d6 vcond=vkchr ELSE !Hubbart + Canuto pour 1.d6< ro < 2.d6 ef=sqrt(1.d0+10.d0**(2.d0*(rol6-vmye_log)/3.d0)) - 1.!A24 gam=22.76d0*10.d0**(rol6/3.d0-tl6)*za !A25 efm=min(1.d0,0.5d0+log10(ef)) !A30 glg=(0.873d0-0.298d0*za_log+(0.333d0-0.168d0*za_log)*efm)* !A29 1 (1.d0-(1.d0+gam)**(-0.85d0)) vkcc=6.753d-8*10.d0**(2.d0*tl6-glg)*zb/(ef**2.5d0*(1.d0+ef))!A28 f=0.5d0*(1.d0-cos(pi*rol6/0.3d0)) !A32 vcond=10.d0**((1.d0-f)*log10(vkchr)+f*log10(vkcc)) !A31 ENDIF ELSE !Canuto pour ro >2.d6 ef=sqrt(1.d0+10.d0**(2.d0*(rol6-vmye_log)/3.d0)) - 1.d0 !A24 gam=22.76d0*10.d0**(rol6/3.d0-tl6)*za !A25 efm=min(1.d0,0.5d0+log10(ef)) !A30 glg=(0.873d0-0.298d0*za_log+(0.333d0-0.168d0*za_log)*efm)* !A29 1 (1.d0-(1.d0+gam)**(-0.85d0)) vkcc=6.753d-8*10.d0**(2.d0*tl6-glg)*zb/(ef**2.5d0*(1.d0+ef)) !A28 vcond=vkcc ENDIF IF(k == 1)THEN kappa_c=vcond ELSEIF(k == 2)THEN t1=t ; dkapdt_c=(vcond-kappa_c)/dstore ELSEIF(k == 3)THEN ro1=ro ; dkapdr_c=(vcond-kappa_c)/dstore ELSEIF(k == 4)THEN xhh=xh(1) ; dkapdx_c=(vcond-kappa_c)/dstore ENDIF ENDDO !k c WRITE(*,2000)kappa,kappa_c,t,ro,xh(1) dkapdr=(kappa**2*dkapdr_c+kappa_c**2*dkapdr)/(kappa+kappa_c)**2 dkapdt=(kappa**2*dkapdt_c+kappa_c**2*dkapdt)/(kappa+kappa_c)**2 dkapdx=(kappa**2*dkapdx_c+kappa_c**2*dkapdx)/(kappa+kappa_c)**2 kappa=kappa*kappa_c/(kappa+kappa_c) RETURN END SUBROUTINE kappa_cond