c******************************************************************** SUBROUTINE static_m13(fait,cx,li,y,be,ae,compt,dt,reprend,ip) c routine private du module mod_static c Calcul des quantités caractéristiques du problème quasi-statique en c lagrangien c variables indépendantes ln P, ln T, r, l, m^1/3 c routine appellée par resout c spline-collocation avec fonction de répartition c avec pression turbulente 7 inconnues c sans pression turbulente 6 inconnues, Ptot=Pgaz c pour avoir un algorithme unique pour la diffusion, la composition c chimique est toujours tabulée en fonction de mu=m^2/3 c l'énergie graviphique, TdS/dt=tds est tabulée en fonction de m^1/3 c Auteur: P.Morel, Département J.D. Cassini, O.C.A., CESAM2k c entrées c fait=1 c -calcule les résidus be(i) des équations au point xcoll i=1,ne c fait=2 c -calcule le "résidu" be(1) de la condition au point limite li c cx : indice du point de collocation c li : numéro de la limite c y : variables au point de collocation xcoll(cx) c compt : compteur du nb. iter. Newton Raphson c dt : pas temporel c ip : indice de couche pour détermination du facteur de répartition c sorties c be : résidus c ae : termes du jacobien pour le point de collocation c reprend=.TRUE. : variation relative de U ou Ro, lnP, lnT trop forte c ou R, L, M < 0, ou acc. centrifuge trop grande c---------------------------------------------------------------- USE mod_donnees, ONLY : alpha, aradia, clight, cpturb, dtmin, d_press, 1 d_temp, fcv, g, ihe4, ini0, Ipg, Krot, kipp, langue, lsol, l_pertm, 2 msol, m_ch, m_ptm, m_tds, nchim, ne, nrot, nucleo, ord_rot, pi, precix, 3 pturb, rsol, secon6, t_inf, m_qs USE mod_atm, ONLY : atm USE mod_etat, ONLY : etat USE mod_evol, ONLY : lmix, lmix_t USE mod_kind USE mod_numerique, ONLY : bsp1dn, entre, no_croiss USE mod_variables, ONLY : bp_t, chim, chim_gram, chim_t, ctem, ctep, 1 cter, ctet, inter, knotc, knotc_t, knotr, knotr_t, knot_tds_t, knot_ptm, 2 knot_t, mc, mct, mct_t, mc_t, mrot, mrot_t, mrott, mrott_t, mstar, m_stat_t, 3 n_ch, n_ch_t, n_ptm, n_qs, n_qs_t, n_rot, n_rot_t, n_tds_t, 4 old_ptm, omega_jpz, q, qt_t, q_t, rota, rota_t, rstar, r_stat_t, 5 sortie, tds_t, wrot, xt_ptm, xt_tds_t, x_ptm, x_tds_t IMPLICIT NONE REAL (kind=dp), INTENT(in), DIMENSION(:,0:) :: y REAL (kind=dp), INTENT(in) :: dt INTEGER, INTENT(in) :: cx, compt, fait, li, ip REAL (kind=dp), INTENT(out), DIMENSION(:,:,0:) :: ae REAL (kind=dp), INTENT(out), DIMENSION(:) :: be LOGICAL, INTENT(out) :: reprend REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: yd REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: xchim0 REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: bes, dern REAL (kind=dp), DIMENSION(nchim) :: depsx, dmu_x, dxchim, 1 dxchim_t, dxchim0, xchim, xchim_t REAL (kind=dp), DIMENSION(ne) :: dfqs, fqs, y_t REAL (kind=dp), DIMENSION(nrot) :: dfrot, frot REAL (kind=dp), DIMENSION(5) :: epsilon REAL (kind=dp), DIMENSION(1) :: dfm, dftds, fm, ftds REAL (kind=dp), SAVE :: ay52, ctev, cte1, cte2, cte3, cte4, 1 cte6, cte11, cte13, cte14, cte14t, cte15, cte21, cte22, 2 cte23, cte25, dd, dmdl0, dmdr0, dpdl0, dpdr0, dptdl0, dptdr0, 3 dtdl0, dtdr0, mext0, pext0, ptext0, text0, unpdd REAL (kind=dp) :: acc_c, alfa, ay3, ay32, ay34, ay4, ay5, ay53, 1 ay55, beta, cp, cp_t, cp0, dcpm, dcpp, dcpp_t, dcpt, dcpt_t, dcpx, 2 dcpx_t, delta, delta_t, deltam, deltap, deltap_t, deltat, deltat_t, 3 deltax, deltax_t, depsm, depsp, depst, dgamlpt, dgamp, 4 dgampt, dgamt, dgamx, dgamm, dgaml, dgamlp, dgamlpp, dgamr, 5 dgam1, dgradadm, dgradadp, dgradadp_t, dgradadt, dgradadt_t, 6 dgradadx, dgradadx_t, gradad0, dgradlp, dgradlpt, dgradp, dgradpt, 7 dgradt, dgradx, dgradm, dgradl, dgradr, dgradlpp, dhpp, dhppt, dhpt, 8 dhpx, dhpr, dhpm, dkapp, dkapt, dkapx, dlnp=0.d0, dlnp0=0.d0, 9 dlnp_tm=0.d0, dlnt=0.d0, dlnt_tm=0.d0, dlpp=0.d0, dlppdp=0.d0, 1 dlppdpt=0.d0, dlnt0=0.d0, dmdl, dmdr, dpn, dpdl, dpdr, dpsisti, 2 dpsitl=0.d0, dpsitm, dpsitp, dpsitr, dpsitt, dptdl, dptdr, 3 dp_tm=0.d0, drom, drop, drop_t, drot, drot_t=0.d0, drox, 4 drox_t=0.d0, dro_tm=0.d0, dstor, dtdl, dtdr, dtdsm=0.d0, 5 dtdsp=0.d0, dtdst=0.d0, dtetal, dtetam, dtetap, 6 dtetar, dtetat, dtn, dt_tm=0.d0, dum, dup, dup_t=0.d0, dut, 7 dut_t=0.d0, dux, dux_t=0.d0, du_tm=0.d0, dvdm, dvdr, 8 dwdm, dwdr, dxy, dyv, dyvdm, dyvdr, dy0, dy1, dy2, dy3, dy4, dy5, 9 dy71=0.d0, dy710, dy71l=0.d0, dy71lp=0.d0, dy71lpt=0.d0, 1 dy71m=0.d0, dy71p=0.d0, dy71pt=0.d0, dy71r=0.d0, dy71t=0.d0, 2 dy72=0.d0, dy720=0.d0, dy72l=0.d0, dy72lp=0.d0, dy72lpt=0.d0, 3 dy72p=0.d0, dy72pt=0.d0, dy72m=0.d0, dy72r=0.d0, dy72t=0.d0, 4 d5tds, d5tdsm, d5tdsp, d5tdst, eps0, eps_rot, gam, gam0, 5 gam1, gamma1, grad, gradad, gradad_t, gradrad, gradrad0, grad0, 6 grad_mu, gravs2, hp, hp0, kap, mext, mk, mu, mue, pext, 7 pgn, prn, psist, psist0, ptext, p_t=0.d0, p_t0, ro, ro0, 8 ro_t=0.d0, ro_t0, rsrst, r_t, sig, stor, stor0, tdst=0.d0, tdst0, 9 teff, teta, teta0, text, trn, t_t=0.d0, 1 t_t0, u, u_t=0.d0, u_t0, u0, x00, v, w, w_t INTEGER, SAVE :: l INTEGER :: i, id, j, jd LOGICAL, SAVE :: der, deriv=.FALSE., init=.TRUE., l_he4 LOGICAL :: radiatif, radia0 CHARACTER (len=7), SAVE, DIMENSION(7) :: variable=(/ 'ln Ptot', 1 'ln T ','R ','L ','M^1/3 ','Psi ','ln Pgaz'/) c------------------------------------------------------------------------ 2000 FORMAT(8es10.3) 2001 FORMAT(5es15.8) c PRINT* ; PRINT*,'entree static_m13',fait,compt,cx,ip,kipp c WRITE(*,2000)dt,xcoll(cx) ; WRITE(*,2000)y(1:ne,0),fac(ip) c WRITE(*,2000)y(1:ne,1) IF(init)THEN init=.FALSE. ctev=aradia*clight/2.d0 cte1=g*msol/rsol**2 cte2=2.d0/3.d0*rsol cte3=2.d0/nucleo(1) cte4=rsol**2/3.d0/secon6 cte11=-3.d0*g/4.d0/pi*(msol/rsol**2)**2 !pour d ln P /d q cte21=cte11*ctep ; cte22=cte11*ctet cte13=msol/rsol*3.d0/4.d0/pi/rsol/rsol !pour d zeta /d q cte23=cte13*cter cte14=3.d0*msol/lsol !pour d lambda/d q cte14t=cte14/secon6 ! " " dt /= 0 cte15=msol/rsol/2.d0/pi !pour la rotation cte25=cte15*ctep cte6=ABS(cpturb)*alpha**2/8.d0 !pour la Pturb der=cpturb < 0.d0 !der=.TRUE. avec dlnPgaz/dlnPtot ALLOCATE(xchim0(nchim)) !pour l'atmosphère c type de calcul du TdS IF(kipp)THEN WRITE(2,10) ; WRITE(*,10) 10 FORMAT('Energie graviphique approx. de Kippenhahn: cp T (dlnT-gradad dlnP)') ELSE WRITE(*,11) ; WRITE(2,11) 11 FORMAT(/,'Energie graviphique forme basique dU + PdV') ENDIF c indice de He4, ihe4=-100 avec pp1, reac_nuc sans hélium explicite, 3 ou 4 sinon l_he4=ihe4 > 0 ENDIF !init c PRINT*,'fait,cx',fait,cx c------------------------------------------------------------------ c ensemble des variables c y(1,0)=ln Ptot y(1,1)=(ln Ptot)' c y(2,0)=ln T y(2,1)=(ln T)' c y(3,0)=r y(3,1)=r' c y(4,0)=l y(4,1)=l' c y(5,0)=m^1/3 y(5,1)=(m^1/3)' c y(6,0)=psi=dQ/dq(=cte) y(6,1)=psi'(=0) c avec Pturb y(Ipg,0)=ln Pgaz y(Ipg,1)=(ln Pgaz)' c xcoll(cx) valeur de q au point de collocation q=1, 2, 3, .. , n_qs c ae(eq,var,der)=ae(ne,ne,0:1)=dérivée de la eq-ieme équation par c rapport à la der-ieme dérivée de la var-ieme variable ae=0.d0 ; be=0.d0 SELECT CASE(fait) c aux points de collocation CASE(1) prn=EXP(y(1,0)) !pression Ptot cgs IF(pturb)THEN !avec pression turbulente pgn=EXP(y(Ipg,0)) !pression Pgaz cgs ELSE !sans pression turbulente pgn=prn ENDIF trn=EXP(y(2,0)) !température K ay3=ABS(y(3,0)) !r=R/Rsol rsrst=ay3*rstar !R/Rstar ay32=ay3**2 !r^2 ay34=ay32**2 !r^4 ay4=y(4,0) !l=L/Lsol ay5=ABS(y(5,0)) !m^1/3 pour éviter des valeurs < 0 ay52=ay5*ay5 !m^2/3 ay53=ay52*ay5 !M/Msol ay55=ay52*ay53 !m^5/3 c l'acc. centrifuge ne doit pas excèder 90% gravité gravs2=cte1*ay53/ay32*0.9d0 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 MAX(mrot(1),MIN(ay52,mrot(n_rot))),l,frot,dfrot) IF(no_croiss)PRINT*,'Pb. en 0 dans static_m13' w=frot(1) CASE(6) CALL omega_jpz(ay3,w) END SELECT c accélération centrifuge acc_c=cte2*ay3*w**2 ; reprend=acc_c > gravs2 IF(reprend)THEN SELECT CASE(langue) CASE('english') WRITE(*,1012)acc_c,gravs2,w,rsrst,ay53 WRITE(2,1012)acc_c,gravs2,w,rsrst,ay53 1012 FORMAT('STOP, in static_m13, the centrifugal acceleration =', 1 es10.3,' > 90% gravity = ',es10.3,/,'angular velicity =',es10.3, 2 ', R/Rstar=',es10.3,', M/Msol=',es10.3) CASE DEFAULT WRITE(*,12)acc_c,gravs2,w,rsrst,ay53 WRITE(2,12)acc_c,gravs2,w,rsrst,ay53 12 FORMAT('ARRET, dans static_m13, acc. centifuge =', 1 es10.3,' > 90% gravité = ',es10.3,/,'vitesse angulaire =',es10.3, 2 ', R/Rstar=',es10.3,', M/Msol=',es10.3) END SELECT CALL sortie('static_m13') ENDIF c si, à cause d'une erreur d'arrondi, ce qui peut arriver au voisinage c du centre r ou m**1/3 soit negatif, on tente de forcer la c convergence en utilisant les ABS(**), les dérivées sont inexactes et c la physique est violée... mais si ça passe les erreurs d'arrondi c disparaissent, parfois, au cours des iterations c pour éviter cette disposition imposer fcv=.FALSE. dans cesam.f ou dans le c fichier de réglage reprend=MIN(y(3,0),y(4,0),y(5,0)) < 0.d0 reprend=.FALSE. IF(reprend)THEN c on tente de converger reprend=.FALSE. IF(ABS(y(4,0)) > precix)THEN SELECT CASE(langue) CASE('english') WRITE(*,1014)y(3,0),y(4,0),y(5,0),ay53,cx 1014 FORMAT(/,'static_m13, r=',es10.3,', l=',es10.3', m^1/3=',es10.3, 1 ', M=',es10.3,/,'n° coll=',i4,', one tries to converge') CASE DEFAULT WRITE(*,14)y(3,0),y(4,0),y(5,0),ay53,cx 14 FORMAT(/,'static_m13, r=',es10.3,', l=',es10.3', m^1/3=',es10.3, 1 ', M=',es10.3,/,'n° coll=',i4,', on tente de converger') END SELECT ENDIF ENDIF c composition chimique au temps t+dt, dxchim: d xchim/dm^2/3 pour thermo CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MIN(ay52,mc(n_ch)),l,xchim,dxchim) IF(no_croiss)PRINT*,'Pb. at 1 in static_m13' WHERE(xchim < 1.d-38) xchim=1.d-38 dxchim=1.d-38 ENDWHERE c pression turbulente IF(pturb .AND. der)THEN !avec pression turbulente 7 inconnues dlpp=y(Ipg,1)/y(1,1) !dlpp=dln Pgaz/dln Ptot dlppdpt=-dlpp/y(1,1) !dérivée dlpp /dln Ptot dlppdp=dlpp/y(Ipg,1) !dérivée dlpp /dln Pgaz ELSE !der=.FALSE. on ne tient pas compte de dln Pgaz/dln Ptot dlpp=1.d0 dlppdpt=0.d0 ; dlppdp=0.d0 ENDIF c pour thermo dxchim est /m^23, xchim et dxchim par mole CALL thermo(prn,pgn,trn,ay53,ay4,ay3,dlpp,xchim,dxchim, 1 ro,drop,drot,drox,u,dup,dut,dux, 2 grad,dgradpt,dgradp,dgradt,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,grad_mu,mu,dmu_x,mue, 8 gradrad,alfa,beta,gamma1,radiatif) c WRITE(*,2000)ay53,prn,pgn,trn,ro,xchim,epsilon(1),depsp,depst,depsx c PRINT*,'couche',cx ; WRITE(*,2000)prn,pgn,trn,ay53,ln,ay3,ro,grad c WRITE(*,2000)gam,epsilon(1),delta,gradad,hp,gradrad c WRITE(*,2000)xchim ; WRITE(*,2000)epsilon c WRITE(*,2000)prn,pgn,dgradpt,dgradp ; PAUSE'thermo' c dXchim/dm^1/3 = 2 m^1/3 dXchim/dm^2/3 dxchim=dxchim*2.d0*ay5 c dérivées des abondances de H et He4 dxy=dxchim(1) IF(l_he4)dxy=dxy-dxchim(ihe4) c dérivées par rapport à : ln Ptot, ln T, zeta, lambda, m^1/3, ln Pgaz drop=drop*pgn !der ro /ln Pgaz drot=drot*trn ; drom=drox*dxy !d ro/m^1/3 dcpp=dcpp*pgn ; dcpt=dcpt*trn ; dcpm=dcpx*dxy deltap=deltap*pgn ; deltat=deltat*trn ; deltam=deltax*dxy dup=dup*pgn ; dut=dut*trn ; dum=dux*dxy dgradadp=dgradadp*pgn !der gradad/ln Pgaz dgradadt=dgradadt*trn !der gradad/ln T dgradadm=dgradadx*dxy !dgradad/m^1/3 dgradpt=dgradpt*prn/grad !der grad /ln Ptot à grad près dgradp=dgradp*pgn/grad !der grad /ln Pgaz à grad près dgradt=dgradt*trn/grad !der grad /ln T à grad près dgradr=dgradr/grad !der grad /r à grad près dgradl=dgradl/grad !der grad /l à grad près c d/dm^1/3 = 3 m^2/3 d/dm = 3 * ay52 d/dm à grad près dgradm=(dgradm*3.d0*ay52+dgradx*dxy)/grad !d grad/m^1/3 dgradlpt=dgradlpp*dlppdpt/grad !der grad /dln Ptot dgradlp =dgradlpp*dlppdp/grad !der grad /dln Pgaz dgampt=dgampt*prn !der gam /ln Ptot dgamp= dgamp*pgn !der gam /ln Pgaz dgamt= dgamt*trn !der gam /ln T c dgamr=dgamr/ay3 !der gam /r c dgaml=dgaml*ay4 !der gam /l dgamm=dgamm*3.d0*ay52+dgamx*dxy !der gam /m^1/3 dgamlpt=dgamlpp*dlppdpt !der gam /dln Ptot dgamlp =dgamlpp*dlppdp !der gam /dln Pgaz dhppt=dhppt*prn/hp !der hp /ln Ptot à hp près dhpp=dhpp*pgn/hp !der hp /ln Pgaz à hp près dhpt=dhpt*trn/hp !der hp /ln T à hp près dhpr=dhpr/hp !der hp /r à hp près dhpm=(dhpm*3.d0*ay52+dhpx*dxy)/hp !der hp/m^1/3 à hp près IF(epsilon(1) > 1.d-25)THEN depsp=depsp*pgn/epsilon(1)!dérivée /ln Pgaz, à 1/epsilon près depst=depst*trn/epsilon(1) !dérivée /ln T, à 1/epsilon près depsm=DOT_PRODUCT(depsx,dxchim)/epsilon(1) ELSE depsp=0.d0 ; depst=0.d0 ; depsm=0.d0 ENDIF c initialisations y_t=0.d0 tdst=0.d0 ; dtdsp=0.d0 ; dtdst=0.d0 ; dtdsm=0.d0 ; dpn=0.d0 ; dtn=0.d0 c sans pression turbulente Ptot=Pgaz c IF(.NOT.pturb)THEN c dgradpt=dgradp ; dgradlpt=dgradlp ; dgampt=dgamp c dgamlpt=dgamlp ; dhppt=dhpp c ENDIF c on ne tient pas compte de la différence Pgaz, Ptot pour l'éner grav c TdS=tds et old_ptm sont tabulés en fonction de m^1/3 c pour avoir un algorithme unique pour la diffusion, la composition c chimique est toujours tabulée en fonction de m^2/3 c d/dm^1/3=2 m^1/3 d/dm^2/3=2 ay5 d/dm^2/3 Idt: IF(dt == 0.d0)THEN tdst=0.d0 ; dtdsp=0.d0 ; dtdst=0.d0 ; dtdsm=0.d0 c interpolation du TdS au temps t ELSE Idt Icp: IF(compt == 0)THEN CALL bsp1dn(1,tds_t,x_tds_t,xt_tds_t,n_tds_t,m_tds, 1 knot_tds_t,.TRUE.,MIN(ay5,x_tds_t(n_tds_t)),l,ftds,dftds) IF(no_croiss)PRINT*,'Pb. at 2 in static_m13' c dtds/dm^1/3 au temps t (tabulation) tdst=ftds(1) ; dtdsm=dftds(1) c WRITE(*,2000)tdst,y(5,0),m_stat_t(n_qs_t) c PAUSE'TdS interpolé' ELSE Icp c perte de masse IF(l_pertm)THEN CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm, 1 .TRUE.,MIN(ay5,x_ptm(n_ptm)),l,fm,dfm) IF(no_croiss)PRINT*,'Pb. at 3 in static_m13' mk=fm(1) !m^1/3 dans static_m13 ELSE mk=ay5 ENDIF CALL inter('m',bp_t,q_t,qt_t,n_qs_t,knot_t,mk,y_t,dfqs,r_stat_t,m_stat_t) IF(pturb)THEN !avec pression turbulente 7 inconnues p_t=EXP(y_t(Ipg)) dp_tm=dfqs(Ipg)*p_t !dp_t / d m^1/3 au temps t dlnp=y(Ipg,0)-y_t(Ipg) dlnp_tm=dfqs(Ipg) ELSE !sans pression turbulente 6 inconnues p_t=EXP(y_t(1)) dp_tm=dfqs(1)*p_t !d p_t / d m^1/3 au temps t dlnp=y(1,0)-y_t(1) dlnp_tm=dfqs(1) !d ln p_t / d m^1/3 au temps t ENDIF t_t=EXP(y_t(2)) dt_tm=dfqs(2)*t_t !d t_t / d m^1/3 au temps t dlnt=y(2,0)-y_t(2) dlnt_tm=dfqs(2) !d ln t_t / d m^1/3 au temps t c limitation des variations temporelles de lnP et lnT IF(compt >= ini0 .AND. trn > 7.d7)reprend=ABS(dlnp) > d_press .OR. 1 ABS(dlnt) > d_temp IRp: IF(reprend)THEN WRITE(*,1)dlnp,d_press,dlnt,d_temp,ip,pgn,p_t,trn,t_t,ay53 1 FORMAT(/,'static_m13, variations temporelles trop importantes',/, 1 'dlnp=',es10.3,' > ',es10.3,', dlnt=',es10.3,' > ',es10.3,', couche:',i4,/, 2 'P=',es10.3,', P_t=',es10.3,', T=',es10.3,', T_t=',es10.3, 3 ', M/Msol=',es10.3) c vérification c i=FLOOR(y_t(6)) c PRINT*,bp_t(5,m_qs*(i-1)+1),mk,y_t(5),bp_t(5,m_qs*i+1),i c PRINT*,bp_t(1,m_qs*(i-1)+1),y(1,0),y_t(1),bp_t(1,m_qs*i+1),y_t(6) c PAUSE'dlnp' RETURN ENDIF Irp c variation locale de l'énergie de rotation SELECT CASE(Krot) CASE(0,1,2,6) eps_rot=0.d0 CASE(3,4,5) c omega_t au temps t CALL bsp1dn(nrot,rota_t,mrot_t,mrott_t,n_rot_t,ord_rot, 1 knotr_t,.TRUE.,mk**2,l,frot,dfrot) !en m^2/3 IF(no_croiss)PRINT*,'Pb. at 3-1 in static_m23' w_t=frot(1) c r_t au temps t r_t=y_t(3) c var. locale de l'énergie cinétique de rotation (ajoutée au TdS) eps_rot=cte4*((ay3*w)**2-(r_t*w_t)**2) END SELECT c Le TdS tdst = -epsilon_g dpn=pgn-p_t ; dtn=trn-t_t IF(kipp)THEN c Approximation de Kippenhahn forme 2: TdS = cp T ( dlnT - grad_ad dlnP )/dt 4-27 tdst=cp*trn*(dlnt-gradad*dlnp)/dt dtdsp=tdst*dcpp/cp-cp*trn*(dgradadp*dlnp+gradad)/dt dtdst=tdst*(dcpt/cp+1.d0)+cp*trn*(1.d0-dgradadt*dlnp)/dt dtdsm=tdst*dcpm/cp+cp*trn*(-dlnt_tm-dgradadm*dlnp+gradad*dlnp_tm)/dt c TdS, forme TdS=dU+PdV=dU-P/ro^2 dro ELSE c composition chimique au temps t-dt mk=mk**2 !m^2/3 CALL bsp1dn(nchim,chim_t,mc_t,mct_t,n_ch_t,m_ch, 1 knotc_t,.TRUE.,MIN(mk,mc_t(n_ch_t)),l,xchim_t,dxchim_t) IF(no_croiss)PRINT*,'Pb. at 4 in static_m13' WHERE(xchim_t < 1.d-38) xchim_t=1.d-38 dxchim_t=1.d-38 ENDWHERE CALL chim_gram(xchim_t,dxchim_t) c WRITE(*,2000)p_t,t_t,xchim_t(1),xcoll(cx) CALL etat(p_t,t_t,xchim_t,.TRUE., !mettre .TRUE. pour d/dX 1 ro_t,drop_t,drot_t,drox_t,u_t,dup_t,dut_t,dux_t, 2 delta_t,deltap_t,deltat_t,deltax_t,cp_t,dcpp_t,dcpt_t,dcpx_t, 3 gradad_t,dgradadp_t,dgradadt_t,dgradadx_t,alfa,beta,gamma1) c dU + P dV dro_tm=drop_t*dp_tm+drot_t*dt_tm+drox_t*cte3*ay5*dxchim_t(1) du_tm =dup_t *dp_tm+dut_t *dt_tm+dux_t *cte3*ay5*dxchim_t(1) tdst=(u-u_t-prn/ro**2*(ro-ro_t))/dt dtdsp=(dup-prn/ro**2*((ro-ro_t)*(1.d0-2.d0*drop/ro)+drop))/dt !d/ln P dtdst=(dut-prn/ro**2*drot*(1.d0-2.d0*(ro-ro_t)/ro))/dt !d/ln T dtdsm=(dum-du_tm-prn/ro**2*(-2.d0*drom/ro*(ro-ro_t)+drom-dro_tm))/dt c PRINT*,'cx/Ptot,Pgaz,T,p_t,t_t,u,u_t,ro,ro_t,tdst',cx c WRITE(*,2000)prn,trn,p_t,t_t,u,u_t,ro,ro_t,tdst ENDIF !kipp c variation de énergie cinétique de rotation (!!!!!incomplet) tdst=tdst-eps_rot/dt ENDIF Icp !compt=0 ENDIF Idt !dt /= 0 c la rotation, w ne dépend que de m c v correspond à l'accélération centrifuge w^2/6 pi R SELECT CASE(Krot) CASE(0,1,2) w=wrot ; dwdm=0.d0 ; dwdr=0.d0 CASE(3,4,5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(ay52,mrot(n_rot))),l,frot,dfrot) IF(no_croiss)PRINT*,'Pb. at 5 in static_m13' w=frot(1) ; dwdm=2.d0*ay5*dfrot(1) ; dwdr=0.d0 CASE(6) CALL omega_jpz(ay3,w,dwdr) ; dwdm=0.d0 END SELECT IF(w == 0.d0)THEN v=0.d0 ; dvdr=0.d0 ; dvdm=0.d0 ; dyv=0.d0 ; dyvdr=0.d0 dyvdm=0.d0 ELSE v=ay52/ay3/prn*w**2 dvdr=-v/ay3+2.d0*v/w*dwdr ; dvdm= v*(1.d0/ay52+2.d0*dwdm/w) dyv=cte25*fac(ip)*v !dérivée/lnP comme celle de dy1: dyvdr=cte25*fac(ip)*dvdr !dyv/dlnp=-dyv dyvdm=cte25*fac(ip)*dvdm ENDIF c la fonction de répartition avec ln Ptot, ln T, ln (r+1), ln (l+1), ln(m^1/3+1) dtetap=0.d0 ; dtetat=0.d0 ; dtetar=0.d0 ; dtetal=0.d0 ; dtetam=0.d0 c teta = ctep dln P/dm^1/3 IF(ctep /= 0.d0)THEN dy1=fac(ip)*ay55/ay34/prn teta=cte11*ctep*dy1 dtetap=-teta !dérivée/ln Ptot dtetar=-teta*4.d0/ay3 !dérivée/r dtetam= teta*5.d0/ay5 !dérivée/m^1/3 ENDIF c teta = teta + rotation IF(Krot > 2)THEN teta=teta+dyv dtetap=dtetap-dyv dtetar=dtetar+dyvdr dtetam=dtetam+dyvdm ENDIF c grad pouvant être < 0 , il faut |grad| pour avoir teta croissant c teta = teta + ctet |grad| dln T/dm^1/3 IF(ctet /= 0.d0)THEN sig=SIGN(1.d0,grad) ; sig=1.d0 dy2=cte11*ctet*dy1*sig*grad teta=teta+dy2 dtetap=dtetap+dy2*(dgradpt-1.d0) dtetat=dy2*dgradt dtetar=dtetar+dy2*(-4.d0/ay3+dgradr) dtetal=dy2*dgradl dtetam=dtetam+dy2*(5.d0/ay5+dgradm) ENDIF c teta = teta + cter d r / dm^1/3 IF(cter > 0.d0)THEN dy3=cte13*cter*fac(ip)/ro*ay52/ay32 teta=teta+dy3 dtetap=dtetap-dy3*drop/ro dtetat=dtetat-dy3*drot/ro dtetar=dtetar-2.d0*dy3/ay3 dtetam=dtetam+2.d0*dy3/ay5 ENDIF c teta = teta + ctem d m^1/3 / d m^1/3 teta=teta+fac(ip)*ctem c psi/teta psist=y(6,0)/teta dpsitp=-psist/teta*dtetap !d (psi/teta) / d ln Ptot dpsitt=-psist/teta*dtetat !d (psi/teta) / d ln T dpsitr=-psist/teta*dtetar !d (psi/teta) / d r dpsitl=-psist/teta*dtetal !d (psi/teta) / d l dpsitm=-psist/teta*dtetam !d (psi/teta) / d m^1/3 dpsisti=1.d0/teta !d (psi/teta) / d psi c équations (Prad est dans Pgas), avec Pturb, y(1,0) est Pgas+Pturb c sans Pturb, y(1,0) est Pgas c dln P/dq=-(3g/4pi Msol^2/ Rsol^4) m^5/3 / r^4 /P c dérivée/ln P comme celle de dy1: dyv/dlnp=-dyv dyv=cte15*v ; dyvdr=cte15*dvdr ; dyvdm=cte15*dvdm dy1=cte11*ay55/ay34/prn ; dy0=dy1+dyv be(1)=y(1,1)-dy0*psist ae(1,1,0)=dy0*(psist-dpsitp) !der./ln Ptot ae(1,1,1)=1.d0 !der./dln Ptot ae(1,2,0)=-dy0*dpsitt !der./dln T ae(1,3,0)=-(dyvdr-dy1*4.d0/ay3)*psist-dy0*dpsitr !dérivée/r ae(1,4,0)=-dy0*dpsitl !derivée /l ae(1,5,0)=-(dyvdm+dy1*5.d0/ay5)*psist-dy0*dpsitm !derivée/m^1/3 ae(1,6,0)=-dy0*dpsisti !dérivée/psi c dln T/dq=dln Ptot/dq grad dy2=dy0*psist*grad ; be(2)=y(2,1)-dy2 ae(2,1,0)=ae(1,1,0)*grad-dy2*dgradpt !dérivée/ln Ptot ae(2,2,0)=ae(1,2,0)*grad-dy2*dgradt !dérivée/ln T ae(2,2,1)=1.d0 !dérivée/dln T ae(2,3,0)=ae(1,3,0)*grad-dy2*dgradr !dérivée/r ae(2,4,0)=ae(1,4,0)*grad-dy2*dgradl !dérivée/l ae(2,5,0)=ae(1,5,0)*grad-dy2*dgradm !dérivée/m^1/3 ae(2,6,0)=ae(1,6,0)*grad !dérivée/psi IF(pturb)THEN !avec pression turb. 7 inc. ae(2,1,1)=-dy2*dgradlpt !dérivée/dln Ptot ae(2,Ipg,0)=-dy2*dgradp !dérivée/ln Pgaz ae(2,Ipg,1)=-dy2*dgradlp !dérivée/dln Pgaz ENDIF c dr/dm^1/3=3/4pi m^2/3 / r^2 ro dy3=cte13/ro*ay52/ay32; be(3)=y(3,1)-dy3*psist IF(pturb)THEN ae(3,1,0)=-dy3*dpsitp !dérivée/ln Ptot ae(3,Ipg,0)= dy3*drop/ro*psist !der. /ln Pgaz ELSE !sans press. turb. 6 inc. ae(3,1,0)=dy3*(drop/ro*psist-dpsitp) !der./ln Ptot ENDIF ae(3,2,0)=dy3*(psist*drot/ro-dpsitt) !dérivée/ln T ae(3,3,0)=dy3*(psist*2.d0/ay3-dpsitr) !dérivée/r ae(3,3,1)=1.d0 !dérivée/dr ae(3,4,0)=-dy3*dpsitl !dérivée/l ae(3,5,0)=dy3*(psist*(-2.d0/ay5+drom/ro)-dpsitm) !dérivée/m^1/3 ae(3,6,0)=-dy3*dpsisti !dérivée/psi c epsilon contribution de epsilon au dl/dmu dy4=cte14*epsilon(1)*ay52 be(4)=y(4,1)-dy4*psist IF(pturb)THEN !avec pression turbulente 7 inconnues ae(4,1,0)=-dy4*dpsitp !dérivée/ln Ptot ae(4,Ipg,0)=-dy4*psist*depsp !dérivée/ln Pgaz ELSE !sans pression turbulente 6 inconnues ae(4,1,0)=-dy4*(psist*depsp+dpsitp) !dérivée/ln Ptot ENDIF ae(4,2,0)=-dy4*(psist*depst+dpsitt) !dérivée/ln T ae(4,3,0)=-dy4*dpsitr !dérivée/r ae(4,4,0)=-dy4*dpsitl !der/l ae(4,4,1)=1.d0 !dérivée/dl ae(4,5,0)=-dy4*(psist*(2.d0/ay5+depsm)+dpsitm) !der/m^1/3 ae(4,6,0)=-dy4*dpsisti !dérivée/psi c contribution du Tds/dt au dl/dmu IF(dt /= 0.d0)THEN dy5=cte14t*ay52 d5tds=dy5*tdst be(4)=be(4)+d5tds*psist IF(compt == 0)THEN d5tdsm=d5tds*2.d0/ay5+dy5*dtdsm !dérivée/m^1/3 si tabulation ae(4,1,0)=ae(4,1,0)+d5tds*dpsitp !dérivée/ln Ptot ae(4,2,0)=ae(4,2,0)+d5tds*dpsitt !dérivée/ln T ae(4,3,0)=ae(4,3,0)+d5tds*dpsitr !dérivée/r ae(4,4,0)=ae(4,4,0)+d5tds*dpsitl !dérivée/l ae(4,5,0)=ae(4,5,0)+d5tds*dpsitm+psist*d5tdsm!dérivée/m^1/3 ae(4,6,0)=ae(4,6,0)+d5tds*dpsisti !dérivée/psi ELSE d5tdsp=dy5*dtdsp d5tdst=dy5*dtdst d5tdsm=dy5*dtdsm+d5tds*2.d0/ay5 IF(pturb)THEN !avec pression turbulente 7 inconnues ae(4,1,0)=ae(4,1,0)+d5tds*dpsitp !dérivée/ln Ptot ae(4,Ipg,0)=ae(4,Ipg,0)+psist*d5tdsp !dérivée/ln Pgaz ELSE !sans pression turbulente 6 inconnues ae(4,1,0)=ae(4,1,0)+d5tds*dpsitp+psist*d5tdsp ENDIF ae(4,2,0)=ae(4,2,0)+d5tds*dpsitt+psist*d5tdst !dérivée/ln T ae(4,3,0)=ae(4,3,0)+d5tds*dpsitr !dérivée/r ae(4,4,0)=ae(4,4,0)+d5tds*dpsitl !dérivée/l ae(4,5,0)=ae(4,5,0)+d5tds*dpsitm+psist*d5tdsm !d/m^1/3 ae(4,6,0)=ae(4,6,0)+d5tds*dpsisti !dérivée/psi ENDIF !compt = 0 ENDIF !dt c pour éviter L < 0 dL/dm=0 c IF(y(4,0) < 1.d-4 .AND. epsilon(1)-tdst < 0.d0)THEN c be(4)=y(4,1) c ae(4,:,:)=0.d0 c ae(4,4,1)=1.d0 c ENDIF c dmu/dq=psi/teta be(5)=y(5,1)-psist !d m^1/3 / dq = psi/teta ae(5,1,0)=-dpsitp !dérivée/ln Ptot ae(5,2,0)=-dpsitt !dérivée/ln T ae(5,3,0)=-dpsitr !dérivée/r ae(5,4,0)=-dpsitl !dérivée/l ae(5,5,0)=-dpsitm !dérivée/m^1/3 ae(5,5,1)= 1.d0 !dérivée/dm^1/3 ae(5,6,0)=-dpsisti !dérivée/psi c dpsi/dq=0 be(6)=y(6,1) !d psi / dq = 0 ae(6,6,1)=1.d0 !dérivée/dpsi c Avec Pturb: c pour les points de collocation dans les zones radiatives, c on resout d ln Ptot - d ln Pgas =0 c pour les points de collocation dans les zones convectives, c on resout -Pturb -Pgaz +Ptot=0 c IF(radiatif)WRITE(*,2000)trn,prn,pgn,gradrad,grad,gradad IF(pturb)THEN !avec pression turbulente 7 inconnues IF(radiatif)THEN be(Ipg)=y(1,1)-y(Ipg,1) !d lnPtot = d lnPgaz ae(Ipg,1,1)= 1.d0 !dérivée /ln Ptot ae(Ipg,Ipg,1)=-1.d0 !dérivée /ln Pgaz ELSE gam1=gam/(gam+1.d0) dgam1=1.d0/(gam+1.d0)**2/gam1 dy71=cte6*gam1*delta*prn*gradad*dlpp dy71pt=dy71*(dgampt*dgam1+1.d0) !der /ln Ptot dy71p= dy71*(dgamp*dgam1+deltap/delta+dgradadp/gradad)!der /ln Pgaz dy71t= dy71*(dgamt*dgam1+deltat/delta+dgradadt/gradad) !der /ln T dy71r= dy71*dgamr*dgamr !der /r dy71l= dy71*dgaml*dgam1 !der /l dy71m= dy71*(dgamm*dgam1+deltam/delta+dgradadm/gradad)!der /m^1/3 dy71lpt=dy71*(dgamlpt*dgam1+dlppdpt/dlpp) !der/dln Ptot dy71lp =dy71*(dgamlp *dgam1+dlppdp /dlpp) !der /dln Pgaz dy72=cte6*gam1*delta*prn*grad dy72pt=dy72*(dgampt*dgam1+1.d0+dgradpt) !der /ln Ptot dy72p= dy72*(dgamp*dgam1+deltap/delta+dgradp)!der /ln Pgaz dy72t= dy72*(dgamt*dgam1+deltat/delta+dgradt) !der /ln T dy72r= dy72*(dgamr*dgam1+dgradr) !der /r dy72l= dy72*(dgaml*dgam1+dgradl) !der /l dy72m= dy72*(dgamm*dgam1+deltam/delta+dgradm)!der /m^1/3 dy72lpt=dy72*(dgamlpt*dgam1+dgradlpt) !der/dln Ptot dy72lp= dy72*(dgamlp *dgam1+dgradlp) !der/dln Pgaz be(Ipg)=dy71-dy72-pgn+prn !-Pturb -Pgaz +Ptot=0 c WRITE(*,2000)be(Ipg),dy71,dy72,prn,pgn,dy71-dy72,prn-pgn ae(Ipg,1,0)=dy71pt-dy72pt+prn !der /ln Ptot ae(Ipg,2,0)=dy71t -dy72t !der /ln T ae(Ipg,3,0)=dy71r -dy72r !der /r ae(Ipg,4,0)=dy71l -dy72l !der /l ae(Ipg,5,0)=dy71m -dy72m !der /m^1/3 ae(Ipg,Ipg,0)=dy71p -dy72p-pgn !der /ln Pgaz ae(Ipg,1,1)=dy71lpt-dy72lpt !der /dln Ptot ae(Ipg,Ipg,1)=dy71lp -dy72lp !der /dln Pgaz ENDIF !radiatif ENDIF !Pturb c PRINT*,'be' ; WRITE(*,2000)be c WRITE(*,2000)cte13,ro,ay5,ay3,dy3,psist,y(3,1) c WRITE(*,2000)cte14,epsilon(1),ay5,ay4,dy4,psist,y(4,1) c PRINT*,'ae' c DO i=1,ne c WRITE(*,2000)ae(i,1:ne,0) c ENDDO c DO i=1,ne c WRITE(*,2000)ae(i,1:ne,1) c ENDDO c PAUSE'fin jacob' c---- test de vérification du jacobien par dérivées numériques----------- c deriv=.TRUE. c deriv=(cx <= 2 .OR. entre(195,200,cx)) .AND. compt > 0 c deriv=cx <= 3 .AND. compt > 0 c deriv=ip <= 3 .AND. compt > 0 c deriv=entre(1,3,cx) .OR. radiatif c deriv=ip <= 3 .OR. radiatif c deriv=entre(179,183,cx) c 1 .OR. entre(2.4d-2,2.6d-2,ay53) c deriv=entre(560,565,cx) .AND. compt > 0 c deriv=entre(340,343,ip) .AND. compt > 0 c 1 .OR. entre(93,95,ip) c 1 .OR. entre(95,98,cx) c 1 .OR. entre(661,665,cx) c deriv=cx <= 4 c deriv=ip == 7 .OR. cx <= 3 c deriv=cx <= 2 .OR. entre(560,565 ,cx) c deriv=entre(232,235,ip) .AND. compt >= 6 c deriv=entre(206,208,ip) .AND. compt >= 6 c deriv=entre(476 ,480 ,cx) c deriv=entre(243 ,245 ,cx) .AND. compt >= 6 c deriv=entre(656 ,661 ,cx) c deriv=entre(1.160d5,1.095d+05,trn) c deriv=entre(508,512,cx) .AND. compt > 0 c deriv=entre(727,729,ip) c 1 .OR. entre(1.44d+16,1.42d+16,prn) c 2 .OR. entre(3.35d+07, 3.37d+07,trn) c deriv=entre(914,916,cx) c 1 .OR. entre(16,19,cx) c 2 .OR. cx == 90 c 3 .OR.radiatif c deriv=cx == 1118 .OR. cx == 344 c deriv=ip == 2240 c deriv=.NOT.radiatif c deriv=radiatif .AND. compt > 0 c deriv=.NOT.radiatif .AND. compt > 0 c deriv=radiatif .AND. compt > 0 c deriv=MAXVAL(ABS(be(1:5))) == 0.d0 c deriv=MAXVAL(ABS(ae(:,:,0))) == 0.d0 c deriv=fac(ip)> 2.d0 c tests de mise au point des dérivées IF(.NOT. deriv)RETURN ALLOCATE(bes(ne),dern(ne),yd(ne,0:1)) ; yd=y ; xchim0=xchim dd=1.d-5 !; dd=1.d-8 IF(ay53 > 0.5d0)THEN unpdd=1.d0-dd ELSE unpdd=1.d0+dd ENDIF IF(radiatif)THEN PRINT*,'RADIATIF' ELSE PRINT*,'CONVECTIF' ENDIF PRINT*,'cx=',cx,', ip=',ip,', lmix:',lmix(ay52) PRINT*,'M=ay53,prn,pgn,trn,ay3,ay4,epsilon,ro,kipp: ',kipp WRITE(*,2000)ay53,prn,pgn,trn,ay3,ay4,epsilon(1),ro PRINT*,'be,xcoll(cx)' ; WRITE(*,2000)be(1:ne),xcoll(cx) PRINT*,'P,T,mstar,dt,cte6,alpha,dconv' WRITE(*,2000)prn,trn,mstar,dt,cte6,alpha,ctev*trn**3*gam/kap/ro**2/cp PRINT*,'kap,dkapp,dkapt,dkapx' WRITE(*,2000)kap,dkapp,dkapt,dkapx PRINT*,'cp,dcpp,dcpt,dcpm,dcpx' WRITE(*,2000)cp,dcpp,dcpt,dcpm,dcpx PRINT*,'ctem,ctep,cter,ctet,fac,Q' WRITE(*,2000)ctem,ctep,cter,ctet,fac(ip), 1 fac(ip)*(y(1,0)*ctep+y(5,0)*ctem+yd(3,0)*cter) PRINT*,'cte14,cte14t,cte21,cte22,cte23,cte25' WRITE(*,2000)cte14,cte14t,cte21,cte22,cte23,cte25 PRINT*,'cx,compt,pturb,der,n_qs,n_qs_t',cx,compt,pturb,der,n_qs,n_qs_t PRINT*,'delta,deltap,deltat,deltam,deltax' WRITE(*,2000)delta,deltap,deltat,deltam,deltax PRINT*,'dlnp,dlnt,dlpp,dlppdpt,dlppdp,dy71,dy72' WRITE(*,2000)dlnp,dlnt,dlpp,dlppdpt,dlppdp,dy71,dy72 PRINT*,'dy71pt,dy71t,dy71r,dy71l,dy71m,dy71p,dy71lpt,dy71lp' WRITE(*,2000)dy71pt,dy71t,dy71r,dy71l,dy71m,dy71p,dy71lpt,dy71lp PRINT*,'dy72pt,dy72t,dy72r,dy72l,dy72m,dy72p,dy72lpt,dy72lp' WRITE(*,2000)dy72pt,dy72t,dy72r,dy72l,dy72m,dy72p,dy72lpt,dy72lp PRINT*,'eps,depsp,depst,depsm/depsx' WRITE(*,2000)epsilon(1),depsp*epsilon(1),depst*epsilon(1), 1 depsm*epsilon(1) WRITE(*,2000)depsx PRINT*,'gam,dgamp,dgamt,dgamr,dgaml,dgamm,dgamx' WRITE(*,2000)gam,dgamp,dgamt,dgamr,dgaml,dgamm,dgamx PRINT*,'grad,dgradp,dgradt,dgradr,dgradl,dgradm,dgradx' WRITE(*,2000)grad,dgradp,dgradt,dgradr,dgradl,dgradm,dgradx PRINT*,'gradad,dgradadp,dgradadt,dgradadm,dgradadx,gradrad' WRITE(*,2000)gradad,dgradadp,dgradadt,dgradadm,dgradadx,gradrad PRINT*,'hp,dhpp,dhpt,dhpr,dhpm,dhpx' WRITE(*,2000)hp,dhpp,dhpt,dhpr,dhpm,dhpx PRINT*,'psist,dpsitp,dpsitt,dpsitr,dpsitl,dpsitm' WRITE(*,2000)psist,dpsitp,dpsitt,dpsitr,dpsitl,dpsitm PRINT*,'ro,drop,drot,drox,drom' WRITE(*,2000)ro,drop,drot,drox,drom PRINT*,'teta,dtetap,dtetat,dtetar,dtetal,dtetam' WRITE(*,2000)teta,dtetap,dtetat,dtetar,dtetal,dtetam PRINT*,'u,dup,dut,dux,dum,tau_rad' WRITE(*,2000)u,dup,dut,dux,dum,pgn*delta/(trn*ro*epsilon(1))/secon6 PRINT*,'tds,dtdsp,dtdst,dtdsm,dp_tm,dlnp_tm,dt_tm,dlnt_tm,compt=',compt WRITE(*,2000)tdst/secon6,dtdsp,dtdst,dtdsm,dp_tm,dlnp_tm,dt_tm,dlnt_tm IF(.NOT. kipp)THEN PRINT*,'ro_t,dro_tm,u_t,du_tm' WRITE(*,2000)ro_t,dro_tm,u_t,du_tm ENDIF PRINT*,'v,dvdr,dvdm,dyv,dyvdr,dyvdm,w' WRITE(*,2000)v,dvdr,dvdm,dyv,dyvdr,dyvdm,w PRINT*,'xchim/dxchim:xchimg' ; WRITE(*,2000)xchim ; WRITE(*,2000)dxchim WRITE(*,2000)xchim*nucleo PRINT*,'y , dy/dq, y_t , y-y_t' ; WRITE(*,2000)y(1:ne,0) WRITE(*,2000)y(1:ne,1) ; WRITE(*,2000)y_t(1:ne-1) WRITE(*,2000)(y(i,0)-y_t(i),i=1,ne-1) PRINT*,'dtn,dpn,trn,dlnt,dlnp,pgn' WRITE(*,2000)dtn,dpn,trn,dlnt,dlnp,pgn c IF(compt > 0 .AND. dt /= 0.d0)PRINT*,'initial fqs(6)=',fqs(6) c valeurs initiales pour les dérivées numériques cp0=cp ; dlnp0=dlnp ; dlnt0=dlnt ; dy710=dy71 ; dy720=dy72 eps0=epsilon(1) ; gam0=gam ; gradad0=gradad ; gradrad0=gradrad grad0=grad ; hp0=hp ; psist0=psist ; p_t0=p_t ; radia0=radiatif ro0=ro ; ro_t0=ro_t ; tdst0=tdst ; teta0=teta ; t_t0=t_t ; u0=u u_t0=u_t ; x00=xchim(1) c test dérivées c DO jd=0,1 !pour tester les dérivées/ y(i,0) et y(i,1) DO jd=0,0 !pour tester les dérivées/ y(i,0) DO id=1,ne stor0=yd(id,jd) ; stor=stor0*unpdd IF(ABS(stor) < dd)stor=SIGN(dd,stor) dstor=stor-stor0 ; yd(id,jd)=stor ay3=yd(3,0) ay32=ay3**2 ay34=ay32**2 ay4=yd(4,0) ay5=yd(5,0) !m^1/3 ay52=ay5**2 !m^2/3 ay53=ay52*ay5 !masse/msol ay55=ay52*ay53 !m^5/3 prn=EXP(yd(1,0)) !Ptot cgs IF(pturb)THEN !avec pression turbulente 7 inconnues pgn=EXP(yd(Ipg,0)) !Pgaz cgs ELSE !sans pression turbulente 6 inconnues pgn=prn ENDIF trn=EXP(yd(2,0)) !temperature K c PAUSE'avant bsp1dn' CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MIN(ay52,mc(n_ch)),l,xchim,dxchim) IF(no_croiss)PRINT*,'Pb. at 6 in static_m13' WHERE(xchim < 1.d-38) xchim=1.d-38 dxchim=1.d-38 ENDWHERE c IF(id == 5)xchim(ihe4-1:nchim)=xchim0(ihe4-1:nchim) c der=t on tient compte de dln Pgaz/dln Ptot IF(pturb .AND. der)THEN dlpp=yd(Ipg,1)/yd(1,1) !dlpp=dln Pgaz/dln Ptot ELSE dlpp=1.d0 !der=f on ne tient pas compte de dln Pgaz/dln Ptot ENDIF CALL thermo(prn,pgn,trn,ay53,ay4,ay3,dlpp,xchim,dxchim, 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 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,grad_mu,mu,dmu_x,mue, 8 gradrad,alfa,beta,gamma1,radiatif) 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 ay52,l,frot,dfrot) IF(no_croiss)PRINT*,'Pb. at 7 in static_m13' w=frot(1) CASE(6) CALL omega_jpz(ay3,w) END SELECT v=ay52/ay3/prn*w**2 !rotation IF(dt == 0.d0)THEN tdst=0.d0 ELSE IF(compt == 0)THEN CALL bsp1dn(1,tds_t,x_tds_t,xt_tds_t,n_tds_t,m_tds, 1 knot_tds_t,.TRUE.,MIN(ay5,x_tds_t(n_tds_t)),l,ftds,dftds) IF(no_croiss)PRINT*,'Pb. at 8 in static_m13' tdst=ftds(1) ELSE IF(l_pertm)THEN CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm, 1 .TRUE.,MIN(ay5,x_ptm(n_ptm)),l,fm,dfm) IF(no_croiss)PRINT*,'Pb. at 9 in static_m13' mk=fm(1) ELSE mk=ay5 ENDIF CALL inter('m',bp_t,q_t,qt_t,n_qs_t,knot_t,mk,fqs,dfqs, 1 r_stat_t,m_stat_t) IF(pturb)THEN !avec pression turbulente 7 inconnues p_t=EXP(fqs(Ipg)) ELSE !sans pression turbulente 6 inconnues p_t=EXP(fqs(1)) ENDIF t_t=EXP(fqs(2)) IF(kipp)THEN !approximation de Kippenhahan dlnp=yd(1,0)-fqs(1) ; dlnt=yd(2,0)-fqs(2) tdst=cp*trn*(dlnt-gradad*dlnp)/dt !4.27 ELSE !TdS=dU+PdV mk=mk**2 !m^2/3 CALL bsp1dn(nchim,chim_t,mc_t,mct_t,n_ch_t,m_ch, 1 knotc_t,.TRUE.,MIN(mk,mc_t(n_ch_t)),l,xchim_t,dxchim_t) IF(no_croiss)PRINT*,'Pb. at 10 in static_m13' WHERE(xchim_t < 1.d-38) xchim_t=1.d-38 dxchim_t=1.d-38 ENDWHERE CALL chim_gram(xchim_t,dxchim_t) CALL etat(p_t,t_t,xchim_t,.FALSE., 1 ro_t,drop_t,drot_t,drox_t,u_t,dup_t,dut_t,dux_t, 2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) tdst=(u-u_t-pgn/ro**2*(ro-ro_t))/dt !+Tds ENDIF !kipp ENDIF !compt=0 ENDIF !dt/=0 c teta teta=ay55/ay34/prn*(cte21+cte11*ctet*grad)+cte25*v 1 +cte13*cter/ro*ay52/ay32+ctem teta=fac(ip)*teta c psi/teta=dm/dq psist=yd(6,0)/teta c équations dy1=cte11*ay55/ay34/prn ; dyv=cte15*v ; dy0=dy1+dyv bes(1)=yd(1,1)-dy0*psist bes(2)=yd(2,1)-dy0*psist*grad !dln t/dq=dln p/dq grad psi/teta dy3=cte13/ro*ay52/ay32 bes(3)=yd(3,1)-dy3*psist !dl/dq=epsilon m13 psi/teta bes(4)=yd(4,1)-cte14*epsilon(1)*ay52*psist IF(dt /= 0.d0)THEN dy5=cte14t*ay52 bes(4)=bes(4)+tdst*dy5*psist !epsilon-tdS ENDIF !dt bes(5)=yd(5,1)-psist bes(6)=yd(6,1) IF(pturb)THEN IF(radiatif)THEN bes(Ipg)=yd(1,1)-yd(Ipg,1) ELSE gam1=gam/(gam+1.d0) ; dy71=cte6*gam1*delta*prn*gradad*dlpp dy72=cte6*gam1*delta*prn*grad ; bes(Ipg)=dy71-dy72-pgn+prn ENDIF !radiatif ENDIF !Pturb dern=(bes-be)/dstor ; PRINT* SELECT CASE(jd) CASE(0) PRINT*,'dérivée / ',variable(id) CASE(1) PRINT*,'dérivée / d',variable(id) END SELECT IF(radiatif)THEN PRINT*,'RADIATIF' ELSE PRINT*,'CONVECTIF' ENDIF IF(radiatif .NEQV. radia0)PRINT*,'chgt de CV: radia0= ',radia0 c IF(compt > 0 .AND. dt /= 0.d0)PRINT*,'deriv fqs(6)=',fqs(6) c WRITE(*,2000)dgradp,dgradp*pgn/grad,pgn,prn c WRITE(*,2000)dgradpt,dgradpt*prn/grad,be(2),bes(2) c WRITE(*,2000)teta,teta0,be(2),bes(2) c PRINT*,teta,teta0,grad,grad0,dstor c IF(id == 5)THEN c PRINT*,ro0,ro c PRINT*,xchim(1:3) c PRINT*,xchim0(1:3) c WRITE(*,2000)(xchim(1)-xchim0(1))/dstor,dxchim(1)*2.d0*ay5, c 1 dxchim(1),ay5,drox,dxchim(1)*ay5*drox*2.d0, c 2 (ro-ro0)/dstor,drom,dstor,(ro-ro0)/(xchim(1)-xchim0(1)), c 3 (ro-ro0)/(xchim(1)-xchim0(1))*dxchim(1)*ay5*2.d0 c ENDIF WRITE(*,2000)(dern(j),ae(j,id,jd),j=1,4) WRITE(*,2000)(dern(j),ae(j,id,jd),j=5,ne) yd(id,jd)=stor0 PRINT*,'grad,teta,eps,psist,ro,ro_t,u,u_t' WRITE(*,2000)(grad-grad0)/dstor/grad0,(teta-teta0)/dstor, 1 (epsilon(1)-eps0)/dstor,(psist-psist0)/dstor, 2 (ro-ro0)/dstor,(ro_t-ro_t0)/dstor,(u-u0)/dstor,(u_t-u_t0)/dstor PRINT*,'hp,tds,p_t,dlnp_t,t_t,dlnt_t' WRITE(*,2000)(hp-hp0)/dstor/hp0,(tdst-tdst0)/dstor,(p_t-p_t0)/dstor, 1 (dlnp-dlnp0)/dstor,(t_t-t_t0)/dstor,(dlnt-dlnt0)/dstor PRINT*,'X,gam,cp,gradrad,gradad,dy71,dy72' WRITE(*,2000)(xchim(1)-x00)/dstor,(gam-gam0)/dstor,(cp-cp0)/dstor, 1 (gradrad-gradrad0)/dstor/gradrad0,(gradad-gradad0)/dstor, 2 (dy71-dy710)/dstor,(dy72-dy720)/dstor PRINT*,'dtds,dgradrad,dgradad,dcp,deps' WRITE(*,2000)tdst-tdst0,gradrad-gradrad0,gradad-gradad0, 1 cp-cp0,epsilon(1)-eps0 ENDDO !jd PRINT*,'cx,ip,n_qs/fac(ip-1:ip+1)',cx,ip,n_qs WRITE(*,2000)fac(MAX(1,ip-1):MIN(n_qs,ip+1)) PRINT*,'psi, dpsi' ; WRITE(*,2000)y(6,0:1) PAUSE'test dérivées' ENDDO !id DEALLOCATE(bes,dern,yd) ; deriv=.FALSE. RETURN c les limites---------------------------------------- CASE(2) c conditions aux limites : résidu be(1) au point limite li c li=1 : limite sur r au centre c li=2 : limite sur l au centre c li=3 : limite sur m au centre c li=4 : limite sur Ptot au raccord c li=5 : limite sur T au raccord c li=6 : limite sur m au raccord c li=Ipg : limite sur Pgaz au raccord avec Pturb c test de dérivation pour les limites c deriv=.TRUE. IF(deriv)THEN ALLOCATE(yd(ne,0:1)) ; yd=y ENDIF c conditions limites SELECT CASE(li) c condition sur r au centre CASE(1) !au centre en r be(1)=y(3,0) !en q=1 r=0 ae(1,3,0)=1.d0 !dérivée/r c PRINT*,'limite au centre',li c WRITE(*,2000)y(1:ne,0),be(1) c condition sur l au centre CASE(2) !au centre en l be(1)=y(4,0) !en q=1 l=0 ae(1,4,0)=1.d0 !dérivée/l c PRINT*,'limite au centre',li c WRITE(*,2000)y(1:ne,0),be(1) c condition sur m au centre CASE(3) !au centre en m be(1)=y(5,0) !en q=1 m^1/3=0 ae(1,5,0)=1.d0 !dérivée/m^1/3 c PRINT*,'limite au centre',li c WRITE(*,2000)(y(i),i=1,ne),be(1) c PAUSE'au centre' c à partir de L et R le modèle d'atmosphère donne P, T, M et leurs dérivées / L et R c condition sur Ptot au raccord (sans Pturb Ptot=Pgaz) CASE(4) ay52=y(5,0)**2 CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MIN(ay52,mc(n_ch)),l,xchim0,dxchim) IF(no_croiss)PRINT*,'Pb. at 11 in static_m13' WHERE(xchim0 < 1.d-38) xchim0=1.d-38 dxchim0=1.d-38 ENDWHERE CALL chim_gram(xchim0,dxchim) ay3=y(3,0) ; ay4=y(4,0) CALL atm(.FALSE.,ay4,ay3,xchim0,ptext0,dptdl0,dptdr0, 1 text0,dtdl0,dtdr0,mext0,dmdl0,dmdr0,pext0,dpdl0,dpdr0,teff) be(1)=y(1,0)-LOG(ptext0) !condition sur Ptot au raccord ae(1,1,0)=1.d0 !dérivée/ln Ptot ae(1,3,0)=-dptdr0/ptext0 !dérivée/r ae(1,4,0)=-dptdl0/ptext0 !dérivée/l c test des dérivées IF(deriv)THEN !test dérivée dd=1.d-5 ; unpdd=1.d0-dd PRINT*,'Pext: y(1,0), Ln(Ptext0)' ; WRITE(*,2000)y(1,0), 1 LOG(ptext0) c WRITE(*,2000)y(1:ne,0) c PRINT*,'exterieur lnp,ln ptext,be(1),ay3,xchim0' c WRITE(*,2000)y(1,0),LOG(ptext),be(1),ay3,ln,xchim0 ptext=ptext0 ; text=text0 ; mext=mext0 ; pext=pext0 DO i=3,4 stor0=yd(i,0) ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 ; yd(i,0)=stor ay3=yd(3,0) ; ay4=yd(4,0) CALL atm(.FALSE.,ay4,ay3,xchim0,ptext,dptdl,dptdr, 1 text,dtdl,dtdr,mext,dmdl,dmdr,pext,dpdl,dpdr,teff) c PRINT*,LOG(ptext),yd(1,0),yd(1,0)-LOG(ptext),be(1) yd(i,0)=stor0 PRINT*,'dérivée pour P, lim=4 / ',variable(i) WRITE(*,2000)(yd(1,0)-LOG(ptext)-be(1))/dstor,ae(1,i,0) yd(i,0)=stor0 ENDDO !i PAUSE'test deriv' ENDIF !deriv c condition sur T au raccord CASE(5) be(1)=y(2,0)-LOG(text0) !condition sur T au raccord ae(1,2,0)=1.d0 !dérivée/ln t ae(1,3,0)=-dtdr0/text0 !dérivée/r ae(1,4,0)=-dtdl0/text0 !dérivée/l IF(deriv)THEN !test de dérivation dd=1.d-5 ; unpdd=1.d0-dd PRINT*,'Text: y(2,0), Ln(Text)' ; WRITE(*,2000)y(2,0),LOG(text) c WRITE(*,2000)y(1:ne,0) c PRINT*,'y(2,0),LOG(text),be(1)' c WRITE(*,2000)y(2,0),LOG(text),be(1) DO i=3,4 stor0=yd(i,0) ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 ; yd(i,0)=stor ay3=yd(3,0) ; ay4=yd(4,0) CALL atm(.FALSE.,ay4,ay3,xchim0,ptext,dptdl,dptdr, 1 text,dtdl,dtdr,mext,dmdl,dmdr,pext,dpdl,dpdr,teff) PRINT*,'dérivée pour T, lim=5 / ',variable(i) WRITE(*,2000)(yd(2,0)-LOG(text)-be(1))/dstor,ae(1,i,0) yd(i,0)=stor0 ENDDO !i PAUSE'test deriv' ENDIF !deriv c condition sur M au raccord CASE(6) be(1)=y(5,0)-mext0**(1.d0/3.d0) !en q=n_qs, m=mext(r,l) ae(1,3,0)=-dmdr0/3.d0/ay52 !dérivée/r ae(1,4,0)=-dmdl0/3.d0/ay52 !dérivée/l ae(1,5,0)=1.d0 !dérivée/m^1/3 IF(deriv)THEN !test de dérivation dd=1.d-3 ; unpdd=1.d0-dd PRINT*,'Mext: y(5,0), Mext0^1/3' WRITE(*,2000)y(5,0),mext0**(1.d0/3.d0) DO i=3,4 stor0=yd(i,0) ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 ; yd(i,0)=stor ay3=yd(3,0) ; ay4=yd(4,0) CALL atm(.FALSE.,ay4,ay3,xchim0,ptext,dptdl,dptdr, 1 text,dtdl,dtdr,mext,dmdl,dmdr,pext,dpdl,dpdr,teff) PRINT*,'dérivée pour M, lim=6 / ',variable(i) PRINT*,'mext0,mext,be(1)',mext0**(1.d0/3.d0), 1 mext**(1.d0/3.d0),yd(5,0)-mext**(1.d0/3.d0)-be(1) WRITE(*,2000)(yd(5,0)-mext**(1.d0/3.d0)-be(1))/dstor,ae(1,i,0) yd(i,0)=stor0 ENDDO !i PAUSE'test deriv' ENDIF !deriv c condition sur Pgaz au raccord CASE(7) be(1)=y(Ipg,0)-LOG(pext0) !en q=n ae(1,3,0)=-dpdr0/pext0 !dérivée/r ae(1,4,0)=-dpdl0/pext0 !dérivée/l ae(1,Ipg,0)=1.d0 !dérivée/ln Pgaz IF(deriv)THEN !test de dérivation PRINT*,'ptgaz ext: y(Ipg,0),pext0' ; WRITE(*,2000)y(Ipg,0),pext0 DO i=3,4 stor0=yd(i,0) ; stor=stor0*unpdd IF(stor == 0.d0)stor=dd dstor=stor-stor0 ; yd(i,0)=stor ay3=yd(3,0) ; ay4=yd(4,0) CALL atm(.FALSE.,ay4,ay3,xchim,ptext,dptdl,dptdr, 1 text,dtdl,dtdr,mext,dmdl,dmdr,pext,dpdl,dpdr,teff) PRINT*,'dérivée pour Pt, lim=7 / ',variable(i) WRITE(*,2000)(yd(Ipg,0)-LOG(pext)-be(1))/dstor,ae(1,i,0) yd(i,0)=stor0 ENDDO !i PAUSE'test deriv' ENDIF !deriv CASE DEFAULT WRITE(*,"('static_m13 erreur sur li =',i3,', ne=',i3)")li,ne PRINT*,'ARRET' STOP END SELECT IF(deriv)THEN DEALLOCATE(yd) ; deriv=.FALSE. ENDIF c PRINT*,'li=',li ; PAUSE'fin des limites' CASE DEFAULT PRINT*,' static_m13 erreur sur fait = 1, ou 2, fait=',fait PRINT*,'ARRET' ; STOP END SELECT !fait reprend=.FALSE. RETURN END SUBROUTINE static_m13