c INCLUDE '../SOURCE/mod_static.f' c****************************************************************** PROGRAM test_thermo c Auteur: P. Morel, Département Cassiopée, O.C.A. c CESAM2k c pour tester les dérivées/Xi il faut garder z constant c dlpp est le rapport des dérivées spatiales dlnPgas / dlnPtot c--------------------------------------------------------------------- USE mod_donnees, ONLY : cpturb, lit_nl, nchim, nom_fich2, 1 pturb USE mod_kind USE mod_nuc, ONLY: nuc c USE mod_static, ONLY : thermo USE mod_variables, ONLY : chim_gram, lim, r_ov, r_zc, wrot IMPLICIT NONE REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: jac REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: xchim, dxchim, ex, ex0 REAL (kind=dp), DIMENSION(5) :: epsilo, epsilo0 REAL (kind=dp), PARAMETER :: dd=1.d-6, unpdd=1.d0+dd REAL (kind=dp) :: be7, b8,et, ero, f17, hh, n13, o15, w_rot REAL (kind=dp) :: pt,p,t,m,l,r,dlpp, 1 ro,drop,drot,drox,u,dup,dut,dux,grad,dgradpt,dgradp, 2 dgradt,dgradx,dgradm,dgradl,dgradr,dgradlpp, 3 gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp, 4 depsp,depst,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm, 8 gradrad,alfa,beta,gamma1,dstor,stor,stor0 REAL (kind=dp) :: ro0,drop0,drot0,drox0,u0,dup0,dut0,dux0,grad0, 1 dgradpt0,dgradp0,dgradt0, 2 dgradx0,dgradm0,dgradl0,dgradr0,dgradlpp0, 3 gam0,dgampt0,dgamp0,dgamt0,dgamx0,dgamm0,dgaml0,dgamr0,dgamlpp0, 4 depsp0,depst0,kap0,dkapp0,dkapt0,dkapx0, 5 delta0,deltap0,deltat0,deltax0,cp0,dcpp0,dcpt0,dcpx0, 6 gradad0,dgradadp0,dgradadt0,dgradadx0, 7 hp0,dhppt0,dhpp0,dhpt0,dhpx0,dhpr0,dhpm0 INTEGER :: ktest, i LOGICAL :: radiatif c------------------------------------------------------------------------ 2000 FORMAT(8es10.3) c lecture des données nom_fich2='test' ; CALL lit_nl(w_rot) c appel d'initialisation pour tabulation des réactions nucléaires c allocations fictives ALLOCATE(xchim(0),dxchim(0),jac(0,0),ex(0)) CALL nuc(1.5d+07,1.5d+02,xchim,dxchim,jac,.FALSE.,0, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c détermination des abondances initiales DEALLOCATE(xchim,dxchim,ex) ALLOCATE(xchim(nchim),dxchim(nchim),ex(nchim),ex0(nchim)) CALL nuc(1.5d+07,1.5d+02,xchim,dxchim,jac,.FALSE.,1, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c Fictif dxchim : dérivée de xchim(nchim) / nu (m^2/3) c non calculé ici dxchim=1.d0 ktest=1 c ktest=8 B2: DO SELECT CASE(ktest) CASE(1) pturb=.FALSE. ; cpturb=1.d0 p=1.754D+05 ; pt=p ; t=6.554D+03 ; m=1.0d0 l=0.7d0 ; r=0.8831d0 ; dlpp=1.d0 CASE(2) pturb=.FALSE. ; cpturb=1.d0 t=10.d0**(3.83644); p=10.d0**(5.34627) ; pt=p ; dlpp=1.d0 m=1.d0 ; l=.7024d0 ; r=.8781*.99911 CASE(3) pturb=.FALSE. ; cpturb=1.d0 t=3.5d6 ; p=828473414360767.4d0 ; pt=p ; dlpp=1.d0 m=.95d0 ; l=1.d0 ; r=.6d0 CASE(4) pturb=.TRUE. ; cpturb=-1.d0 t=6.554D+03 ; p=1.754D+05 ; pt=p*1.1 ; dlpp=0.9D-01 m=1.0d0 ; l=0.7d0 ; r=0.8831d0 CASE(5) pturb=.FALSE. ; cpturb=1.d0 t=6.184E+06 ; p=2.711E+14 ; pt=p ; dlpp=1.d0 m=0.17731*2.830d0 ; l=71.08*62.52d0 ; r=4.511E-02*11.074d0 CASE(6) pturb=.FALSE. ; cpturb=1.d0 t=4.050E+04 ; p=1.567E+05 ; pt=1.567E+05 ; dlpp=1.d0 m=2.850E+00 ; l=9.945E-01*4.358E+01 ; r=9.947E-01*2.5119E+00 CASE(7) pturb=.FALSE. ; cpturb=1.d0 t=8.994E+06 ; p=4.498E+16 ; pt=4.498E+16 ; dlpp=1.d0 m=3.053E-01 ; l=6.405E-01 ; r=2.016E-01 CASE(8) pturb=.FALSE. ; cpturb=1.d0 t=8.962E+06 ; p=4.434E+16 ; pt=4.434E+16 ; dlpp=1.d0 m=2.987E-01 ; l=6.193E-01 ; r=2.018E-01 c CASE() c pturb=.FALSE. ; cpturb=1.d0 c t= ; p= ; pt= ; dlpp= c m= ; l= ; r= CASE DEFAULT EXIT B2 END SELECT ktest=ktest+1 CALL thermo(pt,p,t,m,l,r,dlpp,xchim,dxchim, 1 ro0,drop0,drot0,drox0,u0,dup0,dut0,dux0,grad0,dgradpt0,dgradp0,dgradt0, 2 dgradx0,dgradm0,dgradl0,dgradr0,dgradlpp0, 3 gam0,dgampt0,dgamp0,dgamt0,dgamx0,dgamm0,dgaml0,dgamr0,dgamlpp0, 4 epsilo0,depsp0,depst0,ex0,kap0,dkapp0,dkapt0,dkapx0, 5 delta0,deltap0,deltat0,deltax0,cp0,dcpp0,dcpt0,dcpx0, 6 gradad0,dgradadp0,dgradadt0,dgradadx0, 7 hp0,dhppt0,dhpp0,dhpt0,dhpx0,dhpr0,dhpm0, 8 gradrad,alfa,beta,gamma1,radiatif) WRITE(*,*)'dérivées ro, u, grad, eps, k, delta, cp / grad(l,m,r)' IF(radiatif)THEN WRITE(*,50)grad0,gradad0 50 FORMAT('zone RADIATIVE',/,'grad=',es10.3,', gradad=',es10.3) ELSE WRITE(*,51)grad0,gradad0 51 FORMAT('zone CONVECTIVE',/,'grad=',es10.3,', gradad=',es10.3) ENDIF WRITE(*,*)'pt,p,t,xchim(1:nbelem),m,l,r,dlpp' WRITE(*,2000)pt,p,t,xchim(1:nchim),m,l,r,dlpp WRITE(*,*)'log10(t),log10(p),log10(ro),log10(kappa)' WRITE(*,2000)log10(t),log10(p),log10(ro0),log10(kap0) WRITE(*,*)'ro,drop,drot,drox,u,dup,dut,dux' WRITE(*,2000)ro0,drop0,drot0,drox0,u0,dup0,dut0,dux0 WRITE(*,*)'grad,dgradp,dgradt,dgradl,dgradr,dgradm,dgradx,dgradlpp' WRITE(*,2000)grad0,dgradp0,dgradt0,dgradl0,dgradr0,dgradm0, 1 dgradx0,dgradlpp0 WRITE(*,*)'gam,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp' WRITE(*,2000)gam0,dgamp0,dgamt0,dgamx0,dgamm0,dgaml0,dgamr0, 1 dgamlpp0 WRITE(*,*)'hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm' WRITE(*,2000)hp0,dhppt0,dhpp0,dhpt0,dhpx0,dhpr0,dhpm0 WRITE(*,*)'epsilon,depsp,depst,depsx' WRITE(*,2000)epsilo0(1),depsp0,depst0,ex0(1) WRITE(*,*)'kap,dkapp,dkapt,dkapx' WRITE(*,2000)kap0,dkapp0,dkapt0,dkapx0 WRITE(*,*)'gradad,dgradadp,dgradadt,dgradadx' WRITE(*,2000)gradad0,dgradadp0,dgradadt0,dgradadx0 WRITE(*,*)'cp,dcpp,dcpt,dcpx' WRITE(*,2000)cp0,dcpp0,dcpt0,dcpx0 WRITE(*,*)'alfa,gradrad,gamma1' WRITE(*,2000)alfa,gradrad,gamma1 IF(.NOT.pturb)THEN WRITE(*,9)dgradpt0,dgradp0,dgampt0,dgamp0,dhppt0,dhpp0 9 FORMAT(/,'sans Pturb on doit avoir :' ,/, 1 'dgradpt0=',es10.3,', dgradpt=',es10.3,', dgampt0=',es10.3, 2 ', dgamp0=',es10.3,/,'dhppt0=',es10.3,', dhpp0=',es10.3,/) ENDIF WRITE(*,1)pt,p,t,xchim(1:nchim),m,l,r,dlpp 1 FORMAT('données',/,'pt,p,t,xchim(1:nchim),m,l,r,dlpp',/, 1 8es10.3,/) IF(pturb)THEN WRITE(*,2)pt,dgradpt0*pt/grad0 2 FORMAT('dérivées par rapport à pt=',es10.3,', dgrad/dlnPt=', 1 es10.3) stor0=pt ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 pt=stor CALL thermo(pt,p,t,m,l,r,dlpp,xchim,dxchim, 1 ro,drop,drot,drox,u,dup,dut,dux,grad,dgradpt,dgradp,dgradt, 2 dgradx,dgradm,dgradl,dgradr,dgradlpp, 3 gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp, 4 epsilo,depsp,depst,ex,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm, 8 gradrad,alfa,beta,gamma1,radiatif) pt=stor0 IF(radiatif)THEN WRITE(*,50)grad0,gradad0 ELSE WRITE(*,51)grad0,gradad0 ENDIF WRITE(*,*)'dgradpt dgamdpt dhppt' WRITE(*,2000)dgradpt0,(grad-grad0)/dstor,dgampt0,(gam-gam0)/dstor, 1 dhppt0,(hp-hp0)/dstor PAUSE'dérivées par rapport à Pt' ENDIF WRITE(*,*) WRITE(*,3)p,dgradp0*p/grad0 3 FORMAT('dérivées par rapport à P=',es10.3,', dgrad/dlnP=',es10.3) stor0=p ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 p=stor CALL thermo(pt,p,t,m,l,r,dlpp,xchim,dxchim, 1 ro,drop,drot,drox,u,dup,dut,dux,grad,dgradpt,dgradp,dgradt, 2 dgradx,dgradm,dgradl,dgradr,dgradlpp, 3 gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp, 4 epsilo,depsp,depst,ex,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm, 8 gradrad,alfa,beta,gamma1,radiatif) p=stor0 IF(radiatif)THEN WRITE(*,50)grad0,gradad0 ELSE WRITE(*,51)grad0,gradad0 ENDIF WRITE(*,*)'drop dup dgradp dgamdp' WRITE(*,2000)drop0,(ro-ro0)/dstor,dup0,(u-u0)/dstor, 1 dgradp0,(grad-grad0)/dstor,dgamp0,(gam-gam0)/dstor WRITE(*,*)'depsp dkapp deltap dcpp' WRITE(*,2000)depsp0,(epsilo(1)-epsilo0(1))/dstor,dkapp0, 1 (kap-kap0)/dstor,deltap0,(delta-delta0)/dstor, 2 dcpp0,(cp-cp0)/dstor WRITE(*,*)'dgradadp dhpp' WRITE(*,2000)dgradadp0,(gradad-gradad0)/dstor,dhpp0,(hp-hp0)/dstor PAUSE'dérivées par rapport à P' WRITE(*,4)t 4 FORMAT('dérivées par rapport à T=',es10.3) stor0=t ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 t=stor CALL thermo(pt,p,t,m,l,r,dlpp,xchim,dxchim, 1 ro,drop,drot,drox,u,dup,dut,dux,grad,dgradpt,dgradp,dgradt, 2 dgradx,dgradm,dgradl,dgradr,dgradlpp, 3 gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp, 4 epsilo,depsp,depst,ex,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm, 8 gradrad,alfa,beta,gamma1,radiatif) t=stor0 IF(radiatif)THEN WRITE(*,50)grad0,gradad0 ELSE WRITE(*,51)grad0,gradad0 ENDIF WRITE(*,*)'drot dut dgradt dgamt' WRITE(*,2000)drot0,(ro-ro0)/dstor,dut0,(u-u0)/dstor, 1 dgradt0,(grad-grad0)/dstor,dgamt0,(gam-gam0)/dstor WRITE(*,*)'depst dkapt deltat dcpt' WRITE(*,2000)depst0,(epsilo(1)-epsilo0(1))/dstor,dkapt0, 1 (kap-kap0)/dstor,deltat0,(delta-delta0)/dstor, 2 dcpt0,(cp-cp0)/dstor WRITE(*,*)'dgradadt dhpt' WRITE(*,2000)dgradadt0,(gradad-gradad0)/dstor,dhpt0,(hp-hp0)/dstor PAUSE'dérivées par rapport à T' B1: DO i=1,nchim WRITE(*,*) WRITE(*,5)xchim(i) 5 FORMAT('dérivées par rapport à xchim(i)',es10.3) stor0=xchim(i) stor=stor0*unpdd IF(stor == 0.)stor=dd dstor=stor-stor0 xchim(i)=stor CALL thermo(pt,p,t,m,l,r,dlpp,xchim,dxchim, 1 ro,drop,drot,drox,u,dup,dut,dux,grad,dgradpt,dgradp,dgradt, 2 dgradx,dgradm,dgradl,dgradr,dgradlpp, 3 gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp, 4 epsilo,depsp,depst,ex,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm, 8 gradrad,alfa,beta,gamma1,radiatif) xchim(i)=stor0 IF(radiatif)THEN WRITE(*,50)grad0,gradad0 ELSE WRITE(*,51)grad0,gradad0 ENDIF IF(i == 1)THEN WRITE(*,*)'drox dux dgradx dgamx' WRITE(*,2000)drox0,(ro-ro0)/dstor,dux0,(u-u0)/dstor, 1 dgradx0,(grad-grad0)/dstor,dgamx0,(gam-gam0)/dstor WRITE(*,*)'depsx dkapx deltax dcpx' WRITE(*,2000)ex0(i),(epsilo(1)-epsilo0(1))/dstor, 1 dkapx0,(kap-kap0)/dstor, 2 deltax0,(delta-delta0)/dstor,dcpx0,(cp-cp0)/dstor WRITE(*,*)'dgradadx dhpx' WRITE(*,2000)dgradadx0,(gradad-gradad0)/dstor,dhpx0, 1 (hp-hp0)/dstor c pause'dérivées par rapport à X' EXIT B1 ELSE WRITE(*,*)'depsx' WRITE(*,2000)ex(i),(epsilo(1)-epsilo0(1))/dstor ENDIF ENDDO B1 PRINT*,'ces dérivées ne peuvent être justes que si Z est fixé' PAUSE'dérivées par rapport aux Xi' WRITE(*,6)m 6 FORMAT('dérivées par rapport à m=',es10.3) stor0=m ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 m=stor CALL thermo(pt,p,t,m,l,r,dlpp,xchim,dxchim, 1 ro,drop,drot,drox,u,dup,dut,dux,grad,dgradpt,dgradp,dgradt, 2 dgradx,dgradm,dgradl,dgradr,dgradlpp, 3 gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp, 4 epsilo,depsp,depst,ex,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm, 8 gradrad,alfa,beta,gamma1,radiatif) m=stor0 IF(radiatif)THEN WRITE(*,50)grad0,gradad0 ELSE WRITE(*,51)grad0,gradad0 ENDIF WRITE(*,*)'dgradm dgamm dhpm' WRITE(*,2000)dgradm0,(grad-grad0)/dstor,dgamm0,(gam-gam0)/dstor, 1 dhpm0,(hp-hp0)/dstor PAUSE'dérivées par rapport à m' WRITE(*,7)r 7 FORMAT('dérivées par rapport à r',es10.3) stor0=r ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 r=stor CALL thermo(pt,p,t,m,l,r,dlpp,xchim,dxchim, 1 ro,drop,drot,drox,u,dup,dut,dux,grad,dgradpt,dgradp,dgradt, 2 dgradx,dgradm,dgradl,dgradr,dgradlpp, 3 gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp, 4 epsilo,depsp,depst,ex,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm, 8 gradrad,alfa,beta,gamma1,radiatif) r=stor0 IF(radiatif)THEN WRITE(*,50)grad0,gradad0 ELSE WRITE(*,51)grad0,gradad0 ENDIF WRITE(*,*)'dgradr dgamr dhpr' WRITE(*,2000)dgradr0,(grad-grad0)/dstor,dgamr0,(gam-gam0)/dstor, 1 dhpr0,(hp-hp0)/dstor PAUSE'dérivées par rapport à r' WRITE(*,8)dlpp 8 FORMAT('dérivées par rapport à dlpp',es10.3) stor0=dlpp ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 dlpp=stor CALL thermo(pt,p,t,m,l,r,dlpp,xchim,dxchim, 1 ro,drop,drot,drox,u,dup,dut,dux,grad,dgradpt,dgradp,dgradt, 2 dgradx,dgradm,dgradl,dgradr,dgradlpp, 3 gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp, 4 epsilo,depsp,depst,ex,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm, 8 gradrad,alfa,beta,gamma1,radiatif) dlpp=stor0 IF(radiatif)THEN WRITE(*,50)grad0,gradad0 ELSE WRITE(*,51)grad0,gradad0 ENDIF WRITE(*,*)'dgradlpp dgamlpp' WRITE(*,2000)dgradlpp0,(grad-grad0)/dstor,dgamlpp0, 1 (gam-gam0)/dstor PAUSE'dérivées par rapport à dlpp' ENDDO B2 STOP END PROGRAM test_thermo c********************************************************************* SUBROUTINE thermo(pt,p,t,m,l,r,dlpp,xchim,dxchim, 1 ro,drop,drot,drox,u,dup,dut,dux,grad,dgradpt,dgradp,dgradt, 2 dgradx,dgradm,dgradl,dgradr,dgradlpp, 3 gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp, 4 epsilon,depsp,depst,depsx,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm, 8 gradrad,alfa,beta,gamma1,radiatif) c subroutine private du module mod_static c calcul de la thermodynamique c Auteur: P.Morel, Département J.D. Cassini, c O.C.A., Observatoire de Nice c CESAM2k c MODIF: c 09 10 96 introduction de la gravité effective c dont on tient compte dans Hp c 10 09 97 correction de dgradr c 23 10 97 : pression turbulente c 09 11 99 : correction critère de Ledoux c 19 11 99 : suppression de nh1, nhe1, nhe2, lamb c 30 07 00 : introduction F95 c avec Purb, pt et p sont deux variables indépendantes c la pression turbulente affecte l'échelle de hauteur de pression c de façon différente que la pression gazeuse, c ces différences n'affectent que les dérivées c sans pression turbulente la pression totale est la pression gazeuse c et leurs dérivées sont identiques c entrées : c pt : pression totale Ptot c p : pression thermique+radiative Pgas c t : température K c xchim : composition chimique ATTENTION en 1/mole : *ah pour avoir X c dxchim : d xchim/d nu (nu=m^2/3) en 1/mole c m : masse/msol c l : luminosité/lsol c r : rayon / rsol c lim: nombre de limites ZR/ZC c mstar: masse au temps du calcul, avec perte de masse c dlpp : d lnPgas /d lnPtot c sortie : (drop : d ro /dp etc..) c ro : densité cgs c u : énergie interne c grad : gradient c epsilon : énergie nucléaire cgs c kap : opacité Rosseland cm2/gr c delta, cp, hp, gradad, gradrad : notations évidentes c radiatif=.true. : zone radiative c gam : gamma de la convection MLT, sert avec la pression turbulente c routines externes : c etat, opa, conv, nuc c----------------------------------------------------------------------- USE mod_conv, ONLY: conv USE mod_donnees, ONLY : alpha, aradia, clight, cpturb, g, grad_ovi, 1 grad_ovs, ihe4, Krot, ledoux, lsol, msol, mu_saha, m_ch, m_rot, 2 nchim, nrot, ne, nucleo, ovshts, ovshti, pi, pturb, rsol, tau_max, zi USE mod_etat, ONLY : etat, saha USE mod_kind USE mod_opa, ONLY : opa USE mod_nuc, ONLY: nuc USE mod_numerique, ONLY : bsp1dn, no_croiss USE mod_variables, ONLY : chim_gram, knotc, knotr, lim, mc, mct, 1 mrot, mrott, n_ch, n_rot, rota, r_ov, r_zc, wrot IMPLICIT NONE REAL (kind=dp), INTENT(in), DIMENSION(nchim) :: xchim, dxchim REAL (kind=dp), INTENT(in) :: pt ,p, t, m, l, r, dlpp REAL (kind=dp), INTENT(out), DIMENSION(nchim) :: depsx REAL (kind=dp), INTENT(out), DIMENSION(5) :: epsilon REAL (kind=dp), INTENT(out) :: alfa, beta, cp, dcpp, dcpt, dcpx, 1 delta, deltap, deltat, deltax, depsp, depst, dgaml, 2 dgamlpp, dgamm, dgamp, dgampt, dgamr, dgamt, dgamx, 3 dgradadp, dgradadt, dgradadx, dgradl, dgradlpp, dgradm, dgradpt, 4 dgradp, dgradr, dgradt, dgradx, dhppt, dhpp, dhpt, 5 dhpx, dhpr, dhpm, dkapp, dkapt, dkapx,drop, drot, drox, dup, 6 dut, dux, gam, gamma1, grad, gradad, gradrad, hp, kap, ro, u LOGICAL, INTENT(out) :: radiatif REAL (kind=dp), DIMENSION (0,0) :: jac REAL (kind=dp), DIMENSION (nchim) :: dxchimm, xchimm, xchim0 REAL (kind=dp), DIMENSION (nrot) :: dfrot, frot REAL (kind=dp), DIMENSION (0) :: bidon REAL (kind=dp), SAVE :: cte1, cte13, cte2, cte8, cte9 REAL (kind=dp) :: be7, b8, depsro, dgamcp, dgamdel, dgamgad, dgamhp, 1 dgamkra, dgamgra, dgamgrad, dgamro, dgamtaur,dgradcp, 2 dgraddpp, dgraddpt, dgraddpx, dgradel, dgradgad, dgradgra, dgradgrad, 3 dgradhp, dgradkra, dgradlpp0, dgradl0, dgradpt0, dgradp0, dgradro, 4 dgradr0, dgradtaur, dgradt0, dgradm0, dgradx0, dkapro, 5 dkradp, dkradt, dkradx, dlnmu, dtaurm, dtaurp, dtaurpt, dtaurr, 6 dtaurt, dtaurx, graddp, dgravm, dgravr, dw, d_grad, f17, gravite, 7 grad_mu, hh, krad, ldx, mu, m23, nel, n13, o15, rot, taur, w, Zbar INTEGER :: i=1 LOGICAL, SAVE :: init=.TRUE. LOGICAL :: ovsht, der c------------------------------------------------------------------------ 2000 FORMAT(8es10.3) 2001 FORMAT(8es11.4) IF(init)THEN !initialisations init=.false. cte1=4.d0/3.d0*aradia*clight cte2=2.d0/3.d0*rsol cte8=lsol/4.d0/pi/rsol/rsol c cte8=cte9*lsol/msol !calcul de gradrad sans grav. eff. cte9=3.d0/16.d0/pi/aradia/clight/g cte13=g*msol/rsol/rsol WRITE(*,2) ; WRITE(2,2) 2 FORMAT(/,'---------THERMO----------',/) IF(ledoux)THEN WRITE(*,4) ; WRITE(2,4) 4 FORMAT('convection: critère de LEDOUX') ELSE WRITE(*,5) ; WRITE(2,5) 5 FORMAT('convection: critère de SCHWARZSCHILD') ENDIF c commentaires der=cpturb < 0.d0 IF(pturb)THEN WRITE(*,1)abs(cpturb) ; WRITE(2,1)abs(cpturb) 1 FORMAT('avec pression turbulente C=',es10.3) IF(der)THEN WRITE(*,*)'on tient compte de dln Pgaz/dln Ptot' WRITE(2,*)'on tient compte de dln Pgaz/dln Ptot' ELSE WRITE(*,*)'on ne tient pas compte de dln Pgaz/dln Ptot' WRITE(2,*)'on ne tient pas compte de dln Pgaz/dln Ptot' ENDIF ELSE WRITE(*,3) ; WRITE(2,3) 3 FORMAT('sans pression turbulente') ENDIF ENDIF !initialisations c composition chimique par gramme xchimm=ABS(xchim) ; dxchimm=dxchim ; CALL chim_gram(xchimm,dxchimm) c WRITE(*,*)'pt,p,t,xchim,nucleo,xchimm' c WRITE(*,2000)pt,p,t,xchim(1),nucleo(1),xchimm(1) c WRITE(*,2000)xchimm(1) c WRITE(*,2000)(xchimm(i),i=1,nchim) c PAUSE'1' m23=m**(2.d0/3.d0) c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w=wrot ; dw=0.d0 CASE(3,4,5) IF(m == 0.d0)THEN w=rota(1,1) ; dw=0.d0 ELSE CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,m_rot,knotr,.TRUE., 1 MIN(m23,mrot(n_rot)),i,frot,dfrot) IF(no_croiss)PRINT*,'Pb. en 1 dans thermo' w=frot(1) ; dw=dfrot(1)*2.d0/3.d0/m23 ENDIF END SELECT c pression gazeuse (Pgas + Prad) CALL etat(p,t,xchimm,.TRUE., 1 ro,drop,drot,drox,u,dup,dut,dux, 2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) c WRITE(*,2001)p,t,xchimm(1),ro,drox ; PAUSE'2' drox=drox*nucleo(1) !par mole dux=dux*nucleo(1) deltax=deltax*nucleo(1) dcpx=dcpx*nucleo(1) dgradadx=dgradadx*nucleo(1) c WRITE(*,*)'ro,drop,drot,drox,u,dup,dut,dux' c WRITE(*,2000)ro,drop,drot,drox,u,dup,dut,dux c PAUSE'4' xchim0=xchim CALL nuc(t,ro,xchim0,bidon,jac,.TRUE.,3, 1 epsilon,depst,depsro,depsx,hh,be7,b8,n13,o15,f17) depsp=depsro*drop depst=depsro*drot+depst depsx(1)=depsro*drox+depsx(1) CALL opa(xchimm,t,ro,kap,dkapt,dkapro,dkapx) dkapp=dkapro*drop dkapt=dkapro*drot+dkapt dkapx=dkapro*drox+dkapx*nucleo(1) !par mole c WRITE(*,*)'delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, c 1 gradad,dgradadp,dgradadt,dgradadx' c WRITE(*,2000)delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, c 1 gradad,dgradadp,dgradadt,dgradadx c PAUSE'5' c WRITE(*,*)'epsilon,depsp,depst,depsx,p,t,r,l,m' c WRITE(*,2000)depsp,depst,p,t,r,l,m,ro c WRITE(*,2000)epsilon ; WRITE(*,2000)depsx c WRITE(*,*)'kap,dkapp,dkapt,dkapx' c WRITE(*,2000)kap,dkapp,dkapt,dkapx c PAUSE'6' c diffusivite radiative krad=cte1/kap/ro*t**3 dkradp=krad*(-drop/ro-dkapp/kap) dkradt=krad*(-drot/ro-dkapt/kap+3.d0/t) dkradx=krad*(-drox/ro-dkapx/kap) IF(m*l*r /= 0.d0)THEN !hors du centre c gravité effective gravite=cte13*m/r**2 ; dgravr=-2.d0*gravite/r ; dgravm= gravite/m rot=cte2*w**2*r !gravite effective avec rotation gravite=gravite-rot !remarque de N. Audard dgravr=dgravr-rot/r ; dgravm=dgravm-2.d0*cte2*w*r*dw c échelle de hauteur de pression IF(pturb)THEN hp=pt/gravite/ro dhppt=hp/pt ; dhpp=-hp*drop/ro ELSE hp=p/gravite/ro dhpp=hp*(1.d0/pt-drop/ro) ; dhppt=dhpp ENDIF dhpt=-hp*drot/ro ; dhpx=-hp*drox/ro dhpr=-hp*dgravr/gravite ; dhpm=-hp*dgravm/gravite c gradient radiatif gradrad=cte8*l*hp/r**2/krad/t !gradient radiatif IF(pturb)THEN dgradpt=gradrad*dhppt/hp ; dgradp=gradrad*(dhpp/hp-dkradp/krad) ELSE dgradp=gradrad*(dhpp/hp-dkradp/krad) ; dgradpt=dgradp ENDIF dgradt=gradrad*(dhpt/hp-dkradt/krad-1.d0/t) dgradx=gradrad*(dhpx/hp-dkradx/krad) ; dgradm=gradrad*dhpm/hp dgradr=gradrad*(dhpr/hp-2.d0/r) ; dgradl=gradrad/l c épaisseur optique de la bulle taur=kap*ro*alpha*hp IF(pturb)THEN dtaurpt=taur*dhppt/hp ; dtaurp=taur*(dkapp/kap+drop/ro+dhpp/hp) ELSE dtaurp=taur*(dkapp/kap+drop/ro+dhpp/hp) ; dtaurpt=dtaurp ENDIF dtaurt=taur*(dkapt/kap+drot/ro+dhpt/hp) dtaurx=taur*(dkapx/kap+drox/ro+dhpx/hp) dtaurr=taur*dhpr/hp ; dtaurm=taur*dhpm/hp c WRITE(*,*)'krad,dkradp,dkradt,dkradx,gravite,dgravr,dgravm, c 1 hp,dhpp,dhpt,dhpx,dhpr,dhpm,gradrad,dgradp,dgradt,dgradx, c 2 dgradm,dgradr,dgradl (en dehors du centre)' c WRITE(*,2000)krad,dkradp,dkradt,dkradx,gravite,dgravr,dgravm, c 1 hp,dhpp,dhpt,dhpx,dhpr,dhpm,gradrad,dgradp,dgradt,dgradx, c 2 dgradm,dgradr,dgradl ELSE !au centre, mais approximativement gradrad=cte9*kap*epsilon(1)*pt/t**4 !au centre l/m ~ epsilon hp=1.d38 ; gravite=0.d0 IF(gradrad /= 0.d0)THEN IF(pturb)THEN dgradpt=gradrad/pt ; dgradp=gradrad*(dkapp/kap+depsp/epsilon(1)) ELSE dgradp=gradrad*(1.d0/p+dkapp/kap+depsp/epsilon(1)) dgradpt=dgradp ENDIF dgradt=gradrad*(dkapt/kap+depst/epsilon(1)-4.d0/t) dgradx=gradrad*(dkapx/kap+depsx(1)/epsilon(1)) ELSE dgradpt=0.d0 ; dgradp=0.d0 ; dgradt=0.d0 ; dgradx=0.d0 ENDIF dgradm=0.d0 ; dgradr=0.d0 ; dgradl=0.d0 ; dgradlpp=0.d0 c WRITE(*,*)'gradrad,dgradp,dgradt,dgradx (au centre)' c WRITE(*,2000)gradrad,dgradp,dgradt,dgradx ENDIF c différence des gradients d_grad=gradrad-gradad*dlpp !; PAUSE'7' c critère de Ledoux pour la convection, cas ou ro(P,T,X) c i.e. mu ne depend que de X, e.g. EOS MHD ou OPAL IF(ledoux .AND. ihe4 > 1)THEN CALL mu_mol(dxchim,hp,m23,r,ro,t,xchim,dlnmu,grad_mu,mu,nel,Zbar) ldx=beta/(4.d0-3.d0*beta)*grad_mu d_grad=d_grad+ldx ENDIF c test de convection radiatif=d_grad <= 0.d0 c overshoot IF(radiatif)THEN !zone radiative grad_rad < grad_ad c mise a 0 de gam et dérivées gam=0.d0 ; dgampt=0.d0 ; dgamp=0.d0 ; dgamt=0.d0 ; dgamx=0.d0 dgamm=0.d0 ; dgaml=0.d0 ; dgamr=0.d0 ; dgamlpp=0.d0 c overshoot, gradient adiabatique si grad_ov*=.TRUE., radiatif sinon ovsht=.false. !.true. : il y aura overshoot DO i=1,lim !est-on dans une zone d'overshoot? c si r_zc a été fixé et que la limite ZR/ZC s'est legerement deplacee c r est considere comme en dehors de la zone d'overshoot c on ajoute/retire arbitrairement 0.01 IF(r_ov(i) >= 0.d0)THEN IF(ovshts > 0.d0)ovsht= grad_ovs .AND. 1 (ovsht .OR. (r <= r_ov(i) .AND. r >= r_zc(i)-0.01d0)) IF(ovshti > 0.d0)ovsht=grad_ovi .AND. 1 (ovsht .OR. (r >= r_ov(i) .AND. r <= r_zc(i)+0.01d0)) ENDIF ENDDO !lim c ovsht=.false. !pour avoir grad=grad_rad dans zone ovsht c WRITE(*,*)ovsht,lim c WRITE(*,2000)d_grad,r,(r_zc(i),r_ov(i),i=1,lim) IF(ovsht)THEN !.true. : il y a penetration convective grad=gradad IF(pturb)THEN dgradpt=0.d0 ; dgradp=dgradadp ELSE dgradp=dgradadp ; dgradpt=dgradp ENDIF dgradt=dgradadt ; dgradx=dgradadx ; dgradm=0.d0 ; dgradr=0.d0 dgradl=0.d0 ; dgradlpp=0.d0 ELSE !pc grad=gradrad !pas de penetration convective dgradlpp=0.d0 !les dérivées de grad sont celles de gradrad ENDIF !pc ELSE !zone convective c IF(m*l*r == 0.d0 .or. t > 5.d5)THEN !ZC a m=0, ou profonde IF(m*l*r == 0.d0)THEN !ZC a m=0 grad=gradad IF(pturb)THEN dgradpt=0.d0 ; dgradp=dgradadp ELSE dgradp=dgradadp ; dgradpt=dgradp ENDIF dgradt=dgradadt dgradx=dgradadx ; dgradm=0.d0 ; dgradr=0.d0 dgradl=0.d0 ; dgradlpp=0.d0 c WRITE(*,*)'grad,dgradp,dgradt,dgradx (ZC au centre)' c WRITE(*,2000)grad,dgradp,dgradt,dgradx c zone convective ELSE !zone convective pour m /= 0 IF(pturb .AND. der)THEN !avec pression turbulente graddp=gradad*dlpp+ldx !gradad est * dlnPgaz/dlnPtot dgraddpp=dgradadp*dlpp ; dgraddpt=dgradadt*dlpp dgraddpx=dgradadx*dlpp ; dgradlpp0=gradad ELSE graddp=gradad+ldx ; dgraddpp=dgradadp ; dgraddpt=dgradadt dgraddpx=dgradadx ; dgradlpp0=0.d0 ENDIF c PRINT*,'krad,gravite,delta,cp,ro,hp,taur,gradrad,graddp' c WRITE(*,2000)krad,gravite,delta,cp,ro,hp,taur,gradrad,graddp c PAUSE'avant conv' CALL conv(r,krad,gravite,delta,cp,ro,hp,taur,gradrad,graddp, 1 .TRUE.,grad,dgradkra,dgradgra,dgradel,dgradcp,dgradro, 2 dgradhp,dgradtaur,dgradgrad,dgradgad, 3 gam,dgamkra,dgamgra,dgamdel,dgamcp,dgamro, 4 dgamhp,dgamtaur,dgamgrad,dgamgad) c PRINT*,'grad,gam' c WRITE(*,2000)grad,gam c PAUSE'apres conv' c pour le gradient de la convection dgradpt0=dgradpt ; dgradt0=dgradt ; dgradr0=dgradr dgradl0=dgradl ; dgradm0=dgradm ; dgradp0=dgradp ; dgradx0=dgradx IF(pturb)THEN dgradpt=dgradhp*dhppt +dgradtaur*dtaurpt +dgradgrad*dgradpt0 dgradp=dgradkra*dkradp + dgradel*deltap + dgradcp *dcpp 1 +dgradro* drop + dgradhp*dhpp +dgradtaur*dtaurp 2 +dgradgrad*dgradp0+dgradgad*dgraddpp ELSE dgradp=dgradkra*dkradp + dgradel*deltap + dgradcp *dcpp 1 +dgradro* drop + dgradhp*dhpp +dgradtaur*dtaurp 2 +dgradgrad*dgradp0+dgradgad*dgraddpp dgradpt=dgradp ENDIF c dgradp=dgradkra*dkradp + dgradel*deltap + dgradcp *dcpp c 1 +dgradro* drop + dgradhp*dhpp +dgradtaur*dtaurp c 2 +dgradgrad*dgradp0+dgradgad*dgraddpp dgradt=dgradkra*dkradt + dgradel*deltat + dgradcp*dcpt 1 +dgradro*drot + dgradhp*dhpt +dgradtaur*dtaurt 2 +dgradgrad*dgradt0+dgradgad*dgraddpt dgradx=dgradkra*dkradx +dgradel*deltax + dgradcp*dcpx 1 +dgradro*drox +dgradhp*dhpx +dgradtaur*dtaurx 2 +dgradgrad*dgradx0+dgradgad*dgraddpx dgradm=dgradgra*dgravm +dgradhp*dhpm +dgradtaur*dtaurm 1 +dgradgrad*dgradm0 dgradr=dgradgra*dgravr +dgradhp*dhpr +dgradtaur*dtaurr 1 +dgradgrad*dgradr0 dgradl=dgradgrad*dgradl0 dgradlpp= dgradgad*dgradlpp0 !dérivée/ dln Pgaz/dln Ptot c PRINT*,'dgradlpp,dgradgad,dgradlpp0' c WRITE(*,2000)dgradlpp,dgradgad,dgradlpp0 c pour le gamma de la convection IF(pturb)THEN dgampt=dgamhp*dhppt +dgamtaur*dtaurpt +dgamgrad*dgradpt0 dgamp=dgamkra*dkradp + dgamdel*deltap + dgamcp *dcpp 1 +dgamro* drop + dgamhp*dhpp +dgamtaur*dtaurp 2 +dgamgrad*dgradp0+dgamgad*dgraddpp ELSE dgamp=dgamkra*dkradp + dgamdel*deltap + dgamcp *dcpp 1 +dgamro* drop + dgamhp*dhpp +dgamtaur*dtaurp 2 +dgamgrad*dgradp0+dgamgad*dgraddpp dgampt=dgamp ENDIF c dgamp=dgamkra*dkradp + dgamdel*deltap + dgamcp *dcpp c 1 +dgamro* drop + dgamhp*dhpp +dgamtaur*dtaurp c 2 +dgamgrad*dgradp0+dgamgad*dgraddpp dgamt=dgamkra*dkradt + dgamdel*deltat + dgamcp*dcpt 1 +dgamro*drot + dgamhp*dhpt +dgamtaur*dtaurt 2 +dgamgrad*dgradt0+dgamgad*dgraddpt dgamx=dgamkra*dkradx +dgamdel*deltax + dgamcp*dcpx 1 +dgamro*drox +dgamhp*dhpx +dgamtaur*dtaurx 2 +dgamgrad*dgradx0+dgamgad*dgraddpx dgamm=dgamgra*dgravm +dgamhp*dhpm +dgamtaur*dtaurm 1 +dgamgrad*dgradm0 dgamr=dgamgra*dgravr +dgamhp*dhpr +dgamtaur*dtaurr 1 +dgamgrad*dgradr0 dgaml=dgamgrad*dgradl0 dgamlpp= dgamgad*dgradlpp0 !dérivée/ dln Pgaz/dln Ptot c WRITE(*,*)'gam,dgamkra,dgamgra,dgamdel,dgamcp,dgamro,dgamhp, c 1 dgamtaur,dgamgrad,dgamgad' c WRITE(*,2000)gam,dgamkra,dgamgra,dgamdel,dgamcp,dgamro, c 1 dgamhp,dgamtaur,dgamgrad,dgamgad c WRITE(*,*)'gam,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr,dgamlpp' c WRITE(*,2000)gam,dgamp,dgamt,dgamx,dgamm,dgaml, c 1 dgamr,dgamlpp c PAUSE'les gam' c WRITE(*,*)'zone convective' c WRITE(*,*)'grad,dgradkra,dgradgra,dgradel,dgradcp,dgradro, c 1 dgradhp,dgradgrad,dgradgad,dgradp,dgradt,dgradx,dgradm, c 2 dgradr,dgradl' c WRITE(*,2000)grad,dgradkra,dgradgra,dgradel,dgradcp,dgradro, c 1 dgradhp,dgradgrad,dgradgad,dgradp,dgradt,dgradx,dgradm, c 2 dgradr,dgradl ENDIF !en m=0 ENDIF c WRITE(*,*)'jpz,lim/xchim(1),d_grad,r,r_zc,r_ov',jpz,lim c WRITE(*,2000)grad,gradad,gradrad,d_grad,r,(r_zc(i), c 1 r_ov(i),i=1,lim) c PAUSE'8' RETURN END SUBROUTINE thermo c************************************************************************* SUBROUTINE mu_mol(dxchim,hp,nu,r,ro,t,xchim,dlnmu,grad_mu,mu,nel,Zbar) c subroutine private du module mod_static c calcul de diverses quantités reliées au poids moléculaire moyen c entrées : c dxchim : d xchim / d nu c hp : échelle hauteur de pression c nu : m^2/3 c r : r/Rsol c ro : densité c t : température c xchim : composition chimique par mole c sorties : c dlnmu : d ln mu / d nu c grad_mu : dln mu /d ln P c mu : poids moléculaire moyen c nel : nombre d'électrons libres c Zbar : cherge moyenne c Auteur: P.Morel, Département Cassiopée, O.C.A., Observatoire de Nice c CESAM2k c----------------------------------------------------------------- USE mod_donnees, ONLY : amu, ihe4, msol, mu_saha, nchim, pi, rsol, zi USE mod_etat, ONLY : saha USE mod_kind Use mod_variables, ONLY : chim_gram IMPLICIT NONE REAL (kind=dp), INTENT(in), DIMENSION(nchim) :: dxchim, xchim REAL (kind=dp), INTENT(in) :: hp, nu, r, ro, t REAL (kind=dp), INTENT(out) :: dlnmu, grad_mu, mu, nel, Zbar REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ioni REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: dxchimg, xchimg, 1 z_bar REAL (kind=dp), SAVE :: cte1 REAL (kind=dp) :: eta LOGICAL, SAVE :: init=.TRUE. c-------------------------------------------------------------------- 2000 FORMAT(8es10.3) IF(init)THEN !initialisations init=.FALSE. cte1=8.d0*pi*rsol**2/3.d0/msol IF(mu_saha)THEN ALLOCATE(ioni(0:NINT(MAXVAL(zi)),nchim),z_bar(nchim)) ELSE ALLOCATE(dxchimg(nchim),xchimg(nchim)) ENDIF ENDIF !initialisations c appel à Saha ou approximation X, Y, Z IF(mu_saha)THEN CALL saha(xchim,t,ro,ioni,z_bar,nel,eta) mu=1.d0/DOT_PRODUCT((1.d0+z_bar),xchim) dlnmu=-mu*DOT_PRODUCT((1.d0+z_bar),dxchim) Zbar=DOT_PRODUCT(z_bar,xchim)/SUM(xchim) ELSE dxchimg=dxchim ; xchimg=xchim ; CALL chim_gram(xchimg,dxchimg) mu=16.d0/(23.d0*xchimg(1)+3.d0*xchimg(ihe4)+9.d0) dlnmu=-mu/16.d0*(23.d0*dxchimg(1)+3.d0*dxchimg(ihe4)) Zbar=DOT_PRODUCT(zi,xchim)/SUM(xchim) nel=DOT_PRODUCT(zi,xchim)*ro/amu ENDIF IF(nu <= 0.d0)THEN grad_mu=0.d0 ELSE grad_mu=-cte1*r**2*ro/SQRT(ABS(nu))*dlnmu*hp ENDIF RETURN END SUBROUTINE mu_mol c END PROGRAM test_thermo