c*********************************************************** REAL (kind=dp) FUNCTION t_convec(r_lim) c routine private du module mod_static c estimation du temps caractérisque de convection au rayon R*= R0 + alpha/2 Hp(R*), c R0 rayon à la base de la ZC, si R* est extérieur à la ZC, R* est pris au milieu c en cas de difficulté R*=r_demi, par convention t_convec < 0 c entrées c r_lim(1:2): rayons inférieur et supérieur de la ZC en r/Rsol c Auteur: P.Morel, laboratoire Lagrange, O.C.A., CESAM2k c----------------------------------------------------------------------- USE mod_donnees, ONLY : alpha, aradia, clight, en_m23, g, 1 Krot, lsol, m_ch, msol, nchim, ne, nrot, nucleo, ord_rot, pi, rsol USE mod_conv, ONLY : conv USE mod_etat, ONLY : etat USE mod_kind USE mod_numerique, ONLY : bsp1dn, entre, no_croiss USE mod_opa, ONLY : opa, dehors USE mod_variables, ONLY: bp, chim, knot, knotc, knotr, inter, 1 q, qt, mc, mct, mrot, mrott, m_stat, n_ch, n_qs, 2 n_rot, omega_jpz, rota, r_stat, sortie, wrot IMPLICIT NONE REAL(kind=dp), INTENT(in), DIMENSION(2) :: r_lim REAL(kind=dp), PARAMETER :: epsi=1.d-3 REAL(kind=dp), DIMENSION(ne) :: dfdq, f REAL(kind=dp), DIMENSION(nchim) :: dxchim, xchim REAL (kind=dp), DIMENSION (nrot) :: dfrot, frot REAL(kind=dp), SAVE :: cte1, cte2, cte3, cte8, cte13 REAL (kind=dp) :: alfa, beta, cp, delta, deltap, deltat, deltax, 1 dcpp, dcpt, dcpx, dgamcp, dgamdel, dgamgad, dgamgra, 2 dgamgrad, dgamhp, dgamkra, dgamro, dgamtaur, dgradcp, 3 dgradadp, dgradadt, dgradadx, dgradel, dgradgad, dgradgra, 4 dgradgrad, dgradhp, dgradkra, dgradro, dgradtaur, dkapt, dkapro, 5 dkapx, dlim, dr, drop, drot, drox, dup, dut, dux, dwdr, 6 gam, gamma1, grad, gradad, gradrad, gravite, hp, kap, krad, l, 7 m_moyen, nu, p, pas, rmoins, ro, rot, rplus, r0, r_demi, r_eff, 8 r_moyen, r_new, t, taur, u, w INTEGER, PARAMETER :: np=10, ntour_max=20 INTEGER, SAVE :: lc=1 INTEGER :: ntour LOGICAL, SAVE :: init=.TRUE. LOGICAL :: dico c------------------------------------------------------------------- 2000 FORMAT(8es10.3) c initialisations IF(init)THEN init=.FALSE. cte1=4.d0/3.d0*aradia*clight cte2=2.d0/3.d0*rsol cte3=(alpha*rsol**2/g/msol)**2/6.d0/aradia/clight cte8=lsol/4.d0/pi/rsol/rsol cte13=g*msol/rsol/rsol ENDIF c initialisation de variables dlim=r_lim(2)-r_lim(1) ; r0=r_lim(1) ; r_demi=SUM(r_lim)/2.d0 ntour=0 ; dico=.FALSE. ; pas=dlim/REAL(np-1,dp) r_moyen=r0+pas ; rplus=r0 ; rmoins=-100.d0 c localisation de R*=r_moyen en se déplaçant de R0+pas à r_lim(2) B1: DO c calcul de R_new = R0 + alpha/2 Hp(r_moyen) r_new=rconv_mean(r0,r_moyen) dr=r_new-r_moyen IF(ABS(dr) <= epsi)THEN r_moyen=(r_new+r_moyen)/2.d0 ; EXIT B1 ELSEIF(dr > 0.d0)then rplus=r_moyen ELSE rmoins=r_moyen ; dico=.TRUE. ENDIF c PRINT*,'dr, r_new, rmoyen, r_lim, ntour,dico',ntour,dico c WRITE(*,2000)dr, r_new, r_moyen, r_lim c dichotomie IF(dico)THEN ntour=ntour+1 IF(ntour > ntour_max)THEN r_moyen=r_demi WRITE(*,1)r_moyen 1 FORMAT('t_convec: la dicho. ne converge pas R*=r_milieu=', 1 es10.3) EXIT B1 ENDIF r_moyen=(rplus+rmoins)/2.d0 CYCLE B1 c poursuite de la localisation de R* ELSE r_moyen=r_moyen+pas IF(r_moyen >= r_lim(2))THEN r_moyen=r_demi WRITE(*,2)r_moyen 2 FORMAT('t_convec < 0 : R* non localisé utilisation de r_milieu=', 1 es10.3) EXIT B1 ELSE CYCLE B1 ENDIF ENDIF ENDDO B1 c PAUSE'B1' c r_moyen doit être dans [r_lim(1) , r_lim(2)] IF(.NOT.entre(r_lim(1),r_lim(2),r_moyen))THEN WRITE(*,3)r_lim(1),r_moyen,r_lim(2) 3 FORMAT('ERREUR, sortie de t_convec r_moyen hors limites',3es10.3) CALL sortie('t_convec 1') ENDIF c calcul de t_conv c r_eff: rayon moyen en m^2/3 ou m^1/3 IF(en_m23)THEN r_eff=r_moyen**2 ELSE r_eff=r_moyen ENDIF c variables en r_moyen (r_eff), m_moyen en m/Msol, nu en m^2/3 CALL inter('r',bp,q,qt,n_qs,knot,r_eff,f,dfdq,r_stat,m_stat) p=EXP(f(1)) ; t=EXP(f(2)) IF(en_m23)THEN nu=ABS(f(5)) m_moyen=SQRT(nu)**3 l=SQRT(f(4))**3 ELSE nu=f(5)**2 m_moyen=ABS(f(5))**3 l=f(4) ENDIF c composition chimique en r_moyen (nu=m^2/3) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(nu,mc(n_ch))),lc,xchim,dxchim) c pression gazeuse (Pgas + Prad) xchim=xchim*nucleo !xchim par gramme CALL etat(p,t,xchim,.FALSE., 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 vérification c beta=alpha*rsol/2.d0/g/msol*p*r_moyen-m_moyen*ro !beta:vt c IF(r0 > 0.d0)THEN c PRINT*,'cas r0 > 0,r_moyen,r_demi,r_lim' c WRITE(*,2000)beta*r_moyen+m_moyen*ro*r0,r_moyen,r_demi,r_lim c ELSE c PRINT*,'cas r0 = 0,r_moyen,r_demi,r_lim' c WRITE(*,2000)beta,r_moyen,r_demi,r_lim c ENDIF c PAUSE'verif' c opacité CALL opa(xchim,t,ro,kap,dkapt,dkapro,dkapx) IF(dehors)THEN PRINT*,'t_convec, appel de thermo, arrêt / T, ro, m / xchimg' WRITE(*,2000)t,ro,m_moyen ; WRITE(*,2000)xchim CALL sortie('t_convec 2') ENDIF c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w=wrot CASE(3,4,5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MIN(nu,mrot(n_rot)),lc,frot,dfrot) IF(no_croiss)PRINT*,'Pb. en 1 dans t_convec' w=frot(1) CASE(6) CALL omega_jpz(r_moyen,w,dwdr) END SELECT c gravité effective gravite=cte13*m_moyen/r_moyen**2 rot=cte2*w**2*r_moyen !gravite effective avec rotation gravite=gravite-rot !remarque de N. Audard c échelle de hauteur de pression hp=p/gravite/ro c diffusivite radiative krad=cte1/kap/ro*t**3 c gradient radiatif gradrad=cte8*l*hp/r_moyen**2/krad/t c gradient radiatif gradrad=cte8*l*hp/r_moyen**2/krad/t c épaisseur optique de la bulle taur=kap*ro*alpha*hp c gam: efficacité de la convection CALL conv(r_moyen,krad,gravite,delta,cp,ro,hp,taur,gradrad,gradad, 1 .FALSE.,grad,dgradkra,dgradgra,dgradel,dgradcp,dgradro, 2 dgradhp,dgradtaur,dgradgrad,dgradgad, 3 gam,dgamkra,dgamgra,dgamdel,dgamcp,dgamro, 4 dgamhp,dgamtaur,dgamgrad,dgamgad) c temps caractéristique de convection en années IF(gam > 0.d0)THEN t_convec=cte3*(p*r_moyen**2/m_moyen)**2*kap*cp/t**3/gam ELSE t_convec=0.d0 ENDIF c PRINT*,'sortie t_convec: r_lim,r_moyen,gam,t_convec' c WRITE(*,2000) r_lim,r_moyen,gam,t_convec ; PAUSE't_convec' c en cas de difficulté R*=r_demi, par convention t_convec < 0 IF(r_moyen == r_demi)t_convec=-t_convec RETURN END FUNCTION t_convec