c*********************************************************************** SUBROUTINE cesam c routine public du module mod_cesam c calcul de modèles de structure interne et d'évolution stellaire c programme principal: gestion générale des calculs, des écritures c méthode collocation de de Boor p. 277 c adapté à la structure interne, avec fonction de répartition c Auteur: P.Morel, Département J.D. Cassini, O.C.A., CESAM2k c 07 10 96 : gravité effective dans Väissälä c 10 10 96 : introduction de la rotation dans thermo c 21 03 97 : bug dans calcul de z_opa (signalé par Sacha) c 08 05 97 : introduction de l'option kipp pour calcul du TdS c 26 05 97 : fichiers cohe, coca, début de combustion helium, carbone c 17 09 06 : coox, cosi, début de combustion néon, oxygène, silicium c 11 06 97 : simplification de la gestion de t_inf c 26 06 97 : remplacement du moment angulaire par la vitesse angulaire c 18 07 97 : correction d'énergie exédentaire de rotation c 18 08 97 : modif nom elem chim c 25 08 97 : mise en place des variables eulériennes c 22 09 97 : abondances externes dans le fichier *.HR c 04 11 97 : abondances centrales dans le fichier *.HR c 18 11 97 : mstar dans le fichier *.HR c 12 11 99 : calcul de mue avec les taux d'ionisation c 19 11 99 : rectification de CALL time, des dates, de thermo, c utilisation de saha pour nh1, nhe1, nhe2, degene c 18 01 00 : correction si e_nuc=0 c 02 02 00 : e_grav < 50% e_nuc ----> série pple (au lieu de 99%) c 18 04 00 ; coeff_diff ---> diffm, difft c 30 07 00 : introduction F95 c avec pression turbulente 7 inconnues c sans pression turbulente 6 inconnues, Ptot=Pgaz c noms des routines génériques de physique c etat: équation d'état c opa: opacité c conv: convection c nuc: réactions thermonucléaires c atm: condition limite externe c tdetau: loi t(tau) c diffm: coefficients de diffusion microscopique c difft: coefficients de diffusion turbulente c diffw: coefficients de diffusion turbulente du moment angulaire c pertm: perte de masse c pertw: perte de moment cinétique c des : dessin du modèle durant l'exécution c ini_ctes : initialisation des constantes de physique c fonctions des principales routines internes numériques c resout : initialisation, formation, résolution du pb. aux limites c evol : gestion de l'évolution de la composition chimique c static_m23 : formation des équations de l'équilibre quasi stat. en m^2/3 c static_m13 : formation des équations de l'équilibre quasi stat. en m^2/3 c lim_zc : détermination des limites ZR/ZC c update : transfert t+dt ---> t ---> t-dt c diffus : calcul de la diffusion c coll_atm, eq_atm : restitution d'atmosphère c le centre est à la couche 1, la limite externe en couche n_qs c le nb. de couches est VARIABLE il est fixé dans le SSP lim_zc c une couche est placée à chaque limite ZR/ZC (à ~5% près) 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, TdS/dt=tds est tabulée en fonction c de m^2/3 en lagrangien et de m en eulérien, c incrément temporel : dt=delta t / 10^6 ans c variable indépendante : q l'indice (real) de couche c variables dépendantes : c ksi=ln P pression cgs c eta=ln T température K c dzeta=(r/rsol)^2 c lambda=(l/lsol)^2/3 c mu=(m/mtot)^2/3 c psi=dQ/dq=cte avec Q(mu,t)=fonction de répartition c composition chimique : xchim(i),i=1,nchim : abondances / mole c en f(m/Mtot^2/3) + W (w : vitesse angulaire) c ordre et indices d'identification des éléments definis dans c reac_nuc, Exemple: c H1 H2 He3 He4 Li7 Be7 C12 C13 N14 N15 O16 O17 W c 1 2 3 4 5 6 7 8 9 10 11 12 13 c avant le début d'un pas temporel: c mc_t,mct_t,nc_t,knotc_t,chim_t: comp. chim. à âge-dt c mc,mct,nc,knotc,chim: comp. chim. à âge c bp,q,qt,knot,p,t,r,l,m: var. pples. à âge c x_tds,xt_tds,n_tds,knot_tds,tds: TdS/dt à âge t c en_m23=.TRUE. variables ln Ptot, ln T, r^2, l^2/3, m^2/3, psi, [ln Pgaz] c en_m23=.FALSE. variables ln Ptot, ln T, r, l, m^1/3, psi, [ln Pgaz] c Fichiers utilisés c les fichiers de sortie ont le nom choisi par l'utilisateur c mon_modele suivit d'un attribut qui les distingue c l'extension "_B" est spécifique aux fichiers binaires c FOR002 : listing de sortie, créé dans cesam fermé et garde dans c cesam, nom : mon_modèle.lis c FOR003 : lecture des données, utilisé temporairement dans resout c nom : mon_modele.don c FOR004 : modèle initial ou repris (binaire), c utilisé temporairement dans cesam c FOR011 : lecture des tables d'opacité c utilisé temporairement dans routines d'opacité c FOR011-->017 : lecture des tables d'équation d'état c utilisé temporairement dans certains SSP d'équation d'état c FOR024 : écriture modèle initial d'âge 0 homogène en binaire, c PMS, ZAMS ouvert et fermé dans cesam c nom : mon_modele_B.hom/pms c FOR025 : dernier modèle, ouvert et fermé dans cesam c nom : mon_modele_B.dat en binaire c FOR026 point de reprise, ouvert et fermé dans cesam c nom : mon_modele_B.rep en binaire c FOR030 liste pour oscillation ouvert et fermé dans osc c nom : mon_modele.osc en ASCII c FOR031 : dernier modèle calculé de l'atmosphère restituée c ouvert et fermé dans lim_atm c nom : mon_modele_B.atm en binaire c FOR053 : fichier pour tracer un diag. HR et les ZC, c ouvert et fermé dans cesam c nom : mon_modele.HR en ASCII c FOR060 : tables d'EOS d'OPAL c~~~~B-splines, interpolation, équations différentielles résolution par c collocation, par éléments finis c pour une interpolation et pour une équation équation différentielle c résolue par éléments finis l'ordre des splines est noté m_*, exemple m_ch c pour une équation équation différentielle résolue par collocation c l'ordre des splines est noté ord_*, exemple : ord_qs c l'ordre est la somme d'un odre "primaire" noté m_*, exemple m_qs, et de c l'ordre de l'équation différentielle, exemple r_qs : ord_qs=m_qs+r_qs c dans CESAM2k les équations différentielles résolues par collocation c sont résolues sous la forme de systèmes du premier ordre ord_qs=1 c ce qui évite le calcul de dérivées d'ordre élevé des coefficients c~~~~ c formation des points de collocation, qui serviront pour rota et poisson c pour un ordre d'interpolation m_rot+r_qs=ord_rot, il y a m_rot points de c collocation dans chaque intervalle, la multiplicité des noeuds étant mrot c il y a ncoll_rot=(n_rot-1)*m_rot points de collocation et r_qs points limite c r_qs ordre des équations différentielles c ord_rot=m_rot+r_qs ordre d'interpolation c avec nd discontinuités il convient d'ajouter nzc*r_qs points limite c soit (nd+1)*r_qs points limite c sans discontinuité interne dimension du vecteur nodal knotr=dim_rot+m_rot+r_qs c avec discontinuités internes knotr=dim_rot+m_rot+r_qs(nzc+1) c dim_rot, knotr, éléments de mod_variables, définis dans base_rota c ord_rot est défini dans base_rota c m_rot, r_qs sont des éléments de mod_donnees c voir ligne 4379 pour écrire un modèle d'initialisation en ASCII c rechercher llist c--------------------------------------------------------------------- USE mod_atm, ONLY: atm, dlpp_atm, m_atm, p_atm, pt_atm, 1 r_atm, tau, tdetau, thermo_atm, t_atm USE mod_donnees, ONLY: ab_ini, agemax, ajuste, 1 all_osc, all_rep, alpha, aradia, arret, baratine, cephe, 2 cpturb, diffusion, dn_fixe, dpsi, 3 dpsim, dpsip, dn_sort, dtlist, dtmax, dtmin, dt0, d_lum, 4 d_press, d_ray, d_temp, en_m23, fcv, f_eos, f_opa, g, gmsol, 5 granr, grav_sol, grille_fixe, he_core, ihe4, iLi7, ini0, Ipg, Kdes_rot, 6 Krot, kipp, langue, lim_ro, lit_nl, lisse, li_ini, ln_Tli, loc_zc, 7 LOG_teff, lsol, l_fac, l_hr, l_lb, l_pertm, l_pertmt, l_tr, methode, 8 msol, mdot, mtot, mu_saha, m_ch, m_ptm, m_qs, m_rot, m_tds, ne, 9 nb_max_modeles, nchim, new_bv, nom_conv, nom_elem, nom_des, 1 nom_etat, nom_fich2, nom_opa, nom_nuc, nom_nuc_cpl, 2 nom_abon, nom_atm, nom_tdetau, nom_pertm, nom_pertw, 3 nom_diffm, nom_difft, nom_diffw, nom_ctes, nom_output, 4 nom_qs, nom_thw, no_discon, npt_lisse, nrot, nucleo, 5 nvth, n_atm, n_max, n_min_ZC, ordre, ord_qs, ord_rot, 6 pi, precic, precision, precit, precix, print_ctes, psi0, pturb, 7 q0, rep_atm, re_nu, rsol, rot_min, rot_solid, 8 r_qs, secon6, sigma, thw, typ, t_stop, unit, version, 9 w_rot, x_tams, x_stop, x0, y0, zi, z0 USE mod_etat, ONLY : etat, saha USE mod_evol, ONLY : ecrit_rota, initialise_rota, pulsation, tab_vth USE mod_kind USE mod_nuc, ONLY: nuc, planetoides, vent USE mod_numerique, ONLY: arb_rom, bsp1ddn, bsp1dn, 1 delete_doubles, entre, inside, noedif, newspl, newton, 2 no_croiss, sum_n, zoning USE mod_static, ONLY: d_zc, mmtI, resout, thermo USE mod_variables, ONLY: age, bp, bp_t, chim, chim_gram, ctel, ctem, 1 ctep, cter, ctet, c_iben, dim_ch, dim_qs, dim_rot, dv_tr, inter, jlim, 2 knot, knotc, knotr, knot_t, knot_tds, lconv, lhe_stop, lim, lr_stop, 3 lt_stop, lx_stop, mc, mct, mc_fixe, model_num, mstar, mrot, mrott, mw_tot, 4 m_zc, m_stat, m_stat_t, nb_modeles, nc_fixe, 5 n_ch, n_ch_t, n_qs, n_qs_t, n_rot, n_tds, omega_jpz, q, 6 qt, qt_t, q_t, rota, rota_t, rstar, r_ov, r_zc, 7 r_stat, r_stat_t, sortie, tds, vth, wrot, x_tds, xt_tds IMPLICIT NONE REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: compchim, 1 compg, dcompchim, dcompg, ener, epsilon, fd, ioni, lnta REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: b, bb, jac, 1 new_esp, var REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: alfa, 1 alfa_atm, anube7, anub8, anun13, anuo15, anupep, anupp, 2 beta, beta_atm, delta, delta_atm, cp, cp_atm, dcapdr, 3 dcapdt, depsdr, depsdt, depsx, dlpp, dmu_x, drdt, dxchim, 4 dxchimg, gam, gamma, gamma_atm, glob, grad, gradad, gradconv, 5 gradrad, grad_mj,grad_mj_a, grad_atm, grad_mu, grada_atm, 6 gradc_atm, gradr_atm, gtsg, hp, hp_atm, kap, k_atm, l, 7 degene, ldcapdr_a, ldcapdt_a, lnpt_a, lnpt_at, lnt_a, m, mu, 8 mu_atm, mue, mue_atm, p, pt, r, ro, 9 ro_atm, t, u, vaissala, vais_atm, w, xchim, 1 xchimg, xchim1, xchim1g, z, zbar REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: comp, dcomp, esp, 1 ex, mt, qb, qbt, r_son, v_son REAL (kind=dp), DIMENSION(nvth) :: dfvth, fvth REAL (kind=dp), DIMENSION(7) :: dfqs, fqs REAL (kind=dp), DIMENSION(8) :: dfrot, frot REAL (kind=dp), DIMENSION(0:5) :: absc, fonc, poly REAL (kind=dp), DIMENSION(5) :: epsilo REAL (kind=dp), DIMENSION(4) :: dfen, fen REAL (kind=dp), DIMENSION(1) :: dfatm, dftds, fatm, ftds REAL (kind=dp), SAVE :: agep=1.d-30, cte2, cte3, cte4, 1 cte5, cte6, cte7, cte8, dt, dtp=-1.d2, dts=0.d0, lteffp=-100.d0 REAL (kind=dp) :: a, age_car, alphap, bid, be7, b8, 1 dcpp, dcpt, dcpx, deltap, deltat, deltax, depsp, depst, 2 df_tau, dgaml, dgamlpp, dgamm, dgamp, dgampt, dgamr, 3 dgamrs, dgamt, dgamtau,dgamx, dgradadp, dgradadt, dgradadx, 4 dgradl, dgradlpp, dgradm, dgradp, dgradpt, dgradr, dgradrs, 5 dgradt, dgradtau, dgradx, dhpm, dhpp, dhppt, dhpt, dhpr, 6 dhpx, dkapp, dkapt, dkapx, dml, dmr, dnu02, dnu13, dpl, 7 dpr, dptdl, dptdr, dro_grav, drop, drot, dro_teff, drox, 8 dtdl, dtdr, dtsdg, dtsdteff, dtsdtau, dup, dut, dux,d2f_tau, 9 d2p, d2ro, ero, et, e_cno, e_gr, e_pp, e_nuc, e_3al, f_tau, 1 f17, gam_atm, grav, gtmax=-1.d0, hh, lext, log_l, lp1, lp2, 2 lteff, lb_p0, mext, mgtmax, msh, mshe, mtotp, mtot_ascii, nel, 3 nu0, n13, period, pertmm, pext, ptext, o15, rext, 4 rgtmax, rp1, rp2, ro_ext, teff, text, tp1, tp2, t_max, 5 w_max, w_rotp, x0_ascii, y0_ascii, zsx INTEGER, PARAMETER, DIMENSION(2) :: Ktest=(/ 3, 4 /) INTEGER, DIMENSION(8) :: values INTEGER, PARAMETER :: ntot_ecrit=50 INTEGER, DIMENSION(1) :: igtmax INTEGER, SAVE :: l0, n_ecrit=0, nb_der_out=0 INTEGER :: dim_qsp, i, iglob, ihe4_ascii, ipn, itot, ivar, 1 Krotp, Krot_ascii, i_cno, i_gr, i_pp, i_3a, 2 j, k, knotb, knote, lq, m_chp, m_qsp, m_rotp, n, 3 nadd, nchimp, nchim_ascii, nep, nrotp, n_init, n_tot, ord_qsp, un23 LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:) :: convec LOGICAL, SAVE :: coca=.FALSE., coca_e=.FALSE., coca_pr=.FALSE., 1 cohe=.FALSE., cohe_e=.FALSE., cohe_pr=.FALSE., 2 cone=.FALSE., cone_e=.FALSE., cone_pr=.FALSE., 3 coox=.FALSE., coox_e=.FALSE., coox_pr=.FALSE., 4 cosi=.FALSE., cosi_e=.FALSE., cosi_pr=.FALSE., 5 der, llist=.FALSE., list_pms=.FALSE., list_rep=.FALSE., 6 list_sort=.FALSE., list_zams=.FALSE., passe=.FALSE., 7 post=.FALSE., post_e=.FALSE., post_pr=.FALSE., 8 zams=.FALSE., zams_e=.FALSE., zams_pr=.FALSE. LOGICAL :: diffusionp, ecrHR, ecrHHe, ecrit, ecr_inter, 1 en_m23p, lim_rop, ok, puls, rad_atm, rot_solidp, sort CHARACTER (len=1) :: oui CHARACTER (len=2) :: precisionp, precisione CHARACTER (len=3) :: text1, text2 CHARACTER (len=4), ALLOCATABLE, DIMENSION(:) :: nom_elemp, 1 nom_elem_ascii CHARACTER (len=4), SAVE :: nom_vwrot='Wrot' CHARACTER (len=4) :: suffix CHARACTER (len=5) :: number, zone CHARACTER (len=8) :: date CHARACTER (len=9), PARAMETER, DIMENSION(12) :: mois=(/ 1 'Janvier ','Février ','Mars ','Avril ', 2 'Mai ','Juin ','Juillet ','Aout ', 3 'Septembre','Octobre ','Novembre ','Décembre ' /) CHARACTER (len=9), PARAMETER, DIMENSION(12) :: month=(/ 1 'January ','February ','March ','April ', 2 'May ','June ','July ','August ', 3 'September','October ','November ','December ' /) CHARACTER (len=10), ALLOCATABLE, SAVE, DIMENSION(:,:) :: tx_ioni CHARACTER (len=10) :: time CHARACTER (len=20), SAVE :: nom_save CHARACTER (len=20) :: nom_ctesp, nom_pertmp, nom_pertwp, 1 nom_tdetaup, nom_atmp, nom_convp, nom_nucp, nom_nuc_cplp, 2 nom_diffmp, nom_difftp, nom_diffwp, nom_etatp CHARACTER (len=30) :: instb CHARACTER (len=44), DIMENSION(10) :: list_TR CHARACTER (len=50), DIMENSION(8) :: f_eosp, f_opap CHARACTER (len=50) :: chaine, nom_opap, nom_fich1 CHARACTER (len=60) :: modelb CHARACTER (len=68) :: list_LB CHARACTER (len=80) :: chain, chain_atm, chain_don, chain_rep, 1 titre CHARACTER (len=440) :: list_TRF NAMELIST/nl_rlg/m_qs,m_ch,m_rot,m_tds,m_ptm,ordre,precic,precix, 1 precit,psi0,loc_zc,dtmax,dtmin,dt0,d_lum, 2 d_press,d_ray,d_temp,age_car,ini0,n_atm,kipp,en_m23,ctel,ctep, 3 ctem,cter,ctet,dn_fixe,dpsi,mu_saha,ajuste,lisse,npt_lisse, 4 q0,l0,new_bv,fcv,des_instb,no_discon,n_min_ZC,l_fac,d_zc NAMELIST/nl_langue/langue NAMELIST/nl_blabla/baratine c-------------------------------------------------------------------------- 2000 FORMAT(8es10.3) 2010 FORMAT(5e15.8) c recherche d'un éventuel fichier langue dans l'environnement pour fixer c une langue différente du français INQUIRE(file='langue',exist=ok) IF(ok)THEN OPEN(unit=30,form='formatted',status='old',file='langue') READ(30,nl_langue) ; CLOSE(unit=30) ELSE WRITE(*,30) 30 FORMAT(/,'--------------------',/, 1 'CESAM speaks a bit of english if you include in the',/, 2 'working directory a file named ''langue'' with the',/, 3 'statements : &NL_LANGUE',/,"langue='english'",/,'/',/, 4 'cf. aide_mem2k, chapter Personnalisation',/, 5 '-------------------') langue='francais' ENDIF c langue='english' c-------------menu d'entrée (boucle B1) début------------------------------- B1: DO SELECT CASE(langue) CASE('english') WRITE(*,29) 29 FORMAT(/,'To stop : enter 0 then click on RETURN',/, 1 'To continue the evolution of a stellar model further on,',/, 2 'enter 1 then click on RETURN',/, 3 'To start the evolution with a homogeneous ZAMS model,',/, 4 'enter 2 then click on RETURN',/, 5 'To start the evolution with a homogeneous PMS model,',/, 6 'enter 3 then click on RETURN') CASE DEFAULT WRITE(*,28) 28 FORMAT(/,'Pour arrêter : taper 0 puis RETURN',/, 1 'Pour poursuivre une évolution : taper 1 puis RETURN',/, 2 'Pour initialiser un modèle de ZAMS : taper 2 puis RETURN',/, 3 'Pour initialiser un modèle de PMS : taper 3 puis RETURN') END SELECT c on entre un23 = 1, 2 ou 3 comme indiqué ci dessus c puis, suivant les cas on donné à un23 les valeurs suivantes: c un23 : 1 poursuite d'une évolution, modèle repris en binaire c 2 modèle initial en binaire, ZAMS c 3 modèle initial en binaire, PMS c -1 poursuite d'une évolution, modèle repris en ASCII c -2 modèle initial en ASCII, ZAMS c -3 modèle initial en ASCII, PMS c pms : modèle de PMS c zams : modèle de ZAMS c post : modèle après la ZAMS c cohe : modèle avec combustion de l'helium c coca : modèle avec combustion du carbone c cone : modèle avec combustion du néon c coox : modèle avec combustion de l'oxygène c cosi : modèle avec combustion du silicium READ*,un23 c gestion des modèles repris SELECT CASE(un23) CASE(0) SELECT CASE(langue) CASE('english') PRINT*,'STOP' CASE DEFAULT PRINT*,'ARRET' END SELECT STOP c------------ poursuite d'un modèle--------------------------- c rep_atm=.TRUE.: si existe, reprise atm. restituée en binaire CASE(1) SELECT CASE(langue) CASE('english') WRITE(*,1031) 1031 FORMAT('Is the model which evolution is to be continued',/, 1'stored in a binary file? enter o for yes; enter n for no') CASE DEFAULT WRITE(*,31) 31 FORMAT('Le modèle à poursuivre est-il donné en binaire ? o/n') END SELECT c on doit entrer o ou n B4: DO READ(*,'(a)')oui SELECT CASE(oui) CASE('n', 'o') EXIT B4 CASE DEFAULT SELECT CASE(langue) CASE('english') PRINT*,'enter o or n, you have enter : ',oui CASE DEFAULT PRINT*,'entrer o ou n, vous avez entré : : ',oui END SELECT END SELECT ENDDO B4 IF(oui /= 'n')THEN SELECT CASE(langue) CASE('english') WRITE(*,1033) 1033 FORMAT('Enter the name of the binary file for the model',/, 1 'which evolution is to be continued.',/, 2 'For instance: my_model_B.rep, my_model_B.hom, my_model_B.pms') CASE DEFAULT WRITE(*,33) 33 FORMAT('nom du fichier binaire à poursuivre',/, 1 'Ex.: mon_modele_B.rep, mon_modele_B.hom, mon_modele_B.pms') END SELECT READ*,nom_fich1 INQUIRE(file=TRIM(nom_fich1),exist=ok) IF(ok)THEN SELECT CASE(langue) CASE('english') WRITE(*,1035)TRIM(nom_fich1) 1035 FORMAT('CESAM will use the model : ',a) CASE DEFAULT WRITE(*,35)TRIM(nom_fich1) 35 FORMAT('CESAM utilisera le modèle : ',a) END SELECT zams=.FALSE. ; post=.FALSE. ; cohe=.FALSE. ; coca=.FALSE. cone=.FALSE. ; coox=.FALSE. ; cosi=.FALSE.; rep_atm=.TRUE. EXIT B1 ELSE SELECT CASE(langue) CASE('english') WRITE(*,1036)TRIM(nom_fich1) 1036 FORMAT('This input file is unknown : ',a,/,'STOP') CASE DEFAULT WRITE(*,36)TRIM(nom_fich1) 36 FORMAT('Fichier d''entrée inconnu : ',a/,'ARRET') END SELECT CALL sortie('cesam 1') ENDIF !ok ELSE SELECT CASE(langue) CASE('english') WRITE(*,1081) 1081 FORMAT('Name of the ASCII file for the model which ',/, 1 ' evolution is to be continued; for instance: my_model.ascii') CASE DEFAULT WRITE(*,81) 81 FORMAT('Nom du fichier ASCII du modèle à poursuivre',/, 1 ' Exemple: mon_modele.ascii') END SELECT READ*,nom_fich1 INQUIRE(file=TRIM(nom_fich1),exist=ok) IF(ok)THEN SELECT CASE(langue) CASE('english') WRITE(*,1082)TRIM(nom_fich1) 1082 FORMAT('CESAM use the model : ',a) CASE DEFAULT WRITE(*,82)TRIM(nom_fich1) 82 FORMAT('CESAM utilise le modèle : ',a) END SELECT rep_atm=.FALSE. zams=.FALSE. ; post=.FALSE. ; coox=.FALSE. ; cosi=.FALSE. cone=.FALSE. ; cohe=.FALSE. ; coca=.FALSE. ; un23=-1 ; EXIT B1 ELSE SELECT CASE(langue) CASE('english') WRITE(*,1085)TRIM(nom_fich1) 1085 FORMAT('STOP, this ASCII file is unknown : ',a) CASE DEFAULT WRITE(*,85)TRIM(nom_fich1) 85 FORMAT('ARRET, fichier ASCII inconnu : ',a) END SELECT CALL sortie('cesam 2') ENDIF ENDIF c-------modèle initial de ZAMS--------------------------- CASE(2) zams=.TRUE. ; post=.FALSE. ; cohe=.FALSE. ; cosi=.FALSE. cone=.FALSE. ; coca=.FALSE. ; coox=.FALSE. ; rep_atm=.FALSE. SELECT CASE(langue) CASE('english') WRITE(*,1084) 1084 FORMAT('Is the input ZAMS model stored in a binary file?',/, 1 'enter o for yes; n for no') CASE DEFAULT WRITE(*,84) 84 FORMAT('Le modèle initial de ZAMS est-il en binaire ? o/n') END SELECT READ*,oui IF(oui /= 'n')THEN SELECT CASE(langue) CASE('english') WRITE(*,1086) 1086 FORMAT('Enter the name of the binary file for the initial',/, 1 'model. For instance : my_model_B.hom') CASE DEFAULT WRITE(*,86) 86 FORMAT('Entrer le nom du fichier binaire du modèle initial',/, 1 'Exemples: mon_modele_B.hom') END SELECT READ*,nom_fich1 ; INQUIRE(file=TRIM(nom_fich1),exist=ok) IF(ok)THEN SELECT CASE(langue) CASE('english') WRITE(*,1087)TRIM(nom_fich1) 1087 FORMAT('CESAM will use the following binary file : ',a) CASE DEFAULT WRITE(*,87)TRIM(nom_fich1) 87 FORMAT('CESAM utilise le modèle binaire : ',a) END SELECT un23=2 !fichier en binaire (intruction redondante) EXIT B1 ELSE SELECT CASE(langue) CASE('english') WRITE(*,1088)TRIM(nom_fich1) 1088 FORMAT('STOP unknown input file : ',a) CASE DEFAULT WRITE(*,88)TRIM(nom_fich1) 88 FORMAT('ARRET fichier d''entrée inconnu : ',a) END SELECT CALL sortie('cesam 3') ENDIF !ok ELSE !oui SELECT CASE(langue) CASE('english') WRITE(*,1089) 1089 FORMAT('Enter the ASCII file name of the initial model',/, 1 'Examples: m010.zams, m020.zams, m050.zams, m150.zams') CASE DEFAULT WRITE(*,89) 89 FORMAT('Entrer le nom du fichier ASCII du modèle initial',/, 1 'Exemples: m010.zams, m020.zams, m050.zams, m150.zams') END SELECT READ*,nom_fich1 INQUIRE(file=TRIM(nom_fich1),exist=ok) IF(ok)THEN SELECT CASE(langue) CASE('english') WRITE(*,1090)TRIM(nom_fich1) 1090 FORMAT('CESAM will use the model : ',a) CASE DEFAULT WRITE(*,90)TRIM(nom_fich1) 90 FORMAT('CESAM utilisera le modèle : ',a) END SELECT PRINT*; un23=-2 !fichier en ASCII EXIT B1 ELSE SELECT CASE(langue) CASE('english') WRITE(*,1091)TRIM(nom_fich1) 1091 FORMAT('STOP, unknown input file : ',a) CASE DEFAULT WRITE(*,91)TRIM(nom_fich1) 91 FORMAT('ARRET, fichier d''entrée inconnu : ',a) END SELECT CALL sortie('cesam 4') ENDIF !ok ENDIF !oui c------------- modèle initial de PMS------------------ CASE(3) zams=.FALSE. ; post=.FALSE.; cohe=.FALSE. ; cosi=.FALSE. cone=.FALSE. ; coca=.FALSE. ; coox=.FALSE. ; rep_atm=.FALSE. SELECT CASE(langue) CASE('english') WRITE(*,1092) 1092 FORMAT('Is the initial PMS model stored in a binary file?',/, 1 'enter o for yes, enter n for no') CASE DEFAULT WRITE(*,92) 92 FORMAT('Le modèle initial de PMS est-il donné en binaire ? o/n') END SELECT READ*,oui IF(oui == 'o')THEN SELECT CASE(langue) CASE('english') WRITE(*,1093) 1093 FORMAT('Enter the name of the binary file for the inital',/, 1 'model; for instance : my_model_B.pms') CASE DEFAULT WRITE(*,93) 93 FORMAT('Entrer le nom du fichier binaire du modèle initial',/, 1 'Exemples: mon_modele_B.pms') END SELECT READ*,nom_fich1 ; INQUIRE(file=TRIM(nom_fich1),exist=ok) IF(ok)THEN SELECT CASE(langue) CASE('english') WRITE(*,1094)TRIM(nom_fich1) 1094 FORMAT('CESAM uses the initial PMS model : ',a) CASE DEFAULT WRITE(*,94)TRIM(nom_fich1) 94 FORMAT('CESAM utilise le modèle initial de PMS : ',a) END SELECT PRINT* un23=3 !fichier en binaire et, pour la PMS, pas de modèle homo. EXIT B1 ELSE SELECT CASE(langue) CASE('english') WRITE(*,1095)TRIM(nom_fich1) 1095 FORMAT('STOP, unknown input file : ',a) CASE DEFAULT WRITE(*,95)TRIM(nom_fich1) 95 FORMAT('ARRET, fichier d''entrée inconnu : ',a) END SELECT CALL sortie('cesam 5') ENDIF !ok ELSE !oui SELECT CASE(langue) CASE('english') WRITE(*,1096) 1096 FORMAT('CESAM uses an input PMS model in an ASCII file',/, 1 'enter the name of the ASCII file for the initial PMS model.', 2 /,'For instance: 2d-2.pms(Tc=0.2MK), 5d-4.pms(Tc=0.5MK)',/, 3 '8d-5.pms(Tc=1.0MK), 4d-2.pms (for M*>10Msol)') CASE DEFAULT WRITE(*,96) 96 FORMAT('CESAM utilise un modèle initial de PMS en ASCII',/, 1 'entrer le nom du fichier ASCII du modèle PMS initial',/, 2 'Exemples: 2d-2.pms(Tc=0.2MK), 5d-4.pms(Tc=0.5MK)',/, 3 ' 8d-5.pms(Tc=1.0MK), 4d-2.pms (pour M*>10Msol)') END SELECT READ*,nom_fich1; INQUIRE(file=TRIM(nom_fich1),exist=ok) IF(ok)THEN SELECT CASE(langue) CASE('english') WRITE(*,1097)TRIM(nom_fich1) 1097 FORMAT('CESAM uses the model in the ASCII file : ',a) CASE DEFAULT WRITE(*,97)TRIM(nom_fich1) 97 FORMAT('CESAM utilise le modèle ASCII : ',a) END SELECT un23=-3 !fichier en ASCII et, pour la PMS, pas de modèle homo. EXIT B1 ELSE SELECT CASE(langue) CASE('english') WRITE(*,1098)TRIM(nom_fich1) 1098 FORMAT('STOP, unknown input file : ',a) CASE DEFAULT WRITE(*,98)TRIM(nom_fich1) 98 FORMAT('ARRET, fichier d''entrée inconnu : ',a) END SELECT CALL sortie('cesam 6') ENDIF !ok ENDIF !oui SELECT CASE(langue) CASE('english') WRITE(*,1099) 1099 FORMAT('Error: enter 0, 1, 2 or 3, it was entered : ',i4) CASE DEFAULT WRITE(*,99)un23 99 FORMAT('ERREUR, entrer 0, 1, 2 ou 3, vous avez tapé : ',i4) END SELECT CYCLE B1 END SELECT !menu d'entrée ENDDO B1 c--------menu d'entrées fin, écriture des conditions du calcul début --- c identificateur des fichiers des données, du point de reprise, c des résultats, du modèle d'âge 0, des oscillations SELECT CASE(langue) CASE('english') WRITE(*,1100) 1100 FORMAT('enter the generic name of the model : ',/, 1 'for instance : my_model') CASE DEFAULT WRITE(*,100) 100 FORMAT('entrer l''identificateur du modèle : ',/, 1 'Exemple : mon_modele') END SELECT READ*,nom_fich2 SELECT CASE(langue) CASE('english') WRITE(*,1103)TRIM(nom_fich2) 1103 FORMAT('generic name of the model files : ',a) CASE DEFAULT WRITE(*,103)TRIM(nom_fich2) 103 FORMAT('identificateur des fichiers du modèle : ',a) END SELECT c ouverture du fichier du listing mon_modele.lis OPEN(unit=2,form='formatted',status='unknown', 1 file=TRIM(nom_fich2)//'.lis') !fichiers du listing SELECT CASE(langue) CASE('english') WRITE(2,1001)version ; WRITE(*,1001)version 1001 FORMAT(/,t10, 1'***************************************************************', 2 //,t11, 3 'MODEL OF INTERNAL STRUCTURE computed with CESAM2k version : ',a, 4 //,t10, 5'***************************************************************', 6/) CASE DEFAULT WRITE(2,1)version ; WRITE(*,1)version 1 FORMAT(/,t10, 1 '***************************************************************', 2 //,t11, 3 'MODELE DE STRUCTURE INTERNE calculé par CESAM2k version : ',a, 4 //,t10, 5 '***************************************************************', 6/) END SELECT CALL date_and_time(date,time,zone,values) SELECT CASE(langue) CASE('english') titre='The calculation started on the : '//date(7:8)//' ' 1 //TRIM(month(values(2)))//' '//date(1:4)//' at ' 2 //time(1:2)//'h'//time(3:4) CASE DEFAULT titre='Début du calcul le : '//date(7:8)//' ' 1 //TRIM(mois(values(2)))//' '//date(1:4)//' à ' 2 //time(1:2)//'h'//time(3:4) END SELECT WRITE(*,3)TRIM(titre) ; WRITE(2,3)TRIM(titre) 3 FORMAT(t10,a,/) c lecture du fichier des données mon_modele.don CALL lit_nl(wrot) c recherche d'un éventuel fichier blabla dans l'environnement qui permet de c détourner une partie des sorties sur le moniteur vers des fichiers ASCII c cf. Aide Mémoire §3.13 Limitation des sorties INQUIRE(file='blabla',exist=ok) IF(ok)THEN OPEN(unit=30,form='formatted',status='old',file='blabla') READ(30,nl_blabla) ; CLOSE(unit=30) IF(.NOT.baratine)THEN SELECT CASE(langue) CASE('english') WRITE(*,1301)TRIM(nom_fich2)//'_static',TRIM(nom_fich2)//'_atm', 1 TRIM(nom_fich2)//'_evol' 1301 FORMAT(/,'Part of the on line outputs rerouted towards files :', 1 /,a,/,a,/,a,/) CASE DEFAULT WRITE(*,301)TRIM(nom_fich2)//'_static',TRIM(nom_fich2)//'_atm', 1 TRIM(nom_fich2)//'_evol' 301 FORMAT(/,'Détournement de sorties on line vers les fichiers :',/, 1 a,/,a,/,a,/) END SELECT ENDIF ENDIF c écriture des constantes initialisées dans lit_nl CALL print_ctes(6) ; CALL print_ctes(2) SELECT CASE(langue) CASE('english') WRITE(*,1005) ; WRITE(2,1005) 1005 FORMAT('Input physics :',/) CASE DEFAULT WRITE(*,5) ; WRITE(2,5) 5 FORMAT(t14,'Physique utilisée :',/) END SELECT WRITE(*,2)nom_etat, nom_opa, nom_conv, nom_nuc, nom_nuc_cpl, 1 nom_abon, nom_atm, nom_tdetau, nom_pertm, nom_pertw, nom_diffm, 2 nom_difft, nom_diffw, nom_ctes WRITE(2,2)nom_etat, nom_opa, nom_conv, nom_nuc, nom_nuc_cpl, 1 nom_abon, nom_atm, nom_tdetau, nom_pertm, nom_pertw, nom_diffm, 2 nom_difft, nom_diffw, nom_ctes 2 FORMAT('EOS : ',a,t30,'Opa. : ',a,t60,'conv. : ',a,/, 1 'réac. nuc. : ',a,t30,'compil. nuc. : ',a,t60,'abon_ini. : ',a,/, 2 'atmos.: ',a,t30,'T(tau): ',a,t60,'perte m. : ',a,/, 3 'perte w. : ',a,t30,'diff. micro. : ',a,t60,'diff. turb. : ',a,/ 4 'diff. rot. : ',a,t30,'ctes. : ',a,/) c écriture du fichier mon_modele.mHHe (limites des shells H et He burning) ecrHHe=INDEX(nom_des,'HHe') /= 0 c --------fin des écritures d'entrée, initialisations, réglages c precision(char*5) est lu dans modele.don c precisonp: précision du modèle relu, precisione pour écriture c precisionp(1:2)=precisione(1:2)=precision(1:2) c suffix(1:4)=precision(3:6) pour les extensions A(ajuste), G(TdS), L(lisse), M(n_max), c N(nuc), O(osc). cf. fichier de la SOURCE "arguments du point.don.f" precisione(1:2)=precision(1:2) ; precisionp=precisione suffix(1:4)=precision(3:6) c quelques constantes physiques cte2=mtot**2*gmsol/lsol/rsol*msol*2.d0/secon6 cte3=4.d0*pi*rsol**3/msol cte4=2.d0/3.d0*rsol**3/g/msol cte5=-4.d0*pi*g/3.d0*rsol**2 cte6=0.5d0/g/pi cte7=rsol*1.d-5/secon6 !Rsol/My --> km/s cte8=SQRT(SQRT(lsol/(4.d0*pi*rsol**2*sigma))) !pour Teff c initialisation des constantes de réglage c correspondant à la précision par défaut : pr "réaliste" c elles seront adaptées pour chaque type de précision c la constante de répartition psi0 est fixée de façon à avoir, d'une couche c à la suivante, des variations de P, T de ~0.1dex, et M~0.05, le nombre de c couches < N_MAX, s'ajustant pour respecter ces critères c avec N_MAX < 0 le dernier modèle sera calculé avec N_MAX couches m_qs=2 !ordre des splines-1 pour les variables quasi-stat. m_ch=2 !ordre des splines pour interpolation de la comp.chim. IF(diffusion)m_ch=3 m_rot=2 !ordre des splines-1 pour la rotation m_tds=2 !ordre des splines pour interpolation du TdS m_ptm=2 !ordre des spl. pour int. de la masse (perte de masse) ordre=4 !ordre du schéma d'int. des réac.nuc. avec rk_imps precic=1.d-4!précision/gramme sur iters; NR évol. comp. chim. precix=1.d-3!précision sur les iters NR équilibre quasi-stat. precit=0.1d0!var max/masse pour évolution temp. de comp. chim. c IF(diffusion)precit=0.03d0 loc_zc=1.d-3 !précision de la localisation des limites ZR/ZC dtmin=365.d0*86400.d0/secon6 !pas temporel minimal 1an IF(mtot <= 0.95d0)THEN !pas temporels init pour les évolutions dt0=10.d0 !à partir d'un modèle de ZAMS dtmax=200.d0 !pas temporel maximum ELSEIF(mtot <= 1.25d0)THEN dt0=5.d0 dtmax=100.d0 !pas temporel maximum ELSEIF(mtot <= 1.75d0)THEN dt0=1.d0 dtmax=50.d0 !pas temporel maximum ELSEIF(mtot <= 2.5d0)THEN dt0=1.0d0 dtmax=10.0d0 !pas temporel maximum ELSEIF(mtot <= 3.d0)THEN dt0=1.0d0 dtmax=1.0d0 !pas temporel maximum ELSEIF(mtot <= 5.d0)THEN dt0=0.075d0 dtmax=0.75d0 !pas temporel maximum ELSEIF(mtot <= 8.d0)THEN dt0=0.05d0 dtmax=0.5d0 !pas temporel maximum ELSEIF(mtot <= 15.d0)THEN dt0=0.01d0 dtmax=0.1d0 !pas temporel maximum ELSE dt0=0.001d0 dtmax=0.03d0 !pas temporel maximum ENDIF d_lum=0.1d0 !variation relative max de Lext sur un pas temporel d_press=0.1d0 !variation max de ln P sur un pas temporel d_ray=0.1d0 !variation relative max de Rext sur un pas temporel d_temp=0.1d0 !variation max de ln T sur un pas temporel age_car=dt0*1.d1!âge auquel les cond.ini. sont supp. relaxées ini0=3 !nb. iter. NR avec réestim. de lim. ZR/ZC & evol n_atm=75 !nb. couches pour l'atmosphère restituée kipp=.TRUE. !dS = cp T ( dlnT - grad_ad dlnP ) / dt dn_fixe=0.05d0 !taux pour modif de grille fixe en comp.chim dpsi=0.05d0 !var. rela. max de psi pour modification de n_qs dpsip=0.03d0 !var. rela. > 0 max de psi pour modification de n_qs dpsim=-0.07d0 !var. rela. < 0 max de psi pour modification de n_qs mu_saha=.TRUE. !utilisation de saha pour le calcul de mu ajuste=INDEX(suffix,"A") /= 0 !suffixe de la NAMELIST precision lisse=INDEX(suffix,"L") /= 0 !lissage des discon. de comp_chim npt_lisse=10 !si lisse lissage sur 2*npt+1 points q0=0.05d0!un point à q0 du centre dans fichiers ASCII de sortie l0=0 !nb de points ajoutés, en sortie, au voisinage des lim. ZR/ZC fcv=.FALSE. !si besoin, on forcera la convergence des_instb=.FALSE.! dessin de la zone d'instabilité n_max=MIN(n_max,5000) !nb. max de couches on ne dépassera pas 5000 n_init=500 !nombre de couches du modèle initial cephe=.FALSE. !sauf pour les cepheides typ=3 !zone d'instabilité de Turner(2001) n_min_ZC=5 !nombre de couches MIN dans ZC no_discon=.TRUE. !sans discontinuités de comp.chim. d_zc=100 !sans ajustement de la grille au voisinage des limites ZT/ZC l_fac=.FALSE. !sans placement d'un point de grille sur chaque limite ZR/ZC c valeurs des coefficients et de la constante de répartition ctel=0.d0 !constante de répartition pour la luminosité ctep=-1.d0 !constante de répartition pour la pression ctem=15.d0 !constante de répartition pour la masse cter=0.d0 !constante de répartition pour le rayon ctet=-1.d0 !constante de répartition pour la température psi0=0.d0 !valeur provisoire c personnalisation des paramètres numériques suivant le type de précision SELECT CASE(langue) CASE('english') WRITE(*,1218) ; WRITE(2,1218) 1218 FORMAT('Numerical parameters which are used : ') CASE DEFAULT WRITE(*,218) ; WRITE(2,218) 218 FORMAT('Paramètres numériques utilisés : ') END SELECT c les cas de précision SELECT CASE(precisione) c pour stades avancés, cépheides, flash de l'hélium CASE('AV', 'CE', 'He') q0=0.d0 !point à q0 du centre dans fichiers ASCII sortie dtmin=3600.d0/secon6 !pas temporel minimal 1heure psi0=0.03d0 c utilisation de r, l, m^1/3 pour les modèles avancés en_m23=.FALSE. SELECT CASE(langue) CASE('english') WRITE(*,1133) ; WRITE(2,1133) 1133 FORMAT('Cepheid or advanced models computed with r, l, m^1/3') CASE DEFAULT WRITE(*,133) ; WRITE(2,133) 133 FORMAT('Céphéïdes ou stades avancés calculés avec r, l, m^1/3') END SELECT c Céphéides on écrit tous les *.osc de la bande d'instabilité cephe=precisione == 'CE' all_cephe=INDEX(suffix,'O') /= 0 IF(cephe)des_instb=.TRUE. c pour le flash de l'hélium pour les masses <= ~2Msol IF(precisione == 'He')THEN dtmin=1.d0/secon6 !pas temporel minimal 1 seconde d_lum=0.3d0 npt_lisse=20 !pour lisser utiliser le suffixe L (cf ligne 810) ctet=0.d0 ENDIF c précision CoRoT CASE('co' , 'CO') m_ch=3 precix=1.d-4 precit=0.05d0 !var max/masse pour évol. temp. de la comp. chim. loc_zc=1.d-4 dtmax=MIN(50.d0,dtmax) !pas temporel maximum n_atm=100 q0=0.01d0 l0=5 no_discon=.FALSE. !on tient compte des discontinuités de comp.chim. l_fac=.TRUE. d_zc=2 psi0=2.d-2 kipp=.FALSE. ! TdS= dU + P dV c en r^2, l^2/3, m^2/3 ou r, l, m^1/3 en_m23= precisione == 'co' IF(en_m23)THEN SELECT CASE(langue) CASE('english') WRITE(*,1230) ; WRITE(2,1230) 1230 FORMAT('model with CoROT precision and r^2, l^2/3, m^2/3') CASE DEFAULT WRITE(*,230) ; WRITE(2,230) 230 FORMAT('modèle avec la précision CoROT et r^2, l^2/3, m^2/3') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1231) ; WRITE(2,1231) 1231 FORMAT('model with CoROT precision and r, l, m^1/3') CASE DEFAULT WRITE(*,231) ; WRITE(2,231) 231 FORMAT('modèle avec la précision CoROT et r, l, m^1/3') END SELECT ENDIF c pour petites masses très évoluées CASE('lm' , 'LM') dtmax=MIN(300.d0,dtmax) mu_saha=.FALSE. !calcul de mu simplifié fcv=.TRUE. !on forcera éventuellement la convergence q0=0.d0 c en r^2, l^2/3, m^2/3 ou r, l, m^1/3 en_m23= precisione == 'lm' IF(en_m23)THEN SELECT CASE(langue) CASE('english') WRITE(*,1306) ; WRITE(2,1306) 1306 FORMAT('highly evolved small mass with r^2, l^2/3, m^2/3') CASE DEFAULT WRITE(*,306) ; WRITE(2,306) 306 FORMAT('pour petites masses très evoluées avec r^2, l^2/3, m^2/3') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1307) ; WRITE(2,1307) 1307 FORMAT('highly evolved models with small masses with r, l, m^1/3') CASE DEFAULT WRITE(*,307) ; WRITE(2,307) 307 FORMAT('pour petites masses très evoluées avec r, l, m^1/3') END SELECT ENDIF c précision normale CASE('np', 'NP') m_qs=1 m_rot=1 ordre=1 precix=5.d-3 loc_zc=5.d-3 dtmax=MIN(300.d0,dtmax) n_atm=50 mu_saha=.FALSE. fcv=.TRUE. !on forcera éventuellement la convergence c en r^2, l^2/3, m^2/3 ou r, l, m^1/3 en_m23= precisione == 'np' IF(en_m23)THEN SELECT CASE(langue) CASE('english') WRITE(*,1216) ; WRITE(2,1216) 1216 FORMAT('model with a normal precision with r^2, l^2/3, m^2/3') CASE DEFAULT WRITE(*,216) ; WRITE(2,216) 216 FORMAT('modèle en précision normale avec r^2, l^2/3, m^2/3') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1224) ; WRITE(2,1224) 1224 FORMAT('model computed with r, l, m^1/3, standard precision') CASE DEFAULT WRITE(*,224) ; WRITE(2,224) 224 FORMAT('modèle calculé avec r, l, m^1/3, précision normale') END SELECT ENDIF c précision réaliste c'est le "défaut" CASE('pr' , 'PR') c en r^2, l^2/3, m^2/3 ou r, l, m^1/3 en_m23=precisione == 'pr' IF(en_m23)THEN SELECT CASE(langue) CASE('english') WRITE(*,1130) ; WRITE(2,1130) 1130 FORMAT('the model with r^2, l^2/3, m^2/3, realist precision') CASE DEFAULT WRITE(*,130) ; WRITE(2,130) 130 FORMAT('modèle avec r^2, l^2/3, m^2/3, précision réaliste') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1316) ; WRITE(2,1316) 1316 FORMAT('Model computed with r, l, m^1/3 and realistic precision') CASE DEFAULT WRITE(*,316) ; WRITE(2,316) 316 FORMAT('Modèle calculé avec r, l, m^1/3 et précision réaliste') END SELECT ENDIF c avec réglages externes (en_masse est défini dans les réglages) CASE('rg') chain=TRIM(nom_fich2)//'.rg' INQUIRE(file=TRIM(chain),exist=ok) IF(ok)THEN OPEN(unit=3,form='formatted',status='old',delim='apostrophe', 1 file=TRIM(chain)) ELSE SELECT CASE(langue) CASE('english') WRITE(*,1024)TRIM(chain) ; WRITE(2,1024)TRIM(chain) 1024 FORMAT('The parameter setting file : ',a,/, 1 'is unknown, research of a file named reglages') CASE DEFAULT WRITE(*,24)TRIM(chain) ; WRITE(2,24)TRIM(chain) 24 FORMAT('le fichier des réglages : ',a,/, 1 'non trouvé, recherche du fichier reglages') END SELECT chain='reglages' INQUIRE(file=TRIM(chain),exist=ok) IF(ok)THEN OPEN(unit=3,form='formatted',status='old',delim='apostrophe', 1 file=TRIM(chain)) ELSE SELECT CASE(langue) CASE('english') WRITE(*,1027)TRIM(chain) ; WRITE(2,1027)TRIM(chain) ; CALL sortie('cesam 7') 1027 FORMAT('STOP, The parameter setting file : ',a,' is unknown',/, 1 'in the directory where the calculations are performed') CASE DEFAULT WRITE(*,27)TRIM(chain) ; WRITE(2,27)TRIM(chain) ; CALL sortie('cesam 8') 27 FORMAT('ARRET, le fichier des réglages : ',a,' est inconnu',/, 1 'dans le directory où sont effectués les calculs') END SELECT ENDIF ENDIF SELECT CASE(langue) CASE('english') WRITE(*,1311)TRIM(chain) ; WRITE(2,1311)TRIM(chain) 1311 FORMAT('Model computed with the parameters set of the file:',a,/, 1 'if Error see the file SOURCE/arguments_des_reglages') CASE DEFAULT WRITE(*,311)TRIM(chain) ; WRITE(2,311)TRIM(chain) 311 FORMAT('Modèle utilisant les réglages externes du fichier : ',a,/, 1 'si Error voir le fichier SOURCE/arguments_des_reglages') END SELECT READ(3,nl_rlg) ; CLOSE(unit=3) WRITE(*,nl_rlg) ; WRITE(2,nl_rlg) WRITE(*,*) ; WRITE(2,*) c Solar accuracy pour le soleil actuel il n'y a pas de discontinuité CASE('sa' , 'SA' ) m_ch=3 precix=1.d-5 precit=0.02d0 loc_zc=1.d-5 dtmax=MIN(50.d0,dtmax) n_atm=MAX(100,n_atm) q0=0.01d0 l0=5 new_bv=.FALSE. !ancienne formulation de BV kipp=.FALSE. ! TdS= dU + P dV d_zc=2 psi0=0.02d0 c en r^2, l^2/3, m^2/3 ou r, l, m^1/3 en_m23=precisione =='sa' IF(en_m23)THEN SELECT CASE(langue) CASE('english') WRITE(*,1300) ; WRITE(2,1300) 1300 FORMAT('the model with solar precision and r^2, l^2/3, m^2/3') CASE DEFAULT WRITE(*,300) ; WRITE(2,300) 300 FORMAT('modèle avec la précision solaire er r^2, l^2/3, m^2/3') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1298) ; WRITE(2,1298) 1298 FORMAT('the model with solar precision and r, l, m^1/3') CASE DEFAULT WRITE(*,298) ; WRITE(2,298) 298 FORMAT('modèle avec la précision solaire et r, l, m^1/3') END SELECT ENDIF c super précision CASE('sp' , 'SP') m_ch=3 precix=1.d-4 precit=0.02d0 loc_zc=1.d-4 kipp=.FALSE. ! TdS= dU + P dV no_discon=.FALSE. !on tient compte des discontinuités de comp.chim. l_fac=.TRUE. ! placement des limites ZR/ZC à d_zc % d'un bord de couche d_zc=5 psi0=0.02d0 q0=0.01d0 c en r^2, l^2/3, m^2/3 ou r, l, m^1/3 en_m23=precisione =='sp' IF(en_m23)THEN SELECT CASE(langue) CASE('english') WRITE(*,1215) ; WRITE(2,1215) 1215 FORMAT('model with super precision and r^2, l^2/3, m^2/3 ') CASE DEFAULT WRITE(*,215) ; WRITE(2,215) 215 FORMAT('Modèle en super précision avec r^2, l^2/3, m^2/3 ') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1315) ; WRITE(2,1315) 1315 FORMAT('model computed with super precision and r, l, m^1/3') CASE DEFAULT WRITE(*,315) ; WRITE(2,315) 315 FORMAT('Modèle calculé en super précision avec r, l, m^1/3') END SELECT ENDIF c problème CASE DEFAULT SELECT CASE(langue) CASE('english') WRITE(*,1312)precisione ; WRITE(2,1312)precisione ; CALL sortie('cesam 9') 1312 FORMAT('STOP, this class of precision is unknown : ',a3, 1 ', defined classes :',/, 2 'AV, CE, He, CO, co, LM, lm, NP, np, PR, pr, rg, SA, sa, SP, sp') CASE DEFAULT WRITE(*,312)precisione ; WRITE(2,312)precisione ; CALL sortie('cesam 10') 312 FORMAT('ARRET, type de précision non prévu : ',a3,', types prévus :',/, 1 'AV, CE, He, CO, co, LM, lm, NP, np, PR, pr, rg, SA, sa, SP, sp') END SELECT END SELECT !type de calcul c PRINT*,l_fac,' ',precisione ; PAUSE'cesam2k l_fac' c------------aménagements divers début------------------------------ c avec diffusion l'intégration est globale precic=precix c l'équ. de diff. est d'ordre 2, la spline doit être 2 fois dérivable IF(diffusion .AND. precisione /= 'rg')THEN precic=precix m_ch=MAX(3,m_ch) ENDIF c avec rotation, on lisse les discontinuités SELECT CASE (Krot) CASE(3,4) PRINT*,'ARRET: diffusion du moment cinétique indisponible' CALL sortie('cesam 11') IF(precisione /= 'rg')THEN precix=1.d-3 loc_zc=1.d-3 mu_saha=.FALSE. ord_rot=m_rot+r_qs ENDIF CASE(5) ord_rot=m_rot+r_qs CASE(6) q0=0.d0 ; l0=0 END SELECT c PRINT*,Krot,dtmax ; PAUSE'dtmax' c avec lissage les discontinuités sont créées puis lissées, l_fac et no_discon sont F IF(lisse)THEN no_discon=.TRUE. l_fac=.FALSE. d_zc=100 ENDIF c si on ignore les discontinuités l_fac et d_zc sont sans objet IF(no_discon)THEN l_fac=.FALSE. d_zc=100 ENDIF c si psi0=0, psi0=1.% par constante de répartition non nulle (epsilo est Vt) IF(psi0 == 0.d0)THEN epsilo(1)=ctel ; epsilo(2)=ctem ; epsilo(3)=ctep ; epsilo(4)=cter epsilo(5)=ctet psi0=COUNT(epsilo /= 0.d0)*1.d-2 ENDIF c utilisation du nombre de couche maximal N_MAX IF(INDEX(suffix,"M") /= 0)THEN psi0=1.d-8 SELECT CASE(langue) CASE('english') WRITE(*,1310)n_max ; WRITE(2,1310)n_max 1310 FORMAT('the model is computed with the maximum shell number=', 1 i5) CASE DEFAULT WRITE(*,310)n_max ; WRITE(2,310)n_max 310 FORMAT('modèle calculé avec le nombre de couches maximal=',i5) END SELECT ENDIF c----------aménagements divers fin--------------------------- c ne : nombre d'inconnues du modèle quasi-statique c Ipg : indice de ln Pgaz avec pression turbulente (1 sans Pturb) c pturb=.TRUE. on tient compte de la pression turbulente c avec pression turbulente 7 inconnues c sans pression turbulente 6 inconnues, Ptot=Pgaz der=cpturb < 0.d0 !on tient compte de dln Pgaz/dln Ptot pturb=ABS(cpturb) /= 0.d0 IF(pturb)THEN PRINT*,'ARRET: pression turbulente indisponible' CALL sortie('cesam 12') ne=7 ; Ipg=7 ELSE ne=6 ; Ipg=1 ENDIF c identification des variables ALLOCATE(nom_qs(ne)) SELECT CASE(ne) CASE(6) IF(en_m23)THEN nom_qs=(/'ln Pgaz',' ln T ',' r^2 ',' l^2/3 ',' m^2/3 ', 1 ' Psi ' /) ELSE nom_qs=(/'ln Pgaz',' ln T ',' r ',' l ',' m^1/3 ', 1 ' Psi '/) ENDIF CASE(7) IF(en_m23)THEN nom_qs=(/ 'ln Ptot',' ln T ',' r^2 ',' l^2/3 ',' m^2/3 ', 1 ' Psi ', 'ln Pgaz' /) ELSE nom_qs=(/ 'ln Ptot',' ln T ',' r ',' l ',' m^1/3 ', 1 ' Psi ', 'ln Pgaz' /) ENDIF END SELECT c écritures SELECT CASE (Krot) CASE(3,4,5) WRITE(*,299)m_qs,m_ch,m_rot,ordre,precix,loc_zc,precit, 1 psi0,n_init,d_zc,ctel,ctem,ctep,cter,ctet,kipp,ajuste,no_discon, 2 lisse,new_bv,en_m23,l_fac,d_lum,d_press,d_ray,d_temp WRITE(2,299)m_qs,m_ch,m_rot,ordre,precix,loc_zc,precit, 1 psi0,n_init,ctel,ctem,ctep,cter,ctet,kipp,ajuste,no_discon, 2 lisse,new_bv,en_m23,l_fac,d_lum,d_press,d_ray,d_temp 299 FORMAT('m_qs=',i2,', m_ch=',i2,', m_rot=',i2, 1 ', ordre=',i2,', precix=',es8.1,', loc_zc=',es8.1,/, 2 'precit=',es8.1,', psi0=',es8.1,', n_init=',i3,', d_zc=',i3,/, 3 'ctel=',es10.3,', ctem=',es10.3,', ctep=',es10.3,', cter=',es10.3, 4 ', ctet=',es10.3,/,'kipp=',l2,', ajuste=',l2,', no_discon=',l2, 5 ', lisse=',l2,', new_bv=',l2,', en_m23=',l2,', l_fac=',l2,/, 6 'd_lum=',es10.3,', d_press=',es10.3,', d_ray=',es10.3, 7 ', d_temp=',es10.3,/) CASE DEFAULT WRITE(*,305)m_qs,m_ch,ordre,precix,loc_zc,precit,psi0,n_init, 1 d_zc,ctel,ctem,ctep,cter,ctet,kipp,ajuste,no_discon, 2 lisse,new_bv,en_m23,l_fac,d_lum,d_press,d_ray,d_temp WRITE(2,305)m_qs,m_ch,ordre,precix,loc_zc,precit,psi0,n_init, 1 d_zc,ctel,ctem,ctep,cter,ctet,kipp,ajuste,no_discon, 2 lisse,new_bv,en_m23,l_fac,d_lum,d_press,d_ray,d_temp 305 FORMAT('m_qs=',i2,', m_ch=',i2,', ordre=',i2,', precix=',es8.1, 1 ', loc_zc=',es8.1,/,'precit=',es8.1,', psi0=',es8.1,', n_init=',i3, 2 ', d_zc=',i3,/,'ctel=',es10.3,', ctem=',es10.3,', ctep=',es10.3, 3 ', cter=',es10.3,', ctet=',es10.3,/,'kipp=',l2,', ajuste=',l2, 4 ', no_discon=',l2,', lisse=',l2,', new_bv=',l2,', en_m23=',l2, 5 ', l_fac=',l2,/,'d_lum=',es10.3,', d_press=',es10.3,', d_ray=',es10.3, 6 ', dtemp=', es10.3,/) END SELECT IF(grille_fixe)THEN SELECT CASE(langue) CASE('english') WRITE(*,1222) ; WRITE(2,1222) 1222 FORMAT('Rigid grid in m^2/3 for the chemical composition') CASE DEFAULT WRITE(*,222) ; WRITE(2,222) 222 FORMAT('grille fixe en m^2/3 pour la composition chimique') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1223) ; WRITE(2,1223) 1223 FORMAT('Adaptative grid in m^2/3 for the chemical composition') CASE DEFAULT WRITE(*,223) ; WRITE(2,223) 223 FORMAT('grille variable en m^2/3 pour la comp. chim.') END SELECT ENDIF IF(no_discon)THEN WRITE(*,150) 150 FORMAT('on ignorera les discontinuités de composition chimique') ELSE WRITE(*,151) 151 FORMAT('on tiendra compte explicitement des discontinuités de comp.chim.') ENDIF IF(diffusion)THEN SELECT CASE(langue) CASE('english') WRITE(*,1240)m_ch ; WRITE(2,1240)m_ch 1240 FORMAT('Diffusion of the chemical comp. in radiative (ZR)', 1 /,'and convective zones (ZC), with m_ch=',i2) CASE DEFAULT WRITE(*,240)m_ch ; WRITE(2,240)m_ch 240 FORMAT('Diff. de la comp. chim. dans ZR et ZC, avec m_ch=',i2) END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1241) ; WRITE(2,1241) 1241 FORMAT('model without diffusion') CASE DEFAULT WRITE(*,241) ; WRITE(2,241) 241 FORMAT('modèle sans diffusion') END SELECT ENDIF IF(mdot == 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1242) ; WRITE(2,1242) 1242 FORMAT('model without mass loss') CASE DEFAULT WRITE(*,242) ; WRITE(2,242) 242 FORMAT('modèle sans perte de masse') END SELECT ELSEIF(mdot < 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1212)mdot ; WRITE(2,1212)mdot 1212 FORMAT('model with mass loss :',es10.3,'Msol/an') CASE DEFAULT WRITE(*,212)mdot ; WRITE(2,212)mdot 212 FORMAT('modèle avec perte de masse :',es10.3,'Msol/an') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1211)ABS(mdot) ; WRITE(2,1211)ABS(mdot) 1211 FORMAT('model with mass increase :',es10.3,'Msol/an') CASE DEFAULT WRITE(*,211)ABS(mdot) ; WRITE(2,211)ABS(mdot) 211 FORMAT('modèle avec gain de masse :',es10.3,'Msol/an') END SELECT ENDIF IF(l_pertmt)THEN SELECT CASE(langue) CASE('english') WRITE(*,1243) ; WRITE(2,1243) 1243 FORMAT('The model takes into account the mass los due to E=mc^2') CASE DEFAULT WRITE(*,243) ; WRITE(2,243) 243 FORMAT('Le modèle tient compte de la perte de masse due à E=mc^2') END SELECT ENDIF c type de rotation SELECT CASE(langue) CASE('english') WRITE(*,1017)thw(Krot) ; WRITE(2,1017)thw(Krot) 1017 FORMAT('Kind of rotation used : ',a) CASE DEFAULT WRITE(*,17)thw(Krot) ; WRITE(2,17)thw(Krot) 17 FORMAT('Type de rotation utilisé : ',a) END SELECT c contrôles du pas temporel SELECT CASE(langue) CASE('english') c WRITE(*,1018)precix,precit,d_lum,d_press,d_ray,d_temp c WRITE(2,1018)precix,precit,d_lum,d_press,d_ray,d_temp WRITE(*,1018)precix,precit,d_press,d_temp WRITE(2,1018)precix,precit,d_press,d_temp 1018 FORMAT('Time step controls:',/, 1 '1 - convergence of quasi-statique equilibrium',es10.3,/, 2 '2 - limitation of changes of chemicals per mass unit',es10.3,/, c 3 '3 - limitation of changes of L respect to time',es10.3,/, 4 '4 - limitation of changes of P with respect to time',es10.3,/, c 5 '5 - limitation of changes of R with respect to time',es10.3,/, 6 '6 - limitation of changes of T with respect to time',es10.3,/) CASE DEFAULT c WRITE(*,18)precix,precit,d_lum,d_press,d_ray,d_temp c WRITE(2,18)precix,precit,d_lum,d_press,d_ray,d_temp WRITE(*,18)precix,precit,d_press,d_temp WRITE(2,18)precix,precit,d_press,d_temp 18 FORMAT(/,'Contrôles du pas temporel:',/, 1 '1 - convergence de la solution équilibre quasi-statique',es10.3,/, 2 '2 - limitation de la var. temporelle des abon/ masse',es10.3,/, c 3 '3 - limitation de la variation temporelle de L',es10.3,/, 4 '4 - limitation de la variation temporelle de P',es10.3,/, c 5 '5 - limitation de la variation temporelle de R',es10.3,/, 6 '6 - limitation de la variation temporelle de T',es10.3,/) END SELECT c formation de la chaine methode WRITE(text1,213)m_qs ; WRITE(text2,213)m_ch 213 FORMAT(i3) methode='CESAM2k version '//version IF(en_m23)THEN methode=TRIM(methode)//' lagr23 colloc'//text1//text2 ELSE methode=TRIM(methode)//' lagr13 colloc'//text1//text2 ENDIF methode=TRIM(methode)//' '//precision IF(diffusion)THEN methode=TRIM(methode)//' diffus' ELSE methode=TRIM(methode)//' no diffus' ENDIF c initialisation des modèles chain_atm=TRIM(nom_fich2)//'_B.atm' chain_don=TRIM(nom_fich2)//'.don' chain_rep=TRIM(nom_fich1) c-----initialisations des conditions de calcul fin, reprise des modèles de départ SELECT CASE(un23) c poursuite d'une évolution modèle repris en binaire, début CASE(1) OPEN(unit=4,form='unformatted',status='old',file=nom_fich1) READ(4)nep,m_qsp,n_qs,knot,nchimp,n_ch,m_chp,knotc,Krotp,nrotp, 1 n_rot,m_rotp,knotr,n_tds,knot_tds,mtotp,alphap,w_rotp,lim_rop, 3 diffusionp,rot_solidp,precisionp,en_m23p,f_eosp,f_opap, 4 nom_ctesp,nom_pertmp,nom_pertwp,nom_tdetaup,nom_atmp,nom_convp, 5 nom_nucp,nom_nuc_cplp,nom_diffmp,nom_difftp,nom_diffwp, 6 nom_etatp,nom_opap c PRINT*,nep,m_qsp,n_qs,knot,nchimp,n_ch c PRINT*,m_chp,knotc,n_tds,knot_tds,mtotp c PRINT*,diffusionp,rot_solidp,precisionp,en_m23p c PRINT*,precisione c PAUSE'case 1' c vérification de concordance c pour la masse sauf s'il y a perte de masse IF(mtotp /= mtot .AND. .NOT.(l_pertm .OR. l_pertmt))THEN SELECT CASE(langue) CASE('english') WRITE(*,1700)mtotp,TRIM(chain_rep),mtot,TRIM(chain_don) WRITE(2,1700)mtotp,TRIM(chain_rep),mtot,TRIM(chain_don) 1700 FORMAT('STOP, the initial mass :',es10.3, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial mass : ',es10.3, 3 ' read the input file : ',a) CASE DEFAULT WRITE(*,700)mtotp,TRIM(chain_rep),mtot,TRIM(chain_don) WRITE(2,700)mtotp,TRIM(chain_rep),mtot,TRIM(chain_don) 700 FORMAT('ARRET, masse init. modèle repris : ',es10.3, 1 ', fichier : ',a,/,'/= masse init. des données : ',es10.3, 2' fichier : ',a) END SELECT CALL sortie('cesam 13') ENDIF IF(alphap /= alpha)THEN SELECT CASE(langue) CASE('english') WRITE(*,1701)alphap,TRIM(chain_rep),alpha,TRIM(chain_don) WRITE(2,1701)alphap,TRIM(chain_rep),alpha,TRIM(chain_don) 1701 FORMAT('STOP, the initial mixing length :',es10.3, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial mixing length : ',es10.3, 3 ' read the input file : ',a,/, 4 'does one keep on computing with the continued model?',/, 5 'enter o for yes; n for no') CASE DEFAULT WRITE(*,701)alphap,TRIM(chain_rep),alpha,TRIM(chain_don) WRITE(2,701)alphap,TRIM(chain_rep),alpha,TRIM(chain_don) 701 FORMAT('ARRET, long. mel. modèle repris : ',es10.3, 1 ', fichier : ',a,/,'/= long. mel. des données : ',es10.3, 2 ', fichier : ',a,/, 3 'Poursuit-on avec la nouvelle valeur? entrer o/n') END SELECT READ*,oui IF(oui == 'n')THEN SELECT CASE(langue) CASE('english') WRITE(*,1702) ; WRITE(*,1702) ; CALL sortie('cesam 14') 1702 FORMAT('STOP') CASE DEFAULT WRITE(*,702) ; WRITE(*,702) ; CALL sortie('cesam 15') 702 FORMAT('ARRET') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1703)chain_don ; WRITE(2,1703)chain_don 1703 FORMAT('The computation is carried on with the data of',/, 1 'the file : ',a) CASE DEFAULT WRITE(*,703)chain_don ; WRITE(2,703)chain_don 703 FORMAT('on poursuit avec les données du fichier : ',a) END SELECT ENDIF ENDIF c il faut les mêmes paramètres de rotation IF(nrotp /= nrot)THEN SELECT CASE(langue) CASE('english') WRITE(*,1129)nrotp,TRIM(chain_rep),nrot,TRIM(chain_don) WRITE(2,1129)nrotp,TRIM(chain_rep),nrot,TRIM(chain_don) 1129 FORMAT('STOP, nb. of var. for rotation :',i2, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the nb. of var. for rotation : ',i2, 3 ' read the input file : ',a) CASE DEFAULT WRITE(*,129)nrotp,TRIM(chain_rep),nrot,TRIM(chain_don) WRITE(2,129)nrotp,TRIM(chain_rep),nrot,TRIM(chain_don) 129 FORMAT('ARRET, nb. coeff. de rotation modèle repris : ',i2, 1 ', fichier : ',a,/,'/= nb. coeff. de rotation des données : ', 2 i2,' fichier : ',a) END SELECT CALL sortie('cesam 16') ENDIF IF(m_rotp /= m_rot)THEN SELECT CASE(langue) CASE('english') WRITE(*,1131)m_rotp,TRIM(chain_rep),m_rot,TRIM(chain_don) WRITE(2,1131)m_rotp,TRIM(chain_rep),m_rot,TRIM(chain_don) 1131 FORMAT('STOP, the spline''s order for the rotation :',i2, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from those : ',i2,' read the input file : ',a) CASE DEFAULT WRITE(*,131)m_rotp,TRIM(chain_rep),m_rot,TRIM(chain_don) WRITE(2,131)m_rotp,TRIM(chain_rep),m_rot,TRIM(chain_don) 131 FORMAT('ARRET, l''ordre des splines :',i2, 1 ' pour la rotation du modèle repris dans le fichier : ',a,/, 2 'diffère de celui des données : ',i2,' fichier : ',a) END SELECT CALL sortie('cesam 17') ENDIF IF(Krotp /= Krot)THEN SELECT CASE(langue) CASE('english') WRITE(*,1132)Krotp,TRIM(chain_rep),Krot,TRIM(chain_don) WRITE(2,1132)Krotp,TRIM(chain_rep),Krot,TRIM(chain_don) 1132 FORMAT('STOP, the flag for rotation Krot :',i2, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial flag for rotation Krot : ',i2, 3 ' read the input file : ',a) CASE DEFAULT WRITE(*,132)Krotp,TRIM(chain_rep),Krot,TRIM(chain_don) WRITE(2,132)Krotp,TRIM(chain_rep),Krot,TRIM(chain_don) 132 FORMAT('ARRET, le flag de rotation du modèle repris Krot : ',i2, 1 ', fichier : ',a,/,'/= flag de rotation des données Krot : ',i2, 2 ' fichier : ',a) END SELECT CALL sortie('cesam 18') ENDIF c si la rotation a été entrée en kms/s w = v / R IF(w_rotp /= w_rot .AND. (unit /= 'kms/s'))THEN SELECT CASE(langue) CASE('english') WRITE(*,1704)w_rotp,TRIM(chain_rep),w_rot,TRIM(chain_don) WRITE(2,1704)w_rotp,TRIM(chain_rep),w_rot,TRIM(chain_don) 1704 FORMAT('STOP, the initial angular velocity :',es10.3, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial angular velocity : ',es10.3, 3 ' read the input file : ',a) CASE DEFAULT WRITE(*,704)w_rotp,TRIM(chain_rep),w_rot,TRIM(chain_don) WRITE(2,704)w_rotp,TRIM(chain_rep),w_rot,TRIM(chain_don) 704 FORMAT('ARRET, vit. ang. ini. modèle repris : ',es10.3, 1 ', fichier : ',a,/,'/= vit. ang. init. des données : ',es10.3, 2 ' fichier : ',a) END SELECT CALL sortie('cesam 19') ENDIF c il convient d'utiliser la même condition limite pour l'atmosphère IF(lim_rop .NEQV. lim_ro)THEN SELECT CASE(langue) CASE('english') WRITE(*,1705)lim_rop,TRIM(chain_rep),lim_ro,TRIM(chain_don) WRITE(2,1705)lim_rop,TRIM(chain_rep),lim_ro,TRIM(chain_don) 1705 FORMAT('WARNING, the boundary conditions lim_ro : ',l1,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the boundary conditions lim_ro : ',l1,/, 3 'read the input file : ',a,/, 4 'The computation is carried on with the data of the input file') CASE DEFAULT WRITE(*,705)lim_rop,TRIM(chain_rep),lim_ro,TRIM(chain_don) WRITE(2,705)lim_rop,TRIM(chain_rep),lim_ro,TRIM(chain_don) 705 FORMAT('ATTENTION, la condition limite lim_ro : ',l1,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffèrent des conditions limites lim_ro : ',l1,/, 3 'du fichier des données: ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c on peut poursuivre sans diffusion un modèle avec diffusion et réciproquement IF(diffusionp .NEQV. diffusion)THEN c IF(.FALSE.)THEN SELECT CASE(langue) CASE('english') WRITE(*,1706)diffusionp,TRIM(chain_rep),diffusion, 1 TRIM(chain_don) WRITE(2,1706)diffusionp,TRIM(chain_rep),diffusion, 1 TRIM(chain_don) 1706 FORMAT('WARNIG, the diffusion parameter diffusion : ',l1,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial diffusion parameter diffusion : ', 3 l1,/,'read the input file : ',a) CASE DEFAULT WRITE(*,706)diffusionp,TRIM(chain_rep),diffusion, 1 TRIM(chain_don) WRITE(2,706)diffusionp,TRIM(chain_rep),diffusion, 1 TRIM(chain_don) 706 FORMAT('ATTENTION, le paramètre de diffusion : ',l1,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de diffusion : ',l1,/, 3 'du fichier des données : ',a) END SELECT c CALL sortie('cesam 21') ENDIF c on ne peut poursuivre en rotation solide par un modèle en rotation c non solide et réciproquement IF(rot_solidp .NEQV. rot_solid)THEN SELECT CASE(langue) CASE('english') WRITE(*,1707)rot_solidp,TRIM(chain_rep),rot_solid, 1 TRIM(chain_don) WRITE(2,1707)rot_solidp,TRIM(chain_rep),rot_solid,TRIM(chain_don) 1707 FORMAT('STOP, the rotation parameter rot_solid : ',l1,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial rotation parameter rot_solid : ', 3 l1,/,'read the input file : ',a) CASE DEFAULT WRITE(*,707)rot_solidp,TRIM(chain_rep),rot_solid, 1 TRIM(chain_don) WRITE(2,707)rot_solidp,TRIM(chain_rep),rot_solid, 1 TRIM(chain_don) 707 FORMAT('ARRET, le paramètre de rotation rot_solid : ',l1,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de rotation rot_solid : ',l1,/, 3 'du fichier des données : ',a) END SELECT CALL sortie('cesam 22') ENDIF c on ne peut changer m23 <==> m13 en cours d'évolution IF(en_m23p .NEQV. en_m23)THEN SELECT CASE(langue) CASE('english') WRITE(*,1708)en_m23p,TRIM(chain_rep),en_m23,TRIM(chain_don) WRITE(2,1708)en_m23p,TRIM(chain_rep),en_m23,TRIM(chain_don) 1708 FORMAT('STOP, the structure parameter en_m23 : ',l1,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial structure parameter en_m23 : ',l1,/, 3 'read the input file : ',a,/,'restart the evolution from initial') CASE DEFAULT WRITE(*,708)en_m23p,TRIM(chain_rep),en_m23,TRIM(chain_don) WRITE(2,708)en_m23p,TRIM(chain_rep),en_m23,TRIM(chain_don) 708 FORMAT('ARRET, le paramètre de structure en_m23 : ',l1,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de structure en_m23 : ',l1,/, 3 'du fichier des données : ',a,/,'reprendre à partir du modèle initial') END SELECT CALL sortie('cesam 23') ENDIF c il faut garder le même type de précision IF(precisionp /= precisione)THEN SELECT CASE(langue) CASE('english') WRITE(*,1709)precisionp,TRIM(chain_rep),precisione, 1 TRIM(chain_don) WRITE(2,1709)precisionp,TRIM(chain_rep),precisione, 1 TRIM(chain_don) 1709 FORMAT('WARNING, the precision parameter precision : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial precision parameter precision : ', 3 a,/,'read the input file : ',a,/) CASE DEFAULT WRITE(*,709)precisionp,TRIM(chain_rep),precisione, 1 TRIM(chain_don) WRITE(2,709)precisionp,TRIM(chain_rep),precisione, 1 TRIM(chain_don) 709 FORMAT('ATTENTION, le paramètre de précision : ',a,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de précision : ',a,/, 3 'du fichier des données : ',a,/) END SELECT ENDIF c il convient de conserver les mêmes valeurs des constantes IF(nom_ctesp /= nom_ctes)THEN SELECT CASE(langue) CASE('english') WRITE(*,1710)nom_ctesp,TRIM(chain_rep),nom_ctes,TRIM(chain_don) WRITE(2,1710)nom_ctesp,TRIM(chain_rep),nom_ctes,TRIM(chain_don) 1710 FORMAT('WARNING, the physical constants parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial physical constants parameter : ', 3 a,/, 'read the input file : ',a,/, 4 'The computation is carried on with the data of the input file') CASE DEFAULT WRITE(*,710)nom_ctesp,TRIM(chain_rep),nom_ctes,TRIM(chain_don) WRITE(2,710)nom_ctesp,TRIM(chain_rep),nom_ctes,TRIM(chain_don) 710 FORMAT('ATTENTION, le paramètre de constantes physiques : ', 1 a,/,'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de constantes physiques : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c il convient de conserver le même type de perte de masse IF(nom_pertmp /= nom_pertm)THEN SELECT CASE(langue) CASE('english') WRITE(*,1711)nom_pertmp,TRIM(chain_rep),nom_pertm, 1 TRIM(chain_don) WRITE(2,1711)nom_pertmp,TRIM(chain_rep),nom_pertm, 1 TRIM(chain_don) 1711 FORMAT('WARNING, The mass loss parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial mass loss parameter : ',a,/, 3 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file') CASE DEFAULT WRITE(*,711)nom_pertmp,TRIM(chain_rep),nom_pertm, 1 TRIM(chain_don) WRITE(2,711)nom_pertmp,TRIM(chain_rep),nom_pertm, 1 TRIM(chain_don) 711 FORMAT('ATTENTION, le paramètre de perte de masse : ',a,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de perte de masse : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c il convient de conserver le même type de perte de moment cinétique IF(entre(3,4,Krot) .AND. nom_pertwp /= nom_pertw)THEN SELECT CASE(langue) CASE('english') WRITE(*,1699)nom_pertwp,TRIM(chain_rep),nom_pertw, 1 TRIM(chain_don) WRITE(2,1699)nom_pertwp,TRIM(chain_rep),nom_pertw, 1 TRIM(chain_don) 1699 FORMAT('WARNING, The ang. momentum loss parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial ang. momentum loss parameter : ',a,/, 3 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file') CASE DEFAULT WRITE(*,699)nom_pertwp,TRIM(chain_rep),nom_pertw, 1 TRIM(chain_don) WRITE(2,699)nom_pertwp,TRIM(chain_rep),nom_pertw, 1 TRIM(chain_don) 699 FORMAT('ATTENTION, le paramètre de perte de moment cin. : ', 1 a,/,'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de perte de moment cin. : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c il conviendrait de conserver la même loi T(tau) IF(nom_tdetaup /= nom_tdetau)THEN SELECT CASE(langue) CASE('english') WRITE(*,1712)nom_tdetaup,TRIM(chain_rep),nom_tdetau, 1 TRIM(chain_don) WRITE(2,1712)nom_tdetaup,TRIM(chain_rep),nom_tdetau, 1 TRIM(chain_don) 1712 FORMAT('WRANING, the T(tau)-law parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial T(tau)-law parameter : ',a,/, 3 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file') CASE DEFAULT WRITE(*,712)nom_tdetaup,TRIM(chain_rep),nom_tdetau, 1 TRIM(chain_don) WRITE(2,712)nom_tdetaup,TRIM(chain_rep),nom_tdetau, 1 TRIM(chain_don) 712 FORMAT('ATTENTION, le paramètre de loi T(tau) : ',a,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de loi T(tau) : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c il conviendrait de conserver la même atmosphère IF(nom_atmp /= nom_atm)THEN SELECT CASE(langue) CASE('english') WRITE(*,1713)nom_atmp,TRIM(chain_rep),nom_atm,TRIM(chain_don) WRITE(2,1713)nom_atmp,TRIM(chain_rep),nom_atm,TRIM(chain_don) 1713 FORMAT('WARNING, the atmosphere parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial atmosphere parameter : ',a,/, 3 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file') CASE DEFAULT WRITE(*,713)nom_atmp,TRIM(chain_rep),nom_atm,TRIM(chain_don) WRITE(2,713)nom_atmp,TRIM(chain_rep),nom_atm,TRIM(chain_don) 713 FORMAT('ATTENTION, le paramètre d''atmosphère : ',a,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre d''atmosphère : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c si on poursuit avec une atmosphère simplifiée il faut n_atm=0 IF(nom_atm /= 'lim_atm')n_atm=0 IF(nom_convp /= nom_conv)THEN SELECT CASE(langue) CASE('english') WRITE(*,1714)nom_convp,TRIM(chain_rep),nom_conv,TRIM(chain_don) WRITE(2,1714)nom_convp,TRIM(chain_rep),nom_conv,TRIM(chain_don) 1714 FORMAT('WARNING, the convection parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial convection parameter : ',a,/, 3 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file') CASE DEFAULT WRITE(*,714)nom_convp,TRIM(chain_rep),nom_conv,TRIM(chain_don) WRITE(2,714)nom_convp,TRIM(chain_rep),nom_conv,TRIM(chain_don) 714 FORMAT('ATTENTION, le paramètre de convection : ',a,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de convection : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c on ne peut poursuivre avec un réseau nucléaire différent IF(nom_nucp /= nom_nuc)THEN SELECT CASE(langue) CASE('english') WRITE(*,1715)nom_nucp,TRIM(chain_rep),nom_nuc,TRIM(chain_don) WRITE(2,1715)nom_nucp,TRIM(chain_rep),nom_nuc,TRIM(chain_don) 1715 FORMAT('The nuclear reactions parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial nuclear reactions parameter : ',a,/, 3 'read in the input file : ',a,/, 4 'The computation is carried on with data of the input file', 5 /,'only if the chemical elements are identical') CASE DEFAULT WRITE(*,715)nom_nucp,TRIM(chain_rep),nom_nuc,TRIM(chain_don) WRITE(2,715)nom_nucp,TRIM(chain_rep),nom_nuc,TRIM(chain_don) 715 FORMAT('Le paramètre de réactions nucléaires : ',a,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de réactions nucléaires : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données que si les' 5 ,/,'éléments chimiques utilisés sont identiques') END SELECT IF((TRIM(nom_nucp) == 'ppcno9' .AND. 1 TRIM(nom_nuc) == 'ppcno3a9').OR. 2 (TRIM(nom_nucp) == 'ppcno3a9' .AND. 3 TRIM(nom_nuc) == 'ppcno9'))THEN SELECT CASE(langue) CASE('english') WRITE(*,1787) ; WRITE(2,1787) 1787 FORMAT('The computation is carried on.') CASE DEFAULT WRITE(*,787) ; WRITE(2,787) 787 FORMAT('Le calcul est poursuivi.') END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1788) ; WRITE(2,1788) 1788 FORMAT('STOP, the chemical elements are not identical.') CASE DEFAULT WRITE(*,788) ; WRITE(2,788) 788 FORMAT('ARRET, les éléments chimiques utilisés diffèrent ') END SELECT CALL sortie('cesam 24') ENDIF ENDIF c il convient de poursuivre avec la même compil. de réactions nucléaires IF(nom_nuc_cplp /= nom_nuc_cpl)THEN SELECT CASE(langue) CASE('english') WRITE(*,1716)nom_nuc_cplp,TRIM(chain_rep),nom_nuc_cpl, 1 TRIM(chain_don) WRITE(2,1716)nom_nuc_cplp,TRIM(chain_rep),nom_nuc_cpl, 1 TRIM(chain_don) 1716 FORMAT('WARNING, the compilation of nuclear reactions : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial compilation of nuclear reactions : ', 3 a,/, 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file', 5 /,'if the chemical elements are identical') CASE DEFAULT WRITE(*,716)nom_nuc_cplp,TRIM(chain_rep),nom_nuc_cpl, 1 TRIM(chain_don) WRITE(2,716)nom_nuc_cplp,TRIM(chain_rep),nom_nuc_cpl, 1 TRIM(chain_don) 716 FORMAT('ATTENTION, la compilation des réactions nucléaires : ' 1 ,a,/,'du modèle repris du fichier : ',a,/, 2 'diffère de la compilation des réactions nucléaires : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données si les',/, 5 'éléments chimiques utilisés sont identiques') END SELECT ENDIF c il convient de poursuivre avec la même diffusion microscopique IF(nom_diffmp /= nom_diffm)THEN SELECT CASE(langue) CASE('english') WRITE(*,1717)nom_diffmp,TRIM(chain_rep),nom_diffm, 1 TRIM(chain_don) WRITE(2,1717)nom_diffmp,TRIM(chain_rep),nom_diffm, 1 TRIM(chain_don) 1717 FORMAT('WARNING, the microscopic diffusion parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial microscopic diffusion parameter : ', 3 a,/, 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file') CASE DEFAULT WRITE(*,717)nom_diffmp,TRIM(chain_rep),nom_diffm, 1 TRIM(chain_don) WRITE(2,717)nom_diffmp,TRIM(chain_rep),nom_diffm, 1 TRIM(chain_don) 717 FORMAT('ATTENTION, le paramètre de diffusion microscopique : ' 1 ,a,/,'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de diffusion microscopique : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c il convient de poursuivre avec la même diffusion du moment cinétique IF(nom_diffwp /= nom_diffw)THEN SELECT CASE(langue) CASE('english') WRITE(*,1724)nom_diffwp,TRIM(chain_rep),nom_diffw, 1 TRIM(chain_don) WRITE(2,1724)nom_diffwp,TRIM(chain_rep),nom_diffw, 1 TRIM(chain_don) 1724 FORMAT('WARNING, the turbulent diffusion parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial turbulent diffusion parameter : ', 3 a,/, 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file') CASE DEFAULT WRITE(*,724)nom_diffwp,TRIM(chain_rep),nom_diffw, 1 TRIM(chain_don) WRITE(2,724)nom_diffwp,TRIM(chain_rep),nom_diffw, 1 TRIM(chain_don) 724 FORMAT('ATTENTION, le paramètre de diffusion turbulente : ' 1 ,a,/,'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre de diffusion turbulente : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c il convient de poursuivre avec la même diffusion turbulente IF(nom_difftp /= nom_difft)THEN SELECT CASE(langue) CASE('english') WRITE(*,1718)nom_difftp,TRIM(chain_rep),nom_difft, 1 TRIM(chain_don) WRITE(2,1718)nom_difftp,TRIM(chain_rep),nom_difft, 1 TRIM(chain_don) 1718 FORMAT('WARNING, the turbulent diffusion routine : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial turbulent diffusion routine : ', 3 a,/, 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file') CASE DEFAULT WRITE(*,718)nom_difftp,TRIM(chain_rep),nom_difft, 1 TRIM(chain_don) WRITE(2,718)nom_difftp,TRIM(chain_rep),nom_difft, 1 TRIM(chain_don) 718 FORMAT('ATTENTION, la routine de diffusion turbulente : ',a,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère de la routine de diffusion turbulente : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données') END SELECT ENDIF c il convient de poursuivre avec la même EOS IF(nom_etatp /= nom_etat)THEN SELECT CASE(langue) CASE('english') WRITE(*,1719)nom_etatp,TRIM(chain_rep),nom_etat,TRIM(chain_don) WRITE(2,1719)nom_etatp,TRIM(chain_rep),nom_etat,TRIM(chain_don) 1719 FORMAT('WARNING, the equation of state parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial equation of state parameter : ', 3 a,/, 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file', 5 /,'The program stops if the EOS file is not correct') CASE DEFAULT WRITE(*,719)nom_etatp,TRIM(chain_rep),nom_etat,TRIM(chain_don) WRITE(2,719)nom_etatp,TRIM(chain_rep),nom_etat,TRIM(chain_don) 719 FORMAT('ATTENTION, le paramètre d''équation d''état : ',a,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre d''équation d''état : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données PLANTAGE ' 5 ,/,'si le nom du fichier d''équation d''état est incorrect') END SELECT ENDIF c il convient de poursuivre avec la même opacité IF(nom_opap /= nom_opa)THEN SELECT CASE(langue) CASE('english') WRITE(*,1720)nom_opap,TRIM(chain_rep),nom_opa,TRIM(chain_don) WRITE(2,1720)nom_opap,TRIM(chain_rep),nom_opa,TRIM(chain_don) 1720 FORMAT('WARNING, the opacity parameter : ',a,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the initial opacity parameter : ', 3 a,/, 'read the input file : ',a,/, 4 'The computation is carried on with data of the input file', 5 /,'The program stops if the opacity file is not correct') CASE DEFAULT WRITE(*,720)nom_opap,TRIM(chain_rep),nom_opa,TRIM(chain_don) WRITE(2,720)nom_opap,TRIM(chain_rep),nom_opa,TRIM(chain_don) 720 FORMAT('ATTENTION, le paramètre d''opacité : ',a,/, 1 'du modèle repris du fichier : ',a,/, 2 'diffère du paramètre d''opacité : ',a,/, 3 'du fichier des données : ',a,/, 4 'on poursuit avec la valeur du fichier des données PLANTAGE ' 5 ,/,'si le nom du fichier d''opacité est incorrect') END SELECT ENDIF c il convient de poursuivre avec le même ordre de splines pour les c variables de structure IF(m_qsp /= m_qs)THEN SELECT CASE(langue) CASE('english') WRITE(*,1004)m_qsp,TRIM(chain_rep),m_qs,TRIM(chain_don) WRITE(2,1004)m_qsp,TRIM(chain_rep),m_qs,TRIM(chain_don) 1004 FORMAT('WARNING, the order of splines for the quasi-static',/, 1 'equilibrium m_qs=',i3,/, 2 'read in the file of the model to be continued : ',a,/, 3 'differs from the order of the spline m_qs=',i3, 4 ', read the input file : ',a,/,'one keeps on computing with',/, 5 'the value read in the file of the model to be continued') CASE DEFAULT WRITE(*,4)m_qsp,TRIM(chain_rep),m_qs,TRIM(chain_don) WRITE(2,4)m_qsp,TRIM(chain_rep),m_qs,TRIM(chain_don) 4 FORMAT('ATTENTION, l''ordre de la spline m_qs=',i3,/, 1 'pour l''équilibre quasi-statique du modèle : ',a,/, 2 'diffère de l''ordre de la spline m_qs=',i3,/, 3 'du modèle à poursuivre : ',a,/, 4 'on garde l''ordre de la spline du modèle repris') END SELECT ENDIF m_qs=m_qsp c il convient de poursuivre avec le même ordre de splines pour la c composition chimique IF(m_chp /= m_ch)THEN SELECT CASE(langue) CASE('english') WRITE(*,1006)m_chp,TRIM(chain_rep),m_ch,TRIM(chain_don) WRITE(2,1006)m_chp,TRIM(chain_rep),m_ch,TRIM(chain_don) 1006 FORMAT('WARNING, the order of spline for the interpolation of' 1 ,/,'the chemical composition, m_ch=',i3,/, 2 'read in the file of the model to be continued : ',a,/, 3 'differs from the order of the spline m_ch=',i3, 4 'read the input file : ',a,/,'one keeps on computing with',/, 5 'the value read in the file of the model to be continued') CASE DEFAULT WRITE(*,6)m_chp,TRIM(chain_rep),m_ch,TRIM(chain_don) WRITE(2,6)m_chp,TRIM(chain_rep),m_ch,TRIM(chain_don) 6 FORMAT('ATTENTION, l''ordre de la spline m_ch=',i3,/, 1 'pour la composition chimique du modèle repris : ',a,/, 2 'diffère de l''ordre de la spline m_ch=',i3,/, 3 'du modèle à poursuivre : ',a,/, 4 'on garde l''ordre de la spline du modèle repris') END SELECT ENDIF m_ch=m_chp c reprise du modèle, on ajuste les dim. à celles du modèle repris ord_qs=m_qs+r_qs ; dim_qs=knot-ord_qs ; dim_ch=knotc-m_ch dim_rot=knotr-ord_rot ; m_tds=knot_tds-n_tds ALLOCATE(bp(nep,dim_qs),q(n_qs),qt(knot),nom_elemp(nchimp), 1 chim(nchimp,dim_ch),mc(n_ch),mct(knotc),mrot(n_rot), 2 mrott(knotr),tds(1,knot_tds-m_tds),x_tds(n_tds), 3 xt_tds(knot_tds),m_stat(n_qs),r_stat(n_qs),rota(nrot,dim_rot)) REWIND(unit=4) READ(4)ne,m_qs,n_qs,knot,nchim,n_ch,m_ch,knotc,Krot,nrot,n_rot, 1 m_rot,knotr,n_tds,knot_tds,mtotp,alphap,w_rotp, 2 lim_rop,diffusionp,rot_solidp,precisionp,en_m23p,f_eosp, 3 f_opap,nom_ctesp,nom_pertmp,nom_pertwp,nom_tdetaup,nom_atmp, 4 nom_convp,nom_nucp,nom_nuc_cplp,nom_diffmp,nom_difftp, 5 nom_diffwp,nom_etatp,nom_opap,nom_elemp, 6 bp,q,qt,chim,mc,mct,rota,mrot,mrott,tds,x_tds,xt_tds,m_stat, 7 r_stat,m_zc,r_zc,r_ov,age,dtp,dts,mstar,rstar,mw_tot,wrot,jlim, 8 lconv,lim,model_num CLOSE(unit=4) model_num=model_num-1 !est augmenté d'une unité ligne ~4155 nb_der_out=model_num-dn_sort+1 !n° dernière sortie c PRINT*,m_ch,n_ch,knotc,nom_fich1 c WRITE(*,2000)mct(1:4),mct(knotc-3:knotc) c PAUSE'READ' c écritures SELECT CASE(langue) CASE('english') WRITE(*,1021)model_num+1,age ; WRITE(2,1021)model_num+1,age 1021 FORMAT('Start from binary model #',i5,' aged',es10.3,/) CASE DEFAULT WRITE(*,21)model_num+1,age ; WRITE(2,21)model_num+1,age 21 FORMAT('Reprise du modèle binaire n°',i5,' d''âge',es10.3,/) END SELECT c si on utilise une grille en m^2/3 fixe pour la composition chimique IF(grille_fixe)THEN ALLOCATE(mc_fixe(n_qs)) nc_fixe=n_qs ; mc_fixe=m_stat ENDIF c appel d'initialisation pour tabulation des réactions nucléaires ALLOCATE(comp(0),dcomp(0),jac(0,0),ex(0)) CALL nuc(1.5d+07,1.5d+02,comp,dcomp,jac,.FALSE.,0, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c la rotation WRITE(*,17)thw(Krot) ; WRITE(2,17)thw(Krot) c vitesse angulaire minimale IF(entre(3,5,Krot))rot_min(1)=wrot*1.d-2 SELECT CASE(langue) CASE('english') WRITE(*,1023)wrot ; WRITE(2,1023)wrot 1023 FORMAT('initial angular velocity =',es10.3,' rd/sec',/) CASE DEFAULT WRITE(*,23)wrot ; WRITE(2,23)wrot 23 FORMAT('vitesse angulaire initiale =',es10.3,' rd/sec',/) END SELECT c détermination des abondances initiales DEALLOCATE(comp) ; ALLOCATE(comp(nchim)) CALL nuc(1.5d+07,1.5d+02,comp,dcomp,jac,.FALSE.,1, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c PRINT*,nom_elem ; PRINT*,nom_elemp ; PAUSE'elem' IF(nchimp /= nchim)THEN SELECT CASE(langue) CASE('english') WRITE(*,1010)nchimp,TRIM(chain_rep),nchim,TRIM(chain_don) 1010 FORMAT('STOP, the initial number of chemical elements : ',i3,/, 1 'read in the file of the model to be continued : ',a,/, 2 'differs from the number of chemical elements : ',i3,/, 3 'read the input file : ',a) CASE DEFAULT WRITE(*,10)nchimp,TRIM(chain_rep),nchim,TRIM(chain_don) 10 FORMAT('ARRET, le nombre d''éléments chimimiques = ',i3,/, 1 'du modèle repris : ',a,/, 2 'diffère du nombre d''éléments chimiques = ',i3,/, 3 'du modèle à poursuivre : ',a) END SELECT CALL sortie('cesam 25') ENDIF ok=.TRUE. DO i=1,nchim ok=ok .AND. nom_elem(i) == nom_elemp(i) ENDDO IF(.NOT.ok)THEN SELECT CASE(langue) CASE('english') WRITE(*,1011)nom_elemp WRITE(*,1012)TRIM(chain_rep),nom_elem WRITE(*,1032)TRIM(chain_don) WRITE(2,1011)nom_elemp WRITE(2,1012)TRIM(chain_rep),nom_elem WRITE(2,1032)TRIM(chain_don) 1011 FORMAT('STOP,the chemical elements : ',30a) 1012 FORMAT('read the input file : ',a,/, 1 'differ from the chemical elements : ',30a) 1032 FORMAT('read the input file : ',a) CASE DEFAULT WRITE(*,11)nom_elemp WRITE(*,12)TRIM(chain_rep),nom_elem WRITE(*,32)TRIM(chain_don) WRITE(2,11)nom_elemp WRITE(2,12)TRIM(chain_rep),nom_elem WRITE(2,32)TRIM(chain_don) 11 FORMAT('ARRET, les élements chimiques : ',30a) 12 FORMAT('du modèle repris : ',a,/,'diffèrent de ceux : ',30a) 32 FORMAT('du modèle à poursuivre : ',a) END SELECT CALL sortie('cesam 26') ENDIF c reprise: on utilise la valeur de psi0 en cours et non celle initialisée par défaut c psi0=bp(6,1) !soleil de ZAMS 8.5d-2 c tout semble OK, let's go c suppression/adjonction de la pression turbulente c si le modèle repris est avec/sans Pturb c et que le modèle à calculer est sans/avec Pturb c il faut enlever/ajouter la variable ln Pgaz c ne/nep : nb. inconnues du modèle à calculer/repris IF(ne /= nep)THEN ALLOCATE(b(nep,dim_qs)) ; b=bp DEALLOCATE(bp) ; ALLOCATE(bp(ne,dim_qs)) c changement du nb. inconnues !suppression de Pgaz on passe de ne= 7 à 6 inconnues IF(ne == 6)THEN SELECT CASE(langue) CASE('english') WRITE(*,1721) ; WRITE(2,721) 1721 FORMAT('Pturb included in the starting model,',/, 1 'but not to be included in the model to be calculated') CASE DEFAULT WRITE(*,721) ; WRITE(2,721) 721 FORMAT('modèle repris avec Pturb,',/, 1 'modèle à calculer sans Pturb') END SELECT bp=b(1:6,:) ELSE !adjonction de Pgaz on passe de 6 à 7 inconnues SELECT CASE(langue) CASE('english') WRITE(*,1722) ; WRITE(2,722) 1722 FORMAT('Pturb not included in the starting model,'/, 1 'but to be included in the model to be calculated') CASE DEFAULT WRITE(*,722) ; WRITE(2,722) 722 FORMAT('modèle repris sans Pturb,',/, 1 'modèle à calculer avec Pturb') END SELECT bp(1:6,:)=b ; bp(Ipg,:)=b(1,:) ENDIF !passage de 6 <---> 7 inconnues SELECT CASE(langue) CASE('english') WRITE(*,1723)nep,ne ; WRITE(2,1723)nep,ne 1723 FORMAT('we switch from ne=',i3,' unknowns to',/, 1 'ne=',i3,' unknowns') CASE DEFAULT WRITE(*,723)nep,ne ; WRITE(2,723)nep,ne 723 FORMAT('on passe de ne=',i3,' inconnues à',/, 1 'ne=',i3,' inconnues') END SELECT c b est désormais inutile DEALLOCATE(b) ENDIF !transfert et changement du nb. inconnues c WRITE(*,2000)age,agemax,dtp,dts,dt ; PAUSE'dt' c PRINT*,dts ; PAUSE'1' c gestion du premier pas temporel IF(dts > 0.d0 .AND. agemax > 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1111)model_num+1,age,agemax,dtp,dts 1111 FORMAT('model #,'i5,'continued, age (in units of 10^6 y)', 1 es10.3,', agemax=',es10.3,/,'previous time step=',es10.3,/, 2 'estimated optimum time step=',es10.3,/, 3 'does one choose the optimum time step?',/, 4 'enter o for yes, enter n for no') CASE DEFAULT WRITE(*,111)model_num+1,age,agemax,dtp,dts 111 FORMAT('modèle n°',i5,' repris, âge (10^6 ans):',es10.3, 1 ', agemax=',es10.3,/,'pas temporel précédent=',es10.3,/, 2 'pas temporel optimal estimé =',es10.3,/, 3 'utilise-t-on le pas temporel optimal ? entrer o/n') END SELECT READ*,oui IF(oui /= 'o')THEN SELECT CASE(langue) CASE('english') WRITE(*,1113) 1113 FORMAT('enter the new time step, in units of 10^6 years') CASE DEFAULT WRITE(*,113) 113 FORMAT('entrer le nouveau pas temporel, unité: 10^6 ans') END SELECT READ*,dts ENDIF !oui /= 'o' c PRINT*,dts,dtmax ; PAUSE'2' IF(dts > dtmax)THEN SELECT CASE(langue) CASE('english') WRITE(*,1112)dtmax 1112 FORMAT('dt > dtmax, the time step is set back to dtmax=', 1 es10.3) CASE DEFAULT WRITE(*,112)dtmax 112 FORMAT('dt > dtmax, le pas temporel est ramené à dtmax=', 1 es10.3) END SELECT ENDIF !dt > dtmax dt=MIN(dts,dtmax) ; dts=dt !on garde le pas temporel précédent c PRINT*,dt,dts,dtmax ; PAUSE'3' IF(age+dt > agemax)THEN SELECT CASE(langue) CASE('english') WRITE(*,1053)age+dt,agemax ; WRITE(2,1053)age+dt,agemax 1053 FORMAT('âge+dt=',es10.3,' > agemax =',es10.3,/, 1 'tuning of dt') CASE DEFAULT WRITE(*,53)age+dt,agemax ; WRITE(2,53)age+dt,agemax 53 FORMAT('âge+dt=',es10.3,' > agemax =',es10.3,/, 1 'ajustement de dt') END SELECT dt=MIN(dtmax,MAX(0.d0,agemax-age)) ; dts=dt c WRITE(*,2000)dtmax,agemax-age,dt ; PAUSE'dt' IF(dt <= 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1054) ; WRITE(2,1054) 1054 FORMAT('STOP, agemax is reached, no further evolution',/, 1 'before exit, write listing of the starting model?',/, 2 'enter o for yes, n for no') CASE DEFAULT WRITE(*,54) ; WRITE(2,54) 54 FORMAT('ARRET, agemax est atteint, pas de poursuite',/, 1 'liste du modèle repris avant de sortir? entrer o/n') END SELECT READ*,oui IF(oui == 'o')THEN list_sort=.TRUE. ELSE CALL sortie('cesam 27') ENDIF ENDIF !dt <= 0.d0 ENDIF !age+dt > agemax c PRINT*,dt,dts,dtmax ; PAUSE'4' SELECT CASE(langue) CASE('english') WRITE(2,1050)dt ; WRITE(*,1050)dt 1050 FORMAT(/,'The evolution is continued with the time step', 1 es10.3,' (million years)',/) CASE DEFAULT WRITE(2,50)dt ; WRITE(*,50)dt 50 FORMAT(/,'Poursuite de l''évolution avec le pas temporel', 1 es10.3,' (million d''années)',/) END SELECT list_rep=.TRUE. !listing du modèle repris ELSEIF(agemax > 0.d0)THEN !dtp=0, on poursuit un modèle initial SELECT CASE(langue) CASE('english') WRITE(*,1114)dt0 1114 FORMAT('does one use dt0=',es10.3,' for the time step?',/, 1 'enter o for yes, enter n for no') CASE DEFAULT WRITE(*,114)dt0 114 FORMAT('utilise-t-on dt0=',es10.3,' pour le pas temporel? o/n') END SELECT READ*,oui IF(oui == 'o')THEN dt=MIN(dt0,dtmax) ; dts=dt ELSE SELECT CASE(langue) CASE('english') WRITE(*,1115) 1115 FORMAT('enter the new time step, in units of 10^6 years') CASE DEFAULT WRITE(*,115) 115 FORMAT('entrer le nouveau pas temporel, unité : 10^6 ans') END SELECT READ*,dt ; dts=dt IF(dt > dtmax)THEN dt=dtmax ; dts=dt ; WRITE(*,112)dtmax ENDIF ENDIF !oui /= 'o' c PRINT*,dt,dts,dtmax ; PAUSE'51' IF(age+dt > agemax)THEN SELECT CASE(langue) CASE('english') WRITE(*,1053)age+dt,agemax CASE DEFAULT WRITE(*,53)age+dt,agemax END SELECT dt=MIN(dtmax,MAX(0.d0,agemax-age)) c PRINT*,dt,dts,dtmax ; PAUSE'5' IF(dt <= 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1054) CASE DEFAULT WRITE(*,54) END SELECT READ*,oui IF(oui == 'o')THEN list_sort=.TRUE. ELSE CALL sortie('cesam 28') ENDIF ENDIF !dt <= 0.d0 SELECT CASE(langue) CASE('english') WRITE(2,1051)dt ; WRITE(*,1051)dt ; list_rep=.TRUE. 1051 FORMAT('Evolution starting at age=0 with dt=',es10.3) CASE DEFAULT WRITE(2,51)dt ; WRITE(*,51)dt ; list_rep=.TRUE. 51 FORMAT('Evolution à partir de l''âge=0 avec dt=',es10.3) END SELECT ENDIF !age+dt > agemax list_rep=.TRUE. !listing du modèle repris ELSE !agemax=0 dtp quelconque SELECT CASE(langue) CASE('english') WRITE(*,1009)age,agemax ; WRITE(2,1009)age,agemax 1009 FORMAT('age of the continued model',es10.3,'> agemax : ', 1 es10.3,/,'one writes a listing of the continued model, exits') CASE DEFAULT WRITE(*,9)age,agemax ; WRITE(2,9)age,agemax 9 FORMAT('âge du modèle repris : ',es10.3,' > agemax:',es10.3,/, 1 'on liste le modèle repris et on sort') END SELECT list_sort=.TRUE. ENDIF !dtp > 0 et agemax > 0 c initialisation des paramètres pour le vent et les chutes de planétoïdes CALL vent ; CALL planetoides c Poursuite d'une évolution modèle repris en ASCII CASE(-1) IF(agemax <= 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1610)agemax ; WRITE(2,1610)agemax 1610 FORMAT('With a model to be continued stored in ASCII,',/, 1 'agemax=',es10.3,' is not consistant, STOP') CASE DEFAULT WRITE(*,610)agemax ; WRITE(2,610)agemax 610 FORMAT('Avec un modèle à poursuivre repris en ASCII, agemax=', 1 es10.3,' est incohérent, ARRET') END SELECT CALL sortie('cesam 29') ENDIF OPEN(unit=4,form='formatted',status='old',file=nom_fich1, 1 delim='apostrophe') READ(4,2004)n,nchim_ascii,ihe4_ascii,Krot_ascii 2004 FORMAT(10i5) ALLOCATE(nom_elem_ascii(nchim_ascii)) READ(4,2006)nom_elem_ascii(1:nchim_ascii) 2006 FORMAT(20a4) READ(4,2005)age,dtp,mtot_ascii,x0_ascii,y0_ascii IF(x0_ascii /= x0 .OR. y0_ascii /= y0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1606)x0_ascii,x0,y0_ascii,y0 WRITE(2,1606)x0_ascii,x0,y0_ascii,y0 1606 FORMAT('STOP, x0_ascii=',es10.3,'is different from x0=',es10.3, 1 'or y0_ascii=',es10.3,'is different of y0',es10.3) CASE DEFAULT WRITE(*,606)x0_ascii,x0,y0_ascii,y0 WRITE(2,606)x0_ascii,x0,y0_ascii,y0 606 FORMAT('ARRET,x0_ascii=',es10.3,' diffère de x0=',es10.3, 1 ' ou y0_ascii=',es10.3,' diffère de y0=',es10.3) END SELECT CALL sortie('cesam 30') ENDIF IF(mtot /= mtot_ascii)THEN SELECT CASE(langue) CASE('english') WRITE(*,1607)mtot_ascii,mtot,TRIM(chain_don) WRITE(2,1607)mtot_ascii,mtot,TRIM(chain_don) 1607 FORMAT('the mass of the continued model in ascii mtot=',es10.3, 1 'is different from mtot=',es10.3,/,'read in the file : ',a) CASE DEFAULT WRITE(*,607)mtot_ascii,mtot,TRIM(chain_don) WRITE(2,607)mtot_ascii,mtot,TRIM(chain_don) 607 FORMAT('STOP, la masse du modèle repris en ascii mtot=',es10.3, 1 ' diffère de mtot=',es10.3,/,'lu dans le fichier : ',a) END SELECT ENDIF c appel d'initialisation pour tabulation des réactions nucléaires c allocations fictives ALLOCATE(comp(0),dcomp(0),jac(0,0),ex(0)) CALL nuc(1.5d+07,1.5d+02,comp,dcomp,jac,.FALSE.,0, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c il faut Krot >= 3 IF(entre(3,5,Krot) .AND. Krot /= Krot_ascii)THEN SELECT CASE(langue) CASE('english') WRITE(*,1609)Krot_ascii,TRIM(chain_rep),Krot,TRIM(chain_don) WRITE(2,1609)Krot_ascii,TRIM(chain_rep),Krot,TRIM(chain_don) 1609 FORMAT('STOP, the option for the rotation Krot=',i3,/, 1 'is different for the model to be continued : ',a,/, 1 'from the option Krot=',i3,/,'of the continued model : ',a) CASE DEFAULT WRITE(*,609)Krot_ascii,TRIM(chain_rep),Krot,TRIM(chain_don) WRITE(2,609)Krot_ascii,TRIM(chain_rep),Krot,TRIM(chain_don) 609 FORMAT('ARRET, le type de rotation Krot=',i3,/, 1 'du modèle repris : ',a,/,'est différent de celui Krot=',i3,/, 2 'du modèle à poursuivre : ',a) END SELECT CALL sortie('cesam 31') ELSE SELECT CASE(langue) CASE('english') WRITE(*,1023)wrot ; WRITE(2,1023)wrot CASE DEFAULT WRITE(*,23)wrot ; WRITE(2,23)wrot END SELECT ENDIF c la rotation, avec modèle en ASCII la rotation avec diffusion du moment c cinétique ou conservation locale du moment cinétique c n'est pas incluse dans le modèle repris (modèle repris en ASCII) WRITE(*,17)thw(Krot) ; WRITE(2,17)thw(Krot) SELECT CASE(Krot) CASE(0,1,2,6) ALLOCATE(rota(nrot,n_rot),mrot(n_rot),mrott(knotr)) CASE(3,4,5) SELECT CASE(langue) CASE('english') WRITE(*,1023)wrot ; WRITE(2,1023)wrot CASE DEFAULT WRITE(*,23)wrot ; WRITE(2,23)wrot END SELECT END SELECT c détermination des abondances initiales (modèle repris en ASCII) DEALLOCATE(comp) ; ALLOCATE(comp(nchim)) CALL nuc(1.5d+07,1.5d+02,comp,dcomp,jac,.FALSE.,1, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c PRINT*,nom_elem ; PRINT*,nom_elemp ; PAUSE'elem' c test de cohérence de la composition chimique IF(nchim_ascii /= nchim)THEN SELECT CASE(langue) CASE('english') WRITE(*,1601)nchim_ascii,TRIM(chain_rep),nchim,TRIM(chain_don) WRITE(2,1601)nchim_ascii,TRIM(chain_rep),nchim,TRIM(chain_don) 1601 FORMAT('STOP, the number of elements nchim= ,',i3,/, 1 'of the continued model in ascii : ,'a,/, 2 'is different nchim= ,',i3,/, 3 'from the one built from the data of the input file : ,'a) CASE DEFAULT WRITE(*,601)nchim_ascii,TRIM(chain_rep),nchim,TRIM(chain_don) WRITE(2,601)nchim_ascii,TRIM(chain_rep),nchim,TRIM(chain_don) 601 FORMAT('ARRET, le nombre d''éléments nchim= ,',i3,/, 1 'du modèle repris en ASCII : ',a,/, 2 'diffère de celui nchim= ,',i3,/, 3 'du modèle à calculer avec les données du fichier : ',a) END SELECT CALL sortie('cesam 32') ELSEIF(nchim_ascii /= nchim)THEN SELECT CASE(langue) CASE('english') WRITE(*,1602)nchim_ascii,TRIM(chain_rep),nchim,TRIM(chain_don) WRITE(2,1602)nchim_ascii,TRIM(chain_rep),nchim,TRIM(chain_don) 1602 FORMAT('STOP, the number of elements nchim= ,',i3,/, 1 'of the continued model in ascii : ,'a,/, 2 'is different nchim= ,',i3,/, 3 'from the one built from the data of the input file : ,'a) CASE DEFAULT WRITE(*,602)nchim_ascii,TRIM(chain_rep),nchim,TRIM(chain_don) WRITE(2,602)nchim_ascii,TRIM(chain_rep),nchim,TRIM(chain_don) 602 FORMAT('ARRET, le nombre d''éléments nchim= ,',i3,/, 1 'du modèle repris en ASCII : ',a,/, 2 'diffère de celui nchim= ,',i3,/, 3 'du modèle à calculer avec les données du fichier : ',a) END SELECT CALL sortie('cesam 33') ELSEIF(ihe4_ascii /= ihe4)THEN SELECT CASE(langue) CASE('english') WRITE(*,1605)ihe4_ascii,TRIM(chain_rep),ihe4,TRIM(chain_don) WRITE(2,1605)ihe4_ascii,TRIM(chain_rep),ihe4,TRIM(chain_don) 1605 FORMAT('STOP, the subscript of He4 ihe4=,',i3,/, 1 'of the continued model in ascii : ,'a,/, 2 'is different ih4= ,',i3,/, 3 'from the one built from the data of the input file : ,'a) CASE DEFAULT WRITE(*,605)ihe4_ascii,TRIM(chain_rep),ihe4,TRIM(chain_don) WRITE(2,605)ihe4_ascii,TRIM(chain_rep),ihe4,TRIM(chain_don) 605 FORMAT('ARRET, l''indice de He4 ihe4= ,',i3,/, 1 'du modèle repris en ASCII : ',a,/, 2 'diffère de celui ihe4= ,',i3,/, 3 'du modèle à calculer avec les données du fichier : ',a) END SELECT CALL sortie('cesam 34') ENDIF ok=.TRUE. DO i=1,nchim ok=ok .AND. nom_elem(i) == nom_elem_ascii(i) ENDDO IF(.NOT.ok)THEN SELECT CASE(langue) CASE('english') WRITE(*,1011)nom_elem_ascii WRITE(*,1603)TRIM(chain_rep),nom_elem WRITE(*,1032)TRIM(chain_don) WRITE(2,1011)nom_elem_ascii WRITE(2,1603)TRIM(chain_rep),nom_elem WRITE(2,1032)TRIM(chain_don) 1603 FORMAT('read the input ASCII file : ',a,/, 1 'differ from the chemical elements : ',30a) CASE DEFAULT WRITE(*,11)nom_elem_ascii WRITE(*,603)TRIM(chain_rep),nom_elem WRITE(*,32)TRIM(chain_don) WRITE(2,11)nom_elem_ascii WRITE(2,603)TRIM(chain_rep),nom_elem WRITE(2,32)TRIM(chain_don) 603 FORMAT('du modèle en ASCII repris : ',a,/, 1 'diffèrent de ceux : ',30a) END SELECT CALL sortie('cesam 35') ENDIF c suite de la lecture du modèle ASCII ALLOCATE(bp(ne,n),tds(1,n),chim(nchim,n)) ; bp=0.d0 DO k=1,n READ(4,2005)bp(5,k),bp(1,k),bp(2,k),bp(3,k),bp(4,k),tds(1,k), 1 chim(1:nchim,k) 2005 FORMAT(5es22.15) ENDDO CLOSE(unit=4) c initialisation de la spline d'inter. de la comp. chim. c la composition chimique est tabulée par mole (modèle repris en ASCII) DO i=1,nchim chim(i,:)=chim(i,:)/nucleo(i) ENDDO 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 m^2/3 ou m^1/3 n_ch=n ALLOCATE(mc(n_ch),mct(n_ch+m_ch)) mc(:)=bp(5,:)**(2.d0/3.d0) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.FALSE.,mc(1),lq, 1 xchim,dxchim) IF(no_croiss)PRINT*,'Pb. at 1 in cesam' dim_ch=knotc-m_ch c si on utilise une grille fixe pour la composition chimique IF(grille_fixe)THEN ALLOCATE(mc_fixe(n_ch)) nc_fixe=n_ch ; mc_fixe=mc ENDIF c on passe en variables d'intégration c m--> m^2/3-1, p--> ln p, t--> ln t, r--> r^2, l-->l^2/3 c ou bien m--> m^1/3, p--> ln p, t--> ln t, r--> r, l--> l bp(1:2,1:n)=LOG(bp(1:2,1:n)) !pour Ptot et T IF(pturb)bp(Ipg,:)=bp(1,:) !Pturb IF(en_m23)THEN bp(3,1:n)=bp(3,1:n)**2 !pour R bp(4:5,1:n)=bp(4:5,1:n)**(2.d0/3.d0) !pour L et M ELSE bp(5,1:n)=bp(5,1:n)**(1.d0/3.d0) !pour M ENDIF bp(6,:)=1.d0 !pour psi (fictif) ord_qsp=2 !initialisation sur snoein c initialisation de la spline d'inter. de TdS sur m^1/3 ou m^2/3 n_tds=n ALLOCATE(x_tds(n_tds),xt_tds(n_tds+m_tds)) ; x_tds(:)=bp(5,:) xt_tds=0.d0 CALL bsp1dn(1,tds,x_tds,xt_tds,n_tds,m_tds,knot_tds,.FALSE., 1 x_tds(1),lq,ftds,dftds) IF(no_croiss)PRINT*,'Pb. at 22 in cesam' c bp(repris) --> spline par bsp1dn sur q --> fonction d'espacement qb --> bb(qb(i)) c + esp(i) --> spline par bsp1dn sur bb (sauf bb(6)) --> bb(6,i) --> spline par bsp1dn c sur bb complet --> base de noedif --> newspline bb--> bp c création de la spline d'interpolation des variables principales ALLOCATE(q(n),qt(n+ord_qsp)) ; q=(/ (i, i=1,n) /) CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.FALSE.,q(1),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 23 in cesam' c recherche de la fonction d'espacement (modèle repris en ASCII) ALLOCATE(esp(n)) DO i=1,n CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.TRUE.,q(i),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 2 in cesam' esp(i)=ctep*fqs(1)+ctet*fqs(2)+ctem*fqs(5) ENDDO c WRITE(*,2000)esp(1:8) ; PAUSE'esp' c on dispose les ncouches pour assurer une répartition approx. c uniforme de la fonction de répartition sur n_qs=n=n_ascii couches n_qs=n ALLOCATE(qb(n_qs)) CALL zoning(n,q,esp,n_qs,qb) !choix des nouveaux q dans qb c esp est désormais inutile (modèle repris en ASCII) DEALLOCATE(esp) c dans bb on se place aux n_qs points qb on calcule la nouvelle fonction de répartition ALLOCATE(bb(ne,n_qs),new_esp(1,n_qs)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.TRUE.,qb(i),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 3 in cesam' bb(1:ne,i)=fqs(1:ne) !; WRITE(*,2000)bb(1:ne,i) new_esp(1,i)=ctep*fqs(1)+ctet*fqs(2)+ctem*fqs(5) ENDDO c PAUSE'new_esp' c spline d'ordre 2 sur les nouveaux qb sur n_qs points pour new_esp qb=(/ (i, i=1,n_qs) /) ALLOCATE(qbt(n_qs+ord_qsp)) CALL bsp1dn(1,new_esp,qb,qbt,n_qs,2,knotb,.FALSE.,qb(1),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 24 in cesam' c dér. de la fonction de répartition dans bb(6,:) (modèle repris en ASCII) DO i=1,n_qs CALL bsp1dn(1,new_esp,qb,qbt,n_qs,2,knotb,.TRUE.,qb(i),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 4 in cesam' bb(6,i)=dfqs(1) ENDDO c on spline bb (modèle repris en ASCII) CALL bsp1dn(ne,bb,qb,qbt,n_qs,2,knotb,.FALSE.,qb(1),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 25 in cesam' c DO i=1,n_qs c CALL bsp1dn(ne,bb,qb,qbt,n_qs,2,knotb,.TRUE.,qb(i),lq,fqs,dfqs) c WRITE(*,2000)qb(i),fqs c ENDDO c PAUSE'bb' c dim. de la base de splines pour l'intégration equi. quasi. stat. dim_qs=(n_qs-1)*m_qs+r_qs ; ord_qs=m_qs+r_qs c vecteur nodal pour l'intégration sur n_qs points DEALLOCATE(q,qt) ALLOCATE(q(n_qs),qt(dim_qs+ord_qs)) ; q=(/ (i, i=1,n_qs) /) CALL noedif(q,n_qs,m_qs,r_qs,qt,knot) IF(no_croiss)THEN PRINT*,'ARRET 6 dans cesam / STOP 6 in cesam' ; CALL sortie('cesam 36') ENDIF c PRINT*,ord_qs,knot,dim_qs ; WRITE(*,2000)qt(1:knot) c PAUSE'noedif' c on place la solution initiale dans la base de q, bb ==> bp DEALLOCATE(bp) ; ALLOCATE(bp(ne,dim_qs)) CALL newspl(ne,qb,qbt,knotb,2,qt,knot,ord_qs,bb,bp) c PAUSE'newspl' c bb, qb, qbt sont désormais inutiles DEALLOCATE(bb,qb,qbt) c au centre R=L=M=0 bp(3:5,1)=0.d0 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,fqs,dfqs) c WRITE(*,2000)q(i),fqs c ENDDO c PAUSE'vérification' c initialisations et allocations diverses mstar=mtot c tabulation de r_stat et m_stat qui accélèrent la recherche dans le ssp. c inter valable pour m_stat=m^2/3 et r_stat=r^2 ou m_stat=m^1/3, r_stat=r ALLOCATE(r_stat(n_qs),m_stat(n_qs)) DO i=1,n_qs c(test) CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,fqs,dfqs) r_stat(i)=bp(3,m_qs*(i-1)+1) ; m_stat(i)=bp(5,m_qs*(i-1)+1) ENDDO c estimation provisoire de rstar (modèle repris en ASCII) rstar=r_stat(n_qs) ; IF(en_m23)rstar=SQRT(rstar) c vérification concernant la vitesse angulaire c pour la couche externe : force centifuge < force de gravité / 2. IF(wrot == 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1083) ; WRITE(2,1083) 1083 FORMAT('without rotation nor loss of angular momentum.') CASE DEFAULT WRITE(*,83) ; WRITE(2,83) 83 FORMAT('sans rotation ni perte de moment. cinétique.') END SELECT ELSE CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(n_qs),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 5 in cesam' IF(en_m23)THEN fqs(5)=SQRT(fqs(5))**3 ; fqs(3)=SQRT(fqs(3)) ELSE fqs(5)=ABS(fqs(5))**3 !m^1/3 ENDIF fqs(5)=fqs(5)*msol ; fqs(3)=fqs(3)*rsol w_max=SQRT(g*fqs(5)/fqs(3)**3)/2.d0 c si la rotation est entrée en kms/s w = v / R (modèle repris en ASCII) IF(unit == 'kms/s')THEN wrot=wrot*1.d5/fqs(3) ; w_rot=wrot ENDIF IF(wrot > w_max)THEN SELECT CASE(langue) CASE('english') WRITE(*,1400)wrot,w_max,TRIM(chain_don) 1400 FORMAT('STOP, the centrifugal force is too large w_rot=', 1 es10.3,' > w_max=',es10.3,/, 2 'decrease w_rot in the input file : ',a) CASE DEFAULT WRITE(*,400)wrot,w_max,TRIM(chain_don) 400 FORMAT('ARRET, la force centrifuge est trop grande w_rot=', 1 es10.3,' > w_max=',es10.3,/, 2 'réduire w_rot dans le fichier de données : ',a) END SELECT CALL sortie('cesam 37') ENDIF c rot_min : valeur négligeable, initialisation de rota IF(entre(3,5,Krot))THEN rot_min(1)=wrot*1.d-2 ; CALL initialise_rota ENDIF ENDIF DEALLOCATE(nom_elem_ascii) c gestion du premier pas temporel SELECT CASE(langue) CASE('english') WRITE(*,1058)dtp 1058 FORMAT('the time step of the continued model in ASCII, dtp=', 1 es10.3,/,'is recommanded for the first time step, OK?',/, 2 'enter o for yes , enter n for no') CASE DEFAULT WRITE(*,58)dtp 58 FORMAT('le pas temporel du modèle repris en ASCII, dtp=',es10.3,/, 1 'est conseillé pour le premier pas temporel, ok? o/n') END SELECT READ*,oui IF(oui == 'o')THEN dt=MIN(dtp,dtmax) ; dts=dt ELSE SELECT CASE(langue) CASE('english') WRITE(*,1059) 1059 FORMAT('enter the new time step, in units of 10^6 years') CASE DEFAULT WRITE(*,59) 59 FORMAT('entrer le nouveau pas temporel, unité : 10^6 ans') END SELECT READ*,dt IF(dt > dtmax)THEN dt=dtmax ; dts=dt c PRINT*,'en 1 dt=',dt,dtmax ; PAUSE'0' SELECT CASE(langue) CASE('english') WRITE(*,1112)dtmax CASE DEFAULT WRITE(*,112)dtmax END SELECT ENDIF ENDIF !oui /= 'o' IF(age+dt > agemax)THEN SELECT CASE(langue) CASE('english') WRITE(*,1053)age+dt,agemax CASE DEFAULT WRITE(*,53)age+dt,agemax END SELECT dt=MIN(dtmax,MAX(0.d0,agemax-age)) IF(dt <= 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1054) ; WRITE(2,1054) CASE DEFAULT WRITE(*,54) ; WRITE(2,54) END SELECT READ*,oui IF(oui == 'o')THEN list_sort=.TRUE. ELSE CALL sortie('cesam 38') ENDIF ENDIF !dt == 0.d0 ELSE SELECT CASE(langue) CASE('english') WRITE(*,1051)dt ; WRITE(2,1051)dt CASE DEFAULT WRITE(*,51)dt ; WRITE(2,51)dt END SELECT list_rep=.TRUE. !liste du modèle ENDIF !age+dt > agemax c initialisation des paramètres pour le vent et les chutes de planétoïdes CALL vent ; CALL planetoides c---------initialisation avec un modèle de ZAMS-------------------- c le modèle de ZAMS repris sera actualisé CASE(-2, 2) c modèle d'initialisation en binaire/ASCII, modèle de ZAMS c pour éviter de prendre un modèle d'atmosphère trop différent d'une c atmosphère de ZAMS, suppression du modèle d'atmosphère s'il existe INQUIRE(file=TRIM(chain_atm),exist=ok) IF(ok)CALL system('rm '//TRIM(chain_atm)) c que le modèle soit repris en binaire ou en ASCII, la spline des variables principales c est crée puis pépartie, on génèrera ensuite la composition chimique initiale I1: IF(un23 > 0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1008)nom_fich1 ; WRITE(2,1008)nom_fich1 1008 FORMAT('One chooses as a binary initial model : ',a) CASE DEFAULT WRITE(*,8)nom_fich1 ; WRITE(2,8)nom_fich1 8 FORMAT('On prend pour modèle initial en binaire : ',a) END SELECT OPEN(unit=4,form='unformatted',status='old',file=nom_fich1) READ(4)nep,m_qsp,n,knot,nchimp c PRINT*,nep,m_qsp,n,knot,nchimp ord_qsp=m_qsp+r_qs ; dim_qsp=knot-ord_qsp ALLOCATE(b(nep,dim_qsp),q(n),qt(knot),nom_elemp(nchimp)) REWIND(unit=4) READ(4)nep,m_qsp,n,knot,nchimp,n_ch,m_chp,knotc,Krotp,nrotp, 1 n_rot,m_rotp,knotr,n_tds,knot_tds,mtotp,alphap,w_rotp,lim_rop, 2 diffusionp,rot_solidp,precisionp,en_m23p,f_eosp,f_opap, 3 nom_ctesp,nom_pertmp,nom_pertwp,nom_tdetaup,nom_atmp,nom_convp, 4 nom_nucp,nom_nuc_cplp,nom_diffmp,nom_difftp,nom_diffwp, 5 nom_etatp,nom_opap,nom_elemp,b,q,qt CLOSE(unit=4) c PRINT*,n,ord_qsp,knot,nep,n_tds,knot_tds c WRITE(*,2000)b(1,1),q(1),qt(knot) ; PAUSE'lec' c DO i=1,n c CALL bsp1dn(nep,b,q,qt,n,ord_qsp,knot,.TRUE.,q(i),lq,fqs,dfqs) c WRITE(*,2000)q(i),fqs c PRINT*,i,fqs(5) c ENDDO c PAUSE'test' c on met modèle repris dans bp en l'adaptant, modèle de ZAMS ALLOCATE(bp(ne,dim_qsp)) c suppression/adjonction de la pression turbulente c si le modèle repris est avec/sans Pturb c et que le modèle à calculer est sans/avec Pturb c il faut enlever/ajouter la variable ln Pgaz c ne/nep : nb. inconnues du modèle à calculer/repris I2: IF(ne == nep)THEN bp=b !on a le même nombre d'inconnues ELSE I2 !changement du nb. inconnues I3: IF(ne == 6)THEN !suppression de Pgaz on passe de ne=7 à 6 inconnues SELECT CASE(langue) CASE('english') WRITE(*,1721) ; WRITE(2,721) CASE DEFAULT WRITE(*,721) ; WRITE(2,721) END SELECT bp=b(1:6,:) ELSE I3 !adjonction de Pgaz on passe de 6 à 7 inconnues SELECT CASE(langue) CASE('english') WRITE(*,1722) ; WRITE(2,722) CASE DEFAULT WRITE(*,722) ; WRITE(2,722) END SELECT bp(1:6,:)=b ; bp(Ipg,:)=b(1,:) ENDIF I3 !passage de 6 <---> 7 inconnues SELECT CASE(langue) CASE('english') WRITE(*,1723)nep ; WRITE(2,1723)ne CASE DEFAULT WRITE(*,723)nep ; WRITE(2,723)ne END SELECT ENDIF I2 !transfert et changement du nb. inconnues c b est désormais inutile DEALLOCATE(b) c on adapte le type variables du modèle repris à celles c du modèle à calculer, modèle de ZAMS ord_qs=ord_qsp I4: IF(en_m23p .NEQV. en_m23)THEN c m23 --> m13 I5: IF(en_m23p)THEN SELECT CASE(langue) CASE('english') WRITE(*,1225) ; WRITE(2,1225) 1225 FORMAT('change of variables: R^2==>R, L^2/3==>L, M^2/3==>M^1/3',/) CASE DEFAULT WRITE(*,225) ; WRITE(2,225) 225 FORMAT('chgt. de variables : R^2==>R, L^2/3==>L, M^2/3==>M^1/3',/) END SELECT bp(3,:)=SQRT(ABS(bp(3,:))) !R^2 ==> R bp(4,:)=SQRT(ABS(bp(4,:)))**3 !L^2/3 ==> L bp(5,:)=SQRT(ABS(bp(5,:))) !M^2/3 ==> M^1/3 c m13 --> m23 ELSE I5 SELECT CASE(langue) CASE('english') WRITE(*,1226) ; WRITE(2,1226) 1226 FORMAT('change of variables: R==>R^2, L==>L^2/3, M^1/3==>M^2/3',/) CASE DEFAULT WRITE(*,226) ; WRITE(2,226) 226 FORMAT('chgt. de variables : R==>R^2, L==>L^2/3, M^1/3==>M^2/3',/) END SELECT bp(3,:)=bp(3,:)**2 !R ==> R^2 bp(4,:)=ABS(bp(4,:))**(2.d0/3.d0) !L ==> L^2/3 bp(5,:)=bp(5,:)**2 !M^1/3 ==> M^2/3 ENDIF I5 c pas de changement de variable ELSE I4 SELECT CASE(langue) CASE('english') WRITE(*,1227)en_m23p,en_m23 ; WRITE(2,1227)en_m23p,en_m23 1227 FORMAT('no change of variables: en_m23p: ',l1,', en_m23: ',l1,/) CASE DEFAULT WRITE(*,227)en_m23p,en_m23 ; WRITE(2,227)en_m23p,en_m23 227 FORMAT('pas de chgt. de variables, en_m23p: ',l1,', en_m23: ',l1,/) END SELECT ENDIF I4 c modèle d'initialisation en ASCII de ZAMS ELSE I1 OPEN(unit=4,form='formatted',status='old',file=nom_fich1) c détermination du nombre de couches n=0 DO READ(4,*,END=15)fqs(1) ; n=n+1 ENDDO 15 REWIND(unit=4) c PRINT*,n ; WRITE(*,2000)mtot ; PAUSE'reprise ASCII' c les masses sont en 1-M/Mtot ALLOCATE(bp(ne,n)) ; bp=0.d0 DO i=1,n READ(4,*)bp(5,i),bp(1:4,i) ; bp(5,i)=(1.d0-bp(5,i))*mtot ENDDO CLOSE(unit=4) ! PAUSE'lect' SELECT CASE(langue) CASE('english') WRITE(*,1016)nom_fich1,n ; WRITE(2,1016)nom_fich1,n 1016 FORMAT('initializing model in ASCII : ',a,/, 1 'number of shells : ',i4) CASE DEFAULT WRITE(*,16)nom_fich1,n ; WRITE(2,16)nom_fich1,n 16 FORMAT('modèle d''initialisation en ASCII : ',a,/, 1 'nombre de couches :',i4) END SELECT c retournement 1 au centre, n à la surface bp(:,1:n:+1)=bp(:,n:1:-1) c on passe en variables d'intégration c m--> m^2/3, p--> ln p, t--> ln t, r--> r^2, l-->l^2/3 c ou bien m--> m^1/3, p--> ln p, t--> ln t, r--> r, l-->l bp(1:2,1:n)=LOG(bp(1:2,1:n)) !pour Ptot et T IF(pturb)bp(Ipg,:)=bp(1,:) !Pturb I6: IF(en_m23)THEN bp(3:5,1:n)=bp(3:5,1:n)**2 !pour R, L et M bp(4:5,1:n)=bp(4:5,1:n)**(1.d0/3.d0) !pour L et M ELSE I6 bp(5,1:n)=bp(5,1:n)**(1.d0/3.d0) !pour M ENDIF I6 bp(6,:)=1.d0 !pour psi (fictif) ord_qsp=2 !initialisation sur snoein ALLOCATE(q(n),qt(n+ord_qsp)) ; q=(/ (i, i=1,n) /) CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.FALSE.,q(1),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 26 in cesam' ENDIF I1 !initialisation binaire/ASCII c on fixe le nb. de couches du modèle d'initialisation à ncouche=300 c 300 étant un bon ordre de grandeur pour un modèle solaire c ce nombre sera ensuite ajusté de facon à ce que dQ/dq = psi0 c on détermine les ncouche points de masse m(ncouche) pour c assurer une répartition approximativement uniforme en c ctep lnP + ctet lnT + cter zeta + ctem mu c recherche de la fonction d'espacement, modèle de ZAMS, on ne ne peut pas c tenir compte de Lnuc la composition chimique étant inconnue à ce stade c le modèle d'initialisation est dans bp(ne,n) si on a repris c en ASCII ou dans b(ne,dim_qsp) si on a repris en binaire, comme c on ne fait qu'interpoler le fait que n puisse être différent de c dim_qsp n'a pas d'importance ALLOCATE(esp(n)) DO i=1,n CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.TRUE.,q(i),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 6 in cesam' esp(i)=ctep*fqs(1)+ctet*fqs(2)+ctem*fqs(5) ENDDO c PRINT*,'ord_qsp,knot,n',ord_qsp,knot,n c WRITE(*,2000)ctel,ctep,ctet,ctem,cter ; WRITE(*,2000)esp c WRITE(*,2000)q ; WRITE(*,2000)qt ; PAUSE'esp' c WRITE(*,2000)esp(1:8) c on dispose les ncouches pour assurer une répartition approx. c uniforme de la fonction de répartition sur n_qs couches n_init=MAX(n_init,NINT((esp(n)-esp(1))/psi0)) n_qs=n_init ; ALLOCATE(qb(n_qs)) c PRINT*,'n_init 1 =',n_init ; PAUSE'n_init' CALL zoning(n,q,esp,n_qs,qb) !choix des nouveaux q dans qb c PRINT*,'anciens q',n ; WRITE(*,2000)q(1:n) ; PRINT*,'esp' c WRITE(*,2000)esp(1:n) ; PRINT*,'new q' ; WRITE(*,2000)qb(1:n_qs) c PAUSE'après zoning' c esp est désormais inutile DEALLOCATE(esp) c dans bb on se place aux n_qs points qb c on calcule la nouvelle fonction de répartition, modèle de ZAMS ALLOCATE(bb(ne,n_qs),new_esp(1,n_qs)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.TRUE.,qb(i),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 7 in cesam' bb(1:ne,i)=fqs(1:ne) !; WRITE(*,2000)bb(1:ne,i) new_esp(1,i)=ctep*fqs(1)+ctet*fqs(2)+ctem*fqs(5) ENDDO c PAUSE'new_esp' c spline d'ordre 2 sur les nouveaux qb sur n_qs points pour new_esp ALLOCATE(qbt(n_qs+ord_qsp)) ; qb=(/ (i, i=1,n_qs) /) CALL bsp1dn(1,new_esp,qb,qbt,n_qs,2,knotb,.FALSE.,qb(1),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 27 in cesam' c dérivée de la fonction de répartition dans bb(6,:) DO i=1,n_qs CALL bsp1dn(1,new_esp,qb,qbt,n_qs,2,knotb,.TRUE.,qb(i),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 8 in cesam' bb(6,i)=dfqs(1) ENDDO c WRITE(*,2000)bb(5,1:4),bb(6,1:4) ; PAUSE'bb' c on spline bb CALL bsp1dn(ne,bb,qb,qbt,n_qs,2,knotb,.FALSE.,qb(1),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 28 in cesam' c vérification c DO i=1,n_qs c CALL bsp1dn(ne,bb,qb,qbt,n_qs,2,knotb,.TRUE.,qb(i),lq,fqs,dfqs) c WRITE(*,2000)qb(i),fqs(1:5) !; PRINT*,qb(i),fqs(5) c ENDDO c PAUSE'bb' c dimension de la base de splines pour l'intégration equi. quasi. stat., c modèle de ZAMS dim_qs=(n_qs-1)*m_qs+r_qs ; ord_qs=m_qs+r_qs c vecteur nodal pour l'intégration sur n_qs points DEALLOCATE(q,qt) ALLOCATE(q(n_qs),qt(dim_qs+ord_qs)) ; q=(/ (i, i=1,n_qs) /) CALL noedif(q,n_qs,m_qs,r_qs,qt,knot) IF(no_croiss)THEN PRINT*,'ARRET 10 dans cesam / STOP 10 in cesam' ; CALL sortie('cesam 39') ENDIF c PRINT*,ord_qs,knot,dim_qs ; WRITE(*,2000)qt(1:knot) ; PAUSE'noedif' c on place la solution initiale dans la base de q, bb ==> bp, modèle de ZAMS DEALLOCATE(bp) ; ALLOCATE(bp(ne,dim_qs)) CALL newspl(ne,qb,qbt,knotb,2,qt,knot,ord_qs,bb,bp) c PAUSE'newspl' c bb, qb, qbt sont désormais inutiles DEALLOCATE(bb,qb,qbt) c au centre R=L=M=0 bp(3,1)=0.d0 ; bp(4,1)=0.d0 ; bp(5,1)=0.d0 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,fqs,dfqs) c WRITE(*,2000)q(i),f(1:ne) c PRINT*,i,fqs(5) c ENDDO c PRINT*,'ord_qs=',ord_qs c PAUSE'vérification' c modèle de ZAMS c initialisation du vecteur de composition chimique c pour la couche externe : force centifuge < force de gravité / 2. IF(wrot == 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1083) ; WRITE(2,1083) CASE DEFAULT WRITE(*,83) ; WRITE(2,83) END SELECT ELSE CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(n_qs),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 9 in cesam' IF(en_m23)THEN fqs(5)=SQRT(fqs(5))**3 ; fqs(3)=SQRT(fqs(3)) ELSE fqs(5)=ABS(fqs(5))**3 ENDIF fqs(5)=fqs(5)*msol ; fqs(3)=fqs(3)*rsol w_max=SQRT(g*fqs(5)/fqs(3)**3)/2.d0 c si la rotation a été entrée en kms/s w = v / R IF(unit == 'kms/s')THEN wrot=wrot*1.d5/fqs(3) ; w_rot=wrot ENDIF IF(wrot > w_max)THEN SELECT CASE(langue) CASE('english') WRITE(*,1400)wrot,w_max,TRIM(chain_don) CASE DEFAULT WRITE(*,400)wrot,w_max,TRIM(chain_don) END SELECT CALL sortie('cesam 40') ENDIF c rot_min : valeur négligeable IF(entre(3,5,Krot))rot_min(1)=wrot*1.d-2 ENDIF c appel d'initialisation pour tabulation des réactions nucléaires c allocations fictives, modèle de ZAMS ALLOCATE(comp(0),dcomp(0),jac(0,0),ex(0)) CALL nuc(1.5d+07,1.5d+02,comp,dcomp,jac,.FALSE.,0, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c détermination des abondances initiales, modèle de ZAMS DEALLOCATE(comp) ; ALLOCATE(comp(nchim)) CALL nuc(1.5d+07,1.5d+02,comp,dcomp,jac,.FALSE.,1, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c une trop grande abondance initiale de Lithium est incompatible c avec un modèle de ZAMS homogène, modèle de ZAMS Bli: DO i=1,nchim IF(nom_elem(i)(1:2) == 'Li')THEN IF(comp(i) > 1.d-15)THEN SELECT CASE(langue) CASE('english') WRITE(*,1014)nom_elem(i),comp(i),EXP(ln_Tli) WRITE(2,1014)nom_elem(i),comp(i),EXP(ln_Tli) 1014 FORMAT('With an initial abundance of ',a4,'=',es10.3,/, 1 'the ZAMS might not be able to converge',/, 2 'does one try without lithium for t >',es10.3, 3 'enter o for yes, enter n for no') CASE DEFAULT WRITE(*,14)nom_elem(i),comp(i),EXP(ln_Tli) WRITE(2,14)nom_elem(i),comp(i),EXP(ln_Tli) 14 FORMAT('Avec une abondance initiale de ',a4,'=',es10.3,/, 1 'le modèle de ZAMS risque de ne pas converger',/, 2 'tente-on en supprimant le lithium pour T >',es10.3,/, 3 'entrer o/n.') END SELECT READ*,oui IF(oui == 'o')THEN SELECT CASE(langue) CASE('english') WRITE(*,1025)EXP(ln_Tli) ; WRITE(2,1025)EXP(ln_Tli) 1025 FORMAT('one tries, without lithium for T > ',es10.3) CASE DEFAULT WRITE(*,25)EXP(ln_Tli) ; WRITE(2,25)EXP(ln_Tli) 25 FORMAT('on tente, en supprimant le lithium pour T > ',es10.3) END SELECT EXIT Bli ELSE SELECT CASE(langue) CASE('english') WRITE(*,1026) ; WRITE(2,1026) 1026 FORMAT('STOP, it is suggested :',/, 1 'either to compute the PMS or to initialize Li7 < 1.d-15') CASE DEFAULT WRITE(*,26) ; WRITE(2,26) 26 FORMAT('ARRET, suggestions :',/, 1 '- calculer la PMS ou encore initialiser Li7 < 1.d-15') END SELECT CALL sortie('cesam 41') ENDIF ENDIF ENDIF ENDDO Bli c la rotation, modèle de ZAMS WRITE(*,17)thw(Krot) ; WRITE(2,17)thw(Krot) SELECT CASE(Krot) CASE(0,1,2,6) nrot=0 ; n_rot=0 ; knotr=0 ALLOCATE(rota(nrot,n_rot),mrot(n_rot),mrott(knotr)) END SELECT SELECT CASE(langue) CASE('english') WRITE(*,1023)wrot ; WRITE(2,1023)wrot CASE DEFAULT WRITE(*,23)wrot ; WRITE(2,23)wrot END SELECT c tabulation de r_stat et m_stat qui accélèrent la recherche dans le ssp. c inter pour variables m_stat=m^2/3 et r_stat=r^2 ou m_stat=m^1/3 et r_stat=r ALLOCATE(r_stat(n_qs),m_stat(n_qs)) DO i=1,n_qs c(test) CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,fqs,dfqs) r_stat(i)=bp(3,m_qs*(i-1)+1) ; m_stat(i)=bp(5,m_qs*(i-1)+1) ENDDO 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 n_ch=n_qs ALLOCATE(mc(n_ch)) IF(en_m23)THEN mc=m_stat ELSE mc=m_stat**2 ENDIF c vérification de la stricte croissance des mc, suppression des doubles c réallocation de mc si nécessaire DO i=1,n_ch-1 IF(mc(i) >= mc(i+1))mc(i+1)=mc(i) ENDDO CALL delete_doubles(n_ch,mc) c initialisation de la composition chimique c sur la ZAMS annulation de l'abondance initiale du lithium ALLOCATE(chim(nchim,n_ch)) DO i=1,n_ch chim(:,i)=comp(:) IF(iLi7 /= 0 .AND. bp(2,m_qs*(i-1)+1) > ln_Tli)chim(iLi7,i)=li_ini ENDDO c initialisation de la spline d'inter. de la comp. chim., modèle de ZAMS ALLOCATE(dxchim(nchim),mct(n_ch+m_ch),xchim(nchim)) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.FALSE.,mc(1),lq, 1 xchim,dxchim) IF(no_croiss)PRINT*,'Pb. at 29 in cesam' dim_ch=knotc-m_ch ; DEALLOCATE(dxchim,xchim) c initialisations et allocations diverses mstar=mtot c initialisation de rota, modèle de ZAMS IF(entre(3,5,Krot))CALL initialise_rota c si on utilise une grille fixe pour la composition chimique IF(grille_fixe)THEN ALLOCATE(mc_fixe(n_qs)) ; nc_fixe=n_qs ; mc_fixe=m_stat IF(.NOT.en_m23)mc_fixe=m_stat**3 ENDIF c estimation provisoire de rstar rstar=r_stat(n_qs) ; IF(en_m23)rstar=SQRT(rstar) c PAUSE'avant call resout' c initialisation de l'âge et du pas temporel pour le modèle ZAMS c pour écriture sur (26), on impose dt=0 (et non pas dt_minimum) age=0.d0 ; dt=0.d0 c calcul du modèle de ZAMS homogène------------------------------ CALL resout(2,dt,dts) ; dt=0.d0 ; dts=0.d0 c WRITE(*,2000)bp(6,1) ; PAUSE'bp(6,1)' c on garde le modèle de ZAMS d'âge 0 sur le fichier xxxx_B.hom model_num=0 c écriture en binaire du fichier mon_modele_B.hom modelb=TRIM(nom_fich2)//'_B.hom' c PRINT*,SIZE(mrott),knotr c WRITE(*,2000)mrott ; PAUSE'mrott' OPEN(unit=26,form='unformatted',status='unknown',file=TRIM(modelb)) WRITE(26)ne,m_qs,n_qs,knot,nchim,n_ch,m_ch,knotc,Krot,nrot,n_rot, 1 m_rot,knotr,n_tds,knot_tds,mtot,alpha,w_rot,lim_ro,diffusion, 2 rot_solid,precisione,en_m23,f_eos,f_opa, 3 nom_ctes,nom_pertm,nom_pertw,nom_tdetau,nom_atm,nom_conv,nom_nuc, 4 nom_nuc_cpl,nom_diffm,nom_difft,nom_diffw,nom_etat,nom_opa, 5 nom_elem,bp,q,qt,chim,mc,mct,rota,mrot,mrott,tds,x_tds,xt_tds, 6 m_stat,r_stat,m_zc,r_zc,r_ov,age,dt,dts,mstar,rstar,mw_tot,wrot, 7 jlim,lconv,lim,model_num CLOSE(unit=26) WRITE(*,46)TRIM(modelb) 46 FORMAT('Formation ou remplacement du fichier binaire: ',1x,a,/) c PRINT*,all_rep,model_num,modelb,n_ch,m_ch,knotc c PAUSE'.hom' c le numéro du modèle sera augmenté d'une unité ligne ~4155 model_num=-1 c pas temporel initial pour évolution à partir de la ZAMS dtp=0.d0 ; dt=MIN(dt0,dtmax) ; dts=dt c pour l'entête du modèle du fichier mon_modèle.lis list_zams=.TRUE. c initialisation des paramètres pour le vent et les chutes de planétoïdes CALL vent ; CALL planetoides c initialisation d'un modèle PMS CASE(-3, 3) c pour éviter de prendre un modèle d'atmosphère trop différent d'une c atmosphère de PMS, suppression du modèle d'atmosphère s'il existe INQUIRE(file=TRIM(chain_atm),exist=ok) IF(ok)CALL system('rm '//TRIM(chain_atm)) c modèle d'initialisation en binaire/ASCII IF(un23 > 0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1008)nom_fich1 ; WRITE(2,1008)nom_fich1 CASE DEFAULT WRITE(*,8)nom_fich1 ; WRITE(2,8)nom_fich1 END SELECT OPEN(unit=4,form='unformatted',status='old',file=nom_fich1) READ(4)nep,m_qsp,n,knot,nchimp c PRINT*,nep,m_qsp,n,knot,nchimp ord_qsp=m_qsp+r_qs ; dim_qsp=knot-ord_qsp ALLOCATE(b(nep,dim_qsp),q(n),qt(knot),nom_elemp(nchimp)) REWIND(unit=4) READ(4)nep,m_qsp,n,knot,nchimp,n_ch,m_chp,knotc,Krotp,nrotp, 1 n_rot,m_rotp,knotr,n_tds,knot_tds,mtotp,alphap,w_rotp,lim_rop, 2 diffusionp,rot_solidp,precisionp,en_m23p,f_eosp,f_opap, 3 nom_ctesp,nom_pertmp,nom_pertwp,nom_tdetaup,nom_atmp,nom_convp, 4 nom_nucp,nom_nuc_cplp,nom_diffmp,nom_difftp,nom_diffwp, 5 nom_etatp,nom_opap,nom_elemp,b,q,qt CLOSE(unit=4) c PRINT*,n,ord_qsp,knot,nep,knot_tds c WRITE(*,2000)b(1,1),q(1),qt(knot) ; PAUSE'lec' c DO i=1,n c CALL bsp1dn(nep,b,q,qt,n,ord_qsp,knot,.TRUE.,q(i),lq,fqs,dfqs) c WRITE(*,2000)q(i),fqs ; PAUSE'1' c ENDDO c PAUSE'test' c on met modèle repris dans bp en l'adaptant ALLOCATE(bp(ne,dim_qsp)) c suppression/adjonction de la pression turbulente c si le modèle repris est avec/sans Pturb c et que le modèle à calculer est sans/avec Pturb c il faut enlever/ajouter la variable ln Pgaz c ne/nep : nb. inconnues du modèle à calculer/repris IF(ne == nep)THEN bp=b !on a le même nombre d'inconnues ELSE !changement du nb. inconnues IF(ne == 6)THEN !suppression de Pgaz on passe de ne=7 --> 6 SELECT CASE(langue) CASE('english') WRITE(*,1721) ; WRITE(2,721) CASE DEFAULT WRITE(*,721) ; WRITE(2,721) END SELECT bp=b(1:6,:) ELSE !adjonction de Pgaz on passe de 6 à 7 inconnues SELECT CASE(langue) CASE('english') WRITE(*,1722) ; WRITE(2,722) CASE DEFAULT WRITE(*,722) ; WRITE(2,722) END SELECT bp(1:6,:)=b ; bp(Ipg,:)=b(1,:) ENDIF !passage de 6 <---> 7 inconnues SELECT CASE(langue) CASE('english') WRITE(*,1723)nep ; WRITE(2,1723)ne CASE DEFAULT WRITE(*,723)nep ; WRITE(2,723)ne END SELECT ENDIF !transfert et changement du nb. inconnues c b est désormais inutile DEALLOCATE(b) c on adapte le type variables du modèle repris à celles du modèle à calculer IF(en_m23p .NEQV. en_m23)THEN IF(en_m23p)THEN SELECT CASE(langue) CASE('english') WRITE(*,1225) ; WRITE(2,1225) CASE DEFAULT WRITE(*,225) ; WRITE(2,225) END SELECT bp(3,:)=SQRT(ABS(bp(3,:))) !R^2 ==> R bp(4,:)=SQRT(ABS(bp(4,:)))**3 !L^2/3 ==> L bp(5,:)=SQRT(ABS(bp(5,:))) !M^2/3 ==> M^1/3 ELSE SELECT CASE(langue) CASE('english') WRITE(*,1226) ; WRITE(2,1226) CASE DEFAULT WRITE(*,226) ; WRITE(2,226) END SELECT bp(3,4:)=bp(3,4:)**2 !R ==> R^2, L ==> L^2 bp(4,:)=bp(4,:)**(1.d0/3.d0) !L^2 ==> L^2/3 bp(5,:)=bp(5,:)**2 !M^1/3 ==> M^2/3 ENDIF ENDIF c DO i=1,n c CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.TRUE.,q(i),lq,fqs,dfqs) c WRITE(*,2000)fqs c ENDDO c PAUSE'reprise' c modèle d'initialisation de PMS en ASCII ELSE OPEN(unit=4,form='formatted',status='old',file=nom_fich1) c détermination du nombre de couches n=0 DO READ(4,*,END=160)fqs(1) ; n=n+1 ENDDO 160 REWIND(unit=4) c PRINT*,n ; PAUSE'nb.couches' c les masses sont en 1-M/Mtot ALLOCATE(bp(ne,n)) ; bp=0.d0 DO i=1,n READ(4,*)bp(5,i),bp(1:4,i) bp(5,i)=(1.d0-bp(5,i))*mtot c WRITE(*,2000)bp(1:5,i) ENDDO CLOSE(unit=4) SELECT CASE(langue) CASE('english') WRITE(*,1016)nom_fich1,n ; WRITE(2,1016)nom_fich1,n CASE DEFAULT WRITE(*,16)nom_fich1,n ; WRITE(2,16)nom_fich1,n END SELECT c retournement 1 au centre, n à la surface bp(:,1:n:+1)=bp(:,n:1:-1) c on passe en variables d'intégration c m--> m^2/3-1, p--> ln p, t--> ln t, r--> r^2, l--> l^2/3 c ou bien m --> m^1/3, p--> ln p, t--> ln t, r--> r, l--> l bp(1:2,1:n)=LOG(bp(1:2,1:n)) !pour Ptot et T IF(pturb)bp(Ipg,:)=bp(1,:) !Pgas avec Pturb IF(en_m23)THEN bp(3,1:n)=bp(3,1:n)**2 !pour R bp(4:5,1:n)=bp(4:5,1:n)**(2.d0/3.d0) !pour L et M ELSE bp(5,1:n)=bp(5,1:n)**(1.d0/3.d0) !pour M ENDIF bp(6,:)=1.d0 !pour psi (fictif) ord_qsp=2 !initialisation sur snoein ALLOCATE(q(n),qt(n+ord_qsp)) ; q=(/ (i, i=1,n) /) CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.FALSE.,q(1),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 30 in cesam' ENDIF !initialisation binaire/ASCII c on fixe le nb de couches du modèle d'initialisation à ncouche=300 c 300 étant un bon ordre de grandeur pour un modèle solaire c ce nombre sera ensuite ajusté de facon à ce que dQ/dq = psi0 c on détermine les ncouche points de masse m(ncouche) pour c assurer une répartition approximativement uniforme en c ctep lnP + ctet lnT + cter zeta + ctem mu c recherche de la fonction d'espacement ALLOCATE(esp(n)) DO i=1,n CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.TRUE.,q(i),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 11 in cesam' esp(i)=ctep*fqs(1)+ctet*fqs(2)+ctem*fqs(5) ENDDO c PRINT*,'ord_qsp,knot,n',ord_qsp,knot,n c WRITE(*,2000)ctel,ctep,ctet,ctem,cter ; WRITE(*,2000)esp c WRITE(*,2000)q ; WRITE(*,2000)qt ; PAUSE'esp' c on dispose les ncouches pour assurer une répartition approx. c uniforme de la fonction de répartition sur n_qs couches n_qs=n_init ; ALLOCATE(qb(n_qs)) CALL zoning(n,q,esp,n_qs,qb) !choix des nouveaux q dans qb c PRINT*,'anciens q',n ; WRITE(*,2000)q(1:n) ; PRINT*,'esp' c WRITE(*,2000)esp(1:n) ; PRINT*,'new q' ; WRITE(*,2000)qb(1:n_qs) c PAUSE'après zoning' c esp est désormais inutile DEALLOCATE(esp) c dans bb on se place aux n_qs points qb c on calcule la nouvelle fonction de répartition ALLOCATE(bb(ne,n_qs),new_esp(1,n_qs)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n,ord_qsp,knot,.TRUE.,qb(i),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 12 in cesam' bb(1:ne,i)=fqs(1:ne) !; WRITE(*,2000)bb(1:ne,i) new_esp(1,i)=ctep*fqs(1)+ctet*fqs(2)+ctem*fqs(5) ENDDO c PAUSE'new_esp' c spline d'ordre 2 sur les nouveaux qb sur n_qs points pour new_esp ALLOCATE(qbt(n_qs+ord_qsp)) ; qb=(/ (i, i=1,n_qs) /) CALL bsp1dn(1,new_esp,qb,qbt,n_qs,2,knotb,.FALSE.,qb(1), 1 lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 31 in cesam' c dérivée de la fonction de répartition dans bb(6,:) DO i=1,n_qs CALL bsp1dn(1,new_esp,qb,qbt,n_qs,2,knotb,.TRUE.,qb(i), 1 lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 13 in cesam' bb(6,i)=dfqs(1) ENDDO c on spline bb CALL bsp1dn(ne,bb,qb,qbt,n_qs,2,knotb,.FALSE.,qb(1),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 32 in cesam' c DO i=1,n_qs c CALL bsp1dn(ne,bb,qb,qbt,n_qs,2,knotb,.TRUE.,qb(i),lq,fqs,dfqs) c WRITE(*,2000)qb(i),fqs c ENDDO c PAUSE'bb' c dimension de la base de splines pour l'intégration equi. quasi. stat. dim_qs=(n_qs-1)*m_qs+r_qs ; ord_qs=m_qs+r_qs c vecteur nodal pour l'intégration sur n_qs points DEALLOCATE(q,qt) ALLOCATE(q(n_qs),qt(dim_qs+ord_qs)) ; q=(/ (i, i=1,n_qs) /) CALL noedif(q,n_qs,m_qs,r_qs,qt,knot) IF(no_croiss)THEN PRINT*,'ARRET 15 dans cesam / STOP 15 in cesam' ; CALL sortie('cesam 42') ENDIF c PRINT*,ord_qs,knot,dim_qs ; WRITE(*,2000)qt(1:knot) c PAUSE'noedif' c on place la solution initiale dans la base de q, bb ==> bp DEALLOCATE(bp) ; ALLOCATE(bp(ne,dim_qs)) CALL newspl(ne,qb,qbt,knotb,2,qt,knot,ord_qs,bb,bp) c PAUSE'newspl' c bb, qb, qbt sont désormais inutiles DEALLOCATE(bb,qb,qbt) c au centre R=L=M=0 bp(3,1)=0.d0 ; bp(4,1)=0.d0 ; bp(5,1)=0.d0 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,fqs,dfqs) c WRITE(*,2000)q(i),fqs c ENDDO c PAUSE'vérification' c initialisation du vecteur de rotation c pour la couche externe : force centifuge < force de gravité / 2. IF(wrot == 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1083) ; WRITE(2,1083) CASE DEFAULT WRITE(*,83) ; WRITE(2,83) END SELECT ELSE CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(n_qs), 1 lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 14 in cesam' IF(en_m23)THEN fqs(5)=SQRT(fqs(5))**3 ; fqs(3)=SQRT(fqs(3)) ELSE fqs(5)=ABS(fqs(5))**3 ENDIF fqs(5)=fqs(5)*msol ; fqs(3)=fqs(3)*rsol w_max=SQRT(g*fqs(5)/fqs(3)**3)/2.d0 c si la rotation est entrée en kms/s w = v / R IF(unit == 'kms/s')THEN wrot=wrot*1.d5/fqs(3) ; w_rot=wrot ENDIF IF(wrot > w_max)THEN SELECT CASE(langue) CASE('english') WRITE(*,1400)wrot,w_max,TRIM(chain_don) CASE DEFAULT WRITE(*,400)wrot,w_max,TRIM(chain_don) END SELECT CALL sortie('cesam 43') ENDIF c rot_min : valeur négligeable IF(entre(3,4,Krot))rot_min(1)=wrot*1.d-2 ENDIF c appel d'initialisation pour tabulation des réactions nucléaires c allocations fictives ALLOCATE(comp(0),dcomp(0),jac(0,0),ex(0)) CALL nuc(1.5d+07,1.5d+02,comp,dcomp,jac,.FALSE.,0, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c détermination des abondances initiales DEALLOCATE(comp) ; ALLOCATE(comp(nchim)) CALL nuc(1.5d+07,1.5d+02,comp,dcomp,jac,.FALSE.,1, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) c la rotation WRITE(*,17)thw(Krot) ; WRITE(2,17)thw(Krot) SELECT CASE(Krot) CASE(0,1,2,6) nrot=0 ; n_rot=0 ; knotr=0 ALLOCATE(rota(nrot,n_rot),mrot(n_rot),mrott(knotr)) END SELECT SELECT CASE(langue) CASE('english') WRITE(*,1023)wrot ; WRITE(2,1023)wrot CASE DEFAULT WRITE(*,23)wrot ; WRITE(2,23)wrot END SELECT 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 lagragien ou en eulérien n_ch=n_qs ALLOCATE(chim(nchim,n_ch),mc(n_ch)) DO i=1,n_qs CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,fqs,dfqs) IF(no_croiss)PRINT*,'Pb. at 15 in cesam' c WRITE(*,2000)f(1:ne) fqs(5)=MAX(fqs(5),0.d0) IF(en_m23)THEN mc(i)=fqs(5) !interpolation en m^2/3 ELSE mc(i)=fqs(5)**2 ENDIF chim(1:nchim,i)=comp(:) ENDDO c initialisation de la spline d'inter. de la comp. chim. ALLOCATE(dxchim(nchim),mct(n_ch+m_ch),xchim(nchim)) CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,knotc,.FALSE.,mc(1),lq, 1 xchim,dxchim) IF(no_croiss)PRINT*,'Pb. at 34 in cesam' dim_ch=knotc-m_ch ; DEALLOCATE(dxchim,xchim) c initialisations et allocations diverses mstar=mtot c tabulation de r_stat et m_stat qui accélèrent la recherche dans le ssp. inter c pour variables m_stat=m^2/3 et r_stat=r^2 ou m_stat=m^1/3 et r_stat=r ALLOCATE(r_stat(n_qs),m_stat(n_qs)) DO i=1,n_qs c(test) CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,fqs,dfqs) r_stat(i)=bp(3,m_qs*(i-1)+1) ; m_stat(i)=bp(5,m_qs*(i-1)+1) ENDDO c initialisation de rota, modèle de PMS IF(entre(3,5,Krot))CALL initialise_rota c si on utilise une grille fixe pour la composition chimique IF(grille_fixe)THEN ALLOCATE(mc_fixe(n_qs)) ; nc_fixe=n_qs ; mc_fixe=m_stat IF(.NOT.en_m23)mc_fixe=m_stat**2 ENDIF c estimation provisoire de rstar rstar=r_stat(n_qs) ; IF(en_m23)rstar=SQRT(rstar) c initialisation de l'âge et du pas temporel pour le modèle initial c et calcul du premier modèle de PMS homogène c le premier pas temporel en PMS étant petit on n'en modifie pas la c si on diffuse le moment cinétique (Krot=3, 4,et 5) age=0.d0 ; dt=0.d0 ; dts=dt c on remplace la routine de réac. thermonucléaire lue dans le c fichier de données mon_modèle.don par iben nom_save=nom_nuc ; nom_nuc='iben' SELECT CASE(langue) CASE('english') WRITE(*,1800) 1800 FORMAT('-----------------------------------',/, 1 'calculation of the initial PMS model') CASE DEFAULT WRITE(*,800) 800 FORMAT('-----------------------------------',/, 1 'calcul du modèle initial de PMS') END SELECT B3: DO !1-iere boucle infinie pour le choix de c_iben SELECT CASE(langue) CASE('english') WRITE(*,1801) 1801 FORMAT('enter the contraction constant c_iben, examples : ',/, 1 'for Tc=100 000 K, entrer c_iben= 2.d-2',/, 2 'for Tc=500 000 K, enter c_iben= 5.d-4',/, 3 'for Tc=1 000 000 K, enter c_iben= 8.d-5',/, 4 'for M*>10 Msun, enter c_iben= 4.d-2') CASE DEFAULT WRITE(*,801) 801 FORMAT('entrer la constante de contraction c_iben, exemples : ',/ 1 'pour Tc=100 000 K, entrer c_iben= 2.d-2',/, 2 'pour Tc=500 000 K, entrer c_iben= 5.d-4',/, 3 'pour Tc=1 000 000 K, entrer c_iben= 8.d-5'/, 4 'pour M*>5~10 Msol ou si échec, entrer c_iben= 4.d-2') END SELECT READ*,c_iben c première résolution du Pb. quasi-static PMS avec c_iben------------ c PAUSE'avant resout' CALL resout(3,dt,dts) B2: DO !2-sde boucle infinie pour le choix de c_iben rp1=bp(3,m_qs*(n_qs-1)+1) ; lp1=bp(4,m_qs*(n_qs-1)+1) IF(en_m23)THEN rp1=SQRT(rp1) ; lp1=SQRT(lp1**3) ENDIF tp1=EXP(bp(2,1)) SELECT CASE(langue) CASE('english') WRITE(*,1055)c_iben,rp1,lp1,tp1 1055 FORMAT('Pre main sequence first model : C=',es10.3,/, 1 'Rext=',es10.3,', Lext=',es10.3,', Tc=',es10.3,/, 2 'ok? enter o for yes , n for no, to stop enter s') CASE DEFAULT WRITE(*,55)c_iben,rp1,lp1,tp1 55 FORMAT('Pré Séquence Principale premier modèle : C=',es10.3,/, 1 'Rext=',es10.3,', Lext=',es10.3,', Tc=',es10.3,/, 2 'ok? entrer o/n, pour arrêter entrer s') END SELECT READ*,oui c PRINT*,'oui ',oui IF(oui == 'o')THEN c_iben=c_iben*1.1d0 SELECT CASE(langue) CASE('english') WRITE(*,1802) 1802 FORMAT('calculation of a new model with 1.1 C') CASE DEFAULT WRITE(*,802) 802 FORMAT('calcul d''un nouveau modèle avec 1.1 C') END SELECT c seconde résolution du modèle quasi statique de PMS------------- CALL resout(-3,dt,dts) EXIT B2 ELSEIF(oui == 's')THEN CALL sortie('cesam 44') ELSE SELECT CASE(langue) CASE('english') WRITE(*,1056)c_iben 1056 FORMAT('c_iben requires a new value, previous value : ',es10.3) CASE DEFAULT WRITE(*,56)c_iben 56 FORMAT('il faut nouvelle valeur pour c_iben,' 1 ' valeur précédente : ',es10.3) END SELECT CYCLE B3 ENDIF ENDDO B2 c le second modèle est-il satisfaisant ? rp2=bp(3,m_qs*(n_qs-1)+1) ; lp2=bp(4,m_qs*(n_qs-1)+1) IF(en_m23)THEN rp2=SQRT(rp2) ; lp2=SQRT(lp2**3) ENDIF tp2=EXP(bp(2,1)) c estimation du dt correspondant dt=cte2*ABS(rp1-rp2)/rp1/rp2/(lp1+lp2) SELECT CASE(langue) CASE('english') WRITE(*,1057)c_iben,rp2,lp2,tp2,dt 1057 FORMAT('Pre main sequence second model, : C=',es10.3,/, 1 'Rext=',es10.3,', Lext=',es10.3,', Tc=',es10.3, 2 ', dt=',es10.3,/,'ok? enter o for yes, n for no,', 3 ' to stop enter s') CASE DEFAULT WRITE(*,57)c_iben,rp2,lp2,tp2,dt 57 FORMAT('Pré Séquence Principale second modèle : C=',es10.3,/, 1 'Rext=',es10.3,', Lext=',es10.3,', Tc=',es10.3, 2 ', dt=',es10.3,/,' ok? entrer o/n, pour arrêter entrer s') END SELECT READ*,oui IF(oui == 'n')THEN !nouvelle cte de contraction SELECT CASE(langue) CASE('english') WRITE(*,1803) 1803 FORMAT('a new value is required for c_iben') CASE DEFAULT WRITE(*,803) 803 FORMAT('il faut une nouvelle valeur pour c_iben') END SELECT dt=0.d0 ; CYCLE B3 ELSEIF(oui == 's')THEN CALL sortie('cesam 45') ELSE !on est content c on remplace la routine de réactions thermonucléaires fictive iben c par celle du fichier de données mon_modèle.don nom_nuc=nom_save c on garde le premier modèle de PMS d'âge 0 dans fichier xxxx_B.pms model_num=0 c écriture en binaire du fichier mon_modele_B.pms modelb=TRIM(nom_fich2)//'_B.pms' c WRITE(*,1)ADJUSTL(modelb) ; WRITE(2,1)ADJUSTL(modelb) OPEN(unit=26,form='unformatted',status='unknown',file=modelb) WRITE(26)ne,m_qs,n_qs,knot,nchim,n_ch,m_ch,knotc,Krot,nrot,n_rot, 1 m_rot,knotr,n_tds,knot_tds,mtot,alpha,w_rot,lim_ro,diffusion, 2 rot_solid,precisione,en_m23,f_eos,f_opa, 3 nom_ctes,nom_pertm,nom_pertw,nom_tdetau,nom_atm,nom_conv,nom_nuc, 4 nom_nuc_cpl,nom_diffm,nom_difft,nom_diffw,nom_etat,nom_opa, 5 nom_elem,bp,q,qt,chim,mc,mct,rota,mrot,mrott,tds,x_tds,xt_tds, 6 m_stat,r_stat,m_zc,r_zc,r_ov,age,dt,dts,mstar,rstar,mw_tot,wrot, 7 jlim,lconv,lim,model_num CLOSE(unit=26) WRITE(*,46)TRIM(modelb) model_num=-1 !est augmenté d'une unité ligne ~4155 c dtp=0.d0 : on poursuit le modèle PMS obtenu dtp=0.d0 ; dts=dt EXIT B3 ENDIF ENDDO B3 c pour l'entête du modèle du fichier mon_modèle.lis list_pms=.TRUE. c initialisation des paramètres pour le vent et les chutes de planétoïdes CALL vent ; CALL planetoides CASE DEFAULT SELECT CASE(langue) CASE('english') WRITE(*,1804)un23 1804 FORMAT('ARRET un23=',i3,', is different from -3, -2, -1, 1, 2, 3') CASE DEFAULT WRITE(*,804)un23 804 FORMAT('ARRET, un23=',i3,', est différent de -3, -2, -1, 1, 2, 3') END SELECT CALL sortie('cesam 47') END SELECT !initialisation des modèles c PAUSE'après initialisations' c----------évolution temporelle-------------- c les modèles de départ, soit de PMS, soit de ZAMS, soit de reprise c sont établis ou repris, on démarre l'évolution c ouverture du fichier du diagramme HR c pour les modèles initiaux, PMS ou ZAMS ABS(un23) = 2 ou 3 on ignore c ce qui existe dans le fichier, pour une reprise ABS(un23) = 1 on c écrit à la suite de ce qui existe en omettant celles du modèle repris SELECT CASE(ABS(un23)) CASE(1) !on écrira à la suite ! fichier *.TR à compléter IF(l_tr)THEN OPEN(unit=56,form='formatted',status='unknown',access='append', 1 file=TRIM(nom_fich2)//'_TR.csv') ELSEIF(l_lb)THEN OPEN(unit=56,form='formatted',status='unknown',access='append', 1 file=TRIM(nom_fich2)//'_LB.csv') ENDIF ! fichier *.HR ou LR (HR réduit) à complèter IF(l_hr)THEN OPEN(unit=53,form='formatted',status='unknown',access='append', 1 file=TRIM(nom_fich2)//'.LR') ELSE OPEN(unit=53,form='formatted',status='unknown',access='append', 1 file=TRIM(nom_fich2)//'.HR') ENDIF c PRINT*,l_hr ; PAUSE'cesam, hr1' ! fichier *.HR (et *.TR) ouverts ecrHR=.FALSE. ! fichier des limites des shells H et He burning IF(ecrHHe)OPEN(unit=54,form='formatted',status='unknown', 1 access='append',file=TRIM(nom_fich2)//'.mHHe') ! fichier des pulsations pour les Céphéïdes IF(cephe)OPEN(unit=55,form='formatted',status='unknown', 1 access='append',file=TRIM(nom_fich2)//'.cephe') CASE(2, 3) ! fichiers *.TR.csv, *.LB.csv à créer IF(l_tr)THEN OPEN(unit=56,form='formatted',status='unknown', 1 file=TRIM(nom_fich2)//'_TR.csv') c en-tête 12345678901234567890123456789012345678901234 len=44 list_TR(1)='model_num,age_Ma,Mstar_sun,Rstar_sun,Teff_K,' list_TR(2)='lum_sun,Pc_cgs,Tc_K,Roc_cgs,C_tr,Mmt_In,lim,' list_TR(3)='lconv1,lconv2,lconv3,lconv4,lconv5,' list_TR(4)='r_zc1,r_zc2,r_zc3,r_zc4,r_zc5,' list_TR(5)='m_zc1,m_zc2,m_zc3,m_zc4,m_zc5,' list_TR(6)='ro_zc1,ro_zc2,ro_zc3,ro_zc4,ro_zc5,' list_TR(7)='dv_zc1,dv_zc2,dv_zc3,dv_zc4,dv_zc5,' list_TR(8)='Mmt_zc1,Mmt_zc2,Mmt_zc3,Mmt_zc4,Mmt_zc5,' list_TR(9)='Mmt_atm,t_zc1,t_zc2,t_zc3,t_zc4,t_zc5,S_ctr,' list_TR(10)='S_ent1,S_ent2,S_ent3,S_ent4,S_ent5,S_atm' list_TRF=TRIM(list_TR(1))//TRIM(list_TR(2))//TRIM(list_TR(3)) 1 //TRIM(list_TR(4))//TRIM(list_TR(5))//TRIM(list_TR(6)) 2 //TRIM(list_TR(7))//TRIM(list_TR(8))//TRIM(list_TR(9)) 3 //TRIM(list_TR(10)) c écriture de l'entête WRITE(56,2009)TRIM(list_TRF) 2009 FORMAT(a) c WRITE(*,2009)TRIM(list_TRF) ; PAUSE'cesam TR' ! fichier *.LB.csv' à créer ELSEIF(l_lb)THEN OPEN(unit=56,form='formatted',status='unknown', 1 file=TRIM(nom_fich2)//'_LB.csv') c en-tête len=68 12345678901234567890123456789012345678901234567890123456789012345678 list_LB='model_num,age_Ma,R_sun,Teff_K,lum_sun,A,dnu02,dnu13,Nu0,P0,X_c,Y,Z/X' c écriture de l'entête WRITE(56,2009)TRIM(list_LB) ENDIF ! création du fichier *.HR ou LR (HR réduit) IF(l_hr)THEN OPEN(unit=53,form='formatted',status='unknown', 1 file=TRIM(nom_fich2)//'.LR') ELSE OPEN(unit=53,form='formatted',status='unknown', 1 file=TRIM(nom_fich2)//'.HR') ENDIF c PRINT*,l_hr ; PAUSE'cesam, hr2' ! le fichier *.HR (*.LR) est créé ecrHR=.TRUE. ! fichier des limites des shells H et He burning IF(ecrHHe)OPEN(unit=54,form='formatted',status='unknown', 1 file=TRIM(nom_fich2)//'.mHHe') ! fichier des pulsations pour les Céphéïdes IF(cephe)OPEN(unit=55,form='formatted',status='unknown', 1 file=TRIM(nom_fich2)//'.cephe') END SELECT c on ignore désormais les entrées binaire/ASCII, ZAMS, PMS, ou reprise un23=1 c allocations indépendantes du nombre de couches du modèle quasi-statique ALLOCATE(alfa_atm(n_atm),beta_atm(n_atm),cp_atm(n_atm), 1 delta_atm(n_atm),depsx(nchim),dmu_x(nchim),dxchim(nchim), 2 dxchimg(nchim),fd(ne,0:ord_qs-1),hp_atm(n_atm), 3 gamma_atm(n_atm),grad_atm(n_atm),grada_atm(n_atm), 4 gradc_atm(n_atm),gradr_atm(n_atm),grad_mj_a(n_atm), 5 ioni(0:NINT(MAXVAL(zi)),nchim),k_atm(n_atm),ldcapdt_a(n_atm), 6 ldcapdr_a(n_atm),lnpt_a(n_atm),lnpt_at(n_atm+m_ch), 7 lnta(1,n_atm),lnt_a(n_atm),mu_atm(n_atm),mue_atm(n_atm), 8 ro_atm(n_atm),vais_atm(n_atm),xchim(nchim),xchimg(nchim), 9 xchim1(nchim),xchim1g(nchim),zbar(nchim)) c-------évolution temporelle : boucle infinie Bev------------ Bev: DO c allocations générales pour les écritures, il faut changer c les dimensions lorsque le nombre de couches a varié IF(ALLOCATED(p))THEN IF(SIZE(p) /= n_qs)THEN DEALLOCATE(alfa,anube7,anub8,anun13, 1 anuo15,anupep,anupp,beta,compchim,compg,convec,cp,dcapdr,dcapdt, 2 delta,depsdr,depsdt,dcompchim,dcompg,dlpp,drdt,ener,epsilon, 3 gam,gamma,grad,gradad,gradconv,gradrad,grad_mj,grad_mu,hp,kap, 4 l,degene,m,mu,mue,p,pt,r,ro,t,tx_ioni,u,vaissala,w,z) ALLOCATE(alfa(n_qs),anube7(n_qs),anub8(n_qs),anun13(n_qs), 1 anuo15(n_qs),anupep(n_qs),anupp(n_qs),beta(n_qs), 2 compchim(nchim,n_qs),compg(nchim,n_qs),convec(n_qs),cp(n_qs), 3 dcapdr(n_qs),dcapdt(n_qs),delta(n_qs), 4 depsdr(n_qs),depsdt(n_qs),dcompchim(nchim,n_qs),dcompg(nchim,n_qs), 5 dlpp(n_qs),drdt(n_qs),ener(4,n_qs),epsilon(5,n_qs), 6 gam(n_qs),gamma(n_qs),grad(n_qs),gradad(n_qs),gradconv(n_qs), 7 gradrad(n_qs),grad_mj(n_qs),grad_mu(n_qs),hp(n_qs),kap(n_qs), 8 l(n_qs),degene(n_qs),m(n_qs),mu(n_qs),mue(n_qs),p(n_qs),pt(n_qs), 9 r(n_qs),ro(n_qs),t(n_qs),tx_ioni(nchim,n_qs), 1 u(n_qs),vaissala(n_qs),w(n_qs),z(n_qs)) ENDIF ELSE ALLOCATE(alfa(n_qs),anube7(n_qs),anub8(n_qs),anun13(n_qs), 1 anuo15(n_qs),anupep(n_qs),anupp(n_qs),beta(n_qs), 2 compchim(nchim,n_qs),compg(nchim,n_qs),convec(n_qs),cp(n_qs), 3 dcapdr(n_qs),dcapdt(n_qs),delta(n_qs), 4 depsdr(n_qs),depsdt(n_qs),dcompchim(nchim,n_qs),dcompg(nchim,n_qs), 5 dlpp(n_qs),drdt(n_qs),ener(4,n_qs),epsilon(5,n_qs), 6 gam(n_qs),gamma(n_qs),grad(n_qs),gradad(n_qs),gradconv(n_qs), 7 gradrad(n_qs),grad_mj(n_qs),grad_mu(n_qs),hp(n_qs),kap(n_qs), 8 l(n_qs),degene(n_qs),m(n_qs),mu(n_qs),mue(n_qs),p(n_qs),pt(n_qs), 9 r(n_qs),ro(n_qs),t(n_qs),tx_ioni(nchim,n_qs), 1 u(n_qs),vaissala(n_qs),w(n_qs),z(n_qs)) ENDIF c pour le test \tilde g / g avec diffusion du moment cinétique IF(COUNT(Krot == Ktest) >= 1)THEN IF(ALLOCATED(gtsg))DEALLOCATE(gtsg) ALLOCATE(gtsg(n_qs)) ; gtsg=0.d0 ENDIF c extraction des valeurs physiques des variables c calcul des énergies PP, CNO, 3alpha, grav pour intégration des énergies IF(ALLOCATED(ex))DEALLOCATE(ex) ; ALLOCATE(mt(n_qs+2),ex(nchim)) msh=0.d0 ; mshe=0.d0 !limites des shells H et He B5: DO i=1,n_qs !extraction p,t,r,l,m CALL bsp1ddn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,fd) IF(no_croiss)PRINT*,'Pb. at 16 in cesam' pt(i)=EXP(fd(1,0)) !Ptot IF(pturb)THEN !avec pression turbulente 7 inconnues p(i)=EXP(fd(Ipg,0)) !Pgas ELSE !sans pression turbulente 6 inconnues p(i)=pt(i) ENDIF IF(pturb .AND. der)THEN dlpp(i)=fd(Ipg,1)/fd(1,1) !dln Pgaz/dln Ptot ELSE dlpp(i)=1.d0 ENDIF t(i)=EXP(fd(2,0)) !variable ln T IF(en_m23)THEN r(i)=SQRT(ABS(fd(3,0))) !rayon/Rsol l(i)=SQRT(ABS(fd(4,0)))**3 !l/Lsol m(i)=SQRT(ABS(fd(5,0)))**3 !m/Msol ELSE r(i)=ABS(fd(3,0)) !rayon/Rsol l(i)=fd(4,0) !l/Lsol m(i)=ABS(fd(5,0))**3 !m/Msol fd(5,0)=m(i)**(2.d0/3.d0) ENDIF IF(r(i) > 0.d0)THEN grav=grav_sol*m(i)/r(i)**2 ELSE grav=0.d0 ENDIF c estimation de la vitesse sauf pour une reprise ou une initialisation c contraction locale quand dR/dt < 0 sinon expension IF(ALLOCATED(bp_t) .AND. dt > 0.d0)THEN CALL inter('m',bp_t,q_t,qt_t,n_qs_t,knot_t, 1 inside(m_stat_t(1),m_stat_t(n_qs_t),fd(5,0)),fqs,dfqs, 2 r_stat_t,m_stat_t) IF(en_m23)fqs(3)=SQRT(ABS(fqs(3))) !rayon/Rsol drdt(i)=(r(i)-fqs(3))/dt ELSE drdt(i)=0.d0 ENDIF c le gradient réel pour MJo grad_mj(i)=fd(2,1)/fd(1,1) c WRITE(*,2000)p(i),t(i),r(i),l(i),m(i),grad_mj(i) c WRITE(*,2000)r(i),r_stat(i),m(i),m_stat(i) c la composition chimique CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(fd(5,0),mc(n_ch))),lq,xchim,dxchim) IF(no_croiss)PRINT*,'Pb. at 17 in cesam' IF(fd(5,0) > 0.d0)THEN dxchimg=dxchim*2.d0/3.d0/SQRT(fd(5,0)) ELSE dxchimg=0.d0 ENDIF xchimg=xchim ; CALL chim_gram(xchimg,dxchimg) compchim(:,i)=xchim(:) ; compg(:,i)=xchimg(:) dcompchim(:,i)=dxchim(:) ; dcompg(:,i)=dxchimg(:) c limites des shells H et He IF(xchimg(1) < 1.d-3)msh=m(i) IF(ihe4 > 0)THEN IF(xchimg(ihe4) < 1.d-3)mshe=m(i) ENDIF c cas particuliers IF(i == 1)THEN lx_stop=lx_stop .OR. xchimg(1) <= x_stop ELSEIF(i == n_qs)THEN !pour l'atmosphère xchim1=xchim ; xchim1g=xchimg ENDIF c éléments lourds IF(ihe4 > 1)THEN c z(i)=1.d0-xchimg(1)-xchimg(ihe4)-xchimg(ihe4-1) z(i)=SUM(xchimg(ihe4+1:nchim)) ELSE z(i)=z0 ENDIF c équation d'état CALL etat(p(i),t(i),xchimg,.FALSE., 1 ro(i),drop,drot,drox,u(i),dup,dut,dux, 2 delta(i),deltap,deltat,deltax,cp(i),dcpp,dcpt,dcpx, 3 gradad(i),dgradadp,dgradadt,dgradadx,alfa(i),beta(i),gamma(i)) c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w(i)=wrot CASE(3,4) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot, 1 knotr,.TRUE.,MAX(mrot(1),MIN(fd(5,0),mrot(n_rot))),lq,frot,dfrot) IF(no_croiss)PRINT*,'Pb. at 18 in cesam' w(i)=frot(1) c MAX \tilde g / g IF(i == 1 .AND. w(1) > 0.d0)THEN gtsg(1)=cte6*w(1)**2/ro(1) ELSE gtsg(i)=cte4*r(i)**3*w(i)**2/m(i) ENDIF CASE(5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot, 1 knotr,.TRUE.,MAX(mrot(1),MIN(fd(5,0),mrot(n_rot))),lq,frot,dfrot) IF(no_croiss)PRINT*,'Pb. at 181 in cesam' w(i)=frot(1) CASE(6) CALL omega_jpz(r(i),w(i)) c IF(i < 6)PRINT*,'cesam',i,r(i),w(i),rstar END SELECT c l'énergie graviphique tds=TdS/dt variable globale du module c mod_static, a été calculée dans resout à l'issue c du calcul de ZAMS ou PMS ou est celle du modèle à poursuivre CALL nuc(t(i),ro(i),xchim,dxchim,jac,.FALSE.,3, 1 epsilo,et,ero,ex,hh,be7,b8,n13,o15,f17) CALL bsp1dn(1,tds,x_tds,xt_tds,n_tds,m_tds,knot_tds,.TRUE., 1 MIN(x_tds(i),x_tds(n_tds)),lq,ftds,dftds) IF(no_croiss)PRINT*,'Pb. at 19 in cesam' c après l'appel à nuc on a epsilo(PP+CNO+(3a+C+O), PP, CNO, 3a+C+O) c correction d'Yveline pour le signe du TdS, ABS(TdS) --> -TdS epsilo(5)=ftds(1)/secon6 epsilon(:,i)=epsilo !PP+CNO+(3a+C+O), PP, CNO, 3a+C+O, grav epsilon(1,i)=epsilon(1,i)-epsilo(5) !PP+CNO+(3a+C+O)-(+TdS) c les énergies locales ener(1,i)=epsilo(2) !PP ener(2,i)=epsilo(3) !CNO ener(3,i)=epsilo(4) !3alpha + C + O ener(4,i)=-epsilo(5) !grav ENDDO B5 c PRINT*,'lum' ; WRITE(*,2000)l(1:64) ; PAUSE'cesam, lum' c dr/dt en Rsol/My --> km/s drdt=drdt*cte7 c tabulation et intégration des énergies CALL bsp1dn(4,ener,m,mt,n_qs,2,knote,.FALSE.,m(1),lq,fen,dfen) IF(no_croiss)PRINT*,'Pb. at 20 in cesam' CALL sum_n(4,ener,mt,2,knote,.FALSE.,m(1),m(n_qs),fen) DEALLOCATE(mt) c les énergies globales en Lsol et proportions en %/ énergie nucléaire totale fen=fen*msol/lsol ; e_nuc=SUM(fen(1:3)) e_pp=fen(1) ; e_cno=fen(2) ; e_3al=fen(3) ; e_gr=fen(4) c PRINT*,'fen,e_nuc,e_gr' ; WRITE(*,2000)fen,e_nuc,e_gr c pourcentages respectifs IF(e_nuc > 0.d0)THEN i_pp=NINT(100.d0*(e_pp/e_nuc)) ; i_cno=NINT(100.d0*(e_cno/e_nuc)) i_3a=NINT(100.d0*(e_3al/e_nuc)) c l'énergie graviphique totale peut être négative i_gr=NINT(MIN(100.d0,100.d0*ABS(e_gr/e_nuc))) ELSE i_pp=0 ; i_cno=0 ; i_3a=0 ; i_gr=100 ENDIF c PRINT*,'énergies',i_pp,i_cno,i_3a,i_gr c WRITE(*,2000)e_nuc,e_pp,e_cno,e_3al,e_gr,dt ; PAUSE'après i_gr' c estimation de la nature du modèle. Source des T caractéristiques wikipedia zams_pr=zams ; post_pr=post ; cohe_pr=cohe ; coca_pr=coca cone_pr=cone ; coox_pr=coox ; cosi_pr=cosi t_max=MAXVAL(t) IF(compg(1,1) > x_tams)THEN IF(e_nuc >= ABS(e_gr))THEN chaine='modèle de série principale ' zams=.TRUE. ; post=.FALSE. ; cohe=.FALSE. ; coca=.FALSE. cone=.FALSE. ; coox=.FALSE. ; cosi=.FALSE. ELSE SELECT CASE(langue) CASE('english') chaine='pre-main sequence model' CASE DEFAULT chaine='modèle de pré-série principale' END SELECT zams=.FALSE. ; post=.FALSE. ; cohe=.FALSE. ; coca=.FALSE. cone=.FALSE. ; coox=.FALSE. ; cosi=.FALSE. ENDIF ELSEIF(t_max < 1.d8)THEN SELECT CASE(langue) CASE('english') chaine='post-main sequence model' CASE DEFAULT chaine='modèle post série principale' END SELECT zams=.FALSE. ; post=.TRUE. ; cohe=.FALSE. ; coca=.FALSE. cone=.FALSE. ; coox=.FALSE. ; cosi=.FALSE. ELSEIF(t_max < 6.d8)THEN SELECT CASE(langue) CASE('english') chaine= 'model with helium burning' CASE DEFAULT chaine='modèle avec combustion hélium' END SELECT zams=.FALSE. ; post=.FALSE. ; cohe=.TRUE. ; coca=.FALSE. cone=.FALSE. ; coox=.FALSE. ; cosi=.FALSE. ELSEIF(t_max < 9.d8)THEN SELECT CASE(langue) CASE('english') chaine='model with carbon burning' CASE DEFAULT chaine= 'modèle avec combustion carbone' END SELECT zams=.FALSE. ; post=.FALSE. ; cohe=.FALSE. ; coca=.TRUE. cone=.FALSE. ; coox=.FALSE. ; cosi=.FALSE. ELSEIF(t_max < 1.d9)THEN SELECT CASE(langue) CASE('english') chaine='model with neon burning' CASE DEFAULT chaine= 'modèle avec combustion néon' END SELECT zams=.FALSE. ; post=.FALSE. ; cohe=.FALSE. ; coca=.FALSE. cone=.TRUE. ; coox=.FALSE. ; cosi=.FALSE. ELSEIF(t_max < 1.2d9)THEN SELECT CASE(langue) CASE('english') chaine='model with oxygen burning' CASE DEFAULT chaine= 'modèle avec combustion oxygène' END SELECT zams=.FALSE. ; post=.FALSE. ; cohe=.FALSE. ; coca=.FALSE. cone=.FALSE. ; coox=.TRUE. ; cosi=.FALSE. ELSE SELECT CASE(langue) CASE('english') chaine='model with silicium burning' CASE DEFAULT chaine= 'modèle avec combustion silicium' END SELECT zams=.FALSE. ; post=.FALSE. ; cohe=.FALSE. ; coca=.FALSE. cone=.FALSE. ; coox=.FALSE. ; cosi=.TRUE. ENDIF c l'énergie graviphique totale peut être négative i_gr=i_gr*(SIGN(1.d0,e_gr)) c PRINT*,i_pp,i_cno,i_3a,i_gr,chaine ; PAUSE'chaine' c avec diffusion du moment cinétique, estimation du maximum de \tilde g / g IF(COUNT(Krot == Ktest) >= 1)THEN gtmax=MAXVAL(gtsg) ; igtmax=MAXLOC(gtsg) rgtmax=r(igtmax(1))/rstar ; mgtmax=m(igtmax(1))/mstar ENDIF c actualisation du nombre de modèles c nb_modeles : nombre de modèles depuis le début du calcul ou la reprise c model_num : nombre de modèles depuis l'initialisation ZAMS ou PMS c nb_der_out : numéro de la dernière formation du fichier output nb_modeles=nb_modeles+1 ; model_num=model_num+1 c pour les "long run" vers les stades avancés c après la séquence principale on impose l'approximation de Kippenhahn c ce qui évite le calcul de dV/dt à travers une discontinuité de c composition chimique en particulier à la limite d'un coeur convectif c sauf pour le premier passage (alors passe=.FALSE.) IF(passe)THEN zams_e= .NOT.zams_pr .AND. zams !on passe PMS--->ZAMS post_e= .NOT.post_pr .AND. post !on passe ZAMS-->POST cohe_e= .NOT.cohe_pr .AND. cohe !on passe POST-->COHE coca_e= .NOT.coca_pr .AND. coca !on passe COHE-->COCA cone_e= .NOT.cone_pr .AND. cone !on passe COCA-->CONE coox_e= .NOT.coox_pr .AND. coox !on passe CONE-->COOX cosi_e= .NOT.cosi_pr .AND. cosi !on passe COOX-->COSI ENDIF passe=.TRUE. c suivant le cas: écriture des fichiers mon_modèle.zams, mon_modèle.post, c mon_modèle.cohe, mon_modèle.coca, mon_modèle.cone, mon_modèle.cone, c mon_modèle.coox ou mon_modèle.cosi IF(zams_e .OR. post_e .OR. cohe_e .OR. coca_e .OR. cone_e 1 .OR. coox_e .OR. cosi_e)THEN IF(zams_e)THEN modelb=TRIM(nom_fich2)//'_B.zams' ELSEIF(post_e)THEN modelb=TRIM(nom_fich2)//'_B.post' ELSEIF(cohe_e)THEN modelb=TRIM(nom_fich2)//'_B.cohe' ELSEIF(coca_e)THEN modelb=TRIM(nom_fich2)//'_B.coca' ELSEIF(cone_e)THEN modelb=TRIM(nom_fich2)//'_B.cone' ELSEIF(coox_e)THEN modelb=TRIM(nom_fich2)//'_B.coox' ELSEIF(cosi_e)THEN modelb=TRIM(nom_fich2)//'_B.cosi' ENDIF OPEN(unit=26,form='unformatted',status='unknown',file=TRIM(modelb)) WRITE(26)ne,m_qs,n_qs,knot,nchim,n_ch,m_ch,knotc,Krot,nrot,n_rot, 1 m_rot,knotr,n_tds,knot_tds,mtot,alpha,w_rot,lim_ro,diffusion, 2 rot_solid,precisione,en_m23,f_eos,f_opa, 3 nom_ctes,nom_pertm,nom_pertw,nom_tdetau,nom_atm,nom_conv,nom_nuc, 4 nom_nuc_cpl,nom_diffm,nom_difft,nom_diffw,nom_etat,nom_opa, 5 nom_elem,bp,q,qt,chim,mc,mct,rota,mrot,mrott,tds,x_tds,xt_tds, 6 m_stat,r_stat,m_zc,r_zc,r_ov,age,dt,dts,mstar,rstar,mw_tot,wrot, 7 jlim,lconv,lim,model_num CLOSE(unit=26) WRITE(*,46)TRIM(modelb) ENDIF c llist est initialisé .FALSE. c mettre llist=.TRUE. pour écrire les modèles sur l'écran c et dans le fichier sortie.dat c par exemple, pour récuperer un modèle en ASCII c servant à l'initialisation PMS ou ZAMS c la différence entre Ptot et Pgaz est ici ignorée c llist=.TRUE. IF(llist)THEN OPEN(unit=4,form='formatted',status='new',file='sortie.dat') WRITE(*,101)age 101 FORMAT('age=',es10.3) WRITE(*,102) 102 FORMAT(t9,'1-m',t24,'p',t36,'t',t49,'r',t60,'l',t71,'X') DO k=n_qs,1,-1 WRITE(*,2001)1.d0-m(k)/mstar,p(k),t(k),r(k),l(k),compg(1,k) WRITE(4,2001)1.d0-m(k)/mstar,p(k),t(k),r(k),l(k),compg(1,k) 2001 FORMAT(es17.10,4es12.5,es10.3) ENDDO c pour la fermeture du dessin TRUE pour appel à pgend CLOSE(unit=4) ; CALL des(.TRUE.,dt,teff) PRINT*,'arrêt après écriture du modèle en ASCII' ; CALL sortie('cesam 48') ENDIF c estimation de la température effective et dessin en ligne rext=r(n_qs) ; lext=l(n_qs) ; teff=cte8*SQRT(SQRT(lext)/rext) IF(age > 0.d0)THEN c WRITE(*,2000)teff,r(n_qs),l(n_qs) ; PAUSE'Teff2' CALL des(.FALSE.,dtp,teff) ELSE CALL des(.FALSE.,0.d0,teff) ENDIF c WRITE(*,2000)teff; PAUSE'des1' ; CALL des(.FALSE.,dt,teff) c PAUSE'des2' c--------------------------------------------------------------------- c E C R I T U R E S c--------------------------------------------------------------------- lteff=LOG10(teff) !teff approche pour le test d'arrêt log_l=LOG10(l(n_qs)) c le modèle est-il dans la zone d'instabilité? CALL pulsation(l(n_qs),teff,typ,mstar,puls,period) IF(puls)THEN instb='Inside the instability strip' ELSE instb='Outside the instability strip' ENDIF c on sort? c sort=(age >= agemax-1.d-5) sort=(age >= agemax) 1 .OR. (nb_modeles >= nb_max_modeles) 2 .OR. lx_stop.OR. lt_stop .OR. lhe_stop .OR. lr_stop 3 .OR. (zams .AND. arret == 'zams') 4 .OR. (post .AND. arret == 'post') 5 .OR. (cohe .AND. arret == 'cohe') 6 .OR. (coca .AND. arret == 'coca') 7 .OR. (cone .AND. arret == 'cone') 8 .OR. (coox .AND. arret == 'coox') 9 .OR. (cosi .AND. arret == 'cosi') 1 .OR. list_sort c formation du fichier de rotation IF(nb_modeles > 0)THEN SELECT CASE(Kdes_rot) CASE(1) IF(sort)CALL ecrit_rota(dt) CASE(2,4) CALL ecrit_rota(dt) END SELECT ENDIF c au cours de l'évolution (lteffp > 0) on sort si: c en se déplacant vers les Teff croissantes Log Teff > LOG_teff c en se déplacant vers les Teff décroissantes Log Teff < -LOG_teff c LOG_teff étant positif ou négatif selon le cas IF(lteffp > 0.d0 .AND. age > age_car)THEN IF(LOG_teff > 0.d0 .AND. lteff > lteffp)THEN sort=sort .OR. lteff > LOG_teff !à gauche dans HR ELSEIF(LOG_teff < 0.d0 .AND. lteff < lteffp)THEN sort=sort .OR. lteff < -LOG_teff !à droite dans HR ENDIF ENDIF lteffp=lteff ; n_ecrit=n_ecrit+1 c si l_hr (listing réduit) le nombre de listing complet n'est pas incrémenté IF(.NOT.l_hr)n_ecrit=n_ecrit+1 ecrit=sort .OR. zams_e .OR. post_e .OR. cohe_e .OR. coca_e 1 .OR. cone_e .OR. coox_e .OR. cosi_e .OR. age-agep >= dtlist 2 .OR. list_rep .OR. list_zams .OR. list_pms 3 .OR. n_ecrit >= ntot_ecrit IF(ecrit)n_ecrit=0 c si l_hr (listing réduit) on forme un listing réduit c ecrit=ecrit .AND. .NOT.l_hr c décompte des écritures pour les fichiesr intermédiaires c ecr_inter=.TRUE. les fichiers intermédiaires ASCII mon_modelexxxx.osc c et/ou binaire xxxx-mon_modele_B.rep seront écrits c nb_der_out: numéro du dernier modèle dont les fichiers intermédiaires c seront écrits si, respectivement, all_osc et/ou all_rep sont .TRUE. ecr_inter=nb_der_out+dn_sort == model_num .OR. model_num == 0 IF(ecr_inter)nb_der_out=model_num c PRINT*,nb_der_out,dn_sort,model_num,nom_output c PRINT*,nb_der_out+dn_sort,model_num,ecr_inter ; PAUSE'en C2' c si écrit: écriture, ensuite si sort: sortie c PRINT*,sort,ecrit,list_rep,list_sort ; PAUSE'list_sort' c-----------------début du calcul des éléments pour les sorties------------ c avant écritures, calcul de N^2, nu0, dnu0, grad,... SI + atm c interpolation de variables thermodynamiques CALL tab_vth c calcul de r_son et vson pour la détermination des fréquences c nombre total de lignes avec atmosphère n_tot=n_qs IF(n_atm > 0)n_tot=n_tot+n_atm-1 c allocations IF(ALLOCATED(r_son))DEALLOCATE(r_son,v_son) ALLOCATE(r_son(n_tot),v_son(n_tot)) c pour la stucture interne DO i=1,n_qs xchim=compchim(:,i) ; dxchim=dcompchim(:,i) xchimg=compg(:,i) ; dxchimg=dcompg(:,i) CALL thermo(pt(i),p(i),t(i),m(i),l(i),r(i),dlpp(i),xchim,dxchim, 1 ro(i),drop,drot,drox,u(i),dup,dut,dux, 2 grad(i),dgradpt,dgradp,dgradt,dgradx,dgradm,dgradl,dgradr, 3 dgradlpp,gam(i),dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr, 4 dgamlpp,epsilo,depsp,depst,depsx,kap(i),dkapp,dkapt,dkapx, 5 delta(i),deltap,deltat,deltax,cp(i),dcpp,dcpt,dcpx, 6 gradad(i),dgradadp,dgradadt,dgradadx, 7 hp(i),dhppt,dhpp,dhpt,dhpx,dhpr,dhpm,grad_mu(i),mu(i), 8 dmu_x,mue(i),gradrad(i),alfa(i),beta(i),gamma(i),convec(i)) c r_son & v_son dans la SI v_son(i)=SQRT(gamma(i)*p(i)/ro(i)) ; r_son(i)=r(i)*rsol c dérivées de l'opacité et de l'énergie nucléaire dcapdr(i)=dkapp/drop ; dcapdt(i)=dkapt-dcapdr(i)*drot depsdr(i)=depsp/drop ; depsdt(i)=depst-depsdr(i)*drot c thermo sort avec radiatif=convec=.TRUE. dans ZR convec(i)=.NOT.convec(i) !thermo sort avec convec=.TRUE. dans ZR IF(convec(i))THEN gradconv(i)=grad(i) ELSE gradconv(i)=0.d0 ENDIF c estimation des taux d'ionisation et calcul des poids c moléculaires moyens, si T < 4500K on suppose les éléments c totalement recombinés (saha ne converge pas) IF(t(i) > 4.5d3)THEN IF(mu_saha)THEN CALL saha(xchim,t(i),ro(i),ioni,zbar,nel,degene(i)) ELSE zbar=zi ; degene(i)=-100.d0 ENDIF ELSE zbar=0.d0 ; degene(i)=-100.d0 ENDIF c les taux d'ionisation c arb_rom : transformation nb. arabes --> nb. romains DO j=1,nchim tx_ioni(j,i)=arb_rom(NINT(zbar(j))+1) ENDDO c PRINT*,tx_ioni(:,i) ; PAUSE'tx_ioni' c PRINT*,i,pturb c WRITE(*,2000)pt(i),p(i),t(i),xchim(1),m(i),l(i),r(i),dlpp(i) c WRITE(*,2000)drop,dkapp,dcapdr(i) c WRITE(*,2000)dcapdt(i),dkapt,dcapdr(i),drot c WRITE(*,2000)depsdr(i),depsp,depsdt(i),depst,depsdr(i) c calcul d'une des formes de Väissälä (Kippenhahan & Weigert p. 42) c 1/gamma1 dlnP/dln r - dln ro/dln r = IF(hp(i) <= 0.d0)THEN !au centre vaissala(i)=0.d0 ELSEIF(new_bv)THEN c utilisation de la tabulation vth c fvth(1)=lnP, fvth(2)=lnT, fvth(5)=ln ro, fvth(6)=cp, fvth(7)=delta, c fvth(8)=gamma1, fvth(9)=ln µ, fvth(10)=ln kap, fvth(11)=degene c en m^2/3 : fvth(3)=r^2, fvth(4)=l^2/3 c en m^1/3 : fvth(3)=r, fvth(4)=l CALL bsp1dn(nvth,vth,mc,mct,n_ch,m_ch,knotc,.TRUE., 1 MAX(mc(1),MIN(m(i)**(2.d0/3.d0),mc(n_ch))),lq,fvth,dfvth) IF(no_croiss)PRINT*,'Pb. at 211 in cesam' vaissala(i)=(dfvth(1)/fvth(8)-dfvth(5))*fvth(3)/dfvth(3) IF(en_m23)vaissala(i)=2.d0*vaissala(i) c ancienne formulation ELSE c vaissala(i)=r(i)*rsol/hp(i)*(alfa(i)*(dlpp(i)-1.d0)+ c 1 delta(i)*(gradad(i)-grad_mj(i) c 2 +beta(i)/(4.d0-3.d0*beta(i))*grad_mu(i))) vaissala(i)=r(i)*rsol/hp(i)*(alfa(i)*(dlpp(i)-1.d0) 1 +delta(i)*(gradad(i)-grad_mj(i)))-cte3*r(i)**3*drox*dxchimg(1) ENDIF c WRITE(*,*)'m,r,drot,t,p,grad,gradad,drox/vai,dchimg,nuc' c WRITE(*,2000)m(i),r(i),drot,t(i),p(i),grad(i),gradad(i),drox c WRITE(*,2000)vaissala(i),dxchimg(1) c production de neutrinos CALL nuc(t(i),ro(i),xchim,dcomp,jac,.FALSE.,4, 1 epsilo,et,ero,ex,anupp(i),anupep(i),anub8(i),anube7(i), 2 anun13(i),anuo15(i)) ENDDO c calcul de P0 lb_p0=p0(r,m,vaissala) c WRITE(*,2000)lb_p0 ; PAUSE'appel à p0' c atmosphère, température effective, rayon extérieur c WRITE(*,*)'wrot,l(n_qs),r(n_qs),xchim en n' c WRITE(*,2000)wrot,l(n_qs),r(n_qs),xchim1g(1:nchim) c WRITE(*,*)rstar CALL atm(.TRUE.,l(n_qs),r(n_qs),xchim1g,ptext,dptdl,dptdr, 1 text,dtdl,dtdr,mext,dml,dmr,pext,dpl,dpr,teff) c WRITE(*,*)rstar ; PAUSE'lim ext 1' c l'atmosphère n'est détaillée que si elle est reconstituée ie. n_atm > 1 IF(n_atm > 1)THEN dxchim=0.d0 DO i=1,n_atm !thermo pour l'atmosphère grav=grav_sol*m_atm(i)/r_atm(i)**2 CALL tdetau(tau(i),teff,grav,bid,dtsdtau,dtsdteff,dtsdg, 1 ro_ext,dro_grav,dro_teff,f_tau,df_tau,d2f_tau) c WRITE(*,2000)tau,f_tau,df_tau,d2f_tau c les gradients CALL thermo_atm(pt_atm(i),p_atm(i),t_atm(i),xchim1g,m_atm(i), 1 l(n_qs),r_atm(i),dlpp_atm(i),tau(i),df_tau,d2f_tau,rstar, 2 ro_atm(i),drop,drot,k_atm(i),dkapp,dkapt, 3 grada_atm(i),dgradadp,dgradadt, 4 grad_atm(i),dgradpt,dgradp,dgradt,dgradr,dgradrs,dgradm, 5 dgradtau,dgradlpp,gam_atm,dgampt,dgamp,dgamt,dgamr,dgamrs,dgamm, 6 dgamtau,dgamlpp,hp_atm(i),dhppt,dhpp,dhpt,dhpr,dhpm, 7 delta_atm(i),deltap,deltat,cp_atm(i),gradr_atm(i), 8 alfa_atm(i),beta_atm(i),gamma_atm(i),rad_atm,.TRUE.) IF(rad_atm)THEN gradc_atm(i)=0.d0 ELSE gradc_atm(i)=grad_atm(i) ENDIF c r_son & v_son dans l'atmosphère c la couche 1 de l'atmosphère étant la couche n_qs de la SI IF(i > 1)THEN v_son(n_qs+i-1)=SQRT(gamma_atm(i)*p_atm(i)/ro_atm(i)) r_son(n_qs+i-1)=r_atm(i)*rsol ENDIF c pour calculer les gradients tabulation de ln T en fonction de ln P lnpt_a(i)=LOG(pt_atm(i)) ; lnt_a(i)=LOG(t_atm(i)) ldcapdr_a(i)=dkapp/drop ; ldcapdt_a(i)=dkapt-ldcapdr_a(i)*drot c PRINT*,'p_atm,t_atm,xchim1,r_atm,l_atm,m_atm' c WRITE(*,2000)p_atm(i),t_atm(i),xchim1(1),r_atm(i),l(n_qs), c 1 m_atm(i),lnpt_a(i),lnt_a(i) c PAUSE'cesam atm list' c dans l'atmosphère pas de gradient de comp. chim c pour Väissälä on prend ou bien le vrai gradient grad_mj c ou celui calculé aux points de grille, X=cte dans l'atmosphère vais_atm(i)=r_atm(i)*rsol/hp_atm(i)*(alfa_atm(i)* 1 (dlpp_atm(i)-1.d0)+delta_atm(i)*(grada_atm(i)-grad_atm(i))) c estimation des taux d'ionisation et calcul des poids c moléculaires moyens, si T < 4500K on suppose les éléments c totalement recombinés (saha ne converge pas) IF(t_atm(i) < 4.5d3)THEN zbar=0.d0 ; mue_atm(i)=1.d30 ELSE CALL saha(xchim1,t_atm(i),ro_atm(i),ioni,zbar,nel,bid) mue_atm(i)=1.d0/DOT_PRODUCT(zbar,xchim1) ENDIF mu_atm(i)=1.d0/DOT_PRODUCT((1.d0+zbar),xchim1) ENDDO c estimation du gradient dans l'atmosphère les abs. devant être croissantes c on prend l'opposé pour tabuler ln T en fct. de ln P IF(n_atm > m_ch)THEN lnpt_a=-lnpt_a ; lnta=RESHAPE(lnt_a,SHAPE(lnta)) !;lnta(1,:)=lnt_a CALL bsp1dn(1,lnta,lnpt_a,lnpt_at,n_atm,m_ch, 1 k,.FALSE.,lnpt_a(1),lq,fatm,dfatm) IF(no_croiss)PRINT*,'Pb. at 212 in cesam' DO i=1,n_atm CALL bsp1dn(1,lnta,lnpt_a,lnpt_at,n_atm,m_ch,k,.TRUE., 1 lnpt_a(i),lq,fatm,dfatm) IF(no_croiss)PRINT*,'Pb. at 21 in cesam' grad_mj_a(i)=fatm(1) ENDDO lnpt_a=-lnpt_a ; grad_mj_a=-grad_mj_a ELSE grad_mj_a=0.4d0 ENDIF ENDIF !fin atm c estimation de Nu0 et dnu02 et dnu13 SI + atm CALL dnunl(r_son,v_son,n_tot,nu0,dnu02,dnu13,a) c calcul de la dérivée seconde + normalisation de la pression au centre d2p=cte5*ro(1)**2*rstar**2/p(1) c calcul de la dérivée seconde + normalisation de la densité au centre fonc(0)=ro(1) ; fonc(1)=ro(2) ; fonc(2)=ro(3) fonc(3)=ro(4) ; fonc(4)=ro(5) ; fonc(5)=ro(6) absc(0)=r(1) ; absc(1)=r(2) ; absc(2)=r(3) absc(3)=r(4) ; absc(4)=r(5) ; absc(5)=r(6) CALL newton(0,m_qs+1,fonc,absc,0.d0,poly,2) d2ro=rstar**2/ro(1)*poly(2) c WRITE(*,2000)fonc ; WRITE(*,2000)absc c PRINT*,'rstar,d2ro,ro,poly(2)' c WRITE(*,2000)rstar,d2ro,ro(1),poly(0:2) c PRINT*,'d2p,d2ro,rstar,p(1),ro(1)' c WRITE(*,2000)d2p,d2ro,rstar,p(1),ro(1) c PAUSE'ro/r/deriv' c------------fin du calcul des éléments pour sorties------------------------ c vérification de la solution: estimation des premiers et 2d membres c des équations d'évolution IF(.FALSE.)THEN c IF(.TRUE.)THEN DO i=2,n_qs-1 WRITE(*,2000)(t(i+1)-t(i-1))/(m(i+1)-m(i-1))/(mtot*msol), 1 -g*m(i)*mtot*msol/4.d0/pi/(rsol*r(i))**4*grad(i)*t(i)/p(i), 2 (pt(i+1)-pt(i-1))/(m(i+1)-m(i-1))/(mtot*msol), 3 -g*m(i)*mtot*msol/4.d0/pi/(rsol*r(i))**4, 4 (r(i+1)-r(i-1))*rsol/(m(i+1)-m(i-1))/(mtot*msol), 5 1.d0/4.d0/pi/(rsol*r(i))**2/ro(i), 6 (l(i+1)-l(i-1))*lsol/(m(i+1)-m(i-1))/(mtot*msol), 7 epsilon(1,i) ENDDO PAUSE'verif' ENDIF c --------------------les sorties début--------------------------- c pour sortir il faut former les fichiers ASCII, il en est de même quand c on ne sort que les fichiers ASCII IF(ecrit .OR. all_cephe .OR. ecr_inter)THEN c écritures de gestion CALL ecrit1(list_pms,list_rep,list_sort,list_zams,model_num) c sortie IF(sort)THEN CALL list(a,alfa,anub8,anube7,anun13,anuo15,anupep,anupp,beta, 1 compg,cp,delta,dcapdr,dcapdt,degene,depsdr,depsdt,dnu02,dnu13, 2 drdt,d2p,d2ro,chaine,convec,sort,epsilon,gam,gamma,gradad,grada_atm, 3 gradconv,gradc_atm,gradrad,gradr_atm,hp,i_cno,i_gr,i_pp,i_3a, 4 kap,l,lb_p0,m,mu,mue,m_atm,nu0,p,pt,pt_atm,p_atm,r,ro,ro_atm, 5 r_atm,t,tau,teff,tx_ioni,t_atm,u,vaissala,w,z) c output ASCII, on n'écrira que le *osc standard c PRINT*,sort ; PAUSE'sort' all_cephe=.FALSE. IF(nom_output /= 'no_output')CALL ecrit_ascii c sorties sur écran~~~~~~~~~~~~~~ c Alerte avec diffusion du moment cinétique si \tilde g \ g > 20% IF(COUNT(Krot == Ktest) >= 1)THEN IF(gtmax > 0.2d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1808)nom_thw,gtmax,igtmax,rgtmax,mgtmax WRITE(2,1808)nom_thw,gtmax,igtmax,rgtmax,mgtmax 1808 FORMAT(/,'WARNING The formalism of ',a,/, 1 'is no longer valuable : MAX acce.centr./gravity=',es10.3, 2 ', shell ',i4,', R/R*=',es10.3,', M/M*=',es10.3) CASE DEFAULT WRITE(*,808)nom_thw,gtmax,igtmax,rgtmax,mgtmax WRITE(2,808)nom_thw,gtmax,igtmax,rgtmax,mgtmax 808 FORMAT(/,'ATTENTION le formalisme de ',a,/, 1 'est injustifié : MAX accé.centr./gravité=',es10.3,', couche ',i4, 2 ', R/R*=',es10.3,', M/M*=',es10.3) END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1809)gtmax,igtmax,rgtmax,mgtmax WRITE(2,1809)gtmax,igtmax,rgtmax,mgtmax 1809 FORMAT(/,'MAX acce.centr./gravity=',es10.3, 1 ', shell ',i4,', R/R*=',es10.3,', M/M*=',es10.3) CASE DEFAULT WRITE(*,809)gtmax,igtmax,rgtmax,mgtmax WRITE(2,809)gtmax,igtmax,rgtmax,mgtmax 809 FORMAT(/,'MAX accé.centr./gravité=',es10.3,', couche ',i4, 1 ', R/R*=',es10.3,', M/M*=',es10.3) END SELECT ENDIF ENDIF c écriture du cartouche c PRINT*,mstar,mtot zsx=z(n_qs)/compg(1,n_qs) pertmm=mstar/mtot-1.d0 IF(ihe4 <= 1)THEN WRITE(*,205)age,lteff,log_l,LOG10(rstar), 1 LOG10(mstar/rstar**2*grav_sol),pt(1),t(1),ro(1),compg(1,1), 2 i_pp,i_cno,i_3a,i_gr,pertmm,mstar,TRIM(chaine),TRIM(instb), 3 period,zsx 205 FORMAT(/,'*********',/,'âge=',es10.3,', LogTeff=',es10.3, 1 ', LogL/Lsol=',es10.3,', LogR/Rsol=',es10.3,/, 2 'Log10 g=',es10.3,', Pc=',es10.3,', Tc=',es10.3, 3 ', Roc=',es10.3,', Xc=',es10.3,/, 4 'ePP/eNUC=',i3,'%, eCNO/eNUC=',i3,'%, e3a+C+O/eNUC=',i3, 5 '%, eGRAV/eNUC=',i5,'%',/,'Var. rel. de masse : ',es10.3, 6 ', M*=',es10.3,'Msol, ',a,/,a30,' (P=',es10.3,'day), Z/X=',es10.3,/, 7 '*********',/) ELSEIF(w(n_qs) > 0.d0)THEN WRITE(*,207)age,lteff,log_l,LOG10(rstar), 1 LOG10(mstar/rstar**2*grav_sol),pt(1),t(1),ro(1),compg(1,1), 2 i_pp,i_cno,i_3a,i_gr,compg(ihe4,1),pertmm,mstar,TRIM(chaine), 3 2.d0*pi/24.d0/3600.d0/w(n_qs),rstar*w(n_qs)*rsol*1.d-5, 4 TRIM(instb),period,zsx 207 FORMAT(/,'*********',/,'âge=',es10.3,', LogTeff=',es10.3, 1 ', LogL/Lsol=',es10.3,', LogR/Rsol=',es10.3,/, 2 'Log g=',es10.3,', Pc=',es10.3,', Tc=',es10.3, 3 ', Roc=',es10.3,', Xc=',es10.3,/, 4 'ePP/eNUC=',i3,'%, eCNO/eNUC=',i3,'%, e3a+C+O/eNUC=',i3, 5 '%, eGRAV/eNUC=',i5,'%',', Yc=',es10.3,/, 6 'Var. rel. de masse : ',es10.3, 7 ', M*=',es10.3,'Msol, ',a,/,'Période de rota.=',es10.3, 8 ' jours,',' vitesse de rotation à la surface=',es10.3, 9 ' km/s',/,a30,' (P=',es10.3,'day), Z/X=',es10.3,/,'*********',/) ELSE WRITE(*,206)age,lteff,log_l,LOG10(rstar), 1 LOG10(mstar/rstar**2*grav_sol),pt(1),t(1),ro(1),compg(1,1), 2 i_pp,i_cno,i_3a,i_gr,compg(ihe4,1),pertmm,mstar,TRIM(chaine), 3 TRIM(instb),period,zsx 206 FORMAT(/,'*********',/,'âge=',es10.3,', LogTeff=',es10.3, 1 ', LogL/Lsol=',es10.3,', LogR/Rsol=',es10.3,/, 2 'Log g=',es10.3,', Pc=',es10.3,', Tc=',es10.3, 3 ', Roc=',es10.3,', Xc=',es10.3,/, 4 'ePP/eNUC=',i3,'%, eCNO/eNUC=',i3,'%, e3a+C+O/eNUC=',i3, 5 '%, eGRAV/eNUC=',i5,'%',', Yc=',es10.3,/, 6 'Var. rel. de masse : ',es10.3,', M*=',es10.3,'Msol, ',a,/,a30, 7 ' (P=',es10.3,'day), Z/X=',es10.3,/,'*********',/) ENDIF c écriture du titre et de la raison de la sortie titre=TRIM(version) SELECT CASE(langue) CASE('english') WRITE(*,1501)titre ; WRITE(2,1501)titre 1501 FORMAT('End of evolution with CESAM2k version : ',a) CASE DEFAULT WRITE(*,501)titre ; WRITE(2,501)titre 501 FORMAT('Fin d''évolution avec CESAM2k version : ',a) END SELECT c IF(age >= agemax-1.d-5)THEN IF(age >= agemax)THEN SELECT CASE(langue) CASE('english') WRITE(*,1409) ; WRITE(2,409) 1409 FORMAT('Stop because agemax has been reached') CASE DEFAULT WRITE(*,409) ; WRITE(2,409) 409 FORMAT('Sortie car agemax atteint') END SELECT ELSEIF(nb_modeles >= nb_max_modeles)THEN SELECT CASE(langue) CASE('english') WRITE(*,1410) ; WRITE(2,410) 1410 FORMAT('Stop : the maximum number of models has been reached') CASE DEFAULT WRITE(*,410) ; WRITE(2,410) 410 FORMAT('Sortie : le nombre maximum de modèles est atteint') END SELECT ELSEIF(lx_stop)THEN SELECT CASE(langue) CASE('english') WRITE(*,1401) ; WRITE(2,1401) 1401 FORMAT('Stop because H < x_stop at the center') CASE DEFAULT WRITE(*,401) ; WRITE(2,401) 401 FORMAT('Sortie car H < x_stop au centre') END SELECT ELSEIF(lr_stop)THEN SELECT CASE(langue) CASE('english') WRITE(*,1413) ; WRITE(2,1413) 1413 FORMAT('Stop because R_stop is reached') CASE DEFAULT WRITE(*,413) ; WRITE(2,413) 413 FORMAT('Sortie car R_stop est atteint') END SELECT ELSEIF(t(1) >= t_stop)THEN SELECT CASE(langue) CASE('english') WRITE(*,1402) ; WRITE(2,1402) 1402 FORMAT('Stop because T > t_stop at the center') CASE DEFAULT WRITE(*,402) ; WRITE(2,402) 402 FORMAT('Sortie car T > t_stop au centre') END SELECT ELSEIF(zams .AND. arret == 'zams')THEN SELECT CASE(langue) CASE('english') WRITE(*,1403) ; WRITE(2,1403) 1403 FORMAT('Stop at the ZAMS stage') CASE DEFAULT WRITE(*,403) ; WRITE(2,403) 403 FORMAT('Sortie sur la ZAMS') END SELECT ELSEIF(post .AND. arret == 'post')THEN SELECT CASE(langue) CASE('english') WRITE(*,1404) ; WRITE(2,1404) 1404 FORMAT('stop at the end of the ZAMS stage (TAMS)') CASE DEFAULT WRITE(*,404) ; WRITE(2,404) 404 FORMAT('Sortie à la fin de la ZAMS (TAMS)') END SELECT ELSEIF(cohe .AND. arret == 'cohe')THEN SELECT CASE(langue) CASE('english') WRITE(*,1405) ; WRITE(2,1405) 1405 FORMAT('Stop at the He burning onset') CASE DEFAULT WRITE(*,405) ; WRITE(2,405) 405 FORMAT('Sortie à la combustion de He') END SELECT ELSEIF(coca .AND. arret == 'coca')THEN SELECT CASE(langue) CASE('english') WRITE(*,1406) ; WRITE(2,1406) 1406 FORMAT('Stop at the C combustion onset') CASE DEFAULT WRITE(*,406) ; WRITE(2,406) 406 FORMAT('Sortie à la combustion de C') END SELECT ELSEIF(cone .AND. arret == 'cone')THEN SELECT CASE(langue) CASE('english') WRITE(*,1414) ; WRITE(2,1414) 1414 FORMAT('Stop at the Ne combustion onset') CASE DEFAULT WRITE(*,414) ; WRITE(2,414) 414 FORMAT('Sortie à la combustion du Ne') END SELECT ELSEIF(coox .AND. arret == 'coox')THEN SELECT CASE(langue) CASE('english') WRITE(*,1412) ; WRITE(2,1412) 1412 FORMAT('Stop at the O combustion onset') CASE DEFAULT WRITE(*,412) ; WRITE(2,412) 412 FORMAT('Sortie à la combustion de O') END SELECT ELSEIF(cosi .AND. arret == 'cosi')THEN SELECT CASE(langue) CASE('english') WRITE(*,1415) ; WRITE(2,1415) 1415 FORMAT('Stop at the Si combustion onset') CASE DEFAULT WRITE(*,415) ; WRITE(2,415) 415 FORMAT('Sortie à la combustion de Si') END SELECT ELSEIF(LOG_teff > 0.d0 .AND. lteff > LOG_teff)THEN SELECT CASE(langue) CASE('english') WRITE(*,1407) ; WRITE(2,1407) 1407 FORMAT('Stop because LOG_teff was crossed') CASE DEFAULT WRITE(*,407) ; WRITE(2,407) 407 FORMAT('Sortie car LOG_teff est dépassé') END SELECT ELSEIF(LOG_teff < 0.d0 .AND. lteff < -LOG_teff)THEN SELECT CASE(langue) CASE('english') WRITE(*,1408) ; WRITE(2,1408) 1408 FORMAT('Stop because, LOG_teff was crossed') CASE DEFAULT WRITE(*,408) ; WRITE(2,408) 408 FORMAT('Sortie car on est en dessous de LOG_teff') END SELECT ELSEIF(lhe_stop)THEN SELECT CASE(langue) CASE('english') WRITE(*,1411)he_core ; WRITE(2,1411)he_core 1411 FORMAT('Stop because the he_core extends until ',es10.3,'Msol') CASE DEFAULT WRITE(*,411)he_core ; WRITE(2,411)he_core 411 FORMAT('Sortie car le noyau d''hélium atteint ',es10.3,'Msol') END SELECT ENDIF CALL date_and_time(date,time,zone,values) SELECT CASE(langue) CASE('english') titre='End of the calculation on '//date(7:8)// 1 ' '//TRIM(month(values(2)))// 2 ' '//date(1:4)//' à '//time(1:2)//'h'//time(3:4) CASE DEFAULT titre='Fin du calcul le '//date(7:8)//' '//TRIM(mois(values(2)))// 1 ' '//date(1:4)//' à '//time(1:2)//'h'//time(3:4) END SELECT WRITE(*,502)titre ; WRITE(2,502)titre 502 FORMAT(a) c la litanie des noms des fichiers de sortie CALL ecrit2 c on garde le modèle final, en fixant le pas temporel extrapolé dts au pas dt dts=dt c écriture du fichier mon_modèle_B.dat modelb=TRIM(nom_fich2)//'_B.dat' OPEN(unit=26,form='unformatted',status='unknown',file=TRIM(modelb)) WRITE(26)ne,m_qs,n_qs,knot,nchim,n_ch,m_ch,knotc,Krot,nrot,n_rot, 1 m_rot,knotr,n_tds,knot_tds,mtot,alpha,w_rot,lim_ro,diffusion, 2 rot_solid,precisione,en_m23,f_eos,f_opa, 3 nom_ctes,nom_pertm,nom_pertw,nom_tdetau,nom_atm,nom_conv,nom_nuc, 4 nom_nuc_cpl,nom_diffm,nom_difft,nom_diffw,nom_etat,nom_opa, 5 nom_elem,bp,q,qt,chim,mc,mct,rota,mrot,mrott,tds,x_tds,xt_tds, 6 m_stat,r_stat,m_zc,r_zc,r_ov,age,dt,dts,mstar,rstar,mw_tot,wrot, 7 jlim,lconv,lim,model_num CLOSE(unit=26) WRITE(*,46)TRIM(modelb) c PRINT*,all_rep,model_num,modelb ; WRITE(*,2000)dt,dts ; PAUSE'.dat' c dernières écritures sur le fichier mon_modele.TR.csv et fermeture IF(l_tr)THEN !fichier mon_modele.TR.csv c écritures dans mon_modele.TR.csv CALL csv_tr c fermeture du fichier mon_modele.TR.csv CLOSE(unit=56) ELSEIF(l_lb)THEN c écritures dans mon_modele.LB.csv CALL csv_lb(a,dnu02,dnu13,nu0,lb_p0,compg) c PRINT*,'écriture1 dans LB.csv, p0=',lb_p0 c fermeture du fichier mon_modele.LB.csv CLOSE(unit=56) ENDIF c écritures dans mon_modele .LR et .HR IF(l_hr)THEN WRITE(53,80)age,lteff,log_l,LOG10(rstar),model_num c WRITE(*,80)age,lteff,log_l ; PAUSE'cesam hr4' ELSE WRITE(53,52)age,nchim,Krot,lim,model_num,lconv(1:lim) WRITE(53,2002)lteff,log_l,LOG10(rstar),mstar, 1 (mstar-m_zc(i),r_zc(i),r_ov(i),i=1,lim) WRITE(53,2003)(nom_elem(i),compg(i,1),compg(i,n_qs),i=1,nchim) WRITE(53,2007)nu0,dnu02,dnu13,a,lb_p0 c PRINT*,'écriture1 dans .HR, p0=',lb_p0 ENDIF c compléments avec rotation IF(Krot >=3)THEN !pour Krot = 3, 4, 5, 6 WRITE(53,2003)nom_vwrot,rstar*rsol*1.d-5*w(n_qs),w(n_qs) ENDIF c si rstar < r_zc(lim)' IF(rstar < r_zc(lim))THEN PRINT*,'ATTENTION!!! dans cesam 5253, rstar < r_zc(lim)' WRITE(*,2002)rstar,r_zc(1:lim) !; PAUSE'cesam 5308' ENDIF c fermeture du fichier mon_modele.HR CLOSE(unit=53) c limites des shells H et He IF(ecrHHe)THEN WRITE(54,9999)model_num,age,msh,mshe CLOSE(unit=54) ENDIF c fichier des Céphéïdes IF(cephe .AND. t(1) > 1.d7)THEN CALL pulsation(l(n_qs),teff,typ,mstar,puls,period) IF(puls)THEN WRITE(55,9998)model_num,age,l(n_qs),mstar,period,rstar,teff CLOSE(unit=55) ENDIF ENDIF c fermetures du fichier *.lis CLOSE(unit=2) ; EXIT Bev ENDIF !on sort c output ASCII éventuel c PAUSE'appel ecrit_ascii' IF(nom_output /= 'no_output')CALL ecrit_ascii ENDIF !ecrit c Alerte avec diffusion du moment cinétique si \tilde g \ g > 20% IF(COUNT(Krot == Ktest) >= 1)THEN IF(gtmax > 0.2d0)THEN SELECT CASE(langue) CASE('english') WRITE(*,1808)nom_thw,gtmax,igtmax,rgtmax,mgtmax WRITE(2,1808)nom_thw,gtmax,igtmax,rgtmax,mgtmax CASE DEFAULT WRITE(*,808)nom_thw,gtmax,igtmax,rgtmax,mgtmax WRITE(2,808)nom_thw,gtmax,igtmax,rgtmax,mgtmax END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(*,1809)gtmax,igtmax,rgtmax,mgtmax WRITE(2,1809)gtmax,igtmax,rgtmax,mgtmax CASE DEFAULT WRITE(*,809)gtmax,igtmax,rgtmax,mgtmax WRITE(2,809)gtmax,igtmax,rgtmax,mgtmax END SELECT ENDIF ENDIF c écriture, totale si ecrit, partielle si .not.ecrit c PRINT*,ecrit ; PAUSE'ecrit' CALL list(a,alfa,anub8,anube7,anun13,anuo15,anupep,anupp,beta, 1 compg,cp,delta,dcapdr,dcapdt,degene,depsdr,depsdt,dnu02,dnu13, 2 drdt,d2p,d2ro,chaine,convec,ecrit,epsilon,gam,gamma,gradad,grada_atm, 3 gradconv,gradc_atm,gradrad,gradr_atm,hp,i_cno,i_gr,i_pp,i_3a, 4 kap,l,lb_p0,m,mu,mue,m_atm,nu0,p,pt,pt_atm,p_atm,r,ro,ro_atm,r_atm,t,tau, 5 teff,tx_ioni,t_atm,u,vaissala,w,z) c sorties sur écran c PRINT*,mstar,mtot zsx=z(n_qs)/compg(1,n_qs) pertmm=mstar/mtot-1.d0 IF(ihe4 <= 1)THEN WRITE(*,205)age,lteff,log_l,LOG10(rstar), 1 LOG10(mstar/rstar**2*grav_sol),pt(1),t(1),ro(1),compg(1,1), 2 i_pp,i_cno,i_3a,i_gr,pertmm,mstar,TRIM(chaine),TRIM(instb),period, 3 zsx ELSEIF(w(n_qs) > 0.d0)THEN WRITE(*,207)age,lteff,log_l,LOG10(rstar), 1 LOG10(mstar/rstar**2*grav_sol),pt(1),t(1),ro(1),compg(1,1), 2 i_pp,i_cno,i_3a,i_gr,compg(ihe4,1),pertmm,mstar,TRIM(chaine), 3 2.d0*pi/24.d0/3600.d0/w(n_qs),rstar*w(n_qs)*rsol*1.d-5, 4 TRIM(instb),period,zsx ELSE WRITE(*,206)age,lteff,log_l,LOG10(rstar), 1 LOG10(mstar/rstar**2*grav_sol),pt(1),t(1),ro(1),compg(1,1), 2 i_pp,i_cno,i_3a,i_gr,compg(ihe4,1),pertmm,mstar,TRIM(chaine), 3 TRIM(instb),period,zsx ENDIF c écritures sur les fichiers mon_modele.HR et mon_modele.LR IF(ecrHR)THEN c écritures dans mon_modele.TR.csv et mon_modele.LB.csv IF(l_tr)THEN CALL csv_tr ELSEIF(l_lb)THEN CALL csv_lb(a,dnu02,dnu13,nu0,lb_p0,compg) c PRINT*,'écriture2 dans LB.csv, p0=',lb_p0 ENDIF c écritures dans mon_modele.HR et .LR IF(l_hr)THEN WRITE(53,80)age,lteff,log_l,LOG10(rstar),model_num 80 FORMAT(es22.15,3es13.6,i5) c WRITE(*,80)age,lteff,log_l ; PAUSE'cesam hr3' ELSE WRITE(53,52)age,nchim,Krot,lim,model_num,lconv(1:lim) 52 FORMAT(es22.15,3i3,i5,20l2) WRITE(53,2002)lteff,log_l,LOG10(rstar),mstar, 1 (mstar-m_zc(i),r_zc(i),r_ov(i),i=1,lim) 2002 FORMAT(6es13.6) WRITE(53,2003)(nom_elem(i),compg(i,1),compg(i,n_qs),i=1,nchim) 2003 FORMAT(a4,2es12.5) WRITE(53,2007)nu0,dnu02,dnu13,a,lb_p0 2007 FORMAT(5es10.3) c PRINT*,'écriture2 dans .HR, p0=',lb_p0 ENDIF c compléments avec rotation IF(Krot >=3)THEN !pour Krot = 3, 4, 5, 6 WRITE(53,2003)nom_vwrot,rstar*rsol*1.d-5*w(n_qs),w(n_qs) ENDIF c si rstar < r_zc(lim) IF(rstar < r_zc(lim))THEN PRINT*,'ATTENTION!!! dans cesam 5417, rstar < r_zc(lim)' WRITE(*,2002)rstar,r_zc(1:lim) !; PAUSE'cesam 5417' ENDIF c on écrira désormais sur le *.HR (et le *.TR) ELSE ecrHR=.TRUE. ENDIF c limites des shells H et He IF(ecrHHe)WRITE(54,9999)model_num,age,msh,mshe 9999 FORMAT(i5,es22.15, 2es10.3) c fichier de Céphéïdes et ASCII si on est dans la bande d'instabilité IF(cephe .AND. t(1) > 1.d7)THEN CALL pulsation(l(n_qs),teff,typ,mstar,puls,period) IF(puls)THEN WRITE(55,9998)model_num,age,l(n_qs),mstar,period,rstar,teff 9998 FORMAT(i5,es22.15,5es10.3) c PRINT*,'pulsation, all_cephe',puls,all_cephe ; PAUSE'cesam puls' IF(all_cephe)CALL ecrit_ascii ENDIF ENDIF c -----------------sorties fin--------------------------------------- c on forme éventuellement le fichier mon_modelexxxx*.osc c PRINT*,all_osc,ecr_inter ; PAUSE'cesam ecrit' IF(all_osc .AND. ecr_inter)CALL ecrit_ascii c on garde éventuellement xxxxxmon_modele_B.rep IF(all_rep .AND. ecr_inter)THEN WRITE(number,70)model_num 70 FORMAT(i5.5) modelb=number//'-'//TRIM(nom_fich2)//'_B.rep' OPEN(unit=26,form='unformatted',status='unknown',file=TRIM(modelb)) WRITE(26)ne,m_qs,n_qs,knot,nchim,n_ch,m_ch,knotc,Krot,nrot,n_rot, 1 m_rot,knotr,n_tds,knot_tds,mtot,alpha,w_rot,lim_ro,diffusion, 2 rot_solid,precisione,en_m23,f_eos,f_opa, 3 nom_ctes,nom_pertm,nom_pertw,nom_tdetau,nom_atm,nom_conv,nom_nuc, 4 nom_nuc_cpl,nom_diffm,nom_difft,nom_diffw,nom_etat,nom_opa, 5 nom_elem,bp,q,qt,chim,mc,mct,rota,mrot,mrott,tds,x_tds,xt_tds, 6 m_stat,r_stat,m_zc,r_zc,r_ov,age,dt,dts,mstar,rstar,mw_tot,wrot, 7 jlim,lconv,lim,model_num CLOSE(unit=26) WRITE(*,46)TRIM(modelb) ENDIF c on garde toujours le modèle mon_modele_B.rep c WRITE(*,2000)age,dt,dts ; PAUSE'age,dt,dts' IF(.NOT.(list_rep .OR. list_pms .OR. list_zams))THEN modelb=TRIM(nom_fich2)//'_B.rep' OPEN(unit=26,form='unformatted',status='unknown',file=TRIM(modelb)) WRITE(26)ne,m_qs,n_qs,knot,nchim,n_ch,m_ch,knotc,Krot,nrot,n_rot, 1 m_rot,knotr,n_tds,knot_tds,mtot,alpha,w_rot,lim_ro,diffusion, 2 rot_solid,precisione,en_m23,f_eos,f_opa, 3 nom_ctes,nom_pertm,nom_pertw,nom_tdetau,nom_atm,nom_conv,nom_nuc, 4 nom_nuc_cpl,nom_diffm,nom_difft,nom_diffw,nom_etat,nom_opa, 5 nom_elem,bp,q,qt,chim,mc,mct,rota,mrot,mrott,tds,x_tds,xt_tds, 6 m_stat,r_stat,m_zc,r_zc,r_ov,age,dt,dts,mstar,rstar,mw_tot,wrot, 7 jlim,lconv,lim,model_num CLOSE(unit=26) WRITE(*,46)TRIM(modelb) c PRINT*,all_rep,model_num,modelb ; WRITE(*,2000)dt,dts c PRINT*,n_ch,knotc,TRIM(modelb) c WRITE(*,2000)mct(1:4),mct(knotc-3:knotc) c PAUSE'modèle .rep' ELSE list_pms=.FALSE. ; list_rep=.FALSE. ; list_zams=.FALSE. ENDIF c évolution avec le pas temporel dt, age--> age+dt------------------ c WRITE(*,2000)dt,dts ; PAUSE'cesam, avant resout' CALL resout(1,dt,dts) ; age=age+dt ; dtp=dt ENDDO Bev c PAUSE'fin de cesam' IF(age > 0.d0)THEN CALL des(.TRUE.,dtp,teff) ELSE CALL des(.TRUE.,0.d0,teff) ENDIF CALL sortie('cesam 49') CONTAINS INCLUDE'ecrit_ascii.f' END SUBROUTINE cesam