c******************************************************************** SUBROUTINE lim_zc(new_rep,new_nqs) c routine private du module mod_static c répartition des couches et détermination des limites ZR/ZC. c Le nombre est modifié si new_n /= n c Auteur: P.Morel, Département J.D. Cassini, O.C.A., CESAM2k c 17 09 96 : correction bug dans déclaration de vt; signalé par D. Cordier c 25 08 97 : mise en place des variables eulériennes c 19 11 99 : suppression de nh1, nhe1, nhe2, lamb c 20 02 00 : écriture si non affinement de la loc. des lim. des ZR/ZC c entrées c new_rep=.TRUE. : il faut redistribuer les couches c new_nqs : nouveau nombre de couches c---------------------------------------------------------------- USE mod_donnees, ONLY: alpha, aradia, clight, cpturb, 1 diffusion, dn_fixe, en_m23, 2 g, grad_ovi, grad_ovs, grille_fixe, ihe4, Ipg, jpz, Krot, 3 langue, loc_zc, ledoux, lsol, l_fac, msol, mtot, m_ch, m_qs, nchim, ne, 4 no_discon, nrot, n_min_ZC, ord_qs, ord_rot, ovsht, ovshti, ovshts, 5 pi, pnzc, precision, psi0, pturb, rsol, r_qs, t_inf USE mod_kind USE mod_numerique, ONLY : bsp1dn, entre, inside, newspl, noedif, 1 no_croiss, zoning USE mod_opa, ONLY : dehors USE mod_variables, ONLY: bp, chim, dim_qs, dv_tr, idm6, 1 id_conv, if_conv, inter, jlim, lconv, knot, knotc, knotr, lim, 2 mc, mct, mc_fixe, mrot, mrott, mstar, m_zc, m_zc23, m_stat, 3 nc_fixe, n_ch, n_qs, n_rot, omega_jpz, q, qt, rota, 4 rstar, r_stat, r_ov, r_zc, r_zc_conv, sortie, tot_conv, tot_rad, wrot IMPLICIT NONE INTEGER, INTENT(in) :: new_nqs LOGICAL, INTENT(inout) :: new_rep REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: new_b REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: dfrot, d_grad, frot REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: new_q, new_qt, esp, new REAL (kind=dp), DIMENSION(2*pnzc) :: all_rap, rap REAL (kind=dp), DIMENSION(nchim) :: comp, dcomp, depsx, dmu_x REAL (kind=dp), DIMENSION(ne) :: dfdq, f REAL (kind=dp), DIMENSION(5) :: epsilon REAL (kind=dp) :: alfa, beta, pt, dlpp, p, t, r, l, m, rn, 1 grp, gr, j_new, j1, j2, j3, kip, drot, ro, dkapt, kap, gradad, 2 drop, dkapp, hp, drox, u, dup, dut, dux, grad, dgradpt, dgradp, 3 dgradt, dgradx, dgradm, dgradl, dgradr, dgradlpp, gam, dgampt, 4 dgamp, dgamt, dgamx, dgamm, dgaml, dgamr, dgamlpp, depsp, 5 depst, dkapx, delta, deltap, deltat, deltax, cp, dcpp, 6 dcpt, dcpx, dgradadp, dgradadt, dgradadx, dhppt, dhpp, dhpt, 7 dhpx, dhpr, dhpm, gradrad, grad_mu, gamma1, mu, mue, w INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: jlimp INTEGER, DIMENSION(1) :: idist_max INTEGER, SAVE :: dlt_n=60, lc=1 INTEGER :: compt, i, j, limp, lq, new_knot LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:) :: lconvp LOGICAL, SAVE :: der, init=.TRUE. LOGICAL :: radiatif c------------------------------------------------------------------- 2000 FORMAT(8es10.3) c-----------------------initialisations début----------------------- IF(init)THEN init=.FALSE. der=cpturb < 0.d0 !der=T. on tient compte de dln Pgaz/dln Ptot SELECT CASE(langue) CASE('english') WRITE(2,1020)loc_zc ; WRITE(usl_static,1020)loc_zc 1020 FORMAT(/,'-----Parameters concerning the convectives zones----',/, 1 'Accuracy of the localization of limits RZ/CZ :',es10.3) CASE DEFAULT WRITE(2,20)loc_zc ; WRITE(usl_static,20)loc_zc 20 FORMAT(/,'----Paramètres relatifs aux zones convectives----',/, 1 'Précision de la localisation des limites ZR/ZC :',es10.3) END SELECT c les limites ZR/ZC de mélange doivent être situées sur un bord de l'intervalle IF(l_fac .AND. d_zc < 50)THEN WRITE(usl_static,35)d_zc 35 FORMAT('Un bord de la zone mélangée sera approché à mieux que ',i2, 1 "% de chaque limite ZR/ZC") c PRINT*,'d_zc=', d_zc ; PAUSE'd_zc ELSE WRITE(usl_static,27) 27 FORMAT("Sans essai de rapprochement d'un point de grille des limites ZR/ZC") ENDIF IF(no_discon)THEN WRITE(6,42) 42 FORMAT('Composition chimique continue aux limites ZR/ZC') ELSE WRITE(6,47) 47 FORMAT('Composition chimique discontinue aux limites ZR/ZC') ENDIF SELECT CASE(Krot) CASE(3,4,5) ALLOCATE(dfrot(nrot),frot(nrot)) END SELECT c dans lit_nl ovshti/s ont perdu leur signes, grad_ovi/s T/F : grad_ad/grad_rad c dans les zones overshootées IF(ovsht)THEN IF(jpz)THEN IF(ovshti > 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1009)INT(ovshti*100) ; WRITE(2,1009)INT(ovshti*100) 1009 FORMAT('overshooting JPZ beneath the CZ : ',i3,'% Hp/Kip',/) CASE DEFAULT WRITE(usl_static,9)INT(ovshti*100) ; WRITE(2,9)INT(ovshti*100) 9 FORMAT('overshooting inférieur de JPZ des ZC :',i3,'% Hp/Kip') END SELECT ENDIF IF(ovshts > 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1006)INT(ovshts*100) ; WRITE(2,1006)INT(ovshts*100) 1006 FORMAT('overshooting JPZ above the CZ : :',i3,'% R noyau',/) CASE DEFAULT WRITE(usl_static,6)INT(ovshts*100) ; WRITE(2,6)INT(ovshts*100) 6 FORMAT('overshooting supérieur de JPZ des ZC :',i3, 1 '% R noyau',/) END SELECT ENDIF ELSE IF(ovshti > 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1007)INT(ovshti*100) ; WRITE(2,1007)INT(ovshti*100) 1007 FORMAT('overshooting beneath the CZ :',i3,'% Hp') CASE DEFAULT WRITE(usl_static,7)INT(ovshti*100) ; WRITE(2,7)INT(ovshti*100) 7 FORMAT('overshooting inférieur des ZC :',i3,'% Hp') END SELECT IF(grad_ovi)THEN WRITE(usl_static,60) ; WRITE(2,60) 60 FORMAT('grad=grad_ad dans les parties overshootées',/) ELSE WRITE(usl_static,61) ; WRITE(2,61) 61 FORMAT('grad=grad_rad dans les parties overshootées',/) ENDIF ENDIF IF(ovshts > 0.d0)THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1008)INT(ovshts*100) ; WRITE(2,1008)INT(ovshts*100) 1008 FORMAT('overshooting above the CZ :',i3,'% Hp') CASE DEFAULT WRITE(usl_static,8)INT(ovshts*100) ; WRITE(2,8)INT(ovshts*100) 8 FORMAT('overshooting supérieur des ZC :',i3,'% Hp') END SELECT IF(grad_ovs)THEN WRITE(usl_static,60) ; WRITE(2,60) ELSE WRITE(usl_static,61) ; WRITE(2,61) ENDIF ENDIF ENDIF ELSE SELECT CASE(langue) CASE('english') WRITE(usl_static,1005) ; WRITE(2,1005) 1005 FORMAT('Model without overshooting') CASE DEFAULT WRITE(usl_static,5) ; WRITE(2,5) 5 FORMAT('Modèle sans overshooting') END SELECT ENDIF c allocations initiales ALLOCATE(fac(n_qs-1),d_grad(n_qs),jlimp(n_qs),lconvp(n_qs)) c PRINT*,l_fac,d_zc ; PAUSE'lim_zc initialisation' ENDIF c----------------fin des initialisations---------------------------- c PAUSE'début recherche limites ZR/ZC' SELECT CASE(langue) CASE('english') WRITE(usl_static,1010) 1010 FORMAT(/,'----- Search of the limits RZ/CZ (begin) --------',/) CASE DEFAULT WRITE(usl_static,10) 10 FORMAT(/,'------ Recherche des limites ZR/ZC (début) --------',/) END SELECT c on redéfinit la répartition avant la recherche des limites ZR/ZC c si le modèle n'est pas estimé "bien réparti" ou si le nb. c de couches du modèle quasi. stat. doit changer c on répartit les couches en tenant compte (éventuellement) de L, pour l'ajustement c mises à 0 all_rap=0.d0 ; rap=0.d0 c 3 cas se présentent: c 1 - début du processus itératif, reprise, ou modification du nombre de couches c appel de resout avec new_rep=.TRUE. ou new_nqs /= n_qs. Le modèle sera redistribué, c pour le TdS on utilise la tabulation c 2 - processus itératif en cours, le nb. de couches est fixé, pour tenir compte des c nouvelles valeurs des variables il faut réinitialiser la répartition, les limites c ZR/ZC et les fac, on entre avec new_rep=.TRUE. c 3 - processus itératif en cours, au delà de ini1, il n'y plus d'appel à lim_zc c PRINT*,'entrée lim_zc',new_rep,new_nqs,n_qs I1: IF(new_rep .OR. new_nqs /= n_qs) THEN c les vecteurs fac, d_grad et jlimp, n'ont plus la même dimension DEALLOCATE(fac,d_grad,jlimp,lconvp) c cas 1, on met dans esp la répartition précédente, cette portion sera active à c l'initialisation du pas temporel si le nombre de couches ou que la répartition c doit être modifiée. ALLOCATE(esp(n_qs)) DO i=1,n_qs esp(i)=espmt(q(i)) ENDDO c on dispose les couches pour assurer une répartition approximativement c uniforme de la fonction de répartition sur new_nqs couches ALLOCATE(new_q(new_nqs)) CALL zoning(n_qs,q,esp,new_nqs,new_q) !new_q(new_nqs) nouveaux q c esp est désormais inutile DEALLOCATE(esp) c dans new_b on se place aux new_nqs points new_q c on calcule la nouvelle fonction de répartition ALLOCATE(new_b(ne,new_nqs)) DO i=1,new_nqs CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,new_q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Problème localisé en 2 dans lim_zc' new_b(1:ne,i)=f(1:ne) ENDDO c on spline new_b d'ordre 2 sur les nouveaux q sur new_nqs points new_q=(/ (i, i=1,new_nqs) /) ALLOCATE(new_qt(new_nqs+ord_qs)) CALL bsp1dn(ne,new_b,new_q,new_qt,new_nqs,ord_qs,new_knot, 1 .FALSE.,new_q(1),lq,f,dfdq) IF(no_croiss)THEN PRINT*,'Arrêt 2 dans lim_zc' ; STOP ENDIF c commentaires IF(n_qs /= new_nqs)THEN SELECT CASE(langue) CASE('english') WRITE(2,1015)n_qs,new_nqs ; WRITE(usl_static,1015)n_qs,new_nqs 1015 FORMAT('Change of the number of shells :',i4,'-->',i4,/) CASE DEFAULT WRITE(2,15)n_qs,new_nqs ; WRITE(usl_static,15)n_qs,new_nqs 15 FORMAT('Modification du nombre de couches :',i4,'-->',i4,/) END SELECT n_qs=new_nqs ENDIF IF(new_rep)THEN SELECT CASE(langue) CASE('english') WRITE(2,1016)n_qs ; WRITE(usl_static,1016)n_qs 1016 FORMAT('Updated distribution in lim_zc, with the number of shells :',i4) CASE DEFAULT WRITE(2,16)n_qs ; WRITE(usl_static,16)n_qs 16 FORMAT('Répartition actualisée dans lim_zc, nombre de couches :',i4) END SELECT ENDIF c dimension de la base splines pour l'intégration équi.quasi.stat. dim_qs=(n_qs-1)*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*,'Arrêt 3 dans lim_zc' ; STOP ENDIF c PRINT*,'lim_zc ord_qs',ord_qs,knot,dim_qs,n_qs,ne c on place la solution initiale dans la base de q, bp_t ==> bp DEALLOCATE(bp) ALLOCATE(bp(ne,dim_qs)) CALL newspl(ne,new_q,new_qt,new_knot,ord_qs,qt,knot,ord_qs,new_b,bp) c PAUSE'lim_zc, newspl' c détermination de psi bp(6,:)=espmt(q(10))-espmt(q(9)) c WRITE(*,2000)bp(6,1:8) ; PAUSE'lim_zc, bp(6)' c new_b, new_q, new_qt, inutiles, r_stat et m_stat sont à modifier DEALLOCATE(new_b,new_q,new_qt,r_stat,m_stat) c adaptation de fac, d_grad, r_stat et m_stat à la nouvelle répartition ALLOCATE(fac(n_qs-1),d_grad(n_qs),jlimp(n_qs),r_stat(n_qs),m_stat(n_qs), 1 lconvp(n_qs)) c tabulation de r_stat et m_stat qui accélèrent la recherche dans le ssp. inter c r_stat= r^2; m_stat=m^2/3 ou r_stat=r, m_stat=m^1/3 DO i=1,n_qs c CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lc,f,dfdq) r_stat(i)=bp(3,m_qs*(i-1)+1) ; m_stat(i)=bp(5,m_qs*(i-1)+1) ENDDO DO i=2,n_qs IF(r_stat(i-1) > r_stat(i) .OR. bp(5,i-1) > bp(5,i))THEN PRINT*,'lim_zc, tabul, inversion de rayon ou de masse, couche n°',i PRINT*,' R(i-1:i), M(i-1:i)',r_stat(i-1),r_stat(i),bp(5,i-1),bp(5,i) ENDIF ENDDO c test pour changement de grille fixe pour la composition chimique c il y a changement de grille fixe quand c new_nqs diffère de nc_fixe de plus dn_fixe (5% environ fixé dans cesam) IF(grille_fixe)THEN IF(ABS(REAL(new_nqs-nc_fixe,dp)/REAL(new_nqs,dp)) > dn_fixe)THEN WRITE(usl_static,30)nc_fixe,n_qs ; WRITE(2,30)nc_fixe,n_qs 30 FORMAT(/,'changement de grille fixe: nc_fixe=',i4,' ==> ',i4,/) DEALLOCATE(mc_fixe) ; ALLOCATE(mc_fixe(n_qs)) nc_fixe=n_qs ; mc_fixe=m_stat IF(.NOT.en_m23)mc_fixe=mc_fixe**2.d0 SELECT CASE(langue) CASE('english') WRITE(usl_static,1050) 1050 FORMAT(/,'Change of the fixed grid') CASE DEFAULT WRITE(usl_static,50) 50 FORMAT(/,'Changement de grille fixe') END SELECT ENDIF ENDIF new_rep=.FALSE. ENDIF I1 c------------------recherche des limites ZR/ZC---------------------------- c au centre R, L, M =0 bp(3:5,1)=0.d0 c recherche les limites zones convectives / zones radiatives c en fonction de q, détermination des facteurs de répartition c overshooting inférieur et supérieur : extension de la zone mélangée c ou overshooting inférieur PC de JPZ avec saut du gradient c modèle totalement convectif: lim=1, jlim(1)=n, lconv(1)=.FALSE. c modèle totalement radiatif: lim=0, jlim(i)=-100, lconv(i)=.FALSE. c lconv(i)=..TRUE. en r_zc(i) on passe de ZC à ZR c initialisations r_ov=-1.d10 ; r_zc=-1.d10 ; m_zc=-1.d10 ; dv_tr=-1.d10 c lim: nb de limites ZR/ZC, jlim: indices des lim. ZR/ZC c lconv(i)=T/F début/fin de ZC en jlim(i), i=1,lim lim=0 ; jlim=-100 ; lconv=.TRUE. c tot_conv/tot_rad modèle totalement convectif/radiatif tot_conv=.FALSE. ; tot_rad=.FALSE. c indice au dessus duquel la masse ne varie plus: 1-m(idm6)/mstar < 1.d-6 idm6=n_qs c calcul de grad-grad c PRINT*,n_ch,m_ch,knotc,n_qs,nchim ; WRITE(*,2000)mc c WRITE(*,2000)mct ; PRINT*,mstar ; PAUSE'début lim_zc' c PRINT*,'n_qs',n_qs B5: DO i=1,n_qs !calcul de grad-grad CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,q(i),lq,f,dfdq) IF(no_croiss)PRINT*,'Problème localisé en 4 dans lim_zc' c WRITE(*,2000)f pt=EXP(f(1)) dlpp=1.d0 IF(pturb)THEN !avec pression turbulente p=EXP(f(Ipg)) IF(der)dlpp=dfdq(Ipg)/dfdq(1) !dlpp=dln Pgaz/dln Ptot ELSE !sans pression turbulente p=pt ENDIF t=EXP(f(2)) IF(en_m23)THEN r=SQRT(ABS(f(3))) ; l=SQRT(ABS(f(4)))**3 ; m=SQRT(ABS(f(5)))**3 ELSE r=ABS(f(3)) ; l=f(4) ; m=ABS(f(5))**3 ; f(5)=f(5)**2 ENDIF c indice au dessus duquel la masse ne varie pratiquement plus IF(1.d0-m/mstar < 1.d-6)idm6=MIN(idm6,i) c composition chimique CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(f(5),mc(n_ch))),lc,comp,dcomp) IF(no_croiss)THEN PRINT*,'Problème localisé en 5 dans lim_zc',i,knotc,n_ch WRITE(*,2000)q(i),mc(1),f(5),mc(n_ch),mct(knotc) !; PAUSE'en 5' ENDIF c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w=wrot CASE(3,4,5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lc,frot,dfrot) IF(no_croiss)PRINT*,'Problème localisé en 6 dans lim_zc' w=frot(1) CASE(6) CALL omega_jpz(r,w) END SELECT c IF(t < 8.d4)WRITE(*,2000)p,t,q(i) c IF(t > 3.d6)WRITE(*,2000)q(i),p,t,m,l,r,comp(1),comp(3) d_grad(i)=dgrad(pt,p,t,dlpp,comp,m,l,r,dcomp,w) IF(dehors)RETURN ENDDO B5 !i c WRITE(*,2000)d_grad ; PAUSE'lim_zc, d_grad' c r à l'extérieur, i.e. dernière valeur obtenue dans la boucle B5 rn=r c------------identification et gestion des ZR et ZC début-------------- c les limites ZR/ZC ayant au moins n_min_ZC couches, on posera d_grad(1)=d_grad(2) d_grad(1)=d_grad(2) c état convectif aux points de grille DO i=1,n_qs lconvp(i)=d_grad(i) > 0.d0 !vrai : ZC, faux : ZR ENDDO c suppression des limites consécutives, on garde les débuts de ZC DO i=2,n_qs-1 IF(lconvp(i-1) .AND. .NOT.lconvp(i) .AND. lconvp(i+1))THEN WRITE(usl_static,29)i-1,i+1,lconvp(i-1:i+1) 29 FORMAT("Suppression de lim. ZR/ZC consécutives intervalle [",i4,',',i4,']', 1 ', lconv:',3l1) lconvp(i)=lconvp(i-1) c PRINT*,'modif(i-1,i+2)',lconvp(i-1:i+2) !; PAUSE'suppression' ENDIF ENDDO c recherche des limites ZR/ZC limp=0 !nombre de limites DO i=2,n_qs IF(lconvp(i) .NEQV. lconvp(i-1))THEN limp=limp+1 jlimp(limp)=i WRITE(usl_static,28)limp,i-1,i,lconvp(i) 28 FORMAT('limite ZR/ZC n°',i3,", détectée dans l'intervalle [",i4,',',i4,']', 1 ', lconv:',l1) ENDIF ENDDO c PRINT* ; PRINT*,'0 limp=',limp,jlimp(1:limp) c suppression de limites trop centrales Ilp: IF(limp > 0)THEN i=1 B2: DO IF(jlimp(i) > n_min_ZC)EXIT B2 IF(jlimp(i) <= n_min_ZC)THEN !limites sur premières couches WRITE(usl_static,*)'suppression d''une limite ZR/ZC trop centrale n°',i limp=limp-1 ; IF(limp <= 0)EXIT B2 DO j=1,limp jlimp(j)=jlimp(j+1) lconvp(j)=lconvp(j+1) ENDDO c PRINT*,'1 limp,jlimp',limp,jlimp(1:limp) ELSE i=i+1 ; IF(i > limp)EXIT B2 ENDIF ENDDO B2 c PRINT*,'2 limp=',limp,jlimp(1:limp) c suppression de limites trop proches i=2 DO WHILE (i <= limp) IF(jlimp(i)-jlimp(i-1) <= n_min_ZC)THEN !limites trop proches limp=limp-2 !suppression de i-1 et i WRITE(usl_static,46)i-1,i 46 FORMAT('limites trop proches, suppression des limites',2i4) c décalage DO j=i-1,limp jlimp(j)=jlimp(j+2) lconvp(j)=lconvp(j+2) ENDDO ELSE i=i+1 ENDIF !pour la suppression de limites trop rapprochées ENDDO ENDIF Ilp c limp --> lim lim=limp jlim(1:lim)=jlimp(1:limp) lconv(1:lim)=lconvp(1:limp) c au plus pnzc limites, (2 limites par ZC) I3: IF(lim > 2*pnzc-1)THEN PRINT*,'nombre de limites ZR/ZC=',lim,', supérieur à',pnzc*2-1 PRINT*,'limites ZR/ZC imprécises, floues ?' PRINT*,'difficulté d''initialisation ? erreur ? ' CALL sortie('lim_zc 1') ENDIF I3 c PRINT* c on précise la nature et les localisations des limites----------- DO i=1,lim gr=d_grad(jlim(i)) ; grp=d_grad(jlim(i)-1) c vérification IF(gr*grp > 0.d0)THEN WRITE(usl_static,43)i,grp,gr,jlim(i)-1,jlim(i) 43 FORMAT('ERREUR, limite ZR/ZC',i3,', grp=',es10.3,', gr =',es10.3, 1 'de même signe',', intervalle [',i4,',',i4,']') CALL sortie('lim_zc 2') ENDIF lconv(i)=grp < 0.d0 rap(i)=ABS(grp/gr) rap(i)=rap(i)/(1.d0+rap(i)) c écritures IF(lconv(i))THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1003)INT(rap(i)*100.d0),jlim(i)-1,jlim(i),i 1003 FORMAT('a CZ begins at',i3,'% between the shells #', 1 i5,' and #',i5,', limit RZ/CZ #',i3) CASE DEFAULT WRITE(usl_static,3)INT(rap(i)*100.d0),jlim(i)-1,jlim(i),i 3 FORMAT('début de ZC localisée à',i3,'% entre les couches n°', 1 i5,' et ',i5,', limite ZR/ZC n°',i3) END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(usl_static,1004)INT(rap(i)*100.),jlim(i)-1,jlim(i),i 1004 FORMAT('a CZ ends at ',i3,'% between the shells #', 1 i5,' and #',i5,', limit RZ/CZ #',i3) CASE DEFAULT WRITE(usl_static,4)INT(rap(i)*100.),jlim(i)-1,jlim(i),i 4 FORMAT('fin de ZC localisée à',i3,'% entre les couches n°', 1 i5,' et ',i5,', limite ZR/ZC n°',i3) END SELECT ENDIF ENDDO c PAUSE'lim' c on garde les rapports de localisations des ZC all_rap(1:lim)=rap(1:lim)*1.d2 c PRINT*,'all_rap' c WRITE(*,2000)rap(1:lim),all_rap(1:lim) c indices REAL des points de grille ALLOCATE(new(n_qs)) ; new=q c------------on précise les localisations des limites----------- c il y a une limite ZR/ZC entre jlim(i)-1 et jlim(i), i=1,lim c on détermine l'indice REAL j_new de la limite ZR/ZC DO i=1,lim c pour les limites trop externes on se place sur le point de grille le plus près à c cause de l'imprécision sur la masse IF(jlim(i) >= idm6 )THEN new(jlim(i))=jlim(i) c sinon on fixe la limite soit à droite soit à gauche c jlim(i) devient l'indice INTEGER le plus proche de la limite ELSE j_new=new(jlim(i)-1)+rap(i)*(new(jlim(i))-new(jlim(i)-1)) c PRINT*,'i:',i,',j_new:',j_new,',new:',new(jlim(i)-1),rap(i),new(jlim(i)) jlim(i)=NINT(j_new) ENDIF ENDDO c PRINT*,'limites ZR/ZC, lim=',lim c PRINT*,(jlim(i),new(jlim(i)),i=1,lim) !; PAUSE'lim_zc2' c-------------identification et gestion des ZR et ZC fin----------------- c initialisation des facteurs de répartition fac=1.d0 c nature complètement convectif/radiatif du modèle c si lim<=0 et que la différence des gradients est positive sur la c couche externe le modèle est complètement convectif Ilim: IF(lim <= 0)THEN IF(d_grad(n_qs) >= 0.d0)THEN tot_conv=.TRUE. lim=1 ; jlim(1)=n_qs lconv(1)=.FALSE. ; m_zc(1)=mstar ; r_zc(1)=rn SELECT CASE(langue) CASE('english') WRITE(usl_static,1040) 1040 FORMAT('Fully convective and mixed model') CASE DEFAULT WRITE(usl_static,40) 40 FORMAT('Modèle totalement convectif et mélangé') END SELECT c pour simplifier les algorithmes dans evol, les couches les plus externes c des modèles totalement radiatif seront mélangées ELSE tot_rad=.TRUE. lim=1 ; lconv(1)=.TRUE. ; jlim=n_qs-n_min_ZC SELECT CASE(langue) CASE('english') WRITE(usl_static,1041) 1041 FORMAT('Fully radiative model') CASE DEFAULT WRITE(usl_static,41) 41 FORMAT('Modèle totalement radiatif') END SELECT ENDIF c on affine, par dichotomie, la position des limites retenues, pour chaque limite, c encadrement entre q=j1 et q=j2 ELSE Ilim c pour chaque limite B4: DO i=1,lim c pour les limites dont la position est localisée en dessous de idm6 IF(jlim(i) >= idm6)THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,45)i 45 FORMAT('the limits ',i3, 1 ', and the following, too external, will be nor improved nor approched') CASE DEFAULT WRITE(usl_static,44)i,idm6 44 FORMAT("les limites d'indice >",i3,' localisées au delà de la couche n°' 1 ,i4,', trop externes, ne seront',/,'ni affinées ni approchées') END SELECT c limitation du nombre de limites ZR/ZC c lim=i EXIT B4 ENDIF c j1 --> dgrad(1) c j1=new(jlim(i)) j1=jlim(i)-1 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,j1,lq,f,dfdq) IF(no_croiss)PRINT*,'Problème localisé en 7 dans lim_zc' c WRITE(*,2000)f pt=EXP(f(1)) ; dlpp=1.d0 IF(pturb)THEN !avec pression turbulente p=EXP(f(Ipg)) IF(der)dlpp=dfdq(Ipg)/dfdq(1) !dlpp=dln Pgaz/dln Ptot ELSE !sans pression turbulente p=pt ENDIF t=EXP(f(2)) IF(en_m23)THEN r=SQRT(ABS(f(3))) ; l=SQRT(ABS(f(4)))**3 ; m=SQRT(ABS(f(5)))**3 ELSE r=ABS(f(3)) ; l=f(4) ; m=ABS(f(5))**3 ; f(5)=f(5)**2 ENDIF c composition chimique CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(f(5),mc(n_ch))),lc,comp,dcomp) IF(no_croiss)PRINT*,'Problème localisé en 8 dans lim_zc' c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w=wrot CASE(3,4,5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lc,frot,dfrot) IF(no_croiss)PRINT*,'Problème localisé en 9 dans lim_zc' w=frot(1) CASE(6) CALL omega_jpz(r,w) END SELECT c PRINT*,'dgrad1-1, j1',j1 d_grad(1)=dgrad(pt,p,t,dlpp,comp,m,l,r,dcomp,w) IF(dehors)RETURN c WRITE(*,2000)d_grad(1) c j2 en jlim(i) ---> dgrad(2) j2=jlim(i) CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,j2,lq,f,dfdq) IF(no_croiss)PRINT*,'Problème localisé en 10 dans lim_zc' c WRITE(*,2000)f pt=EXP(f(1)) dlpp=1.d0 IF(pturb)THEN !avec pression turbulente p=EXP(f(Ipg)) IF(der) dlpp=dfdq(Ipg)/dfdq(1) !dlpp=dln Pgaz/dln Ptot ELSE !sans pression turbulente p=pt ENDIF t=EXP(f(2)) IF(en_m23)THEN r=SQRT(ABS(f(3))) ; l=SQRT(ABS(f(4)))**3 m=SQRT(ABS(f(5)))**3 ELSE r=ABS(f(3)) ; l=f(4) ; m=ABS(f(5))**3 ; f(5)=f(5)**2 ENDIF CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(f(5),mc(n_ch))),lc,comp,dcomp) IF(no_croiss)PRINT*,'Problème localisé en 11 dans lim_zc' c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w=wrot CASE(3,4,5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lc,frot,dfrot) IF(no_croiss)PRINT*,'Problème localisé en 12 dans lim_zc' w=frot(1) CASE(6) CALL omega_jpz(r,w) END SELECT c PRINT*,'dgrad2-1, j2',j2 d_grad(2)=dgrad(pt,p,t,dlpp,comp,m,l,r,dcomp,w) IF(dehors)RETURN c WRITE(*,2000)d_grad(2) c pas d'inversion entre j1 et j2=jlim(i) --> j2 en jlim(i)-1 IF(d_grad(1)*d_grad(2) > 0.d0)THEN j2=jlim(i)-1 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,j2,lq,f,dfdq) IF(no_croiss)PRINT*,'Problème localisé en 13 dans lim_zc' c WRITE(*,2000)f pt=EXP(f(1)) ; dlpp=1.d0 IF(pturb)THEN !avec pression turbulente p=EXP(f(Ipg)) IF(der)dlpp=dfdq(Ipg)/dfdq(1) !dlpp=dln Pgaz/dln Ptot ELSE !sans pression turbulente p=pt ENDIF t=EXP(f(2)) IF(en_m23)THEN r=SQRT(ABS(f(3))) ; l=SQRT(ABS(f(4)))**3 m=SQRT(ABS(f(5)))**3 ELSE r=ABS(f(3)) ; l=f(4) ; m=ABS(f(5))**3 ; f(5)=f(5)**2 ENDIF CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(f(5),mc(n_ch))),lc,comp,dcomp) IF(no_croiss)PRINT*,'Problème localisé en 14 dans lim_zc' c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w=wrot CASE(3,4,5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lc,frot,dfrot) IF(no_croiss)PRINT*,'Problème localisé en 15 dans lim_zc' w=frot(1) CASE(6) CALL omega_jpz(r,w) END SELECT c PRINT*,'dgrad2-2,j2',j2 d_grad(2)=dgrad(pt,p,t,dlpp,comp,m,l,r,dcomp,w) IF(dehors)RETURN c WRITE(*,2000)d_grad(2) c pas d'inversion entre j1 et j2=jlim(i)-1 ---> j2 en jlim(i)+1 IF(d_grad(1)*d_grad(2) > 0.d0)THEN j2=jlim(i)+1 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,j2,lq,f,dfdq) IF(no_croiss)PRINT*,'Problème localisé en 16 dans lim_zc' c WRITE(*,2000)f pt=EXP(f(1)) ; dlpp=1.d0 IF(pturb)THEN !avec pression turbulente p=EXP(f(Ipg)) IF(der)dlpp=dfdq(Ipg)/dfdq(1) !dlpp=dln Pgaz/dln Ptot ELSE !sans pression turbulente p=pt ENDIF t=EXP(f(2)) IF(en_m23)THEN r=SQRT(ABS(f(3))) ; l=SQRT(ABS(f(4)))**3 m=SQRT(ABS(f(5)))**3 ELSE r=ABS(f(3)) ; l=f(4) ; m=ABS(f(5))**3 ; f(5)=f(5)**2 ENDIF CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(f(5),mc(n_ch))),lc,comp,dcomp) IF(no_croiss)PRINT*,'Problème localisé en 17 dans lim_zc' c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w=wrot CASE(3,4,5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lc,frot,dfrot) IF(no_croiss)PRINT*,'Problème localisé en 18 dans lim_zc' w=frot(1) CASE(6) CALL omega_jpz(r,w) END SELECT c PRINT*,'dgrad2-3,j2',j2 d_grad(2)=dgrad(pt,p,t,dlpp,comp,m,l,r,dcomp,w) IF(dehors)RETURN c WRITE(*,2000)d_grad(2) IF(d_grad(1)*d_grad(2) > 0.d0)THEN PRINT*,'lim_zc : erreur dans la dichotomie, limite ZR/ZC',i PRINT*,'j1, jlim(i)',j1,jlim(i) ; PRINT*,'on arrète' c PAUSE'on arrète' STOP ENDIF !jlim(i)+1 ENDIF !jlim(i)-1 ENDIF !on a encadré la limite entre j1 et j2 c limite encadrée entre j1 et j2, dichotomie avec j3 compt=0 j3=(j1+j2)/2.d0 B3: DO WHILE(ABS(j1-j2) > loc_zc) CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,j3,lq,f,dfdq) IF(no_croiss)PRINT*,'Problème localisé en 19 dans lim_zc' c WRITE(*,2000)f pt=EXP(f(1)) ; dlpp=1.d0 IF(pturb)THEN !avec pression turbulente p=EXP(f(Ipg)) IF(der) dlpp=dfdq(Ipg)/dfdq(1) !dlpp=dln Pgaz/dln Ptot ELSE !sans pression turbulente p=pt ENDIF t=EXP(f(2)) IF(en_m23)THEN r=SQRT(ABS(f(3))) ; l=SQRT(ABS(f(4)))**3 ;m=SQRT(ABS(f(5)))**3 ELSE r=ABS(f(3)) ; l=f(4) ; m=ABS(f(5))**3 ; f(5)=f(5)**2 ENDIF CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(f(5),mc(n_ch))),lc,comp,dcomp) IF(no_croiss)PRINT*,'Problème localisé en 20 dans lim_zc' c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w=wrot CASE(3,4,5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lc,frot,dfrot) IF(no_croiss)PRINT*,'Problème localisé en 21 dans lim_zc' w=frot(1) CASE(6) CALL omega_jpz(r,w) END SELECT c PRINT*,'dgrad3-1,j3',j3 d_grad(3)=dgrad(pt,p,t,dlpp,comp,m,l,r,dcomp,w) IF(dehors)RETURN c WRITE(*,2000)d_grad(3) c PRINT*,j1,j2,j3,compt ; WRITE(*,2000)d_grad(1:3) c limite entre j2 et j3 IF(d_grad(2)*d_grad(3) < 0.d0)THEN j1=j3 ; d_grad(1)=d_grad(3) c limite entre j1 et j3 ELSEIF(d_grad(1)*d_grad(3) <= 0.d0)THEN j2=j3 ; d_grad(2)=d_grad(3) ENDIF c PRINT*,'j1,j2,j3',j1,j2,j3,(j1+j2)/2.d0 j3=(j1+j2)/2.d0 ; compt=compt+1 IF(compt > dlt_n)THEN PRINT*,'dichotomie dans lim_zc, compt>',compt PRINT*,'on force j1,d_grad(1),j2,d_grad(2)',j1,d_grad(1),j2,d_grad(2) EXIT B3 ENDIF ENDDO B3 c on se place coté convectif pour faire plaisir à Maurice IF(d_grad(1) >= 0.d0)THEN j3=j1 !gradrad >= gradad: convectif ELSE j3=j2 !gradrad < gradad: radiatif ENDIF jlim(i)=NINT(j3) ; new(jlim(i))=j3 c PRINT*,'fin dic, i:',i,'jlim(i):',jlim(i),'new:',new(jlim(i)),new(jlim(i)+1) ENDDO B4 !i c---------écritures------------------------------- IF(jlim(1) <= idm6)THEN PRINT* SELECT CASE(langue) CASE('english') WRITE(usl_static,1017)idm6 1017 FORMAT('Limits RZ/CZ after the refinement of locations, i_lim < idm6=',i4) CASE DEFAULT WRITE(usl_static,17)idm6 17 FORMAT("Localisations affinées des ZR/ZC d'indice < ",i4) END SELECT ENDIF c masses et rayons aux limites Bzc: DO i=1,lim c test de localisation CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,new(jlim(i)),lq,f,dfdq) IF(no_croiss)PRINT*,'Problème localisé en 22 dans lim_zc' IF(en_m23)THEN r_zc(i)=SQRT(ABS(f(3))) ; m_zc(i)=SQRT(ABS(f(5)))**3 ELSE r_zc(i)=ABS(f(3)) ; m_zc(i)=ABS(f(5))**3 ENDIF c sauf pour les limites ZR/ZC trop externes IF(jlim(i) > idm6)EXIT Bzc c PRINT*,'n°,lq,jlim(i),new(jlim(i))',i,lq,qt(lq),jlim(i),new(jlim(i)) c PRINT*,'n_qs,ord_qs,knot,size',n_qs,ord_qs,knot,SIZE(bp,2) c PRINT*,f(5),bp(5,lq),bp(5,lq+1),bp(5,lq+2),bp(5,SIZE(bp,2)) IF(1.d0-m_zc(i)/mstar <= 0.d0)THEN c PRINT*,'n°,lq,jlim(i),new(jlim(i))',i,lq,qt(lq),jlim(i),new(jlim(i)) WRITE(usl_static,34)i,1.d0-m_zc(i)/mstar,m_zc(i),mstar 34 FORMAT('Erreur, lim_zc, limite n°', i3,', 1.d0-m_zc(i)/mstar=',es10.3, 1 ', m_zc(i)=',e15.8,/,'mstar=',2es15.8) CALL sortie('lim_zc 3') ENDIF c remplacement des localisations des limites par leurs valeurs affinées all_rap(i)=(new(jlim(i))-REAL(INT(new(jlim(i))),dp))*100.d0 c PRINT*,'i ',NINT(all_rap(i)),INT(new(jlim(i))) c PRINT*,all_rap(i),new(jlim(i)),REAL(INT(new(jlim(i))),dp) IF(lconv(i))THEN !début de ZC SELECT CASE(langue) CASE('english') WRITE(usl_static,1024)i,NINT(all_rap(i)),INT(new(jlim(i))), 2 INT(new(jlim(i)))+1,r_zc(i)/rn,m_zc(i),1.d0-m_zc(i)/mstar 1024 FORMAT('for the limit RZ/CZ #',i4,' at ',i3, 1 '% between shell',i5,' and',i5,', beginning of CZ',/, 2 'R_zc/Rtot=',es10.3,', M_zc/Msol=',es10.3,', 1-M_zc/Mstar=',es10.3) CASE DEFAULT WRITE(usl_static,24)i,NINT(all_rap(i)),INT(new(jlim(i))), 2 INT(new(jlim(i)))+1,r_zc(i)/rn,m_zc(i),1.d0-m_zc(i)/mstar 24 FORMAT('pour la limite affinée ZR/ZC n°',i4,' à ',i3, 1 '% entre couches',i5,' et',i5,', début de ZC',/, 2 'R_zc/Rtot=',es10.3,', M_zc/Msol=',es10.3,', 1-M_zc/Mstar=',es10.3) END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(usl_static,1025)i,NINT(all_rap(i)),INT(new(jlim(i))), 2 INT(new(jlim(i)))+1,r_zc(i)/rn,m_zc(i),1.d0-m_zc(i)/mstar 1025 FORMAT('for the limit RZ/CZ #',i4,' at ',i3, 1 '% between shell',i5,' and #',i5,', end of CZ',/, 2 'R_zc/Rtot=',es10.3,', M_zc/Msol=',es10.3,', 1-M_zc/Mstar=',es10.3) CASE DEFAULT WRITE(usl_static,25)i,NINT(all_rap(i)),INT(new(jlim(i))), 2 INT(new(jlim(i)))+1,r_zc(i)/rn,m_zc(i),1.d0-m_zc(i)/mstar 25 FORMAT('pour limite ZR/ZC n°',i4,' à ',i3,'% entre les couches n°', 1 i5,' et ',i5,', fin de ZC',/,'R_zc/Rtot=',es10.3, 2 ', M_zc/Msol=',es10.3,', 1-M_zc/Mstar=',es10.3) END SELECT ENDIF ENDDO Bzc !i c PRINT*,lconv(1:lim) ; PAUSE'lconv' c---------------------overshoot------------------------------------ c calcul (en Rsol) de r_zc et de r_ov si overshooting ou PC Iov: IF(ovsht)THEN !s'il y a overshooting ou PC Blim: DO i=1,lim c on repositionne r_zc aux new(jlim(i)) r_zc(i)=-100d0 ; r_ov(i)=-100d0 CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,new(jlim(i)), 1 lq,f,dfdq) c j3=REAL(jlim(i),dp) c CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,j3,lq,f,dfdq) IF(no_croiss)PRINT*,'Problème localisé en 24 dans lim_zc' pt=EXP(f(1)) ; dlpp=1.d0 IF(pturb)THEN !avec pression turbulente p=EXP(f(Ipg)) IF(der)dlpp=dfdq(Ipg)/dfdq(1) !dlpp=dln Pgaz/dln Ptot ELSE !sans pression turbulente p=pt ENDIF t=EXP(f(2)) IF(en_m23)THEN l=SQRT(ABS(f(4)))**3 ; r_zc(i)=SQRT(ABS(f(3))) m_zc(i)=SQRT(ABS(f(5)))**3 ELSE l=f(4) ; r_zc(i)=ABS(f(3)) ; m_zc(i)=ABS(f(5))**3 ; f(5)=f(5)**2 ENDIF CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch, 1 knotc,.TRUE.,MAX(mc(1),MIN(f(5),mc(n_ch))),lc,comp,dcomp) IF(no_croiss)PRINT*,'Problème localisé en 25 dans lim_zc' c PRINT*,'m_zc,r_zc,new(jlim(i))',lconv(i),i c WRITE(*,2000)m_zc(i),r_zc(i),new(jlim(i)) c WRITE(*,2000)p,t,r_zc(i),l,m_zc(i),comp(1) c vitesse angulaire SELECT CASE(Krot) CASE(0,1,2) w=wrot CASE(3,4,5) CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE., 1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lc,frot,dfrot) IF(no_croiss)PRINT*,'Problème localisé en 26 dans lim_zc' w=frot(1) CASE(6) CALL omega_jpz(r_zc(i),w) END SELECT CALL thermo(pt,p,t,m_zc(i),l,r_zc(i),dlpp,comp,dcomp, 1 ro,drop,drot,drox,u,dup,dut,dux, 2 grad,dgradpt,dgradp,dgradt,dgradx,dgradm,dgradl,dgradr, 3 dgradlpp,gam,dgampt,dgamp,dgamt,dgamx,dgamm,dgaml,dgamr, 4 dgamlpp,epsilon,depsp,depst,depsx,kap,dkapp,dkapt,dkapx, 5 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx, 6 gradad,dgradadp,dgradadt,dgradadx, 7 hp,dhppt,dhpp,dhpt,dhpx,dhpr,dhpm,grad_mu,mu,dmu_x,mue, 8 gradrad,alfa,beta,gamma1,radiatif) c overshoot: ovshti(s)*Hp ou ovshti*kip de JPZ si T > 5d5 à la limite It5: IF(t > 5.d5)THEN Ilc: IF(lconv(i))THEN !début de ZC, ext. vers le bas IF(ovshti > 0.d0)THEN IF(jpz)THEN kip=(3.d0-t*(drot/ro+dkapt/kap))*gradad-p*(drop/ro+dkapp/kap) r_ov(i)=MAX(r_zc(i)-hp*ovshti/rsol/kip,0.d0) c PRINT*,'i/r_zc(i),hp*ovshti/rsol/kip,hp,ovshti,rsol,kip',i c WRITE(*,2000)r_zc(i),hp*ovshti/rsol/kip,hp,ovshti,rsol,kip c WRITE(*,2000)r_ov(i) ELSE !ext. vers le bas: 0.r_ov>r_zc IF(ovshts > 0.d0)THEN !il y a ovsht sup. IF(jpz)THEN r_ov(i)=r_zc(i)*(1.d0+ovshts) !overshoot sup.= RcoreXovshts ELSE r_ov(i)=MIN(r_zc(i)+hp*ovshts/rsol,r_zc(i)*(1.d0+ovshts)) ENDIF !jpz r_ov(i)=MIN(r_ov(i),rn) c PRINT*,'ovshts' ; WRITE(*,2000)r_ov(i),r_zc(i),hp,rn ENDIF !ovshts ENDIF Ilc !lconv c PAUSE'après if sur lconv' c masse en Rsol d'overshooting ou de pénétration (r_ov en Rsol) IF(r_ov(i) >= 0.d0)THEN !r_ov=-100: no overshoot en r_zc(i) IF(en_m23)THEN CALL inter('r',bp,q,qt,n_qs,knot,r_ov(i)**2,f,dfdq,r_stat,m_stat) m_zc(i)=SQRT(ABS(f(5)))**3 ELSE CALL inter('r',bp,q,qt,n_qs,knot,r_ov(i),f,dfdq,r_stat,m_stat) m_zc(i)=ABS(f(5))**3 ENDIF c indice REAL en q de m_zc j_new=f(6) c WRITE(*,2000)f(1:ne) c l'ancien new(jlim(i)) est désormais un point ordinaire c qui redevient l'indice entier jlim(i) c la limite est désormais en jlim(i)= + proche entier de j_new c PRINT*,'jlim(i),new(jlim(i)),j_new',jlim(i),new(jlim(i)),j_new jlim(i)=NINT(j_new) new(jlim(i))=j_new !nouvel indice pour les limites ZR/ZC c PRINT*,'jlim(i),new(jlim(i)),j_new',jlim(i),new(jlim(i)),j_new c PRINT*,'j_new/M_ov,R_ov',j_new ; WRITE(*,2000)m_zc(i),r_ov(i) c PAUSE'après j_new' ENDIF !r_ov /= 0 ELSE It5 r_ov(i)=r_zc(i) c PAUSE't > 1.d5' ENDIF It5 !t > 1.d5 ENDDO Blim !i sur lim c il ne doit pas y avoir de chevauchement après les extensions c PRINT*,lim,ovsht,(jlim(i),i=1,lim),(lconv(i),i=1,lim) c WRITE(*,2000)(r_ov(j),j=1,lim) c WRITE(*,2000)(r_zc(j),j=1,lim) c WRITE(*,2000)(m_zc(j),j=1,lim) IF(lconv(1))THEN !début de ZC à la première limite i=2 ELSE i=1 ENDIF DO WHILE (i < lim) IF(jlim(i) >= jlim(i+1) .OR. (jlim(i+1)-jlim(i) <= n_min_ZC))THEN c IF(jlim(i) >= jlim(i+1))THEN DO j=i,lim-2 !chevauchement: suppression d'une ZR décalage jlim(j)=jlim(j+2) ; m_zc(j)=m_zc(j+2) ; lconv(j)=lconv(j+2) r_zc(j)=r_zc(j+2) ; r_ov(j)=r_ov(j+2) ENDDO lim=lim-2 ELSE i=i+2 ENDIF ENDDO !WHILE c PRINT*,lim,jlim(1:lim),lconv(1:lim) ; PAUSE'lim,jlim,lconv' c à cause d'un overshoot, le modèle devient totalement convectif IF(lim == 1 .AND. jlim(1) == 1 .AND. lconv(1))THEN jlim(1)=n_qs ; lconv(1)=.FALSE. ; r_zc(1)=rn tot_conv=.TRUE. ELSE tot_conv=.FALSE. c pas de limite ZR/ZC aux extrémités IF(jlim(1) == 1)THEN !en 1 DO i=1,lim-1 jlim(i)=jlim(i+1) ; lconv(i)=lconv(i+1) ; m_zc(i)=m_zc(i+1) r_zc(i)=r_zc(i+1) ; r_ov(i)=r_ov(i+1) ENDDO lim=lim-1 ENDIF IF(jlim(lim) == n_qs)lim=lim-1 !en n ENDIF PRINT* SELECT CASE(langue) CASE('english') WRITE(usl_static,1014) 1014 FORMAT('limits RZ/CZ after examination of overshoot :') CASE DEFAULT WRITE(usl_static,14) 14 FORMAT('limites ZR/ZC après examen de l''overshoot :') END SELECT IF(tot_conv)THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1033) 1033 FORMAT('Model fully mixed',/) CASE DEFAULT WRITE(usl_static,33) 33 FORMAT('Modèle complètement mélangé',/) END SELECT ELSE DO j=1,lim IF(lconv(j))THEN !début de ZC IF(ovshti > 0.d0)THEN IF(jpz)THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1031)j,jlim(j),r_zc(j)/rn,r_ov(j)/rn,m_zc(j)/mstar 1031 FORMAT('Convective penetration of JPZ limit CZ/RZ #',i4, 1 ', shell #',i5,/,3x,'radius CZ/Rtot=',es10.3, 2 ', reduced radius/Rtot=',es10.3,/, 3 'm/Mstar at the limit of the convective penetration=',es10.3,/) CASE DEFAULT WRITE(usl_static,31)j,jlim(j),r_zc(j)/rn,r_ov(j)/rn,m_zc(j)/mstar 31 FORMAT('Pénétration Convective JPZ limite ZC/ZR n°',i4, 1 ', couche n°',i5,/,3x,'rayon ZC/Rtot=',es10.3, 2 ', rayon reduit/Rtot=',es10.3,/, 3 'm/Mstar à la limite de la pénétration convective=',es10.3,/) END SELECT ELSE !jpz SELECT CASE(langue) CASE('english') WRITE(usl_static,1002)j,jlim(j),r_zc(j)/rn,r_ov(j)/rn,m_zc(j)/mstar 1002 FORMAT('negative overshooting limit CZ/RZ #',i4, 1 ', shell #',i5,/,'rayon ZC/Rtot=',es10.3, 2 ', reduced raius/Rtot=',es10.3,/, 3 'm/Mstar at the limit of overshooting=',es10.3,/) CASE DEFAULT WRITE(usl_static,2)j,jlim(j),r_zc(j)/rn,r_ov(j)/rn,m_zc(j)/mstar 2 FORMAT('Overshooting inférieur limite ZC/ZR n°',i4, 1 ', couche n°',i5,/,'rayon ZC/Rtot=',es10.3, 2 ', rayon reduit/Rtot=',es10.3,/, 3 'm/Mstar à la limite de l''overshooting=',es10.3,/) END SELECT ENDIF!jpz ELSE !sans overshoot SELECT CASE(langue) CASE('english') WRITE(usl_static,1011)j,jlim(j),r_zc(j)/rn,m_zc(j)/mstar 1011 FORMAT('Limit without overshooting nor PC, limit RZ/CZ #:', 1 i4,', shell #',i5,/,'radius CZ/Rtot=',es10.3, 2 ', m/Mstar at the end of the CZ=',es10.3,/) CASE DEFAULT WRITE(usl_static,11)j,jlim(j),r_zc(j)/rn,m_zc(j)/mstar 11 FORMAT('Limite sans overshooting ni PC, limite ZR/ZC n°:', 1 i4,', couche n°',i5,/,'rayon ZC/Rtot=',es10.3, 2 ', m/Mstar à la fin de la ZC=',es10.3,/) END SELECT ENDIF !ovshti ELSE !fin de ZC IF(ovshts > 0.d0)THEN IF(jpz)THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1012)j,jlim(j),r_zc(j)/rn,r_ov(j)/rn,m_zc(j)/mstar 1012 FORMAT('positive overshoot * R core, limit CZ/RZ #',i4, 1 ', shell #',i5,/,3x,'radius CZ/Rtot=',es10.3, 2 ', extented radius/Rtot=',es10.3,/, 3 'm/Mstar at the limit of overshoot=',es10.3,/) CASE DEFAULT WRITE(usl_static,12)j,jlim(j),r_zc(j)/rn,r_ov(j)/rn,m_zc(j)/mstar 12 FORMAT('Overshoot supérieur * R noyau, limite ZC/ZR n°',i4, 1 ', couche n°',i5,/,3x,'rayon ZC/Rtot=',es10.3, 2 ', rayon etendu/Rtot=',es10.3,/, 3 'm/Mstar à la limite de l''overshoot=',es10.3,/) END SELECT ELSE !jpz IF(r_ov(j) == r_zc(j))THEN SELECT CASE(langue) CASE('english') WRITE(usl_static,1013)j,jlim(j),r_zc(j)/rn,m_zc(j)/mstar 1013 FORMAT('Limit CZ/RZ without positive overshoot #',i4, 1 ', shell #',i5,/,'because T < 5.d5, radius CZ/Rtot=', 2 es10.3,', m/Mstar on the limit =',es10.3,/) CASE DEFAULT WRITE(usl_static,13)j,jlim(j),r_zc(j)/rn,m_zc(j)/mstar 13 FORMAT('Limite ZC/ZR sans overshooting supérieur n°',i4, 1 ', couche n°',i5,/,'car T < 5.d5, rayon ZC/Rtot=',es10.3, 2 ', m/Mstar à la limite =',es10.3,/) END SELECT ELSE SELECT CASE(langue) CASE('english') WRITE(usl_static,1001)j,jlim(j),r_zc(j)/rn,r_ov(j)/rn,m_zc(j)/mstar 1001 FORMAT('Overshooting positive * Hp, limit CZ/RZ #',i4, 1 ', shell #',i5,/,'radius CZ/Rtot=',es10.3, 2 ', extended radius/Rtot=',es10.3,/, 3 'm/Mstar at the limit of overshooting=', 1 es10.3,/) CASE DEFAULT WRITE(usl_static,1)j,jlim(j),r_zc(j)/rn,r_ov(j)/rn,m_zc(j)/mstar 1 FORMAT('Overshooting supérieur * Hp, limite ZC/ZR n°',i4, 1 ', couche n°',i5,/,'rayon ZC/Rtot=',es10.3, 2 ', rayon etendu/Rtot=',es10.3,/, 3 'm/Mstar à la limite de l''overshooting=',es10.3,/) END SELECT ENDIF ENDIF!jpz ELSE !pas d'overshoot SELECT CASE(langue) CASE('english') WRITE(usl_static,1021)j,jlim(j),r_zc(j)/rn,m_zc(j)/mstar 1021 FORMAT('Limit without overshooting nor PC, limit RZ/CZ #', 1 i4,', shell #',i5,/,'radius CZ/Rtot=',es10.3, 2 ', m/Mstar at the beginning of the CZ=',es10.3,/) CASE DEFAULT WRITE(usl_static,21)j,jlim(j),r_zc(j)/rn,m_zc(j)/mstar 21 FORMAT('Limite sans overshooting ni PC, limite ZR/ZC n°',i4, 1 ', couche n°',i5,/,'rayon ZC/Rtot=',es10.3, 2 ', m/Mstar au début de la ZC=',es10.3,/) END SELECT ENDIF !ovshts ENDIF !lconv ENDDO !j limites ZR/ZC ENDIF !tot_conv ENDIF Iov !ovsht c PAUSE'après les overshoot' c-------------------------zone externe---------------------- c pour le mélange la partie supérieure de l'enveloppe doit être conv. si c tel n'est pas le cas on décide abitrairement que les dernières n_min_ZC c couches le sont c PRINT*,'lim,jlim(1:lim),n_qs',lim,jlim(1:lim),n_qs IF(tot_rad)THEN lim=1 ; lconv(1)=.TRUE. jlim(1)=n_qs-n_min_ZC c ajout d'une limite arbitraire ELSEIF(.NOT.tot_conv .AND. .NOT.lconv(lim))THEN IF(jlim(lim) < n_qs-2*n_min_ZC)THEN lim=lim+1 ; lconv(lim)=.TRUE. jlim(lim)=n_qs-n_min_ZC ENDIF c retrait ou suppression de la dernière limite ELSEIF(.NOT.tot_conv .AND. jlim(lim) > n_qs-n_min_ZC)THEN IF(jlim(lim-1) > n_qs-2*n_min_ZC)THEN lim=lim-1 ELSE jlim(lim)=n_qs-n_min_ZC ENDIF ENDIF c positionnements, en new(jlim(i)), de r_zc et m_zc en Rsol, Msol, m_zc23 en Msol^2/3 DO i=1,lim CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,new(jlim(i)),lq,f,dfdq) IF(no_croiss)PRINT*,'Pb. at 27 in lim_zc' IF(en_m23)THEN r_zc(i)=SQRT(ABS(f(3))) ; m_zc(i)=SQRT(ABS(f(5)))**3 ; m_zc23(i)=f(5) ELSE r_zc(i)=ABS(f(3)) ; m_zc(i)=ABS(f(5))**3 ; m_zc23(i)=f(5)**2 ENDIF ENDDO c PRINT,'lim_zc mzc23',m_zc23(1:lim),m_zc(1:lim) c facteurs d'espacement fac /= 1 au voisinage des limites ZR/ZC c l'utilisation de ces positions relatives des points de grille et des limites c s'est avérée moins efficace que l'utilisation de la fonction de répartition c limite ZR/ZC entre i et i+1 fac*(esp(i+1)-esp(i)) = psi Ilf: IF(l_fac)THEN c WRITE(*,36)idm6 36 FORMAT('déplacements de bords de couche vers les limites ZR/ZC i < ',i4) c PRINT*,'l_fac',l_fac ; PAUSE'l_fac' c les facteurs d'espacement fac ALLOCATE(esp(n_qs)) DO i=1,n_qs esp(i)=espmt(new(i)) ENDDO DO i=1,n_qs-1 fac(i)=bp(6,1)/(esp(i+1)-esp(i)) ENDDO DEALLOCATE(esp) ENDIF Ilf c avec un coeur convectif on prend l=Rcoeur-r, ce qui sera c réalisé en prenant r_zc_conv(id_conv)=-r_zc(1) IF(.NOT.(tot_conv .OR. tot_rad))THEN IF(lconv(1))THEN id_conv=1 ; r_zc_conv(id_conv)=r_zc(1) !coeur non convectif ELSE id_conv=0 ; r_zc_conv(id_conv)=-r_zc(1) !coeur convectif ENDIF DO i=1,lim r_zc_conv(i+id_conv)=r_zc(i) ENDDO IF(lconv(lim))THEN if_conv=lim+1 !enveloppe totalement convective r_zc_conv(if_conv)=MAX(rstar,r_zc(lim)) ELSE if_conv=lim !enveloppe non totalement convective ENDIF ENDIF ENDIF Ilim !lim<=0 c new est encore alloué DEALLOCATE(new) c cas sans limite ZR/ZC, sortie IF(tot_conv .OR. tot_rad)THEN id_conv=0 ; r_zc_conv(id_conv)=0.d0 if_conv=1 ; r_zc_conv(if_conv)=rstar l_dmax=.TRUE. ; dist_max=0.d0 SELECT CASE(langue) CASE('english') WRITE(usl_static,1026) CASE DEFAULT WRITE(usl_static,26) END SELECT RETURN ENDIF c valeur maximale des localisations des limites ZR/ZC < idm6 DO i=1,lim IF(jlim(i) > idm6)all_rap(i)=0.d0 ENDDO WHERE(all_rap >= 50.d0)all_rap=100.d0-all_rap idist_max=MAXLOC(all_rap) ; i=idist_max(1) dist_max=NINT(MAXVAL(all_rap)) c PRINT*,'all_rap:',all_rap(1:lim),i,dist_max c l_dmax les limites ZR/ZC sont à mieux que d_zc% du bord des couches IF(l_fac .AND. d_zc < 50)THEN l_dmax=dist_max <= d_zc IF(l_dmax)THEN WRITE(usl_static,48) 48 FORMAT('un point de grille a été approché de chaque lim. ZR/ZC affinée') ELSE WRITE(usl_static,49)d_zc 49 FORMAT("approche d'un point de grille à mieux que d_zc=", 1 i2,'% de chaque limite ZR/ZC affinée') ENDIF ENDIF c écritures SELECT CASE(langue) CASE('english') IF(l_fac .AND. d_zc < 50)THEN WRITE(usl_static,1022)i,dist_max,jlim(i) 1022 FORMAT(/,'-->--> the farest limit ZR/ZC ',i2,' is located at ',i3, 1 '% from the shell limits',i5,'<--<--<--') ELSE WRITE(usl_static,1027) 1027 FORMAT('without bringing closer limits ZR/ZC and grid points') l_dmax=.TRUE. ENDIF WRITE(usl_static,1026) 1026 FORMAT(/,'---------Search of the limits ZR/ZC(end)-----------') CASE DEFAULT II2: IF(l_fac .AND. d_zc < 50)THEN II1: IF(l_dmax)THEN WRITE(usl_static,32)d_zc 32 FORMAT('les limites ZR/ZC affinées sont à moins de ',i2, 1 '% des points de grille') ELSE II1 II0: IF(rap(i) <= 0.5d0)THEN j=MIN(n_qs,jlim(i)+1) WRITE(usl_static,22)i,jlim(i),NINT(q(j)),dist_max 22 FORMAT('la limite ZR/ZC n°',i2,", de l'intervalle [",i4,',',i4, 1 "] reste éloignée de ",i2,"% d'un bord ") ELSE II0 j=MAX(1,jlim(i))-1 WRITE(usl_static,22)i,NINT(q(j)),jlim(i),dist_max c23 FORMAT('la limite ZR/ZC n°',i2,", dans l'intervalle [",i4,',',i4, c 1 "] reste éloignée d'un des bords de ",i2,'%') ENDIF II0 ENDIF II1 ELSE II2 WRITE(usl_static,27) l_dmax=.TRUE. ENDIF II2 IF(no_discon)THEN WRITE(6,42) ELSE WRITE(6,47) ENDIF WRITE(usl_static,26) 26 FORMAT(/,'---------Recherche des limites ZR/ZC(fin)-----------') END SELECT c CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,1.211d0,lq,f,dfdq) c WRITE(*,2000)f c PAUSE'sortie lim_zc' RETURN END SUBROUTINE lim_zc