c************************************************************* SUBROUTINE coeff_rota4(dt,nui,y,frl,dfrl,coeff) c routine private du module mod_evol c les coefficients pour la rotation avec Krot=4 c entrées : c dt : pas temporel c nu : abscisse en m^2/3 c y : variables de la rotation c sorties : c frl, dfrl : coefficients de rotation et dérivées c deff, dh, dv : coefficients de diffusion du moment cinétique c coeff : OPTIONAL, sortie des coefficients pour dessins c Auteur: A. Moya Institut d'Astrophysique d'Andalousie, c adaptation cesam2k : P. Morel, Département Cassiopée, O.C.A., CESAM2k c NOTATIONS (hélas incohérentes) pour les développements sur B-splines c n_ch : nombre VARIABLE de points élément de mod_variables c nch : nombre FIXE de fonctions élément de mod_donnees c m_ch : ordre FIXE des splines élément de mod_donnees c mch(n_ch) : abscisses VARIABLES élément de mod_variables c dans une ZC sont seuls utiles c frl(25, 26, 27, 28, 29, 30, 31, 32, 33, 34), Deff, Dv, Dh c--------------------------------------------------------------------- USE mod_donnees, ONLY : amu, aradia, clight, diffusion, d_conv, 1 en_masse, g, ihe4, lsol, m_ch, m_ptm, m_rot, msol, mu_saha, ne, nchim, 2 nom_rot, nucleo, nrl, nrot, ord_qs, phi, pi, rsol, secon6, x_tams, 3 zi, z0 USE mod_etat, ONLY: etat, saha USE mod_kind USE mod_nuc, ONLY : lw_planet, mzc_ext, nuc, nuzc_ext, planetoides USE mod_numerique, ONLY : bsp1dn, no_croiss USE mod_opa, ONLY : opa USE mod_variables, ONLY : bp, bp_t, chim, chim_t, chim_gram, 1 inter, knot, knotc, knotc_t, knotr, knotr_t, knot_ptm, knot_t, 2 mc, mct, mc_t, mct_t, mrot, mrot_t, mrott, mrott_t, m23, m23_t, 3 n_ch, n_ch_t, n_ptm, n_qs, n_qs_t, n_rot, n_rot_t, 4 old_ptm, q, qt, qt_t, q_t, rota, rota_t, r2, r2_t, x_ptm, xt_ptm IMPLICIT NONE REAL (kind=dp), INTENT(in), DIMENSION(nrot,0:1) :: y REAL (kind=dp), INTENT(in) :: dt, nui REAL (kind=dp), INTENT(out), OPTIONAL, DIMENSION(4,2,0:1) :: dfrl REAL (kind=dp), INTENT(out), DIMENSION(nrl) :: frl REAL (kind=dp), INTENT(out), OPTIONAL, DIMENSION(ncoeff+nchim) :: coeff REAL (kind=dp), DIMENSION(2,0:1) :: ys REAL (kind=dp), DIMENSION(2,2,0:1) :: dhv REAL (kind=dp), DIMENSION(0:NINT(MAXVAL(zi)),nchim) :: ioni REAL (kind=dp), DIMENSION(0,0) :: jac REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: cte_0 REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: hve, hved REAL (kind=dp), DIMENSION(nchim) :: depsx, dxchim, dxchim_t, dxchimg, 1 dxchimg_t, xchim, xchim_t, xchimg, xchimg_t, z_bar REAL (kind=dp), DIMENSION(ne) :: dfqs, fqs REAL (kind=dp), DIMENSION(nrot) :: dfrot, y_t REAL (kind=dp), SAVE, DIMENSION(13) :: cts_0 REAL (kind=dp), DIMENSION(13) :: cts REAL (kind=dp), DIMENSION(5) :: epsilon REAL (kind=dp), DIMENSION(1) :: dfm, fm REAL (kind=dp), DIMENSION(0) :: dcomp REAL (kind=dp), PARAMETER :: dx=1.d-5, unpdx=1.d0+dx REAL (kind=dp), SAVE :: cte1_0, cte2_0, cte3_0, cte4_0, cte5_0, 1 cte6_0, ln13 REAL (kind=dp) :: alfa, Abar, beta, be7e, b8e, chi, chi_mu, chit, chi_t, 1 cp, dcpp, dcpt, dcpx, deff, delta, 2 deltap, deltat, deltax, depsro, depst, dgrad, dgradadp, dgradadt, 3 dgradadx, dh, dkapr, dkapt, dkapx, dlngrav, dlldlm, dlnmu, dlnmu_t, 4 dlnro_nu, dlnro_t, drop, drot, drox, dstor, dup, dut, dux, 5 dv, eps_mu, eps_nuc, eps_T, eta, f17e, gamma1, grad, gradad, 6 grad_mu_shp, grand_k, grav, hhe, hp, ht, kap, lnP, lnP_t, ln_lambda, 7 lum, mw_planet, mr, mu, mu_s, mu_t, mw_dot, nel, nu, nu12, n13e, 8 nu32, nu52, n2mu, n2t, o15e, p, p_phi, p_t, rr_t, r, ro, ro_t, 9 stor, stor0, t, t_t, u, vis_cin, Zbar INTEGER, SAVE :: l=1 INTEGER :: i, id, j LOGICAL, SAVE :: deriv=.FALSE., init=.TRUE. LOGICAL :: zc c------------------------------------------------------------------------ 2000 FORMAT(8es10.3) c définition des constantes IF(init)THEN init=.FALSE. c pour le calcul de grad_mu, chi, grav, d ln grav / d ln R cte1_0=8.d0*pi*rsol**2/3.d0/msol cte2_0=4.d0*aradia*clight/3.d0 cte3_0=g*msol/rsol**2 cte4_0=3.d0*g*msol**2/8.d0/pi/rsol**4 cte5_0=4.d0*rsol**3/msol cte6_0=2.d-15 ln13=LOG(1.3d4) ca.moya ALLOCATE(cte_0(nrl)) cte_0(1)=-rsol/3.d0/pi/g cte_0(2)=-16.d0*rsol**4/27.d0/msol/g cte_0(3)=2.d0*rsol/3.d0 cte_0(4)=-8.d0*rsol**3/9.d0/g/msol cte_0(5)=8.d0*rsol**2/3.d0/g/msol cte_0(6)=-1.d0/pi/g/rsol cte_0(7)=128.d0*pi*rsol**5/27.d0/msol**2/g cte_0(8)=-16.d0*rsol**2/9.d0/msol/g cte_0(9)=2.d0/rsol cte_0(10)=-16.d0*pi*rsol**2/3.d0/msol cte_0(11)=-2.d0/3.d0 cte_0(12)=-8.d0*pi*rsol**2/3.d0/msol cte_0(13:14)=-1.d0 cte_0(15)=-msol/lsol/secon6 cte_0(16)=msol/lsol cte_0(17)=16.d0*pi*rsol**2/9.d0/msol cte_0(18)=-2.d0/3.d0 cte_0(19)=16.d0*pi*rsol**4/9.d0/msol cte_0(20)=-1.d0 cte_0(21)=1.d0 cte_0(22)=-1.d0 cte_0(23)=1.d0 cte_0(24)=-1.d0/secon6 cte_0(25)=8.d0*pi*rsol**3/3.d0/msol cte_0(26)=1.d0/secon6 cte_0(27)=1.d0 cte_0(28)=8.d0*pi*g*rsol**3/3.d0 c cte_0(28)=0.d0 !TEST!coeff de Omega dOmega =0 cte_0(29)=4.d0*pi*g*rsol**3/3.d0 cte_0(30)=1.d0 cte_0(31)=8.d0*pi*rsol**2/15.d0/msol cte_0(32)=64.d0*pi**2*rsol**4/9.d0/msol**2 cte_0(33)=40.d0*pi*rsol**2/3.d0/msol cte_0(34)=1.d0/secon6 cte_0(35)=8.d0*rsol**3/3.d0/g/msol cte_0(36)=-4.d0*rsol**3/3.d0/pi/g**2/msol cte_0(37)=-64.d0*rsol**6/27.d0/g**2/msol**2 cts_0(1)=msol/lsol cts_0(2)=1.d0 cts_0(3)=-msol/lsol/secon6 cts_0(4)=-3.d0*msol/2.d0/pi/rsol**4 cts_0(5)=msol/lsol cts_0(6)=-1.d0 cts_0(7)=-msol/lsol cts_0(8)=6.d0/rsol**2 cts_0(9)=-1.d0/secon6 cts_0(10)=9.d0*msol/4.d0/pi/rsol**3 cts_0(11)=4.d0*pi*g*rsol cts_0(12)=1.d0/secon6 cts_0(13)=-1.d0/rsol**2 ENDIF c charges moyennes dans l'hypothèse totalement ionisé IF(.NOT.mu_saha)z_bar=zi c mises à 0 cts=0.d0 ; frl=0.d0 ; mw_dot=0.d0 IF(PRESENT(dfrl))dfrl=0.d0 c nu12 etc... nu=MAX(nui,nu_min) ; nu12=SQRT(nu) ; nu32=nu12**3 ; nu52=nu*nu32 c est-on dans une ZC? zc=lmix(nu) c~~~~~~~~~~~~~~~ au temps t ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF(en_masse)THEN mr=nu ELSE mr=nu32 ENDIF mr=MAX(x_ptm(1),MIN(mr,x_ptm(n_ptm))) !perte de masse CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,.TRUE., 1 mr,l,fm,dfm) !masse au temps t : fm(1) IF(no_croiss)PRINT*,'Pb en 2 dans coeff_rota4' IF(.NOT.en_masse)fm(1)=fm(1)**(2.d0/3.d0) mr=MAX(mrot_t(1),MIN(fm(1),mrot_t(n_rot_t))) CALL bsp1dn(nrot,rota_t,mrot_t,mrott_t,n_rot_t,m_rot, 1 knotr_t,.TRUE.,mr,l,y_t,dfrot) IF(no_croiss)PRINT*,'Pb en 3 dans coeff_rota4' mr=MAX(m23_t(1),MIN(fm(1),m23_t(n_qs_t))) CALL inter('m23',bp_t,q_t,qt_t,n_qs_t,knot_t,mr,fqs,dfqs,r2_t,m23_t) IF(no_croiss)PRINT*,'Pb. en 4 dans coeff_rota4' lnP_t=fqs(1) ; p_t=EXP(lnP_t) ; t_t=EXP(fqs(2)) IF(en_masse)THEN rr_t=ABS(fqs(3)) ELSE rr_t=fqs(3)**2 ENDIF c composition chimique mr=MAX(mc_t(1),MIN(nu,mc_t(n_ch_t))) CALL bsp1dn(nchim,chim_t,mc_t,mct_t,n_ch_t,m_ch,knotc_t,.TRUE.,mr,l, 1 xchim_t,dxchim_t) IF(no_croiss)PRINT*,'Pb en 4 dans coeff_rota4' c composition chimique par gramme, équation d'état xchimg_t=ABS(xchim_t) ; dxchimg_t=dxchim_t CALL chim_gram(xchimg_t,dxchimg_t) CALL etat(p_t,t_t,xchimg_t,.FALSE.,ro_t,drop,drot,drox,u,dup,dut,dux, 1 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 2 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) c poids moléculaire moyen, appel à saha pour les charges moyennes IF(mu_saha)CALL saha(xchim_t,t_t,ro_t,ioni,z_bar,nel,eta) mu_t=1.d0/DOT_PRODUCT((1.d0+z_bar),xchim_t) c mu_t=16.d0/(23.d0*xchimg_t(1)+3.d0*xchimg_t(ihe4)+9.d0) IF(zc)THEN dlnmu_t=0.d0 ELSE dlnmu_t=-mu_t*DOT_PRODUCT((1.d0+z_bar),dxchim_t) c dlnmu_t=-mu_t/16.d0*(23.d0*dxchimg_t(1)+3.d0*dxchimg_t(ihe4)) ENDIF c ~~~~~~~~~~~~~~~~au temps t+dt ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c ATTENTION, dérivées / m^2/3 ou m effectuées dans inter CALL inter('m23',bp,q,qt,n_qs,knot,nu,fqs,dfqs,r2,m23) IF(no_croiss)PRINT*,'Pb en 0 dans coeff_rota4' lnP=fqs(1) ; p=EXP(fqs(1)) ; t=EXP(fqs(2)) IF(en_masse)THEN r=SQRT(ABS(fqs(3))) ; lum=SQRT(ABS(fqs(4)))**3 ELSE r=ABS(fqs(3)) ; lum=fqs(4) ENDIF dlldlm=dfqs(4)*fqs(5)/fqs(4) !der/m^2/3 dans inter c grad = d ln T / d ln P grad=dfqs(2)/dfqs(1) c hp, ht hp=-rsol*dfqs(3)/dfqs(1) ; IF(en_masse)hp=0.5d0*hp/r ht=hp/grad c composition chimique mr=MAX(mc(1),MIN(nu,mc(n_ch))) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,mr,l,xchim,dxchim) IF(no_croiss)PRINT*,'Pb en 1 dans coeff_rota4' c composition chimique par gramme, équation d'état xchimg=ABS(xchim) ; dxchimg=dxchim CALL chim_gram(xchimg,dxchimg) CALL etat(p,t,xchimg,.FALSE.,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) c gravité, d ln grav / d ln R grav=cte3_0*nu32/r**2 ; dlngrav=cte5_0*r**3*ro/nu32-2.d0 c poids moléculaire moyen c appel à saha pour les charges moyennes IF(mu_saha)CALL saha(xchim,t,ro,ioni,z_bar,nel,eta) mu=1.d0/DOT_PRODUCT((1.d0+z_bar),xchim) c mu=16.d0/(23.d0*xchimg(1)+3.d0*xchimg(ihe4)+9.d0) c--------------- c on dispose de toutes les quantités pour les coefficients dans les ZC IF(zc)THEN c d ln ro / d nu et / dt c phi = d ln ro / d ln mu phi mal connu, on a posé phi=1.d0 dans mod_donnees dlnro_nu=-cte4_0*(alfa-delta*grad)*nu**2/r**4/p dlnro_t=(alfa-delta*grad)*(lnP-lnP_t)+phi*LOG(mu/mu_t) c coefficients de diffusion deff=d_conv ; dh=d_conv ; dv=d_conv c perte de moment cinétique si on est dans la ZC externe IF(nu >= nuzc_ext)THEN IF(lw_perte)CALL pertw(r,mw_dot) IF(lw_planet)THEN CALL planetoides(mw_planet=mw_planet) mw_dot=mw_dot+mw_planet ENDIF ENDIF c les coefficients pour les ZC cts(10)=nu12/r**3/ro*cts_0(10) cts(11)=r*ro/grav*dlnro_nu*cts_0(11) cts(12)=rr_t*y_t(1)/dt*cts_0(12) !rr_t est r_t**2 cts(13)=mw_dot*cts_0(13) frl(25)=r**3*ro/nu12*cte_0(25) frl(26)=r**2*nu12/dt*cte_0(26) frl(27)=(cts(10)+cts(11))*cte_0(27) frl(28)=r**3*ro/grav*cte_0(28) frl(29)=r**3*ro/grav*dlnro_nu*cte_0(29) frl(30)=(cts(12)+cts(13))*nu12*cte_0(30) frl(31)=r**4*ro*cte_0(31) frl(32)=r**6*ro**2*dv/nu12*cte_0(32) frl(33)=r**2*ro*dv/nu12*cte_0(33) frl(34)=nu12/dt*(r**2*(2.d0+dlnro_t)-rr_t)*cte_0(34) c arguments optionnels IF(PRESENT(coeff))THEN coeff(1)=r !R coeff(2)=nu32 !M coeff(3)=y(1,0) !Omega coeff(4)=0.d0 !U coeff(5)=0.d0 !Psi coeff(6)=0.d0 !Lambda coeff(7)=y(5,0) !Phi coeff(8)=frl(32)*y(1,1) !Flux coeff(9)=deff !Deff coeff(10)=dh !Dh coeff(11)=dv !Dv coeff(12)=t !T coeff(13)=ro !ro coeff(14)=0.d0 !grad_mu c appel à nuc pour les epsilon on suppose tout ionisé dans les ZC CALL nuc(t,ro,xchim,dcomp,jac,.TRUE.,3, 1 epsilon,depst,depsro,depsx,hhe,be7e,b8e,n13e,o15e,f17e) eps_nuc=epsilon(1) coeff(15)=cte_0(16)*eps_nuc*nu32/lum/dlldlm !f_epsilon IF(eps_nuc > 0.d0)THEN z_bar=zi mu=1.d0/DOT_PRODUCT((1.d0+z_bar),xchim) eps_T=t/eps_nuc*depst eps_mu=-SUM(depsx(1:nchim)/(1.d0+z_bar(1:nchim)))/mu/eps_nuc ELSE eps_mu=0.d0 ; eps_T=0.d0 ENDIF coeff(16)=eps_T !epsilon_T coeff(17)=eps_nuc !epsilon coeff(18)=dlldlm !dlnL/dlnM coeff(19)=eps_mu !epsilon_mu coeff(20)=0.d0 !chi coeff(21)=0.d0 !chi_t coeff(22)=0.d0 !chi_mu coeff(23)=0.d0 !K coeff(24)=dlnro_nu !d ln ro / dnu coeff(25)=dlnro_t ! dln ro / dt coeff(26)=y(1,1) !d Omega / dnu coeff(27)=y(2,1) !d U / d nu coeff(28)=0.d0 !n2t coeff(29)=0.d0 !n2mu coeff(ncoeff+1:ncoeff+nchim)=xchimg !Xchim WHERE(ABS(coeff) < 1.d-30)coeff=0.d0 ENDIF c----------- c quantités supplémentaires pour les coefficients dans les ZR ELSE c d Ln mu dlnmu=-mu*DOT_PRODUCT((1.d0+z_bar),dxchim) c dlnmu=-mu/16.d0*(23.d0*dxchimg(1)+3.d0*dxchimg(ihe4)) c Grad µ / Hp grad_mu_shp=-cte1_0*r**2*ro/nu12*dlnmu ctest~~~~~~~~~ c grad_mu_shp=0.d0 c d ln ro / d nu et / dt c phi = d ln ro / d ln mu, phi mal connu, on a posé phi=1.d0 dans mod_donnees dlnro_nu=-cte4_0*(alfa-delta*grad)*nu**2/r**4/p+phi*dlnmu dlnro_t=(alfa-delta*grad)*(lnP-lnP_t)+phi*LOG(mu/mu_t) c estimation de vis_cin, la viscosité cinématique, Lang 3.15, 3.25 et 3.28 Abar=DOT_PRODUCT(xchim,nucleo)/SUM(xchim) Zbar=DOT_PRODUCT(z_bar,xchim)/SUM(xchim) IF(.NOT.mu_saha)nel=DOT_PRODUCT(z_bar,xchim)*ro/amu ln_lambda=ln13+1.5d0*LOG(t)-LOG(nel)/2.d0 vis_cin=cte6_0*SQRT(t)**5*SQRT(Abar)/ro/Zbar**4/ln_lambda ca.moya, p_phi = fonction thermodynamique (6.8) de Maeder & Zahn A&A 98 c n'est valable que si Z=cte, i.e. sans diffusion et sur la MS c pour cohérence il faut utiliser mu simplifié et Z initial IF(diffusion .OR. xchim(1) < x_tams)THEN p_phi=0.d0 ELSE mu_s=16.d0/(20.d0*xchimg(1)+12.d0-3.d0*z0) p_phi=-mu_s*(12.d0*(z0-4.d0)*LOG(4.d0*(16.d0+3.d0*mu_s*(z0-4.d0))) 1 +(32.d0-23.d0*z0)*LOG(-16.d0+mu_s*(32.d0-23.d0*z0)) 2 +5.d0*z0*LOG(5.d0*z0) 3 +2.d0*(8.d0+3.d0*z0)*LOG(2.d0*(16.d0+mu_s*(8.d0+3.d0*z0))))/(80.d0*cp) ENDIF c opacité CALL opa(xchimg,t,ro,kap,dkapt,dkapr,dkapx) ca.moya conductivité radiative chi, dérivées et diffusivité thermique K chi=cte2_0*t**3/ro/kap chit=3.d0-t/kap*dkapt+delta chi_mu=-nucleo(1)/mu/kap/(1.d0+z_bar(1))*dkapx-phi grand_k=chi/ro/cp c chi_t=chit /=0.d0 fait diverger chi_t=0.d0 c vaïssälä dgrad=ABS(gradad-grad) n2t=grav*delta/hp*dgrad ; n2mu=grav*phi*grad_mu_shp c appel à nuc pour les epsilon CALL nuc(t,ro,xchim,dcomp,jac,.TRUE.,3, 1 epsilon,depst,depsro,depsx,hhe,be7e,b8e,n13e,o15e,f17e) eps_nuc=epsilon(1) IF(eps_nuc > 0.d0)THEN eps_mu=-SUM(depsx(1:nchim)/(1.d0+z_bar(1:nchim)))/mu/eps_nuc c eps_mu=-16.d0/mu/eps_nuc*(depsx(1)/23.d0/nucleo(1)+ c 1 depsx(ihe4)/3.d0/nucleo(ihe4)) eps_T=t/eps_nuc*depst ELSE eps_mu=0.d0 ; eps_T=0.d0 ENDIF c coefficients de diffusion du moment cinétique y --> ys c on utilise Omega et U et leurs dérivées au pas temporel précédent c dérivées dhv sont trop instables IF(PRESENT(coeff))THEN ys(1:2,:)=y(1:2,:) ELSE ys(1:2,0)=y_t(1:2) ; ys(1:2,1)=dfrot(1:2) c ys(1:2,:)=y(1:2,:) !totalement implicite, ne marche pas ENDIF CALL diffw(dlnro_nu,grand_k,nu,vis_cin,n2mu,n2t,r,ro,ys,deff,dh,dv,dhv) dhv=0.d0 !dérivées nulles c ~~~~~~~~~~~~~~~~~les coefficients~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c les dérivées c dhv(1,1:2,0:1)= d dh (dv) / d Omega : / (d U : d dU) c dfrl(1,1:2,0:1) : d frl(13) / (d Omega : d dOmega) : / (d U : d dU) c dfrl(2,1:2,0:1) : d frl(23) / (d Omega : d dOmega) : / (d U : d dU) c dfrl(3,1:2,0:1) : d frl(32) / (d Omega : d dOmega) : / (d U : d dU) c dfrl(4,1:2,0:1) : d frl(33) / (d Omega : d dOmega) : / (d U : d dU) ca.moya les F* cts(1)=eps_nuc*nu12*(eps_T-delta)/lum cts(2)=dlldlm/nu*(2.d0*delta-1.d0-chi_t) cts(3)=nu12*cp*T/lum/dt c dfrl(1,1:2,0:1) : d frl(13) / (d Omega : d dOmega) : / (d U : d dU) cts(4)=ht*nu12/r**4/ro*(1.d0+Dh/grand_k) IF(PRESENT(dfrl))dfrl(1,:,:)=ht*nu12/r**4/ro/grand_k*dhv(1,:,:) 1 *cts_0(4)*cte_0(13) cts(5)=eps_nuc*nu12*(eps_mu+phi)/lum cts(6)=dlldlm/nu*(2.d0*phi+chi_mu) cts(7)=nu12*cp*T*p_phi*(dlnmu-dlnmu_t)/lum/dt c dfrl(2,1,0:1) : d frl(23) / (d Omega : d dOmega) : / (d U : d dU) cts(8)=1.d0/r**2 IF(PRESENT(dfrl))dfrl(2,:,:)=dhv(1,:,:)*cts(8)*cts_0(8)*cte_0(23) cts(8)=cts(8)*dh*cte_0(23) cts(9)=(dlnmu-dlnmu_t-1.d0)/dt cts(10)=nu12/r**3/ro cts(11)=r*ro/grav*dlnro_nu cts(12)=rr_t*y_t(1)/dt !rr_t est r_t**2 cts(13)=mw_dot cts=cts*cts_0 c les F frl(1)=r/ro/nu/grav*(dlngrav-2.d0) frl(2)=r**4/nu32/grav*(dlngrav-2.d0) frl(3)=r/nu/grav*(1.d0-dlldlm)*(dlngrav-2.d0) frl(4)=r**3/nu32 frl(5)=r**2/grav/nu32 frl(6)=1.d0/grav/ro/nu/r*dlngrav frl(7)=r**5*ro/grav/nu**2 frl(8)=r**2/grav/nu32*dlngrav frl(9)=(1.d0-dlldlm)/grav/nu/r*dlngrav frl(10)=r**2*ro/grav/nu32*(1.d0-dlldlm) frl(11)=phi+chi_mu frl(12)=r**2*ro*ht/nu32*dlldlm frl(13)=cts(1)+cts(2)+cts(3)+cts(4) frl(14)=cts(5)+cts(6)+cts(7) frl(15)=nu12*cp*T*y_t(3)/lum/dt frl(16)=grav*ro*cp*T*nu12/lum/P*dgrad frl(17)=r**2*ro*ht/nu12 frl(18)=1.d0-delta+chi_t !chi_t non nul ici fait diverger frl(19)=r**4*ro/grav/nu12 frl(20)=phi frl(21)=delta frl(22)=grad_mu_shp frl(23)=cts(8)+cts(9) frl(24)=y_t(4)/dt frl(25)=r**3*ro/nu12 frl(26)=r**2*nu12/dt frl(27)=cts(10)+cts(11) frl(28)=r**3*ro/grav frl(29)=r**3*ro/grav*dlnro_nu frl(30)=(cts(12)+cts(13))*nu12 frl(31)=r**4*ro c dfrl(3,1:2,0:1) : d frl(32) / (d Omega : d dOmega) : / (d U : d dU) frl(32)=r**6*ro**2/nu12 IF(PRESENT(dfrl))dfrl(3,:,0:1)=frl(32)*dhv(2,:,:)*cte_0(32) frl(32)=frl(32)*dv c dfrl(4,1:2,0:1) : d frl(33) / (d Omega : d dOmega) : / (d U : d dU) frl(33)=r**2*ro/nu12 IF(PRESENT(dfrl))dfrl(4,:,:)=frl(33)*dhv(2,:,:)*cte_0(33) frl(33)=frl(33)*dv frl(34)=nu12/dt*(r**2*(2.d0+dlnro_t)-rr_t) frl(35)=r**3/nu52*(1.d0-dlldlm) frl(36)=r**3/ro/nu52 frl(37)=r**6/nu**3 frl=frl*cte_0 c arguments optionnels IF(PRESENT(coeff))THEN coeff(1)=r !R coeff(2)=nu32 !M coeff(3)=y(1,0) !Omega coeff(4)=y(2,0) !U coeff(5)=y(3,0) !Psi coeff(6)=y(4,0) !Lambda coeff(7)=y(5,0) !Phi coeff(8)=frl(31)*y(1,0)*y(2,0)+frl(32)*y(1,1) !flux coeff(9)=deff !Deff coeff(10)=dh !Dh coeff(11)=dv !Dv coeff(12)=t !T coeff(13)=ro !ro coeff(14)=grad_mu_shp*hp !grad_mu coeff(15)=cte_0(16)*eps_nuc*nu32/lum/dlldlm !f_epsilon coeff(16)=eps_T !epsilon_T coeff(17)=eps_nuc !epsilon coeff(18)=dlldlm !dlnL/dlnM coeff(19)=eps_mu !epsilon_mu coeff(20)=chi !chi coeff(21)=chit !chi_t coeff(22)=chi_mu !chi_mu coeff(23)=grand_k !K coeff(24)=dlnro_nu !d ln ro / dnu coeff(25)=dlnro_t ! dln ro / dt coeff(26)=y(1,1) !d Omega / dnu coeff(27)=y(2,1) !d U / d nu coeff(28)=n2t !n2t coeff(29)=n2mu !n2mu coeff(ncoeff+1:ncoeff+nchim)=xchimg !Xchim WHERE(ABS(coeff) < 1.d-30)coeff=0.d0 ENDIF ENDIF !zc c tests de dérivation pour les frl(13,23,32,33) c deriv=nu >= 6.E-01 .AND. nu <= 6.3E-01 c deriv=nu >= 6.2E-01 .AND. nu <= 6.25E-01 c deriv=dh > 1.d2 IF(deriv)THEN PRINT*,'nu' ; WRITE(*,2000)nu PRINT*,'y' ; WRITE(*,2000)y(1:2,0) PRINT*,'dy' ; WRITE(*,2000)y(1:2,1) ALLOCATE(hve(4),hved(4)) hve(1)=frl(13) ; hve(2)=frl(23) ; hve(3)=frl(32) ; hve(4)=frl(33) DO id=0,1 B1: DO i=1,2 ys(1:2,:)=y(1:2,:) stor0=ys(i,id) ; stor=stor0*unpdx ys(i,id)=stor ; dstor=stor-stor0 CALL diffw(dlnro_nu,grand_k,nu,vis_cin,n2mu,n2t,r,ro,ys,deff, 1 dh,dv,dhv) c frl(13) cts(4)=ht*nu12/r**4/ro*(1.d0+dh/grand_k)*cts_0(4) hved(1)=(cts(1)+cts(2)+cts(3)+cts(4))*cte_0(13) c frl(23) cts(8)=1.d0/r**2*dh*cts_0(8)*cte_0(23) hved(2)=cts(8)+cts(9) c frl(32) hved(3)=r**6*ro**2/nu12*dv*cte_0(32) c frl(33) hved(4)=r**2*ro/nu12*dv*cte_0(33) IF(id == 0)THEN PRINT*,'frl(13,23,32,33), dérivées ',nom_rot(i) ELSE PRINT*,'frl(13,23,32,33), dérivées d',nom_rot(i) ENDIF WRITE(*,2000)(dfrl(j,i,id),(hved(j)-hve(j))/dstor,j=1,4) PAUSE ENDDO B1 ENDDO DEALLOCATE(hve,hved) ; deriv=.FALSE. ENDIF !deriv RETURN END SUBROUTINE coeff_rota4