c******************************************************************** SUBROUTINE static_m23(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 lagrangien, 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/Msol)^2/3 c que ce soit en lagrangien ou en eulérien c l'énergie graviphique, Eg=-TdS/dt=(dU+PdV)/dt=tds est tabulée en fonction c de m^2/3 ou de m 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 trop forte ou R, L, M < 0 c ou acc. centrifuge trop grande c---------------------------------------------------------------- USE mod_donnees, ONLY : alpha, cpturb, dtmin, d_press, d_temp, 1 fcv, g, ihe4, ini0, Ipg, Krot, kipp, langue, lsol, l_pertm, msol, 2 m_ch, m_ptm, m_tds, nchim, ne, nrot, nucleo, ord_rot, pi, 3 pturb, rsol, secon6, t_inf 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, bp_t, chim, chim_gram, chim_t, ctem, 1 ctep, cter, ctet, inter, knotc, knotc_t, knotr, knotr_t, knot_tds_t, 2 knot_ptm, knot_t, mc, mct, mct_t, mc_t, mrot, mrot_t, mrott, mrott_t, 3 mstar, m_stat_t, n_ch, n_ch_t, n_ptm, n_qs, n_qs_t, n_rot, n_rot_t, 4 n_tds_t, old_ptm, omega_jpz, 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, 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), PARAMETER :: dd=1.d-5 REAL (kind=dp), SAVE :: cte1, cte2, cte3, cte6, cte11, cte13, cte14, 1 cte14t, cte15, cte21, cte22, cte23, cte25, dmdl0, dmdr0, dpdl0, 2 dpdr0, dptdl0, dptdr0, dtdl0, dtdr0, dt_tds, l13, mext0, m13, 3 pext0, ptext0, ray, text0, unpdd REAL (kind=dp) :: acc_c, alfa, ay3, ay4, ay5, beta, cp, 1 dcpm, dcpp, dcpt, dcpx, delta, deltam, deltap, deltat, deltax, 2 depsm, depsp, depst, dgamlpt, dgamp, 3 dgampt, dgamt, dgamx, dgamm, dgaml, dgamlp, dgamlpp, dgamr, 4 dgam1, gradad, dgradadm, dgradadp, dgradadt, dgradadx, 5 dgradlp, dgradlpt, dgradp, dgradpt, dgradt, dgradx, dgradm, dgradl, 6 dgradr, dgradlpp, dhpp, dhppt, dhpt, dhpx, dhpr, dhpm, dkapp, dkapt, 7 dkapx, dlnp, dlnp_tm, dlnt, dlnt_tm, dlpp, dlppdp, dlppdpt, dmdl, 8 dmdr, dpdl, dpdr, dpsitl, dpsitm, dpsitp, dpsitr, dpsitt, dptdl, 9 dptdr, dp_tm=0.d0, drom, drop, drop_t, drot, drot_t=0.d0, drox, 1 drox_t=0.d0, dro_tm=0.d0, dstor, dtdl, dtdr, dtdsl=0.d0, dtdsm=0.d0, 2 dtdsp=0.d0, dtdst=0.d0, dtdst_m=0.d0, 3 dtetal, dtetam, dtetap, dtetar, dtetat, dt_tm=0.d0, 4 dum, dup, dup_t=0.d0, dut, dut_t=0.d0, dux, dux_t=0.d0, 5 du_tm=0.d0, dvdm, dvdr, dwdm, dwdr, dxy, dyv, dyvdm, dyvdr, 6 dy0, dy1, dy2, dy3, dy4, dy5, dy71=0.d0, dy710, dy71l=0.d0, 7 dy71lp=0.d0, dy71lpt=0.d0, dy71m=0.d0, dy71p=0.d0, dy71pt=0.d0, 8 dy71r=0.d0, dy71t=0.d0, dy72=0.d0, dy720=0.d0, dy72l=0.d0, 9 dy72lp=0.d0, dy72lpt=0.d0, dy72p=0.d0, dy72pt=0.d0, dy72m=0.d0, 1 dy72r=0.d0, dy72t=0.d0, d5tds, d5tdsm, d5tdsp, d5tdst, eps0, 2 eps_rot, gam, gam0, gam1, gamma1, grad, grad0 , gradrad, grad_mu, 3 gravs2, hp, kap, ln, mext, mk, mn, msmst,mu, mue, pext, pgn, prn, 4 psist, psist0, ptext, p_t=0.d0, p_t0, ro, ro0, ro_t=0.d0, ro_t0, 5 rsrst, r_t, stor, stor0, tdst=0, tds0, teff, teta, teta0, 6 text, trn, t_t=0.d0, 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 CHARACTER (len=7), SAVE, DIMENSION(7) :: variable=(/ 'ln Ptot', 1 'ln T ','R**2 ','L**2/3 ','M**23 ','Psi ','ln Pgaz'/) c------------------------------------------------------------------------ 2000 FORMAT(8es10.3) c PRINT* ; PRINT*,'entree static_m23',fait,compt,cx,ip 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. cte1=g*msol/rsol**2 cte2=2.d0/3.d0*rsol cte3=rsol**2/3.d0/secon6 cte11=-3.d0*g/8.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=msol/lsol !pour d lambda/d q cte14t=cte14/secon6 ! " " dt /= 0 cte15=msol/rsol/4.d0/pi !pour la rotation cte25=cte15*ctep cte6=ABS(cpturb)*alpha**2/8.d0 !pour la Pturb der=cpturb < 0.d0 !der=.TRUE. on tient compte de dlnPgaz/dlnPtot dt_tds=10.d0*dtmin !dt < dt_tds on utilise la tabulation TdS ALLOCATE(xchim0(nchim)) !pour l'atmosphère 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**2 y(3,1)=(r**2)' c y(4,0)=l**2/3 y(4,1)=(l**2/3)' c y(5,0)=m**2/3 y(5,1)=(m**2/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) CASE(1) !le point courant 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=MAX(y(3,0),1.d-8) !r**2 pour éviter des valeurs < 0 à r**2 c ay3=ABS(y(3,0)) !r**2 ray=SQRT(ay3) !rayon en Rsol rsrst=ray*rstar !rayon en Rstar c ay4=MAX(y(4,0),1.d-8) !l**2/3 pour éviter des valeurs < 0 à l**2/3 ay4=y(4,0) !l**2/3 l13=SIGN(SQRT(ABS(ay4)),ay4) !l**1/3 ln=l13**3 !l/lsol ay5=MAX(y(5,0),1.d-8) !m**2/3 pour éviter des valeurs < 0 à m**2/3 c ay5=ABS(y(5,0)) !m**2/3 m13=SQRT(ay5) !m**1/3 mn=m13**3 !masse en Msol msmst=mn*mstar !masse en Mstar c PRINT*,'mn=',mn c l'acc. centrifuge ne doit pas excèder 90% gravité gravs2=cte1*mn/ay3*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(ay5,mrot(n_rot))),l,frot,dfrot) IF(no_croiss)PRINT*,'Pb. en 0 dans static_m23' w=frot(1) CASE(6) CALL omega_jpz(ray,w) END SELECT c accélération centrifuge acc_c=cte2*ray*w**2 ; reprend=acc_c > gravs2 IF(reprend)THEN SELECT CASE(langue) CASE('english') WRITE(*,1012)acc_c,gravs2,w,rsrst,msmst WRITE(2,1012)acc_c,gravs2,w,rsrst,msmst 1012 FORMAT('STOP, in static_m23, the centrifugal acceleration =', 1 es10.3,' > 90% gravity = ',es10.3,/,'angular velicity =',es10.3, 2 ', R/Rstar=',es10.3,', M/Mstar=',es10.3) CASE DEFAULT WRITE(*,12)acc_c,gravs2,w,rsrst,msmst WRITE(2,12)acc_c,gravs2,w,rsrst,msmst 12 FORMAT('ARRET, dans static_m23, acc. centifuge =', 1 es10.3,' > 90% gravité = ',es10.3,/,'vitesse angulaire =',es10.3, 2 ', R/Rstar=',es10.3,', M/Mstar=',es10.3) END SELECT CALL sortie('static_m23') ENDIF c si, à cause d'une erreur d'arrondi, ce qui peut arriver au voisinage c du centre, r**2, l**2/3 ou m**23 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 Si la difficulté subsiste il convient d'utiliser static_m13 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 IF(reprend)THEN c on tente de converger IF(fcv)THEN reprend=.FALSE. SELECT CASE(langue) CASE('english') WRITE(*,1014)y(3:5,0),SQRT(ABS(y(3,0))), 1 SQRT(ABS(y(4,0)))**3,SQRT(ABS(y(5,0)))**3,cx 1014 FORMAT(/,'static_m23, r2=',es10.3,', l23=',es10.3,', m23=',es10.3,/, 1 'R=',es10.3,', L=',es10.3,', M=',es10.3,/, 2 'cx=',i4,', one try to converge, better: use static_m13') CASE DEFAULT WRITE(*,14)y(3:5,0),SQRT(ABS(y(3,0))), 1 SQRT(ABS(y(4,0)))**3,SQRT(ABS(y(5,0)))**3,cx 14 FORMAT(/,'static_m23, r2=',es10.3,', l23=',es10.3,', m23=',es10.3,/, 1 'R=',es10.3,', L=',es10.3,', M=',es10.3,/, 2 'cx=',i4,', on tente de converger, mieux: utiliser static_m13') END SELECT c dt --> dt/2 ELSE SELECT CASE(langue) CASE('english') WRITE(*,1016)y(3:5,0),SQRT(ABS(y(3,0))), 1 SQRT(ABS(y(4,0)))**3,SQRT(ABS(y(5,0)))**3,cx 1016 FORMAT(/,'static_m23, r2=',es10.3,', l23=',es10.3,', m23=',es10.3,/, 1 'R=',es10.3,', L=',es10.3,', M=',es10.3,/, 2 'cx=',i4,', dt --> dt/2') CASE DEFAULT WRITE(*,16)y(3:5,0),SQRT(ABS(y(3,0))), 1 SQRT(ABS(y(4,0)))**3,SQRT(ABS(y(5,0)))**3,cx 16 FORMAT(/,'static_m23, r2=',es10.3,', l23=',es10.3,', m23=',es10.3,/, 1 'R=',es10.3,', L=',es10.3,', M=',es10.3,/, 2 'cx=',i4,', dt --> dt/2') END SELECT RETURN ENDIF ENDIF c composition chimique au temps t+dt CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MIN(ay5,mc(n_ch)),l,xchim,dxchim) IF(no_croiss)PRINT*,'Pb. at 1 in static_m23' WHERE(xchim< 1.d-38)xchim=1.d-38 c IF(SUM(xchim*nucleo) > 1.001d0)THEN c WRITE(*,2000)SUM(xchim*nucleo) ; WRITE(*,2000)xchim ; PAUSE'static_m23' c ENDIF c PRINT*,'cx,nchim/prn,pgn,trn,mn,ray / xchim / y /dy',cx,nchim c WRITE(*,2000)prn,pgn,trn,mn,ray,xchim(1),dxchim(1) c WRITE(*,2000)xchim ; WRITE(*,2000)dxchim c WRITE(*,2000)y(1:ne,0) ; WRITE(*,2000)y(1:ne,1) c PRINT*,' ' 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,mn,ln,ray,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)mn,prn,pgn,trn,ro,xchim,epsilon(1),depsp,depst,depsx c PRINT*,'couche',cx ; WRITE(*,2000)prn,pgn,trn,mn,ln,ray,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 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, mu, ln Pgaz drop=drop*pgn !der ro /ln Pgaz drot=drot*trn ; drom=drox*dxy 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/gradad !der gradad/ln Pgaz gradad près dgradadt=dgradadt*trn/gradad !der gradad/ln T à gradad près dgradadm=dgradadx*dxy/gradad !dgradad/m**2/3 gradad près 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 dgradr=dgradr/2.d0/ray/grad !der grad /r**2 dgradl=dgradl*3.d0/2.d0*l13/grad !der grad /l**2/3 dgradm=(dgradm*3.d0/2.d0*m13+dgradx*dxy)/grad !d grad/m**2/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 dgamr=dgamr/2.d0/ray !der gam /r**2 dgaml=dgaml*3.d0/2.d0*l13 !der gam /l**2/3 dgamm=dgamm*3.d0/2.d0*m13+dgamx*dxy !der gam /m**2/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/2.d0/ray/hp !der hp /r**2 à hp près dhpm=(dhpm*3.d0/2.d0*m13+dhpx*dxy)/hp !der hp/m**2/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) c depsm=SUM(depsx(1:nchim)*dxchim(1:nchim))/epsilon(1) 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 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 c l'énergie graviphique TdS=tds et old_ptm sont tabules en fonction de m^2/3 c pour avoir un algorithme unique pour la diffusion, la composition c chimique est toujours tabulée en fonction de m^2/3 Idt: IF(dt == 0.d0)THEN tdst=0.d0 ; dtdsp=0.d0 ; dtdst=0.d0 ; dtdsm=0.d0 ELSE Idt Icp: IF(compt == 0 .OR. dt <= dt_tds)THEN !interp. du TdS au temps t 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_m23' c dtds/dm^2/3 au temps t (tabulation) tdst=ftds(1) ; dtdst_m=dftds(1) ; dtdsp=0.d0 ; dtdst=0.d0 c WRITE(*,2000)tdst,y(5,0),m_stat_t(n_qs_t) c PAUSE'TdS interpolé' ELSE Icp 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_m23' IF(l_pertm)THEN CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,.TRUE., 1 MIN(ay5,x_ptm(n_ptm)),l,fm,dfm) IF(no_croiss)PRINT*,'Pb. at 3 in static_m23' mk=fm(1) 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 !d p_t / d mu 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 mu 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 mu 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,msmst 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,', M/M*=',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,knotr_t, 1 .TRUE.,mk,l,frot,dfrot) IF(no_croiss)PRINT*,'Pb. at 3-1 in static_m23' w_t=frot(1) c r_t au temps t r_t=SQRT(fqs(3)) c variation locale de l'énergie cinétique de rotation (ajoutée au TdS) eps_rot=cte3*((ray*w)**2-(r_t*w_t)**2) END SELECT c le TdS IF(kipp)THEN !approximation de Kippenhahan 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 ELSE !TdS=dU+PdV 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_m23' WHERE(xchim_t < 1.d-38)xchim_t=1.d-38 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,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1) tdst=(u-u_t-prn/ro**2*(ro-ro_t))/dt du_tm =dup_t *dp_tm+dut_t *dt_tm+dux_t *dxchim_t(1) dro_tm=drop_t*dp_tm+drot_t*dt_tm+drox_t*dxchim_t(1) dtdsp=(dup-pgn/ro**2*((ro-ro_t)*(1.d0-2.d0*drop/ro)+drop))/dt dtdst=(dut-pgn/ro**2*drot*(1.d0-2.d0*(ro-ro_t)/ro))/dt dtdsm= (dum-du_tm-pgn/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(ay5,mrot(n_rot))),l,frot,dfrot) IF(no_croiss)PRINT*,'Pb. at 5 in static_m23' w=frot(1) ; dwdm=dfrot(1) ; dwdr=0.d0 CASE(6) CALL omega_jpz(ray,w,dwdr) ; dwdm=0.d0 ; dwdr=dwdr/2.d0/ray 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=SQRT(ay5/ay3)/prn*w**2 dvdr=-v/2.d0/ay3+2.d0*v/w*dwdr ; dvdm= v/2.d0/ay5+2.d0*v/w*dwdm 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 sur ln Ptot, ln T, l^2/3 et m^2/3 dtetap=0.d0 ; dtetat=0.d0 ; dtetar=0.d0 ; dtetal=0.d0 ; dtetam=0.d0 c teta = ctep dksi/dmu IF(ctep /= 0.d0)THEN dy1=fac(ip)*(ay5/ay3)**2/prn teta=cte11*ctep*dy1 dtetap=-teta !dérivée/ln Ptot dtetar=-teta*2.d0/ay3 !dérivée/r**2 dtetam=teta*2.d0/ay5 !dérivée/m**2/3 ENDIF c teta = teta + rotation IF(Krot > 2)THEN teta=teta+dyv dtetap=dtetap-dyv dtetar=dtetar+dyvdr dtetam=dtetam+dyvdm ENDIF c teta = teta + ctet dksi/dmu grad IF(ctet /= 0.d0)THEN dy2=cte11*ctet*dy1*grad teta=teta+dy2 dtetap=dtetap+dy2*(dgradpt-1.d0) dtetat=dy2*dgradt dtetar=dtetar+dy2*(-2.d0/ay3+dgradr) dtetal=dy2*dgradl dtetam=dtetam+dy2*(2.d0/ay5+dgradm) ENDIF c teta = teta + cter d zeta/dmu IF(cter /= 0.d0)THEN dy3=cte13*cter*fac(ip)/ro*m13/ray teta=teta+dy3 dtetap=dtetap-dy3*drop/ro dtetat=dtetat-dy3*drot/ro dtetar=dtetar-dy3/ay3/2.d0 dtetam=dtetam+dy3/ay5/2.d0 ENDIF c teta = teta + ctem d mu/ d mu 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**2 dpsitl=-psist/teta*dtetal !d (psi/teta) / d l**2/3 dpsitm=-psist/teta*dtetam !d (psi/teta) / d m**2/3 c équations (Prad est dans Pgas), avec Pturb, y(1,0) est Pgas+Pturb c sans Pturb, y(1,0) est Pgas c dérivée/ln P comme celle de dy1: dyv/dlnp=-dyv dyv=cte15*v ; dyvdr=cte15*dvdr ; dyvdm=cte15*dvdm dy1=cte11*(ay5/ay3)**2/prn ; dy0=dy1+dyv be(1)=y(1,1)-dy0*psist !dln Ptot/dmu=3g/8 pi (mu+1/zeta+1)**2/P 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*2.d0/ay3)*psist-dy0*dpsitr !dérivée/r**2 ae(1,4,0)=-dy0*dpsitl !der. /l**2/3 ae(1,5,0)=-(dyvdm+dy1*2.d0/ay5)*psist-dy0*dpsitm !der. /m**2/3 ae(1,6,0)=-dy0/teta !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**2 ae(2,4,0)=ae(1,4,0)*grad-dy2*dgradl !dérivée/l**2/3 ae(2,5,0)=ae(1,5,0)*grad-dy2*dgradm !dérivée/m**2/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**2/dmu=3/4pi /ro SQRT(mu/zeta) dy3=cte13/ro*m13/ray; 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/ay3/2.d0-dpsitr) !dérivée/r**2 ae(3,3,1)=1.d0 !dérivée/dr**2 ae(3,4,0)=-dy3*dpsitl !dérivée/l**2/3 ae(3,5,0)=dy3*(psist*(-0.5d0/ay5+drom/ro)-dpsitm) !dérivée/m**2/3 ae(3,6,0)=-dy3/teta !dérivée/psi c dl/dmu=epsilon m13/l13 +Tds/dt contribution de epsilon dy4=cte14*epsilon(1)*m13/l13 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**2 ae(4,4,0)= dy4*(psist/2.d0/ay4-dpsitl) !der/l**2/3 ae(4,4,1)=1.d0 !dérivée/dl**2/3 ae(4,5,0)=-dy4*(psist*(0.5d0/ay5+depsm)+dpsitm)!der/m**2/3 ae(4,6,0)=-dy4/teta !dérivée/psi c contribution du Tds/dt IF(dt /= 0.d0)THEN dy5=cte14t*m13/l13 d5tds=dy5*tdst be(4)=be(4)+d5tds*psist !epsilon-(tdS+d Iw2)/dt dtdsl=-d5tds/2.d0/ay4 !dérivée/l**2/3 dtdsm= d5tds/2.d0/ay5 !dérivée/m**2/3 IF(compt == 0 .OR. dt <= dt_tds)THEN dtdsm= dtdsm+dy5*dtdst_m !dérivée/m**2/3 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**2 ae(4,4,0)=ae(4,4,0)+d5tds*dpsitl+psist*dtdsl !dérivée/l**2/3 ae(4,5,0)=ae(4,5,0)+d5tds*dpsitm+psist*dtdsm !dérivée/m**2/3 ae(4,6,0)=ae(4,6,0)+d5tds/teta !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**2 ae(4,4,0)=ae(4,4,0)+d5tds*dpsitl+psist*dtdsl !dérivée/l**2/3 ae(4,5,0)=ae(4,5,0)+d5tds*dpsitm+psist*d5tdsm !dérivée/m**2/3 ae(4,6,0)=ae(4,6,0)+d5tds/teta !dérivée/psi ENDIF !compt = 0 ENDIF !dt c dmu/dq=psi/teta be(5)=y(5,1)-psist 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**2 ae(5,4,0)=-dpsitl !dérivée/l**2/3 ae(5,5,0)=-dpsitm !dérivée/m**2/3 ae(5,5,1)= 1.d0 !dérivée/dm**2/3) ae(5,6,0)=-1.d0/teta !dérivée/psi c dpsi/dq=0 be(6)=y(6,1) ; 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) !der /ln Pgaz dy71t= dy71*(dgamt*dgam1+deltat/delta+dgradadt) !der /ln T dy71r= dy71*dgamr*dgam1 !der /r**2 dy71l= dy71*dgaml*dgam1 !der /l**2/3 dy71m= dy71*(dgamm*dgam1+deltam/delta+dgradadm) !der /m**2/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**2 dy72l= dy72*(dgaml*dgam1+dgradl) !der /l**2/3 dy72m= dy72*(dgamm*dgam1+deltam/delta+dgradm) !der /m**2/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**2 ae(Ipg,4,0)=dy71l -dy72l !der /l**2/3 ae(Ipg,5,0)=dy71m -dy72m !der /m**2/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),m13,l13,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 c deriv=cx <= 3 c deriv=cx == 1 c deriv=cx == 90 c deriv=entre(705,708,cx) c deriv=cx >= 867 .AND. cx <= 870 c deriv=cx >= 925 .AND. cx <= 930 c deriv=cx >= 345 .AND. cx <= 349 c 1 .OR. cx == 30 c 2 .OR. cx == 90 c 2 .OR. cx == 189 c 2 .OR. cx == 1090 c 2 .OR. cx >= 188 .AND. cx <= 190 c 2 .OR. cx >= 199 .AND. cx <= 202 c 3 .OR. (cx > 74 .AND. cx < 76) c 3 .OR. (cx > 43 .AND. cx < 48) c 4 .OR. (cx > 88 .AND. cx < 92) c 5 .OR. cx == 220 c 6 .OR. cx == 103 c 7 .OR. cx. ge. 148 c 8 .OR. cx. ge. 149 c 9 .OR. (cx. ge. 148 .AND. cx <150) c 1 .OR. cx >= 380 c 2 .OR. ABS(fac(ip)-1.d0) > 1.d-5 c 2 .OR. cx >= 358 .AND. cx <= 362 c 3 .OR.radiatif c deriv=cx == 100 c deriv=.NOT.radiatif c deriv= cx >= 299 .AND. cx <= 303 IF(.NOT. deriv)RETURN PRINT* PRINT*,'tests de mise au point pour les dérivées' ALLOCATE(bes(ne),dern(ne),yd(ne,0:1)) ; yd=y IF(mn > 0.5d0)THEN unpdd=1.d0-dd ELSE unpdd=1.d0+dd ENDIF IF(radiatif)THEN PRINT*,'RADIATIF' ELSE PRINT*,'CONVECTIF' ENDIF PRINT*,'cx,compt,radiatif,pturb,der,n_qs_t', 1 cx,compt,pturb,der,n_qs PRINT*,'ctem,ctep,cter,ctet,fac,Q' WRITE(*,2000)ctem,ctep,cter,ctet,fac(ip), 1 fac(ip)*(yd(1,0)*ctep+yd(5,0)*ctem+yd(3,0)*cter) PRINT*,'cte21,cte22,cte23,cte25' WRITE(*,2000)cte21,cte22,cte23,cte25 PRINT*,'mstar,dt,cte6,alpha' ; WRITE(*,2000)mstar,dt,cte6,alpha PRINT*,'be,xcoll(cx)' ; WRITE(*,2000)be(1:ne),xcoll(cx) PRINT*,'y / dy' ; WRITE(*,2000)yd(1:ne,0) ; WRITE(*,2000)yd(1:ne,1) PRINT*,'mn,prn,pgn,trn,ray,ln,epsilon,ro' WRITE(*,2000)mn,prn,pgn,trn,ray,ln,epsilon(1),ro PRINT*,'xchim(i)' ; WRITE(*,2000)xchim PRINT*,'dxchim(i)' ; WRITE(*,2000)dxchim PRINT*,'despx(i)' ; WRITE(*,2000)depsx PRINT*,'kap,dkapp,dkapt,dkapx' WRITE(*,2000)kap,dkapp,dkapt,dkapx PRINT*,'grad,dgradx,gradrad,gradad' WRITE(*,2000)grad,dgradx,gradrad,gradad PRINT*,'dgradpt,dgradt,dgradr,dgradl,dgradm,dgradp,dgradlpt,dgradlp' WRITE(*,2000)dgradpt,dgradt,dgradr,dgradl,dgradm,dgradp,dgradlpt,dgradlp PRINT*,'dgradlpp,gradad,dgradadp,dgradadt,dgradadx,t_t,p_t' WRITE(*,2000)dgradlpp,gradad,dgradadp,dgradadt,dgradadx,t_t,p_t PRINT*,'psist,teta,dpsitp,dpsitt,dpsitr,dpsitl,dpsitm,dpsist' WRITE(*,2000)psist,teta,dpsitp,dpsitt,dpsitr,dpsitl,dpsitm,1.d0/teta PRINT*,'dtetap,dtetat,dtetar,dtetal,dtetam,dp_tm,dt_tm' WRITE(*,2000)dtetap,dtetat,dtetar,dtetal,dtetam,dp_tm,dt_tm PRINT*,'depsp,depst,depsm,tdst,dtdsp,dtdst,dtdsl,dtdsm' WRITE(*,2000)depsp*epsilon(1),depst*epsilon(1),depsm*epsilon(1),tdst, 1 dtdsp,dtdst,dtdsl,dtdsm PRINT*,'tdst,dlpp,dlppdpt,dlppdp,dy71,dy72' WRITE(*,2000)tdst,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*,'u,dup,dut,dux,dum,u_t,du_tm,dtdst_m' WRITE(*,2000)u,dup,dut,dux,dum,u_t,du_tm,dtdst_m PRINT*,'ro,drop,drot,drox,drom,ro_t,dro_tm' WRITE(*,2000)ro,drop,drot,drox,drom,ro_t,dro_tm PRINT*,'gam,dgamx,dgamlpp' ; WRITE(*,2000)gam,dgamx,dgamlpp PRINT*,'dgampt,dgamt,dgamr,dgaml,dgamm,dgamp,dgamlpt,dgamlp' WRITE(*,2000)dgampt,dgamt,dgamr,dgaml,dgamm,dgamp,dgamlpt,dgamlp 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) psist0=psist ; teta0=teta ; grad0=grad ; eps0=epsilon(1) ro0=ro ; u0=u ; gam0=gam ; x00=xchim(1) ; dy710=dy71 ; dy720=dy72 p_t0=p_t ; t_t0=t_t ; ro_t0=ro_t ; u_t0=u_t ; tds0=tdst 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=ABS(yd(3,0)) ray=SQRT(ay3) !rayon ay4=ABS(yd(4,0)) l13=SQRT(ay4) !l**1/3 ln=l13**3 !l/lsol ay5=ABS(yd(5,0)) m13=SQRT(ay5) !m**1/3 mn=m13**3 !masse/mtot 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(ay5,mc(n_ch)),l,xchim,dxchim) IF(no_croiss)PRINT*,'Pb. at 6 in static_m23' WHERE(xchim < 1.d-38)xchim=1.d-38 IF(pturb .AND. der)THEN !der=t on tient compte de dln Pgaz/dln Ptot 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,mn,ln,ray,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) 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 ay5,l,frot,dfrot) IF(no_croiss)PRINT*,'Pb. at 7 in static_m23' w=frot(1) CASE(6) CALL omega_jpz(ray,w) END SELECT v=SQRT(ay5/ay3)/prn*w**2 !rotation teta=fac(ip)*((ay5/ay3)**2/prn*(cte21+cte22*grad)+cte25*v 1 +cte23/ro*m13/ray+ctem) psist=yd(6,0)/teta !psi/teta IF(dt /= 0.d0)THEN c PRINT*,dt IF(compt == 0)THEN !interpolation du TdS au temps t 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_m23' tdst=ftds(1) ; dtdst_m=dftds(1) c PRINT*,'tdst,mn,l',tdst,dtdst_m,mn,l ELSE !TdS au temps t+dt IF(l_pertm)THEN CALL bsp1dn(1,old_ptm,x_ptm,xt_ptm,n_ptm,m_ptm,knot_ptm,.TRUE., 1 MIN(ay5,x_ptm(n_ptm)),l,fm,dfm) IF(no_croiss)PRINT*,'Pb. at 9 in static_m23' 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)) c le TdS IF(kipp)THEN !approximation de Kippenhahan c Approximation de Kippenhahn forme 2 TdS = cp T ( dlnT - grad_ad dlnP ) / dt 4-27 dlnp=yd(1,0)-fqs(1) ; dlnt=yd(2,0)-fqs(2) tdst=cp*trn*(dlnt-gradad*dlnp)/dt ELSE !TdS=dU+PdV 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_m23' WHERE(xchim_t < 1.d-38)xchim_t=1.d-38 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, 4 alfa,beta,gamma1) tdst=(u-u_t-pgn/ro**2*(ro-ro_t))/dt !+Tds ENDIF !kipp c variation de énergie cinétique de rotation (!!!!!incomplet) tdst=tdst-eps_rot/dt ENDIF !compt=0 ENDIF !dt/=0 dy1=cte11*(ay5/ay3)**2/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*m13/ray bes(3)=yd(3,1)-dy3*psist !dl/dq=epsilon m13 psi/teta bes(4)=yd(4,1)-cte14*epsilon(1)*m13/l13*psist IF(dt /= 0.d0)THEN dy5=cte14t*m13/l13 tdst=dy5*tdst bes(4)=bes(4)+tdst*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 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 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*,'tds,p_t,t_t,X,gam,dy71,dy72' WRITE(*,2000)(tdst-tds0)/dstor, 1 (p_t-p_t0)/dstor,(t_t-t_t0)/dstor,(xchim(1)-x00)/dstor, 2 (gam-gam0)/dstor,(dy71-dy710)/dstor,(dy72-dy720)/dstor ENDDO !jd PRINT*,'cx=',cx,ip,fac(ip) PAUSE'test dérivées' ENDDO !id DEALLOCATE(bes,dern,yd) ; deriv=.FALSE. RETURN CASE(2) !les limites 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 deriv=.TRUE. IF(deriv)THEN ALLOCATE(yd(ne,0:1)) ; yd=y ; unpdd=1.d0-dd ENDIF SELECT CASE(li) c condition sur r au centre CASE(1) !au centre en r be(1)=y(3,0) !en q=1 r**2=0 ae(1,3,0)=1.d0 !dérivée/r**2 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**2/3=0 ae(1,4,0)=1.d0 !dérivée/l**2/3 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**2/3=0 ae(1,5,0)=1.d0 !dérivée/m**2/3 c PRINT*,'limite au centre',li c WRITE(*,2000)(y(i),i=1,ne),be(1) c PAUSE'au centre' c condition sur Ptot au raccord CASE(4) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MIN(y(5,0),mc(n_ch)),l,xchim0,dxchim) IF(no_croiss)PRINT*,'Pb. at 11 in static_m23' WHERE(xchim0 < 1.d-38)xchim0=1.d-38 CALL chim_gram(xchim0,dxchim) ray=SQRT(y(3,0)) ; l13=SQRT(y(4,0)) ; ln=l13**3 CALL atm(.FALSE.,ln,ray,xchim0,ptext0,dptdl0,dptdr0,text0,dtdl0,dtdr0, 1 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/2.d0/ray !dérivée/r**2 ae(1,4,0)=-dptdl0/ptext0*3.d0/2.d0*l13 !dérivée/l**2/3 IF(deriv)THEN !test dérivée PRINT*,'limite pext',li ; WRITE(*,2000)y(1,0),LOG(ptext0) c WRITE(*,2000)y(1:ne,0) c PRINT*,'exterieur lnp,ln ptext,be(1),ray,xchim0' c WRITE(*,2000)y(1,0),LOG(ptext),be(1),ray,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 ray=SQRT(yd(3,0)) ; ln=SQRT(yd(4,0))**3 CALL atm(.FALSE.,ln,ray,xchim0,ptext,dptdl,dptdr,text,dtdl,dtdr, 1 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 en li=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/2.d0/ray !dérivée/r**2 ae(1,4,0)=-dtdl0/text0*3.d0/2.d0*l13 !dérivée/l**2/3 IF(deriv)THEN !test de dérivation PRINT*,'limite Text',li ; 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 ray=SQRT(yd(3,0)) ; ln=SQRT(yd(4,0))**3 CALL atm(.FALSE.,ln,ray,xchim0,ptext,dptdl,dptdr,text,dtdl,dtdr, 1 mext,dmdl,dmdr,pext,dpdl,dpdr,teff) PRINT*,'dérivée pour li=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**(2.d0/3.d0) !en q=n_qs, m=mext(r,l) ae(1,3,0)=-dmdr0/3.d0/ray/m13 !dérivée/r**2 ae(1,4,0)=-dmdl0*l13/m13 !dérivée/l**2/3 ae(1,5,0)=1.d0 !dérivée/m**2/3 c WRITE(*,2000)y(5,0) IF(deriv)THEN !test de dérivation PRINT*,'limite mext',li WRITE(*,2000)y(5,0),mext0**(2.d0/3.d0) WRITE(*,2000)y(1:ne,0) PRINT*,'limite y(5),mext',li PRINT*,'y(5,0),mext,be(1)' WRITE(*,2000)y(5,0),mext0,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 ray=SQRT(yd(3,0)) ; ln=SQRT(yd(4,0))**3 CALL atm(.FALSE.,ln,ray,xchim0,ptext,dptdl,dptdr,text,dtdl,dtdr, 1 mext,dmdl,dmdr,pext,dpdl,dpdr,teff) PRINT*,'dérivée pour li=6 / ',variable(i) WRITE(*,2000)(yd(5,0)-mext**(2.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/2.d0/ray !dérivée/r**2 ae(1,4,0)=-dpdl0/pext0*3.d0/2.d0*l13 !dérivée/l**2/3 ae(1,Ipg,0)=1.d0 !dérivée/ln Pgaz IF(deriv)THEN !test de dérivation PRINT*,'limite pext, uint',li WRITE(*,2000)y(Ipg,0),pext0 ; WRITE(*,2000)y(1:ne,0) PRINT*,'limite y(Ipg)',li ; PRINT*,'y(Ipg,0),pext,be(1)' WRITE(*,2000)y(Ipg,0),pext0,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 ray=SQRT(yd(3,0)) ; ln=SQRT(yd(4,0))**3 CALL atm(.FALSE.,ln,ray,xchim,ptext,dptdl,dptdr,text,dtdl,dtdr, 1 mext,dmdl,dmdr,pext,dpdl,dpdr,teff) PRINT*,'dérivée pour li=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_m23 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_m23 erreur sur fait = 1, ou 2, fait=',fait PRINT*,'ARRET' ; STOP END SELECT !fait reprend=.FALSE. RETURN END SUBROUTINE static_m23