c******************************************************************** SUBROUTINE resout(un23,dt,dts) c routine public du module mod_static c gestion de la résolution des équations avec fonction de répartition c Auteur: P.Morel, Département J.D. Cassini, O.C.A., CESAM2k c on distingue par SELECT CASE(un23) les 3 cas: c un23=1 : poursuite d'une évolution c un23=2 : calcul d'un modèle initial de ZAMS homogène c un23=3 : calcul d'un modèle initial de PMS homogène c Pour le problème aux limites: c résolution du système d'équations différentielles non linéaires c de la structure interne par développement sur B-splines c par itération Newton avec dérivées analytiques c (ou numériques pour tests) c les équations de la SI sont écrites pour tout point dans c static_m23/m13, le jacobien de la méthode NR est formé et résolu c dans coll_qs dans l'espace des splines c pour le modèle quasi-statique, au temps àge+dt c paramètres indépendants du temps c m_qs : ordre des splines pour l'intégration par collocation c ord_qs=m_qs+r_qs=m_qs+1 : ordre de la représentation spline c ne : nombre d'inconnues (6 ou 7) c paramètres dépendants du temps c nombre de couches au temps t : n_qs c knot : nombre de points de table, c dim_qs=knot-ord_qs : dimension de l'espace des splines c bp(ne,dim_qs) : coefficients des splines c q(n_qs), qt(knot) : abscisses, et vecteur nodal c pour la composition chimique, au temps àge+dt c paramètres indépendants du temps c m_ch : ordre de la représentation spline c nchim : nombre d'éléments chimiques c paramètres dépendants du temps c nombre de couches au temps t : n_ch c knotc : nombre de points de table (nodaux), c dim_ch=knotc-m_ch : dimension de l'espace des splines c chim(nchim,dim_ch) : coefficients des splines c q(n_ch), qt(knotc) : abscisses, et vecteur nodal c les quantités au temps t-dt ont l'extension _t c Exemple : knot_t est nombre de points de table de la c représentation spline des variables du modèle quasi-statique c au temps t=age-dt c Pour le problème de valeurs initiales appel à evol c en_m23 = .TRUE. variables m_stat=m^2/3, r_stat=r^2 c en_m23 = .FALSE. variables m_stat=m^1/3, r_stat=r c pour avoir un algorithme unique pour la diffusion, la composition c chimique est toujours tabulée en fonction de (m/Msol)^2/3 c l'énergie graviphique, TdS/dt=tds est tabulée en fonction c de m_stat=m^2/3 ou de m_stat=m^1/3 c entrées c un23 : 1 poursuite d'une évolution c 2 calcul d'un modèle initial de ZAMS c 3 calcul d'un modèle initial de PMS c dt : pas temporel c sortie : c dts : estimation du pas temporel suivant c variables globales de module utilisées c àge : àge du modèle c mc,mct,nc,knotc,chim: comp. chim. à t+dt c n_qs,bp,q,qt,knot,p,t,r,l,m : var. pples. à t+dt c lim : nombre de limites ZR/ZC c jlim : abscisses des limites ZR/ZC c r_zc, r_ov : rayons de limites ZR/ZC et d'overshoot au temps t+dt c m_zc : masses des limites des ZR/ZC au temps t+dt c lconv =.TRUE. : début de ZC au temps t+dt c mstar: masse avec perte de masse c x_tds,xt_tds,m_tds,n_tds,knot_tds,tds: éléments pour interpoler c le TdS au temps t c m_zct : masses des limites des ZR/ZC au temps t c lconv_t =.TRUE. : debut de ZC au temps t c le réseau de points de raccord q, posséde n_qs points c------------------------------------------------------------------------- USE mod_donnees, ONLY: ajuste, agemax, ah, aradia, baratine, clight, 1 dpsi, dpsim, dpsip, dtmin, dx_tams, en_m23, fcv, g, he_ajuste, he_core, 2 hhe_core, ini0, Kdes_stat, Krot, langue, lnt_stop, loc_zc, lsol, l_fac, 3 l_pertm, l_pertmt, l_tr, msol, m_ch, m_ptm, m_qs, m_tds, nb_max_modeles, 4 nchim, nc_max, ne, nom_fich2, nom_output, nrot, nucleo, n_max, n_min, 5 ord_qs, ord_rot, ovshti, ovshts, pi, pnzc, precision, precix, 6 psi0, pturb, rsol, r_ajuste, r_qs, r_stop, secon6, sigma, t_ajuste, 7 t_stop, w_rot, x_ajuste, x_stop, x_tams USE mod_evol, ONLY: evol, nzc USE mod_etat, ONLY: calc_entropie, dege_elec, etat USE mod_kind USE mod_nuc, ONLY: nuc USE mod_numerique, ONLY: bsp1dn, bvald, entre, noein, no_croiss, sum_n USE mod_opa, ONLY: dehors, opa USE mod_variables, ONLY: age, bp, bp_t, chim, chim_gram, 1 chim_t, inter, jlim, knot, knotc, knotc_t, knotr, knot_t, 2 knot_ptm, knot_tds, lconv, lhe_stop, lim, lr_stop, 3 lt_stop, lx_stop, mc, mc_t, mct, mct_t, MI_atm, model_num, mrot, 4 mrott, mstar, mw_tot, m_zc, M_iner, m_stat, m_stat_t, 5 nb_modeles, n_ch, n_ch_t, n_ptm, n_qs, n_qs_t, 6 n_rot, n_tds, n_tds_t, old_ptm, q, qt, q_t, qt_t, 7 rota, rstar, rstar_t, r_ov, r_zc, r_stat, r_stat_t, sortie, 8 S_entro, S_0atm, tds, tot_conv, tot_rad, t_conv, x_ptm, 9 xt_ptm, x_tds, xt_tds_t, xt_tds, wrot IMPLICIT NONE INTEGER, INTENT(in) :: un23 REAL (kind=dp), INTENT(inout) :: dt, dts REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: vtr, jac REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: r_lim REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: dfrot, frot, l, m, 1 p, pt, r, t, vtm, vtmt REAL (kind=dp), DIMENSION(nchim) :: dxchim, depsx, xchim, xchimg, xchim_t REAL (kind=dp), DIMENSION(ne) :: dfdq, f, f_t REAL (kind=dp), DIMENSION(5) :: epsilon REAL (kind=dp), PARAMETER :: drstar=5.d-3 REAL (kind=dp), SAVE :: cte1, mcin_totp, mcin_tot0, preci01, rstarp, 1 x_stopm, x_tamsm REAL (kind=dp) :: alfa, beta, be7, b8, corr, cp, 1 dcpp, dcpt, dcpx, delta, deltap, deltat, deltax, depst, depsro, 2 dgradadp, dgradadt, dgradadx, dkapro, dkapt, dkapx, 3 drop, drot, drox, dup, dut, dux, d1, d2, err, f17, gamma1, grad, 4 gradad, gradrad, kap, hh, mcin_tot, nu, n13, o15, pp, ro, tp, u INTEGER, PARAMETER :: iter_max0=100 INTEGER, SAVE :: ini1, ini2, iter_max, lq=1 INTEGER :: compt, i, idep, j, jinf, jsup, knotw, new_nqs, nqsp, ntest, 1 qmax, vare LOGICAL, SAVE :: cmax=.FALSE., init=.TRUE., ini_mcin=.TRUE.,ovsht, pass1=.TRUE. LOGICAL :: reprend, logic, logic1, l_rstar, modif, no_last, new_rep c------------------------------------------------------------------------- 2000 FORMAT(8es10.3) c PAUSE'entrée resout' c -----------------Début des initialisations------------- Iini: IF(init)THEN init=.FALSE. c constantes cte1=3.d0*lsol/16.d0/pi/aradia/clight/g/msol c baratine=.FALSE. permet de dérouter sur le fichier mon_modele_101 les c informations concernant le déroulement des calculs pour l'équilibre c quasi-statique IF(baratine)THEN usl_static=6 ELSE usl_static=101 OPEN(unit=101,form='formatted',status='unknown',!access='append', 1 file=TRIM(nom_fich2)//'_static') ENDIF c écritures SELECT CASE(langue) CASE('english') WRITE(usl_static,1001) ; WRITE(2,1001) 1001 FORMAT(/,'---- Parameters for the quasi-static model ------',/) CASE DEFAULT WRITE(usl_static,1) ; WRITE(2,1) 1 FORMAT('----- paramètres pour le modèle quasi-statique ----',/) END SELECT c initialisations diverses c WRITE(*,2000)x_stopm,x_stop,ah !arret des que X(centre) <= x_stop c PAUSE'ah' c on ne dépassera pas iter_max IF(l_fac)THEN iter_max=ini0+15 ELSE iter_max=ini0+8 ENDIF c au delà de ini1, si les lim ZR/ZC sont bien positionnées pas d'appel à lim_zc, evol ini1=ini0+5 c on impose au moins ini2 itérations en cas de convergence forte ini2=3 c il y aura ajustement du nombre de couches et une nouvelle répartition c si la cte. de répartition a varié de plus de dpsi SELECT CASE(langue) CASE('english') WRITE(usl_static,1214)psi0,NINT(dpsim*100.d0),NINT(dpsip*100.d0) WRITE(2,1214)psi0,NINT(dpsim*100.d0),NINT(dpsip*100.d0) 1214 FORMAT('Change of the repartition function per shell = ',es10.3, 1 /,'The number of shells changes if the repartition constant',/, 2 'varies by more than [',i2'%,+',i2,'%]',/) CASE DEFAULT WRITE(usl_static,214)psi0,NINT(dpsim*100.d0),NINT(dpsip*100.d0) WRITE(2,214)psi0,NINT(dpsim*100.d0),NINT(dpsip*100.d0) 214 FORMAT('saut de la fonction de répartition par couche =',es10.3, 1 /,'modification du nombre de couches si la constante de',/, 2 'répartition varie de plus de [',i2'%,+',i2,'%]',/) END SELECT c vérification c DO i=1,n_qs c CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) c WRITE(*,2000)q(i),f(1:ne) c ENDDO c PAUSE'entrée resout' c si ovsht, il y a des overshoot ovsht=MAX(ABS(ovshts),ABS(ovshti)) > 0.d0 c le pas temporel dt est estimé dans update, à partir de sa c valeur précédente et d'une estimation dts issue de evol c basée sur la variation des abondances au pas temporel précédent c au premier appel, dans le cas d'une reprise d'évolution (un23=1) c dts est inconnu, on l'initialise à dt c dans les autres cas (un23=1,3) dt=0 et dts=0 sont inutilisés c x_stop par mole x_stopm=x_stop/ah c x_tamsm x_tams par mole x_tamsm=x_tams/ah c rstarp: valeur de rstar à l'itération précédente de NR globale rstarp=rstar c rayons des limites des ZC (Tristan) IF(l_tr)ALLOCATE(M_iner(pnzc),r_lim(2),S_entro(pnzc), 1 S_0atm(0:1),t_conv(pnzc+1)) c limitation de la variation du nombre de couches r_add=75 c si err < preci01 convergence forcée si d'autres critères ne sont pas atisfaits preci01=precix/100.d0 ENDIF Iini !init c ---------------Fin des initialisations---------- calcul des modèles c on différencie les cas suivant un23: c un23 = 1 : poursuite d'une évolution c un23 = 2 : modèle quasi-statique de ZAMS c un23 = 3 : modèle quasi-statique de PMS c PRINT*,'un23=',un23,', n_qs=',n_qs ; WRITE(*,2000)dt,dts ; PAUSE'resout entrée' SELECT CASE(un23) c évolution d'un modèle CASE(1) SELECT CASE(langue) CASE('english') WRITE(usl_static,1002) 1002 FORMAT('**** Integration of the quasi-static model (begin) ****',/) CASE DEFAULT WRITE(usl_static,2) 2 FORMAT('**** Intégration du modèle quasi-statique (début) ****',/) END SELECT c détermination du nombre de couches suivant la précision requise c pour le dernier modèle ie. agemax-age < dt des modèles en "sa" c (solar accuracy), "co' (COROT) c on augmente le nombre de couches au maximum NMAX en fixant psi0 petit c on augmente aussi la précision de l'intégration globale et c celle de la localisation des limites ZR/ZC CALL fcmax(cmax) c PRINT*,cmax,n_qs; WRITE(*,2000)age,dt ; PAUSE'resout, cmax' c on poursuit une évolution modèle àge --> modèle àge+dt c initialisation du nouveau dt CALL update(.TRUE.,dt,dts) SELECT CASE(langue) CASE('english') WRITE(usl_static,1020)model_num+1,dt ; WRITE(2,1020)model_num+1,dt 1020 FORMAT('The computations start for the next model #',i5, 1 ', estimated time step, dt=',es10.3) CASE DEFAULT WRITE(usl_static,20)model_num+1,dt ; WRITE(2,20)model_num+1,dt 20 FORMAT('Début du calcul du modèle n°',i5,', pas temporel estimé,dt=', 1 es10.3) END SELECT c Première boucle infinie Bnr1 gestion de la résolution du système c d'équations dans sa globalité c si compt=0 : utilisation du TdS/dt du pas précédent c PAUSE'resout avant boucle' compt=0 ; err=1.d20 ; reprend=.FALSE. c au début du calcul d'un pas temporel, dans lim_zc, la répartition sera c actualisée si le nombre de couches change new_nqs=actu_nqs(cmax) c PRINT*,'resout ',new_nqs c no_last=.TRUE. : ce n'est pas le dernier modèle à calculer no_last=.TRUE. c Séquence globale infinie Bnr1: DO c initialisation de l_dmax pour imposer l'appel à lim_zc IF(compt == 0)l_dmax=.FALSE. c Seconde séquence infinie Bnr10: chaque pas constitue une itération c globale: détermination des limites ZR/ZC, évolution temporelle de la c composition chimique et de la rotation, résolution des équations c d'équilibre quasi-statique. En cas d'échec de l'une des étapes c (Ex. non convergence) la séquence est reprise. c reprend=.FALSE. !initialisation Bnr10: DO c lrstar= variation relative de R* d'une iter. à l'autre > drstar(0.005) c on poursuit les itérations si Rext varie trop l_rstar=ABS(1.d0-rstar/rstarp) > drstar .OR. r_zc(lim) > rstar c valeur de Rstar de l'itération précédente rstarp=rstar c tant que compt <= ini0 ou que la distance des limites ZR/ZC d'une couche est > 5%, c il y a détermination de la c perte de masse, des limites ZR/ZC et évolution de la comp. chimique logic1=(compt < ini0) .OR. (compt < ini1 .OR. l_rstar .OR. .NOT.l_dmax) IR1: IF(logic1)THEN c PRINT*,'resout logic1', logic1 c détermination la new masse mstar en tenant compte de la perte de masse, c avec perte_tot détermination des diminutions de masse c dues à E=mc^2 d'où les abscisses (en m_stat ou m^1/3) old_ptm c qui, à l'àge t, correspondaient aux abscisses x_ptm à l'àge t+dt IF(l_pertm .OR. l_pertmt)CALL pertm(dt) c s'il y a reprise, réactualisation de la répartition new_rep= compt == 0 .OR. reprend c détermination des limites ZR/ZC et des facteurs de répartition c PRINT*,new_rep,cmax,new_nqs ; PAUSE'avant lim_zc' CALL lim_zc(new_rep,new_nqs) IF(dehors)THEN !il y a sortie de table d'opacité reprend=.TRUE. c évolution temporelle de la composition chimique et de la rotation ELSE CALL evol(compt,dt,dts,reprend) ENDIF ELSE IR1 WRITE(*,120)compt,l_rstar,l_dmax,logic1 120 FORMAT('Sans actu. des lim. ZR/ZC et de comp. chim., compt:', 1 i2,', l_rstar:',l1,', l_dmax:',l1,', logic1:',l1) ENDIF IR1 c il n'y a pas eu de pb. dans evol (CV dans rk_imps, diff_rota, diff_chim) c résolution des équations de la structure interne, le nombre d'itérations c maximal est d'autant plus grand que corr est voisin de 1 IF(.NOT.reprend)CALL coll_qs(dt,compt,reprend,err,vare,qmax,corr) c il peut y avoir eu un problème dans coll_qs (cv ou TdS varie trop) c new_rep=.FALSE.: pas de redistribution de la solution IR38: IF(reprend)THEN dt=dt*0.75d0 IR39: IF(dt < dtmin)THEN WRITE(usl_static,5)dt,dtmin ; WRITE(2,5)dt,dtmin CALL sortie('resout 1') 5 FORMAT('dt =',es10.3,' < dt min.=',es10.3,', abandon') ELSE IR39 !reprise avec pas moitié WRITE(usl_static,6)dt,dtmin ; WRITE(2,6)dt,dtmin 6 FORMAT('réduction du pas temporel dt =',es10.3, 1 ', dtmin =',es10.3,/) new_nqs=n_qs ; CALL update(.FALSE.,dt,dts) IF(pass1)new_nqs=n_qs compt=0 ; err=1.d20 ; CYCLE Bnr10 ENDIF IR39 ELSE IR38 EXIT Bnr10 ENDIF IR38 ENDDO Bnr10 c décompte du nombre d'itérations compt=compt+1 c écriture des variables quasi-statiques au max d'erreur, indice qmax WRITE(usl_static,*) WRITE(usl_static,100)compt,err,nom_qsr(vare),qmax,1.d0/corr 100 FORMAT('iter. globale:',i3,', err. max:',es8.1, 1 ', var:',a,', couche: ',i5,', corr:',es8.1) CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(qmax),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 4 dans resout' pp=EXP(f(1)) ; tp=EXP(f(2)) IF(en_m23)THEN IF(pturb)THEN WRITE(usl_static,101)pp,tp,SQRT(ABS(f(3))), 1 (SQRT(ABS(f(i)))**3,i=4,5),EXP(f(7)) 101 FORMAT('variables Ptot, T, R, L, M, Pgaz :',6es10.3) ELSE WRITE(usl_static,103)pp,tp,SQRT(ABS(f(3))), 1 (SQRT(ABS(f(i)))**3,i=4,5) 103 FORMAT('variables Pgaz, T, R, L, M :',5es10.3) ENDIF ELSE !en m^1/3 IF(pturb)THEN WRITE(usl_static,101)pp,tp,(f(i),i=3,4),f(5)**3,EXP(f(7)) ELSE WRITE(usl_static,103)pp,tp,(f(i),i=3,4),f(5)**3 ENDIF f(5)=f(5)**2 !m^2/3 ENDIF c WRITE(*,2000)bp(6,1) ; PAUSE'cv' c PAUSE'après écritures (EVOLUTION)' c écriture de la comp. chim au max d'erreur CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(f(5),mc(n_ch))),lq,xchim,dxchim) IF(no_croiss)PRINT*,'Pb. en 5 dans resout' xchimg=xchim*nucleo WRITE(usl_static,102)(xchimg(i),i=1,MIN(nchim,6)) 102 FORMAT('abon. -iers elem. :',6es10.3) c EOS au maximum d'erreur CALL etat(pp,tp,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 opacité au maximum d'erreur CALL opa(xchimg,tp,ro,kap,dkapt,dkapro,dkapx) c gradients au maximum d'erreur, grad_rad=0 si err. max. au centre IF(qmax == 1)THEN gradrad=0.d0 ELSE gradrad=cte1*kap*pp/tp**4 ; grad=dfdq(2)/dfdq(1) IF(en_m23)THEN gradrad=gradrad*SQRT(f(4)/f(5))**3 ELSE gradrad=gradrad*f(4)/SQRT(f(5))**3 !f(5) est m^2/3 ENDIF ENDIF grad=dfdq(2)/dfdq(1) c écritures au maximum d'erreur WRITE(usl_static, 104)dege_elec(xchim,tp,ro),ro,gradrad,grad 104 FORMAT('dégéné. élec.:',es10.3,', densité:',es10.3, 1 ', grad_rad=',es10.3,', grad=',es10.3) c kdes_stat=1: écriture/plot des var. quasi-stat. (compt a déjà été incrémenté) IF(Kdes_stat == 1)CALL ecrit_static(compt-1,dt) c pour écrire (debug) la solution intermédiaire IR37: IF(.FALSE.)THEN c IR37: IF(.TRUE.)THEN c IR37: IF(err > 1.d3)THEN PAUSE'écriture' PRINT*,'p(i),t(i),r(i),l(i),m(i)' ALLOCATE(l(n_qs),m(n_qs),p(n_qs),pt(n_qs),r(n_qs),t(n_qs)) Bt: DO i=1,n_qs c IF(i < n_qs-2)CYCLE Bt !pour limiter les écritures c IF(i > 3)CYCLE Bt !pour limiter les écritures IF(.NOT. entre(qmax-10,qmax+10,i))CYCLE Bt CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 6 dans resout' c variables au temps àge, m_stat=m^2/3, r_stat=r^2 ou m_stat=m^1/3, r_stat=r CALL inter('m',bp_t,q_t,qt_t,n_qs_t,knot_t,f(5),f_t, 1 dfdq,r_stat_t,m_stat_t) f_t(1)=EXP(f_t(1)) !variable ln Ptot au temps t pt(i)=EXP(f(1)) !variable ln Ptot f_t(2)=EXP(f_t(2)) !variable ln Tau temps t t(i)=EXP(f(2)) !variable ln T IF(pturb)THEN !avec pression turbulente 8 inconnues p(i)=EXP(f(7)) !variable ln Pgaz ELSE !sans pression turbulente 7 inconnues p(i)=pt(i) ENDIF IF(en_m23)THEN r(i)=SQRT(ABS(f(3))) !rayon/Rsol f_t(3)=SQRT(ABS(f_t(3))) !rayon/Rsol au temps t l(i)=SQRT(ABS(f(4)))**3 !l/Lsol f_t(4)=SQRT(ABS(f_t(4))) !l/Lsol au temps t m(i)=SQRT(ABS(f(5)))**3 !m/Msol ELSE r(i)=ABS(f(3)) !rayon/Rsol l(i)=f(4) !l/Lsol m(i)=f(5)**3 !m/Msol f_t(5)=f_t(5)**3 !m/Msol au temps t f(5)=f(5)**2 !m^2/3 ENDIF c IF(t(i) > 8.d4)CYCLE Bt !pour limiter les écritures PRINT*,i c WRITE(*,2000)pt(i),p(i),t(i),r(i),l(i),m(i),f(6) c WRITE(*,2000)p(i),t(i),r(i),l(i),m(i),f(6) WRITE(usl_static,2000)p(i),f_t(1),t(i),f_t(2),r(i),l(i), 1 f_t(4),m(i) c WRITE(*,2000)f(1:ne) c CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, c 1 knotc,.TRUE.,MIN(f(5),mc(n_ch)),lq,xchim,dxchim) c IF(no_croiss)PRINT*,'Pb. en 7 dans resout' c WRITE(*,2000)xchim(1:8) ; WRITE(*,2000)xchim(9:nchim) ENDDO Bt DEALLOCATE(l,m,p,pt,r,t) PAUSE'solution intermédiaire' ENDIF IR37 !écriture de la solution intermédiaire c contrôle des var. de P, T IF(compt >= 2 .AND. .NOT.(l_pertm .OR. l_pertmt))THEN CALL control_dPT(reprend) IF(reprend)THEN !variation trop importante de L, P, R ou T dt=dt*0.75d0 IF(dt < dtmin)THEN WRITE(usl_static,5)dt,dtmin ; WRITE(2,5)dt,dtmin CALL sortie('resout 2') ELSE !reprise pas moitié WRITE(usl_static,6)dt,dtmin ; WRITE(2,6)dt,dtmin new_nqs=n_qs ; CALL update(.FALSE.,dt,dts) IF(pass1)new_nqs=n_qs compt=0 ; err=1.d20 ; CYCLE Bnr1 c WRITE(usl_static,10)dt,dtmin ; WRITE(2,6)dt,dtmin 10 FORMAT('poursuite avec un pas temporel réduit de 75%, dt =',es10.3, 1 ', dtmin =',es10.3,/) ENDIF ENDIF ENDIF c la réactualisation de m_stat et r_stat qui accélérent la recherche c dans le ssp. inter est effectuée dans coll_qs' c gestion des itérations IF(l_fac)THEN logic=(err < precix .AND. (l_dmax .OR. compt >= iter_max-1)) 1 .OR. (err < preci01 .AND. compt >= ini0) c PRINT*,logic,l_dmax,compt logic=logic .OR. (err < preci01 .AND. compt >= ini2) logic=logic .AND. compt >= ini0 ELSE logic=err < precix .AND. compt >= ini0 ENDIF c PRINT*,'resout logic',logic,compt,ini0,l_fac,l_dmax,err,precix IR7: IF(logic)THEN c le codage des ajustements sur t_stop, r_stop sont introduits par un INCLUDE INCLUDE 'ajustements.f' c tests de convergence, convergence forcée IR13: IF(err > precix)THEN dts=dts/2.d0 SELECT CASE(langue) CASE('english') WRITE(usl_static,1004)model_num,dts WRITE(2,1004)model_num+1,dts 1004 FORMAT(/,'----------WARNING-----------',/, 1 'Forced convergency for the model #',i5, 2 ', estimated next time step reduced, dts=',es10.3,/, 3 '-----------WARNING------------') CASE DEFAULT WRITE(usl_static,4)model_num,dts ; WRITE(2,4)model_num+1,dts 4 FORMAT(/,'----------ATTENTION-----------',/, 1 'Convergence forcée pour le modèle #',i5, 2 ', pas temporel estimé diminué, dts=',es10.3,/, 3 '----------ATTENTION-----------') END SELECT ELSE IR13 SELECT CASE(langue) CASE('english') WRITE(usl_static,1015)model_num+1 1015 FORMAT(/,'Convergency for the model #',i5) CASE DEFAULT WRITE(usl_static,15)model_num+1 15 FORMAT(/,'Convergence du modèle #',i5) END SELECT ENDIF IR13 c on se satisfait de la convergence atteinte, tabulation du TdS CALL tabul_TdS(dt,modif) c problème dans tabul TdS ou sortie IR16: IF(modif)THEN dt=dt*0.75d0 IF(dt < dtmin)THEN WRITE(usl_static,5)dt,dtmin ; WRITE(2,5)dt,dtmin CALL sortie('resout 3') ELSE WRITE(usl_static,6)dt,dtmin ; WRITE(2,6)dt,dtmin new_nqs=n_qs ; CALL update(.FALSE.,dt,dts) IF(pass1)new_nqs=n_qs compt=0 ; err=1.d20 ; new_rep=.FALSE. ; CYCLE Bnr1 ENDIF ENDIF IR16 c on sort de Bnr1 EXIT Bnr1 ELSE IR7 c la convergence n'est pas atteinte on ne peut dépasser iter_max itérations IR5: IF(compt >= iter_max)THEN WRITE(usl_static,12) ; WRITE(2,12) ; modif=.TRUE. !PAUSE'ici1' 12 FORMAT(/,'********',/,'pas de convergence dans resout',/, 1 'diminution du pas temporel, réinitialisation',/,'********',/) c poursuite des itérations ELSE IR5 modif=.FALSE. ENDIF IR5 c reprise avec pas temporel moitié ou sortie IR6: IF(modif)THEN dt=dt*0.75d0 IF(dt < dtmin)THEN WRITE(usl_static,5)dt,dtmin ; WRITE(2,5)dt,dtmin CALL sortie('resout 4') ELSE WRITE(usl_static,6)dt,dtmin ; WRITE(2,6)dt,dtmin new_nqs=n_qs ; CALL update(.FALSE.,dt,dts) IF(pass1)new_nqs=n_qs compt=0 ; err=1.d20 ; new_rep=.TRUE. ; CYCLE Bnr1 ENDIF ENDIF IR6 ENDIF IR7 !compt >= 2 ENDDO Bnr1 c PAUSE'sortie de Bnr1' c le modèle quasi-statique a convergé, détermination de Teff, c de la situation convective, du TdS/dt, c du moment angulaire, puis retour vers cesam SELECT CASE(langue) CASE('english') WRITE(2,1058)n_qs,n_ch 1058 FORMAT(/,'Shell number for structure variables : ',i5,/, 1 'For chemicals : ',i5) CASE DEFAULT WRITE(2,58)n_qs,n_ch 58 FORMAT(/,'Nombre de couches pour les variables de structure : ' 1 ,i4,/,'Pour la composition chimique : ',i4) END SELECT c écritures du dt utilisé SELECT CASE(langue) CASE('english') WRITE(usl_static,1070)dt ; WRITE(2,1070)dt 1070 FORMAT(/,'Time step used, dt=',es10.3) CASE DEFAULT WRITE(usl_static,70)dt ; WRITE(2,70)dt 70 FORMAT(/,'Pas temporel utilisé, dt=',es10.3) END SELECT c modèle totalement convectif: lim=1,jlim(1)=n,lconv(1)=.FALSE. c modèle totalement radiatif: lim=0,jlim(i)=-100,lconv(i)=.FALSE. IR8: IF(lim > 0)THEN IR9: IF(lim == 1 .AND. jlim(1) == 1 .AND. .NOT.lconv(1))THEN WRITE(2,7) 7 FORMAT('modèle complétement convectif',/) ELSE IR9 WRITE(2,*) ; WRITE(2,9) 9 FORMAT('limites Zones radiatives/convectives en Msol et Rstar') DO j=1,lim IR10: IF(ovsht)THEN IR11: IF(lconv(j))THEN !debut de ZC IF(ovshts > 0.d0)THEN !overshoot supérieur pas d'extension WRITE(2,1221)jlim(j),m_zc(j),j,r_zc(j)/rstar 1221 FORMAT('couche:',i5,', masse:',es10.3, 1 ', rayon du début de la ZC n°',i2,'=',es10.3) ELSE !overshoot inférieur extension WRITE(2,122)jlim(j),m_zc(j),r_ov(j)/rstar 122 FORMAT('couche:',i5,', masse:',es10.3, 1 ', rayon du début de l''overshoot inférieur=',es10.3) ENDIF ELSE IR11 !fin de ZC IF(ovshts > 0.d0)THEN !overshoot supérieur extension WRITE(2,222)jlim(j),m_zc(j),r_ov(j)/rstar 222 FORMAT('couche:',i5,', masse:',es10.3, 1 ', rayon de la fin de l''overshoot supérieur=',es10.3) ELSE !overshoot inférieur pas d'extension WRITE(2,2221)jlim(j),m_zc(j),j,r_zc(j)/rstar 2221 FORMAT('couche:',i5,', masse:',es10.3, 1 ', rayon de la fin de la ZC n°',i2,'=',es10.3) ENDIF ENDIF IR11 ELSE IR10 !sans overshoot IF(lconv(j))THEN WRITE(2,1221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ELSE WRITE(2,2221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ENDIF ENDIF IR10 ENDDO WRITE(2,*) ENDIF IR9 ELSE IR8 WRITE(2,8) 8 FORMAT('modèle complétement radiatif',/) ENDIF IR8 c moment d'inertie Rext**2 *Mext INT_0^M r^2 dm ALLOCATE(vtr(1,n_qs),vtm(n_qs),vtmt(n_qs+ord_qs)) DO i=1,n_qs c extraction R et M IF(en_m23)THEN vtr(1,i)=bp(3,m_qs*(i-1)+1) !r^2 vtm(i)=SQRT(bp(5,m_qs*(i-1)+1))**3 !m ELSE vtr(1,i)=bp(3,m_qs*(i-1)+1)**2 !r^2 vtm(i)=bp(5,m_qs*(i-1)+1)**3 !m ENDIF ENDDO c spline r^2(m) CALL bsp1dn(1,vtr,vtm,vtmt,n_qs,m_ch,knotw,.FALSE.,vtm(1), 1 lq,f,dfdq) c intégrale, moment d'inertie CALL sum_n(1,vtr,vtmt,m_ch,knotw,.FALSE.,vtm(1),vtm(n_qs),f) mmtI=f(1)*2.d0/3.d0 !finalement c---------------l_tr évolution début--------------------- c moments d'inertie des ZR et ZC pour le csv de Tristan c tot_conv=.TRUE. lim=1 jlim(1)=n_qs m_zc(1)=mstar r_zc(1)=rn IF(l_tr)THEN M_iner=0.d0 ; S_entro=0.d0 c entropie au centre pp=EXP(bp(1,1)) ; tp=EXP(bp(2,1)) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,mc(1), 1 lq,xchim,dxchim) xchimg=xchim*nucleo S_0atm(0)=calc_entropie(pp,tp,xchimg) c moment d'inertie cas totalement ZC ou ZR est calculé 11 lignes plus haut IF(tot_conv .OR. tot_rad)THEN MI_atm=f(1) ELSE c moments d'inertie de la jlim-(i-1)éme à la jlim-iéme limite DO i=1,lim IF(i == 1)THEN CALL sum_n(1,vtr,vtmt,m_ch,knotw,.TRUE.,vtm(1), 1 vtm(jlim(i)),f) ELSE CALL sum_n(1,vtr,vtmt,m_ch,knotw,.TRUE.,vtm(jlim(i-1)), 1 vtm(jlim(i)),f) ENDIF M_iner(i)=f(1) c entropie(i) à la jlim-iéme limite pp=EXP(bp(1,m_qs*(jlim(i)-1)+1)) tp=EXP(bp(2,m_qs*(jlim(i)-1)+1)) nu=bp(5,m_qs*(jlim(i)-1)+1) IF(.NOT.en_m23)nu=nu**2 !m^2/3 CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,nu, 1 lq,xchim,dxchim) xchimg=xchim*nucleo S_entro(i)=calc_entropie(pp,tp,xchimg) ENDDO c moments d'inertie de la lim-iéme (derniére) limite à la surface CALL sum_n(1,vtr,vtmt,m_ch,knotw,.TRUE.,vtm(jlim(lim)), 1 vtm(n_qs),f) MI_atm=f(1) ENDIF !tot_conv c facteur 2/3 M_iner=2.d0/3.d0*M_iner ; MI_atm=2.d0/3.d0*MI_atm c entropie à la surface pp=EXP(bp(1,m_qs*(n_qs-1)+1)) tp=EXP(bp(2,m_qs*(n_qs-1)+1)) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,mc(n_ch), 1 lq,xchim,dxchim) xchimg=xchim*nucleo S_0atm(1)=calc_entropie(pp,tp,xchimg) c calcul des temps caractérisques des ZC t_conv=0.d0 c ntest devra être égal à nzc (sauf pour modèles initiaux ZAMS et PMS) ntest=0 IRev: IF(tot_conv)THEN ntest=ntest+1 jinf=1 ; jsup=n_qs r_lim(1)=bp(3,m_qs*(jinf-1)+1) r_lim(2)=bp(3,m_qs*(jsup-1)+1) IF(en_m23)r_lim=SQRT(r_lim) c temps caractéristique de convection t_conv(ntest)=t_convec(r_lim) c il ya des ZC ELSEIF(.NOT.tot_rad)THEN IRev c si lim est pair ZC au fond IF(MOD(lim,2) == 0)THEN idep=2 ntest=ntest+1 jinf=1 ; jsup=jlim(1) r_lim(1)=bp(3,m_qs*(jinf-1)+1) ; r_lim(2)=bp(3,m_qs*(jsup-1)+1) IF(en_m23)r_lim=SQRT(r_lim) c temps caractéristique de convection t_conv(ntest)=t_convec(r_lim) c si lim est impair ZR au fond ELSE idep=1 ENDIF c les rayons des limites des ZC DO i=idep,lim-1,2 ntest=ntest+1 jinf=jlim(i) ; jsup=jlim(i+1) r_lim(1)=bp(3,m_qs*(jinf-1)+1) ; r_lim(2)=bp(3,m_qs*(jsup-1)+1) IF(en_m23)r_lim=SQRT(r_lim) c temps caractéristique de convection t_conv(ntest)=t_convec(r_lim) ENDDO c pour la ZC externe ntest=ntest+1 jinf=jlim(lim) ; jsup=n_qs r_lim(1)=bp(3,m_qs*(jinf-1)+1) ; r_lim(2)=bp(3,m_qs*(jsup-1)+1) IF(en_m23)r_lim=SQRT(r_lim) c temps caractéristique de convection t_conv(ntest)=t_convec(r_lim) ENDIF IRev c test c PRINT*,'resout',ntest,nzc,lim c WRITE(*,2000)M_iner(1:lim),MI_atm c WRITE(*,2000)t_conv(1:ntest) c WRITE(*,2000)S_entro(1:lim),S_0atm c PAUSE'en évolution' ENDIF !l_tr c-----------------l_tr évolution fin--------------------------- c suppression des tableaux pour l'intégration du moment d'inertie DEALLOCATE(vtr,vtm,vtmt) c rotation SELECT CASE(Krot) c calcul de la nouvelle valeur de la vitesse angulaire c cas de la rotation solide avec conservation globale CASE(2) ALLOCATE(vtr(1,n_qs),vtm(n_qs),vtmt(n_qs+m_ch)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 8 dans resout' IF(en_m23)THEN vtr(1,i)=f(3) !r^2 vtm(i)=SQRT(f(5))**3 !m ELSE vtr(1,i)=f(3)**2 !r^2 vtm(i)=f(5)**3 !m ENDIF ENDDO vtm(1)=0.d0 ; vtr(1,1)=0.d0 CALL bsp1dn(1,vtr,vtm,vtmt,n_qs,m_ch,knotw,.FALSE., 1 vtm(1),lq,f,dfdq) IF(no_croiss)THEN PRINT*,'arrêt 1 dans resout' ; CALL sortie('resout 5') ENDIF CALL sum_n(1,vtr,vtmt,m_ch,knotw,.FALSE.,vtm(1),vtm(n_qs),f) c nouvelle valeur de la vitesse angulaire wrot=mw_tot/f(1) SELECT CASE(langue) CASE('english') WRITE(usl_static,1382)wrot,mw_tot*2.d0/3.d0 1382 FORMAT('angular velocity=',es10.3, 1 ' rad/sec, total angular momentum=',es10.3,' sol. unit',/) CASE DEFAULT WRITE(usl_static,382)wrot,mw_tot*2.d0/3.d0 382 FORMAT('Vitesse angulaire=',es10.3, 1 ' rad/sec, Moment angulaire total=',es10.3,' unite sol.',/) END SELECT DEALLOCATE(vtr,vtm,vtmt) c cas de la rotation non solide, calcul du moment cinétique total CASE(3,4,5) ALLOCATE(dfrot(nrot),frot(nrot),vtr(1,n_qs),vtm(n_qs), 1 vtmt(n_qs+m_ch)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 8.1 dans resout' IF(en_m23)THEN vtr(1,i)=f(3) !r^2 vtm(i)=SQRT(f(5))**3 !m ELSE vtr(1,i)=f(3)**2 !r^2 vtm(i)=f(5)**3 !m f(5)=f(5)**2 !m^2/3 ENDIF CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lq,frot,dfrot) IF(no_croiss)PRINT*,'Pb. en 8.2 dans resout' vtr(1,i)=vtr(1,i)*frot(1) ENDDO vtm(1)=0.d0 ; vtr(1,1)=0.d0 CALL bsp1dn(1,vtr,vtm,vtmt,n_qs,m_ch,knotw,.FALSE., 1 vtm(1),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. at 03 in resout' CALL sum_n(1,vtr,vtmt,m_ch,knotw,.FALSE.,vtm(1),vtm(n_qs),f) mcin_tot=f(1) IF(ini_mcin)THEN mcin_tot0=f(1) ; mcin_totp=f(1) ; ini_mcin=.FALSE. ENDIF SELECT CASE(langue) CASE('english') WRITE(usl_static,1383)mcin_tot, 1 NINT((mcin_tot/mcin_totp-1.d0)*100.d0), 2 NINT((mcin_tot/mcin_tot0-1.d0)*100.d0) WRITE(2,1383)mcin_tot,NINT((mcin_tot/mcin_totp-1.d0)*100.d0), 1 NINT((mcin_tot/mcin_tot0-1.d0)*100.d0) 1383 FORMAT(/,'Total angular momentum (solar units)=',es10.3,/, 1 'Change on the time step=',i2,'%, Total change=',i2,'%'/) CASE DEFAULT WRITE(usl_static,383)mcin_tot, 1 NINT((mcin_tot/mcin_totp-1.d0)*100.d0), 2 NINT((mcin_tot/mcin_tot0-1.d0)*100.d0) WRITE(2,383)mcin_tot,NINT((mcin_tot/mcin_totp-1.d0)*100.d0), 1 NINT((mcin_tot/mcin_tot0-1.d0)*100.d0) 383 FORMAT(/,'Moment cinétique total (unités sol)=',es10.3,/, 1 'Variation sur le pas temporel=',i3,'%, Variation tot.=',i3,'%'/) END SELECT mcin_totp=f(1) DEALLOCATE(dfrot,frot,vtr,vtm,vtmt) END SELECT c kdes_stat=1: écriture/plot des variables quasi-statiques IF(Kdes_stat == 2)CALL ecrit_static c l'évolution entre àge et àge+dt est terminée, retour vers cesam SELECT CASE(langue) CASE('english') WRITE(usl_static,1003) 1003 FORMAT(/,'--- Integration of the quasi-static model (end)----') CASE DEFAULT WRITE(usl_static,3) 3 FORMAT(/,'---- Intégration du modèle quasi-statique (fin)----') END SELECT c PAUSE'resout solution' c arrêt si t_stop est dépassé au centre lt_stop=lt_stop .OR. bp_t(2,1) > lnt_stop c ------------modèle quasi-statique de ZAMS------------------------- CASE(2) SELECT CASE(langue) CASE('english') WRITE(usl_static,1041) ; WRITE(2,1041) 1041 FORMAT(/,'**** Integration of the ZAMS quasi-static model (begin) ****') CASE DEFAULT WRITE(usl_static,41) ; WRITE(2,41) 41 FORMAT(/,'**** Intégration du modèle quasi-statique de ZAMS(début) ****') END SELECT c détermination du nombre de couches suivant la précision requise c WRITE(*,2000)bp(6,1),psi0, bp(6,1)-psi0, dpsi c PAUSE'psi0' c nombre de couches et répartition new_nqs=actu_nqs(.FALSE.) new_rep=.TRUE. c PAUSE'resout ZAMS' c initialisation, limitation de la variation du nombre de couches nqsp=-100 c détermination des limites ZR/ZC et des facteurs de répartition c précédée, éventuellement, de l'optimisation de la répartition c PRINT*,new_rep,new_nqs CALL lim_zc(new_rep,new_nqs) c Intégration modèle quasi-statique boucle Bnr2 infinie de NR compt=0 ; err=1.d20 Bnr2: DO CALL coll_qs(dt,compt,reprend,err,vare,qmax,corr) compt=compt+1 c écriture des variables quasi-statiques au max d'erreur WRITE(usl_static,*) WRITE(usl_static,100)compt,err,nom_qsr(vare),qmax,1.d0/corr CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(qmax),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 12 dans resout' IF(en_m23)THEN IF(pturb)THEN WRITE(usl_static,101)(EXP(f(i)),i=1,2),SQRT(ABS(f(3))), 1 (SQRT(ABS(f(i)))**3,i=4,5),EXP(f(7)) ELSE WRITE(usl_static,103)(EXP(f(i)),i=1,2),SQRT(ABS(f(3))), 1 (SQRT(ABS(f(i)))**3,i=4,5) ENDIF ELSE !en m^1/3 IF(pturb)THEN WRITE(usl_static,101)(EXP(f(i)),i=1,2),(f(i),i=3,4),f(5)**3, 1 EXP(f(7)) ELSE WRITE(usl_static,103)(EXP(f(i)),i=1,2),(f(i),i=3,4),f(5)**3 ENDIF f(5)=f(5)**2 !m^2/3 ENDIF c WRITE(*,2000)bp(6,1) ; PAUSE'bp(6,1)' c PAUSE'après écritures (ZAMS)' c la réactualisation de m_stat et r_stat qui accélérent la recherche c dans le ssp. inter est effectuée dans coll_qs c pour debug, écriture de la solution intermédiaire IF(.FALSE.)THEN c IF(.TRUE.)THEN PAUSE'écriture' ALLOCATE(l(n_qs),m(n_qs),p(n_qs),pt(n_qs),r(n_qs),t(n_qs)) DO i=1,n_qs !convergence totale si dt>0 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 13 dans resout' pt(i)=EXP(f(1)) !variable ln Ptot t(i)=EXP(f(2)) !variable ln T IF(pturb)THEN !avec pression turbulente 8 inconnues p(i)=EXP(f(7)) !variable ln Pgaz ELSE !sans pression turbulente 7 inconnues p(i)=pt(i) ENDIF IF(en_m23)THEN r(i)=SQRT(ABS(f(3))) !rayon/Rsol l(i)=SQRT(ABS(f(4)))**3 !l/Lsol m(i)=SQRT(ABS(f(5)))**3 !m/Msol ELSE r(i)=ABS(f(3)) !rayon/Rsol l(i)=f(4) !l/Lsol m(i)=f(5)**3 !m/Msol f(5)=f(5)**2 !m^2/3 ENDIF PRINT*,i ; WRITE(usl_static,2000)pt(i),p(i),t(i),r(i),l(i), 1 m(i),f(6) ENDDO DEALLOCATE(l,m,p,pt,r,t) ; PAUSE'solution intermédiaire' ENDIF !écriture de la solution intermédiaire c sur la ZAMS on impose la cohérence de psi0 et de n_qs IS0: IF(compt >= 2 .AND. err < precix .AND. nqsp == n_qs)THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1215) 1215 FORMAT('the model of ZAMS is obtained') CASE DEFAULT WRITE(usl_static,215) 215 FORMAT('le modèle de ZAMS a convergé') END SELECT WRITE(2,58)n_qs,n_ch c modèle totalement convectif: lim=1,jlim(1)=n,lconv(1)=.FALSE. c modèle totalement radiatif: lim=0,jlim(i)=-100,lconv(i)=.FALSE. IF(lim > 0)THEN IF(lim == 1 .AND. jlim(1) == 1 .AND. .NOT. lconv(1))THEN WRITE(2,7) ELSE WRITE(2,*) ; WRITE(2,9) DO j=1,lim IF(ovsht)THEN IF(lconv(j))THEN IF(ovshts > 0.d0)THEN !overshoot supérieur pas d'extension WRITE(2,1221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ELSE !overshoot inférieur extension WRITE(2,122)jlim(j),m_zc(j)/mstar,r_ov(j)/rstar ENDIF ELSE !fin de ZC IF(ovshts > 0.d0)THEN !overshoot supérieur extension WRITE(2,222)jlim(j),m_zc(j)/mstar,r_ov(j)/rstar ELSE !overshoot inférieur pas d'extension WRITE(2,2221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ENDIF ENDIF ELSE !sans overshoot IF(lconv(j))THEN WRITE(2,1221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ELSE WRITE(2,2221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ENDIF ENDIF ENDDO ENDIF ELSE WRITE(2,8) ENDIF PRINT* c moment d'inertie Rsol**2 *Msol INT_0^M r dm ALLOCATE(vtr(1,n_qs),vtm(n_qs),vtmt(n_qs+ord_qs)) DO i=1,n_qs c extraction R et M IF(en_m23)THEN vtr(1,i)=bp(3,m_qs*(i-1)+1) !r^2 vtm(i)=SQRT(bp(5,m_qs*(i-1)+1))**3 !m ELSE vtr(1,i)=bp(3,m_qs*(i-1)+1)**2 !r^2 vtm(i)=bp(5,m_qs*(i-1)+1)**3 !m ENDIF ENDDO c spline r^2(m) CALL bsp1dn(1,vtr,vtm,vtmt,n_qs,m_ch,knotw,.FALSE., 1 vtm(1),lq,f,dfdq) c intégrale, moment d'inertie CALL sum_n(1,vtr,vtmt,m_ch,knotw,.FALSE.,vtm(1),vtm(n_qs),f) mmtI=f(1)*2.d0/3.d0 c----------------l_tr ZAMS début----------------------- c moments d'inertie des ZR et ZC pour le csv de Tristan c tot_conv=.TRUE. lim=1 jlim(1)=n_qs m_zc(1)=mstar r_zc(1)=rn IF(l_tr)THEN M_iner=0.d0 ; S_entro=0.d0 c entropie au centre pp=EXP(bp(1,1)) ; tp=EXP(bp(2,1)) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,mc(1), 1 lq,xchim,dxchim) xchimg=xchim*nucleo S_0atm(0)=calc_entropie(pp,tp,xchimg) c moment d'inertie cas totalement ZC ou ZR calculé précédemment IF(tot_conv .OR. tot_rad)THEN MI_atm=f(1) ELSE c moments d'inertie de la jlim-(i-1)éme à la jlim-iéme limite DO i=1,lim IF(i == 1)THEN CALL sum_n(1,vtr,vtmt,m_ch,knotw,.TRUE.,vtm(1), 1 vtm(jlim(i)),f) ELSE CALL sum_n(1,vtr,vtmt,m_ch,knotw,.TRUE.,vtm(jlim(i-1)), 1 vtm(jlim(i)),f) ENDIF M_iner(i)=f(1) c entropie(i) à la jlim-iéme limite pp=EXP(bp(1,m_qs*(jlim(i)-1)+1)) tp=EXP(bp(2,m_qs*(jlim(i)-1)+1)) nu=bp(5,m_qs*(jlim(i)-1)+1) IF(.NOT.en_m23)nu=nu**2 !m^2/3 CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,nu, 1 lq,xchim,dxchim) xchimg=xchim*nucleo S_entro(i)=calc_entropie(pp,tp,xchimg) ENDDO c moments d'inertie de la lim-iéme (derniére) limite à la surface CALL sum_n(1,vtr,vtmt,m_ch,knotw,.TRUE.,vtm(jlim(lim)), 1 vtm(n_qs),f) MI_atm=f(1) ENDIF !tot_conv c facteur 2/3 M_iner=2.d0/3.d0*M_iner ; MI_atm=2.d0/3.d0*MI_atm c entropie à la surface pp=EXP(bp(1,m_qs*(n_qs-1)+1)) tp=EXP(bp(2,m_qs*(n_qs-1)+1)) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,mc(n_ch), 1 lq,xchim,dxchim) xchimg=xchim*nucleo S_0atm(1)=calc_entropie(pp,tp,xchimg) c calcul des temps carac. des ZC, on distingue les cas lim pair/impair t_conv=0.d0 c ntest devra être égal à nzc (sauf pour les modèles init de ZAMS et de PMS) ntest=0 IRzc: IF(tot_conv)THEN ntest=ntest+1 ; jinf=1 ; jsup=n_qs r_lim(1)=bp(3,m_qs*(jinf-1)+1) ; r_lim(2)=bp(3,m_qs*(jsup-1)+1) IF(en_m23)r_lim=SQRT(r_lim) c temps caractéristique de convection t_conv(ntest)=t_convec(r_lim) c il ya des ZC ELSEIF(.NOT.tot_rad)THEN IRzc c avec lim pair: ZC au fond IF(MOD(lim,2) == 0)THEN idep=2 ntest=ntest+1 ; jinf=1 ; jsup=jlim(1) r_lim(1)=bp(3,m_qs*(jinf-1)+1) r_lim(2)=bp(3,m_qs*(jsup-1)+1) IF(en_m23)r_lim=SQRT(r_lim) c temps caractéristique de convection t_conv(ntest)=t_convec(r_lim) c lim impair: ZR au fond ELSE idep=1 ENDIF c les rayons des limites des ZC DO i=idep,lim-1,2 ntest=ntest+1 jinf=jlim(i) ; jsup=jlim(i+1) r_lim(1)=bp(3,m_qs*(jinf-1)+1) r_lim(2)=bp(3,m_qs*(jsup-1)+1) IF(en_m23)r_lim=SQRT(r_lim) c temps caractéristique de convection t_conv(ntest)=t_convec(r_lim) ENDDO c pour la ZC externe ntest=ntest+1 jinf=jlim(lim) ; jsup=n_qs r_lim(1)=bp(3,m_qs*(jinf-1)+1) ; r_lim(2)=bp(3,m_qs*(jsup-1)+1) IF(en_m23)r_lim=SQRT(r_lim) c temps caractéristique de convection t_conv(ntest)=t_convec(r_lim) nzc=ntest !pour écriture dans le fichier csv ENDIF IRzc c test c PRINT*,ntest ; WRITE(*,2000)t_conv(1:ntest) c WRITE(*,2000)S_entro(1:lim),S_0atm c PAUSE'en zams' ENDIF !l_tr c-------------l_tr ZAMS fin ----------------------------------- c suppression des tableaux pour l'intégration du moment d'inertie DEALLOCATE(vtr,vtm,vtmt) c rotation SELECT CASE(Krot) c calcul du moment angulaire total du modèle de ZAMS CASE(2) ALLOCATE(vtr(1,n_qs),vtm(n_qs),vtmt(n_qs+m_ch)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 14 dans resout' IF(en_m23)THEN vtr(1,i)=f(3) !r^2 vtm(i)=SQRT(f(5))**3 !m ELSE vtr(1,i)=f(3)**2 !r^2 vtm(i)=f(5)**3 !m ENDIF ENDDO vtm(1)=0.d0 ; vtr(1,1)=0.d0 CALL bsp1dn(1,vtr,vtm,vtmt,n_qs,m_ch,knotw,.FALSE.,vtm(1), 1 lq,f,dfdq) IF(no_croiss)THEN PRINT*,'arrêt 3 dans resout' ; CALL sortie('resout 6') ENDIF CALL sum_n(1,vtr,vtmt,m_ch,knotw,.FALSE.,vtm(1),vtm(n_qs),f) mw_tot=f(1)*wrot SELECT CASE(langue) CASE('english') WRITE(usl_static,1038)wrot,mw_tot*2.d0/3.d0 1038 FORMAT('For the homogeneous ZAMS model:',/, 1 'angular velocity=',es10.3, 2 ' rad/sec, total angular momentum =',es10.3,' sol. unit',/) CASE DEFAULT WRITE(usl_static,38)wrot,mw_tot*2.d0/3.d0 38 FORMAT('Pour le modèle de ZAMS homogène:',/, 1 'Vitesse angulaire=',es10.3, 2 ' rad/sec, Moment angulaire total=',es10.3,' unité sol.',/) END SELECT DEALLOCATE(vtr,vtm,vtmt) c cas de la rotation non solide, calcul du moment cinétique total CASE(3,4,5) ALLOCATE(dfrot(nrot),frot(nrot),vtr(1,n_qs),vtm(n_qs), 1 vtmt(n_qs+m_ch)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 8.2 dans resout' vtr(1,i)=ABS(f(3)) ; vtm(i)=ABS(f(5)) IF(en_m23)THEN vtr(1,i)=ABS(f(3)) !r^2 vtm(i)=SQRT(f(5))**3 !m ELSE vtr(1,i)=f(3)**2 !r^2 vtm(i)=ABS(f(5))**3 !m f(5)=f(5)**2 !m^2/3 ENDIF CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lq,frot,dfrot) IF(no_croiss)PRINT*,'Pb. at 05 in resout' vtr(1,i)=vtr(1,i)*frot(1) ENDDO vtm(1)=0.d0 ; vtr(1,1)=0.d0 CALL bsp1dn(1,vtr,vtm,vtmt,n_qs,m_ch,knotw,.FALSE.,vtm(1), 1 lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. at 06 in resout' CALL sum_n(1,vtr,vtmt,m_ch,knotw,.FALSE.,vtm(1),vtm(n_qs),f) mcin_tot=f(1) ; mcin_tot0=f(1) ; mcin_totp=f(1) ; ini_mcin=.FALSE. DEALLOCATE(dfrot,frot,vtr,vtm,vtmt) END SELECT c initialisation de l'interpolation de TdS/dt c on ne tient pas compte de la différence Pgaz, Ptot pour c l'énergie graviphique, TdS/dt=tds qui est tabulée en fonction c de m_stat=m^2/3 ou de m_stat=m^1/3 ; m_stat a été formé dans coll_qs c TdS/dt~0 sur la ZAMS n_tds=n_qs IF(ALLOCATED(tds))DEALLOCATE(tds,x_tds,xt_tds) ALLOCATE(tds(1,n_tds),x_tds(n_tds),xt_tds(n_tds+m_tds)) tds=0.d0 ; x_tds=ABS(m_stat) c vérif. de la stricte croissance des abscisses pour les 40 dernières couches j=n_tds Btds1: DO IF(j <= n_qs-40)EXIT Btds1 IF(x_tds(j) <= x_tds(j-1))n_tds=j-1 j=j-1 ENDDO Btds1 CALL noein(x_tds,xt_tds,n_tds,m_tds,knot_tds) IF(no_croiss)THEN PRINT*,'arrêt 6 dans resout' ; CALL sortie('resout 7') ENDIF c modèle convergé, sortie de brn2, puis retour vers cesam EXIT Bnr2 ELSEIF(compt >= iter_max0)THEN IS0 WRITE(usl_static,220) ; WRITE(2,220) ; CALL sortie('resout 8') 220 FORMAT('no conv. du modèle quasi-statique de ZAMS, ARRET') ELSE IS0 c plus d'itérations sont nécessaires, établissement de n_qs en cohérence avec psi0 nqsp=n_qs new_nqs=actu_nqs(.FALSE.) CALL lim_zc(new_rep,new_nqs) ENDIF IS0 !compt >= 2 ENDDO Bnr2 c le modèle quasi-static de ZAMS a convergé, retour vers cesam SELECT CASE(langue) CASE('english') WRITE(usl_static,1051) ; WRITE(2,1051) 1051 FORMAT('**** Integration of the ZAMS quasi-static model (end) ****',/) CASE DEFAULT WRITE(usl_static,51) ; WRITE(2,51) 51 FORMAT('**** Intégration modèle quasi-statique de ZAMS(fin) ****',/) END SELECT c ------- modèle quasi-statique de ZAMS fin, de PMS début -------------- c un23=3 pour le premier modèle de PMS calculé avec c_iben c un23=-3 pour le second modèle de PMS calculé avec c_iben*1.1 c le modèle est totalement convectif CASE(-3, 3) SELECT CASE(langue) CASE('english') WRITE(usl_static,1061) 1061 FORMAT(/,'**** Integration of the PMS quasi-static model (begin) ****',/) CASE DEFAULT WRITE(usl_static,61) 61 FORMAT(/,'**** Intégration du modèle quasi-statique de PMS (début) ****',/) END SELECT new_rep=.TRUE. c un23=3 pour le premier modèle de PMS calculé avec c_iben IF(un23 == 3)THEN c au début du calcul d'un pas temporel, dans lim_zc, la répartition sera c actualisée si le nombre de couches change ou (disjonctif) si les facteurs de c répartition sont modifiés c new_nqs=actu_nqs(.FALSE.) ; new_rep=.TRUE. new_nqs=n_qs ; new_rep=.FALSE. c PRINT*,new_nqs,n_qs ; PAUSE'resout esp PMS' c détermination des limites ZR/ZC et des facteurs de répartition c précédée, éventuellement, de l'optimisation de la répartition CALL lim_zc(new_rep,new_nqs) ENDIF c pour tous les modèles initiaux de PMS, un23=3 et -3 c Intégration du modèle quasi-statique boucle Bnr3 infinie de NR compt=0 ; err=1.d20 Bnr3: DO CALL coll_qs(dt,compt,reprend,err,vare,qmax,corr) compt=compt+1 c écriture des variables quasi-statiques au max d'erreur WRITE(usl_static,*) WRITE(usl_static,100)compt,err,nom_qsr(vare),qmax,1.d0/corr CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(qmax),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 15 dans resout' IF(en_m23)THEN IF(pturb)THEN WRITE(usl_static,101)(EXP(f(i)),i=1,2),SQRT(ABS(f(3))), 1 (SQRT(ABS(f(i)))**3,i=4,5),EXP(f(7)) ELSE WRITE(usl_static,103)(EXP(f(i)),i=1,2),SQRT(ABS(f(3))), 1 (SQRT(ABS(f(i)))**3,i=4,5) ENDIF ELSE !en rayon IF(pturb)THEN WRITE(usl_static,101)(EXP(f(i)),i=1,2),(f(i),i=3,4),f(5)**3, 1 EXP(f(7)) ELSE WRITE(usl_static,103)(EXP(f(i)),i=1,2),(f(i),i=3,4),f(5)**3 ENDIF f(5)=f(5)**2 !m^2/3 ENDIF c PAUSE'après écritures (PMS)' c la réactualisation de m_stat et r_stat qui accélérent la recherche dans c le ssp. inter valable pour variables m_stat=m^2/3, r_stat=r^2 ou c m_stat=m^1/3, r_stat=r est effectuée dans coll_qs c pour écrire (debug) la solution intermédiaire IF(.FALSE.)THEN c IF(.TRUE.)THEN PRINT*,'un23=',un23 PAUSE'écriture' ALLOCATE(l(n_qs),m(n_qs),p(n_qs),pt(n_qs),r(n_qs),t(n_qs)) DO i=1,n_qs !convergence totale si dt>0 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 1 dans resout' pt(i)=EXP(f(1)) !variable ln Ptot t(i)=EXP(f(2)) !variable ln T IF(pturb)THEN !avec pression turbulente 8 inconnues p(i)=EXP(f(7)) !variable ln Pgaz ELSE !sans pression turbulente 7 inconnues p(i)=pt(i) ENDIF IF(en_m23)THEN r(i)=SQRT(ABS(f(3))) !rayon/Rsol l(i)=SQRT(ABS(f(4)))**3 !l/Lsol m(i)=SQRT(ABS(f(5)))**3 !m/Msol ELSE r(i)=ABS(f(3)) !rayon/Rsol l(i)=f(4) !l/Lsol m(i)=f(5)**3 !m/Msol f(5)=f(5)**2 !m^2/3 ENDIF PRINT*,i WRITE(usl_static,2000)pt(i),p(i),t(i),r(i),l(i),m(i),f(6) ENDDO DEALLOCATE(l,m,p,pt,r,t) PAUSE'solution intermédiaire' ENDIF !écriture de la solution intermédiaire c gestion des itérations Icv: IF(compt >= 2 .AND. err < precix)THEN WRITE(usl_static,*) WRITE(usl_static,*)'modèle de PMS convergé' c avec c_iben =3, 1ier modèle PMS tot.conv. convergé, sortie de brn3, retour vers cesam c acquisition d'une nouvelle valeur ou une petite variation de c_iben c reprise avec un23=-3 si c'est Ok, 3 sinon IF(un23 == 3)EXIT Bnr3 !vers ligne 1694 WRITE(2,58)n_qs,n_ch c traitement du modèle initial de PMS et retour vers CESAM c modèle totalement convectif: lim=1,jlim(1)=n,lconv(1)=.FALSE. c modèle totalement radiatif: lim=0,jlim(i)=-100,lconv(i)=.FALSE. IF(lim > 0)THEN IF(lim == 1 .AND. jlim(1) == 1 .AND. .NOT. lconv(1))THEN WRITE(2,7) ELSE WRITE(2,*) ; WRITE(2,9) DO j=1,lim IF(ovsht)THEN IF(lconv(j))THEN IF(ovshts > 0.d0)THEN !overshoot sup. pas d'extension WRITE(2,1221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ELSE !overshoot inférieur extension WRITE(2,122)jlim(j),m_zc(j)/mstar,r_ov(j)/rstar ENDIF ELSE !fin de ZC IF(ovshts > 0.d0)THEN !overshoot supérieur extension WRITE(2,222)jlim(j),m_zc(j)/mstar,r_ov(j)/rstar ELSE !overshoot inférieur pas d'extension WRITE(2,2221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ENDIF ENDIF ELSE !sans overshoot IF(lconv(j))THEN WRITE(2,1221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ELSE WRITE(2,2221)jlim(j),m_zc(j)/mstar,j,r_zc(j)/rstar ENDIF ENDIF ENDDO ENDIF ELSE WRITE(2,8) ENDIF PRINT* c moment d'inertie Rsol**2 *Msol INT_0^M r dm ALLOCATE(vtr(1,n_qs),vtm(n_qs),vtmt(n_qs+ord_qs)) DO i=1,n_qs c extraction R et M IF(en_m23)THEN vtr(1,i)=bp(3,m_qs*(i-1)+1) !r^2 vtm(i)=SQRT(bp(5,m_qs*(i-1)+1))**3 !m ELSE vtr(1,i)=bp(3,m_qs*(i-1)+1)**2 !r^2 vtm(i)=bp(5,m_qs*(i-1)+1)**3 !m ENDIF ENDDO c spline r^2(m) CALL bsp1dn(1,vtr,vtm,vtmt,n_qs,m_ch,knotw,.FALSE.,vtm(1), 1 lq,f,dfdq) c intégrale, moment d'inertie CALL sum_n(1,vtr,vtmt,m_ch,knotw,.FALSE.,vtm(1),vtm(n_qs),f) mmtI=f(1)*2.d0/3.d0 c----------------l_tr PMS début----------------------- c moments d'inertie des ZR et ZC pour le csv de Tristan c le modèle initial de PMS est totalement convectif c tot_conv=.TRUE. lim=1 jlim(1)=n_qs m_zc(1)=mstar r_zc(1)=rn IF(l_tr)THEN M_iner=0.d0 ; S_entro=0.d0 MI_atm=2.d0/3.d0*f(1) c entropie au centre pp=EXP(bp(1,1)) ; tp=EXP(bp(2,1)) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,mc(1), 1 lq,xchim,dxchim) xchimg=xchim*nucleo S_0atm(0)=calc_entropie(pp,tp,xchimg) c entropie à la surface pp=EXP(bp(1,m_qs*(n_qs-1)+1)) tp=EXP(bp(2,m_qs*(n_qs-1)+1)) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.TRUE.,mc(n_ch), 1 lq,xchim,dxchim) xchimg=xchim*nucleo S_0atm(1)=calc_entropie(pp,tp,xchimg) c modèle totalement convectif, calcul du temps caractérisque de la ZC t_conv=0.d0 ; ntest=1 jinf=1 ; jsup=n_qs r_lim(1)=bp(3,m_qs*(jinf-1)+1) r_lim(2)=bp(3,m_qs*(jsup-1)+1) IF(en_m23)r_lim=SQRT(r_lim) c temps caractéristique de convection t_conv(1)=t_convec(r_lim) nzc=ntest !pour écritures dans le ficier csv c test c PRINT*,ntest c WRITE(*,2000)t_conv(1:ntest) c WRITE(*,2000)S_entro(1:2) c PAUSE'en PMS' ENDIF c-------------l_tr PMS fin ----------------------------------- c suppression des tableaux pour l'intégration du moment d'inertie DEALLOCATE(vtr,vtm,vtmt) c rotation SELECT CASE(Krot) c calcul du moment angulaire total du modèle de PMS CASE(2) ALLOCATE(vtr(1,n_qs),vtm(n_qs),vtmt(n_qs+m_ch)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 2 dans resout' IF(en_m23)THEN vtr(1,i)=f(3) !r^2 vtm(i)=SQRT(f(5))**3 !m ELSE vtr(1,i)=f(3)**2 !r^2 vtm(i)=f(5)**3 !m ENDIF ENDDO vtm(1)=0.d0 ; vtr(1,1)=0.d0 CALL bsp1dn(1,vtr,vtm,vtmt,n_qs,m_ch,knotw,.FALSE.,vtm(1),lq,f,dfdq) IF(no_croiss)THEN PRINT*,'arrêt 4 dans resout' ; CALL sortie('resout 9') ENDIF CALL sum_n(1,vtr,vtmt,m_ch,knotw,.FALSE.,vtm(1),vtm(n_qs),f) mw_tot=f(1)*wrot SELECT CASE(langue) CASE('english') WRITE(usl_static,1381)wrot,mw_tot*2.d0/3.d0 1381 FORMAT('For the homogeneous PMS model:',/, 1 'angular velocity=',es10.3, 2 ' rad/sec, total angular momentum =',es10.3,' sol. unit',/) CASE DEFAULT WRITE(usl_static,381)wrot,mw_tot*2.d0/3.d0 381 FORMAT('Pour le modèle de PMS homogène:',/, 1 'Vitesse angulaire=',es10.3, 2 ' rad/sec, Moment angulaire total=',es10.3,' unite sol.'/) END SELECT DEALLOCATE(vtr,vtm,vtmt) c cas de la rotation non solide, calcul du moment cinétique total CASE(3,4,5) ALLOCATE(dfrot(nrot),frot(nrot),vtr(1,n_qs),vtm(n_qs), 1 vtmt(n_qs+m_ch)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. en 8.3 dans resout' vtr(1,i)=f(3) ; vtm(i)=f(5) IF(en_m23)THEN vtr(1,i)=f(3) !r^2 vtm(i)=SQRT(f(5))**3 !m ELSE vtr(1,i)=f(3)**2 !r^2 vtm(i)=f(5)**3 !m f(5)=f(5)**2 !m^2/3 ENDIF CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MIN(f(5),mrot(n_rot)),lq,frot,dfrot) IF(no_croiss)PRINT*,'Pb. at 07 in resout' vtr(1,i)=vtr(1,i)*frot(1) ENDDO vtm(1)=0.d0 ; vtr(1,1)=0.d0 CALL bsp1dn(1,vtr,vtm,vtmt,n_qs,m_ch,knotw,.FALSE.,vtm(1),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. at 08 in resout' CALL sum_n(1,vtr,vtmt,m_ch,knotw,.FALSE.,vtm(1),vtm(n_qs),f) mcin_tot=f(1) ; mcin_tot0=f(1) ; mcin_totp=f(1) ; ini_mcin=.FALSE. DEALLOCATE(dfrot,frot,vtr,vtm,vtmt) END SELECT c initialisation de l'interpolation de TdS/dt c l'énergie graviphique, TdS/dt=tds est tabulée en fonction c de m_stat=m^2/3 en ou de m_stat=m^1/3 c pour le premier modèle de PMS avec c_iben, TdS/dt = 0 c on a pose un23=3 c pour le second modèle de PMS avec c_iben*1.1, TdS/dt \= 0 c on a pose un23=-3 c le tableau tds peut avoir été alloué lors du premier passage n_tds=n_qs IF(un23 < 0)THEN IF(ALLOCATED(tds))DEALLOCATE(tds,x_tds,xt_tds) ALLOCATE(tds(1,n_tds),x_tds(n_tds),xt_tds(n_tds+m_tds), 1 jac(nchim,nchim)) c la composition chimique étant uniforme xchim(1)=chim(1,1), etc. xchim(:)=chim(:,1) ; dxchim=0.d0 xchimg=xchim*nucleo DO i=1,n_tds CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. 3 en dans resout' IF(pturb)THEN !avec pression turbulente 8 inconnues pp=EXP(f(7)) !variable ln Pgaz ELSE !sans pression turbulente 7 inconnues pp=EXP(f(1)) ENDIF tp=EXP(f(2)) !variable ln T x_tds(i)=f(5) !m^2/3 ou m^1/3 CALL etat(pp,tp,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) CALL nuc(tp,ro,xchim,dxchim,jac,.FALSE.,3, 1 epsilon,depst,depsro,depsx,hh,be7,b8,n13,o15,f17) tds(1,i)=-epsilon(1)*secon6 !cT sort dans epsilon(1) ENDDO CALL bsp1dn(1,tds,x_tds,xt_tds,n_tds,m_tds,knot_tds,.FALSE., 1 x_tds(1),lq,f,dfdq) IF(no_croiss)THEN PRINT*,'arrêt 5 dans resout' ; CALL sortie('resout 10') ENDIF DEALLOCATE(jac) c on garde dans ***_t le modèle de PMS obtenu c CALL update(.TRUE.,dt,dts) ENDIF c modèle PMS convergé, sortie de brn3, puis retour vers cesam EXIT Bnr3 ELSEIF(compt >= iter_max0)THEN Icv PRINT* ; PRINT* PRINT*,'no conv. du modèle quasi-statique de PMS, ARRET' WRITE(2,*) ; WRITE(2,*) ; CALL sortie('resout 11') WRITE(2,*)'no conv. du modèle quasi-statique de PMS, ARRET' ENDIF Icv !compt >= 2 ENDDO Bnr3 c PRINT*,un23 c CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(1),lq,f,dfdq) c WRITE(*,2000)f(1:5) c CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(n_qs),lq,f,dfdq) c WRITE(*,2000)f(1:5) c PAUSE'sortie pms' c le modèle quasi-statique de PMS a convergé, retour vers cesam SELECT CASE(langue) CASE('english') WRITE(usl_static,1071) 1071 FORMAT('Integration of the PMS quasi-static model (end)',/) CASE DEFAULT WRITE(usl_static,71) 71 FORMAT('Intégration du modèle quasi-statique de PMS(fin)',/) END SELECT c fin de la gestion des un23 CASE DEFAULT PRINT*,'erreur |un23| est différent de 1, 2 ou 3, un23=',un23 PRINT*,'ARRET dans resout' ; CALL sortie('resout 12') END SELECT !pour un23 c le premier modèle d'initialisation ou de reprise a été calculé pass1=.FALSE. c retour général vers cesam (pour toutes les valeurs de un23) RETURN CONTAINS INCLUDE 'fcmax.f' END SUBROUTINE resout