c************************************************************************ PROGRAM calib2k_pms c calibration zsx pour le soleil avec une éventuelle variation de la c composition en X et Y dans la ZC externe au cours de l'évolution c en tenant compte de la perte de masse et de la rotation c Auteur: P.Morel, Département J.D. Cassini, O.C.A. c CESAM2k c le Z/X est fixé dans ce programme: zsx=0.0245 par exemple c cf. A.Noel & N.Grevesse, "Proceedings of the first SOHO c workshop, Anapolis, Maryland, USA, c 25-28 August 1992 (ESA SP-348, november 1992)". c la calibration donne X et Y d'ou Z=1-X-Y qu'il faut ajuster à zsx c glob: variables globales c glob(1)=mstar*msol c glob(2)=rtot*rsol c glob(3)=ltot*lsol c glob(4)=z0 c glob(5)=x0 c glob(6)=alpha c glob(7)=9.d0/4.d0 paramètre de convection arbitraire c glob(8)=1./162. paramètre de convection arbitraire c glob(9)=X dans ZC c glob(10)=Y dans ZC c glob(11)=d2p c glob(12)=d2ro c glob(13)=age c glob(14)=wrot vitesse de rotation globale c glob(15)=w_rot vitesse de rotation initiale c var: variables c nvar=22 pour oscillations adiabatiques c var(1,i)=r*rsol c var(2,i)=log(m/mstar) -1.d38 au centre c var(3,i)=t c var(4,i)=Ptot c var(5,i)=ro c var(6,i)=gradient reel d ln T / d ln P c var(7,i)=l c var(8,i)=kap c var(9,i)=énergie thermo+gravifique c var(10,i)=grand Gamma1 c var(11,i)=gradient adiabatique c var(12,i)=delta c var(13,i)=cp c var(14,i)=mu elec. c var(15,i)=vaissala, 0 au centre c var(16,i)=vitesse angulaire, radian/sec c var(17,i)=d ln kappa / d ln T c var(18,i)=d ln kappa / d ln ro c var(19,i)=d epsilon(nuc) / d ln T c var(20,i)=d epsilon(nuc) / d ln ro c var(21,i)=Ptot / Pgas c var(22,i)=gradient radiatif c nvar=25 pour inversion c var(23,i)=d Gamma1 / d log P c var(24,i)=d Gamma1 / d log T c var(25,i)=d Gamma1 / dY = d Gamma1 / dZ c nvar=44 pour oscillations non adiabatiques c var(26,i)=dP / dro (TX) c var(27,i)=dP / dT (roX) c var(28,i)=dP / dX (Tro) c var(29,i)=du / dro (TX) c var(30,i)=du / dT (roX) c var(31,i)=du / dx(Tro) c var(32,i)=énergie interne c var(33,i)=d^2P / dro^2 (TX) c var(34,i)=d^2P / dro dT (X) c var(35,i)=d^2P / dT^2(roX) c var(36,i)=d^2U / dro^2 (TX) c var(37,i)=d^2U / dro dT (X) c var(38,i)=d^2U / dT^2 (X) c var(39,i)=dK / dX c var(40,i)=d^2K / dT^2 c var(41,i)=d epsi / dX c var(42,i)=dX / dR c var(43,i)=J-B c var(44,i)=Edding. facteur c vecteur de composition chimique c var(nvar+j,i)=xchim(j)*nucleo(j), j=1,nchim c------------------------------------------------------------------ USE mod_donnees, ONLY : alpha, garde_xish, Krot, lsol, msol, mtot, 1 nchim, nom_elem, nom_fich2, nom_output, pi, rsol, unit, 2 w_rot, x0, y0, zsx0, zsx_sol, z0, nom_pertm USE mod_exploit, ONLY : lit_nl, write_nl USE mod_kind USE mod_numerique, ONLY : gauss_band USE mod_variables, ONLY : wrot IMPLICIT NONE REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: var REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: glob INTEGER, PARAMETER :: ppara=4 !nombre de paramètres à fixer REAL (kind=dp), DIMENSION(ppara,ppara) :: a REAL (kind=dp), DIMENSION(1,ppara) :: b REAL (kind=dp), DIMENSION(0:ppara) :: al, beta, betaf, lf, rf, 1 x, xf, y, yf, w, wf REAL (kind=dp), DIMENSION(ppara) :: d REAL (kind=dp), PARAMETER :: cal=1.d-4, p_sol=28.d0 REAL (kind=dp) :: cal_max, l, m, r, w_sol, x_zc, y_zc, zsx_zc, z_zc INTEGER, DIMENSION(ppara) :: indpc=1 INTEGER :: i, itot, j, nglob, nvar, un LOGICAL :: err=.FALSE., inversible, ok CHARACTER (len=80), DIMENSION(4) :: abid CHARACTER (len=80) :: chaine c------------------------------------------------------------------- 2000 FORMAT(8es10.3) WRITE(*,20)cal 20 FORMAT(/,'********************************',/, 1 'calibration en Z/X au niveau',es10.3,' d''un modèle solaire',/, 2 'Entrer le nom générique du modèle à calibrer, Ex: soleil') READ(5,'(a)')nom_fich2 WRITE(*,1)TRIM(nom_fich2) 1 FORMAT('Calibration du modèle: ',a) c lecture du fichier de données CALL lit_nl(wrot) c pour une calibration on conserve les rapports métal/Z et non métal/X c de la mixture garde_xish=.FALSE. c période en radian/s w_sol=tr_period(p_sol,.TRUE.) WRITE(*,12)rsol,lsol,zsx_sol,p_sol,w_sol 12 FORMAT('Rsol, Lsol, Z/Xsol, M/Msol et Psol à obtenir',/, 1 'Rsol=',es12.5,', Lsol=',es12.5,', Z/Xsol=',es12.5,/, 2 'Msol=1, Psol=',es9.2,' jours =',es12.5,' rad/sec',/) IF(nom_output == 'osc_adia')THEN chaine=TRIM(nom_fich2)//'-ad.osc' ELSEIF(nom_output == 'osc_nadia')THEN chaine=TRIM(nom_fich2)//'-nad.osc' ELSEIF(nom_output == 'osc_invers')THEN chaine=TRIM(nom_fich2)//'-inv.osc' ELSE WRITE(*,3)TRIM(chaine) ; WRITE(*,4) ; STOP 4 FORMAT('Pour une calibration nom_output doit être', 1 'osc_adia ou osc_nadia ou osc_invers , Arrêt') ENDIF INQUIRE(file=TRIM(chaine),exist=ok) IF(ok)THEN OPEN(unit=30,form='formatted',status='old',file=TRIM(chaine)) WRITE(*,5)TRIM(chaine) 5 FORMAT('fichier ASCII utilisé: ',a) ELSE WRITE(*,3)TRIM(chaine) ; STOP 3 FORMAT('Arrêt, fichier ASCII non trouvé: ',a) ENDIF c on lit les global et la première couche si rotation READ(30,'(a)')abid(1) !m 100X707a173_$.dat du 18/ 8/1993 etc... c PRINT*,abid(1) READ(30,'(a)')abid(2) !fichier pour calcul des oscillations etc.. c PRINT*,abid(2) READ(30,'(a)')abid(3) !methode: CESAM 3.0.0.0 colloc. etc... c PRINT*,abid(3) READ(30,'(a)')abid(4) c PRINT*,abid(4) READ(30,89)nchim ! 9 H He3 He4 etc... 89 FORMAT(i3) ALLOCATE(nom_elem(nchim)) BACKSPACE (unit=30) READ(30,90)nchim,nom_elem(1:nchim) 90 FORMAT(i3,14(1x,a4)) READ(30,137)itot,nglob,nvar,nchim,Krot 137 FORMAT(5i10) c PRINT*,itot,nglob,nvar,nchim,Krot c CALL pause('itot,nglob,nvar,nchim,Krot') ALLOCATE(glob(nglob)) READ(30,138)glob 138 FORMAT(5es19.12) m=glob(1) ; r=glob(2) ; l=glob(3) ; zsx0=z0/x0 x_zc=glob(7) ; y_zc=glob(8) ; z_zc=1.d0-x_zc-y_zc zsx_zc=z_zc/x_zc ; wrot=glob(12) IF(Krot > 0)THEN ALLOCATE(var(nvar+nchim,1)) READ(30,138)var(:,1) w_rot=var(16,1) ENDIF CLOSE(unit=30) WRITE(*,10)x0,y0,z0,zsx0,alpha,mtot,tr_period(wrot,.FALSE.),wrot 10 FORMAT(/,'modèle initial:',/, 1 'X0=',es12.5,', Y0=',es12.5,', Z0=',es12.5,/, 2 'Z0/X0=',es12.5,', alpha=',es12.5,', Mtot=',es12.5,/, 3 'période=',es12.5,' jours =',es12.5' rad/sec',/) WRITE(*,11)r,l,zsx_zc,m/msol,x_zc,y_zc,z_zc, 1 sign(tr_period(w_rot,.FALSE.),wrot),sign(w_rot,wrot) 11 FORMAT('modèle évolué:',/, 1 'R=',es12.5,', L=',es12.5,', Z/X=',es12.5,', Mstar=',es10.3,/, 3 'Xzc=',es12.5,', Yzc=',es12.5,', Zzc=',es12.5,/, 4 'période=',es10.3,' jours =',es12.5,' rad/sec',/) c Initialisation du jacobien de la calibration c valeurs avant évolution pour le calcul des coeff. de calibration al(0)=1.922d0 ; al(1)=1.85d0 ; al(2)=1.922d0 ; al(3)=1.922d0 ; al(4)=1.922d0 x(0)=0.6973d0 ; x(1)=0.6973d0 ; x(2)=0.72d0 ; x(3)=0.6973d0 ; x(4)=0.6973d0 y(0)=0.28285d0 ; y(1)=0.28285d0 ; y(2)=0.26236d0 ; y(3)=0.286d0 ; y(4)=0.28285d0 w(0)=6.0d-9 ; w(1)=6.0d-9 ; w(2)=6.0d-9 ; w(3)=6.0d-9 ; w(4)=6.2d-9 c valeurs après évolution pour le calcul des coeff. de calibration rf(0)=9.9961d-01 !Rfinal/Rsol rf(1)=1.0061d+00 ; rf(2)=9.7205d-01 ; rf(3)=1.0319d+00 ; rf(4)=9.9961d-01 lf(0)=-3.8550d-04 !log10(Lfinal/Lsol) lf(1)=-1.0608d-03 ; lf(2)=-4.6988E-02 ; lf(3)=6.1537d-02 ; lf(4)=-3.8589d-04 xf(0)=7.283d-01 !Xfinal xf(1)=7.289d-01 ; xf(2)=7.470d-01 ; xf(3)=7.329d-01 ; xf(4)=7.283d-01 yf(0)=2.536d-01 !Yfinal yf(1)=2.531d-01 ; yf(2)=2.365d-01 ; yf(3)=2.517d-01 ; yf(4)=2.536d-01 wf(0)=2.222d-06 !Wfinal wf(1)=2.260d-6 ; wf(2)=2.260d-06 ; wf(3)=2.070d-06 ; wf(4)=2.296d-06 DO i=0,ppara lf(i)=10.d0**lf(i) !Lfinal/Lsol beta(i)=(1.d0-x(i)-y(i))/x(i) !beta =Z/X betaf(i)=(1.d0-xf(i)-yf(i))/xf(i) ENDDO c PRINT*,'beta/beta*' c WRITE(*,2000)(beta(i),i=0,3) c WRITE(*,2000)(betaf(i),i=0,3) c WRITE(*,2000) a=0.d0 DO i=1,ppara a(i,i)=1.d0 ENDDO a(1,1)=(rf(1)-rf(0))/(al(1)-al(0)) !dR / d alpha a(1,2)=(rf(2)-rf(0))/(x(2)-x(0)) !dR / d X a(1,3)=(rf(3)-rf(0))/(beta(3)-beta(0)) !dR / d beta IF(w_rot /= 0.d0)a(1,4)=(rf(4)-rf(0))/(w(4)-w(0)) !dR / d w a(2,1)=(lf(1)-lf(0))/(al(1)-al(0)) !dL / d alpha a(2,2)=(lf(2)-lf(0))/(x(2)-x(0)) !dL / d X a(2,3)=(lf(3)-lf(0))/(beta(3)-beta(0)) !dL / d beta IF(w_rot /= 0.d0)a(2,4)=(lf(4)-lf(0))/(w(4)-w(0)) !dL / d w a(3,1)=(betaf(1)-betaf(0))/(al(1)-al(0)) !d beta* / d alpha a(3,2)=(betaf(2)-betaf(0))/(x(2)-x(0)) !d beta* / d X a(3,3)=(betaf(3)-betaf(0))/(beta(3)-beta(0)) !d beta* / d beta IF(w_rot /= 0.d0)a(4,3)=(betaf(4)-betaf(0))/(w(4)-w(0)) !d beta / d w IF(w_rot /= 0.d0)THEN a(4,1)=(wf(1)-wf(0))/(al(1)-al(0)) !d w* / d alpha a(4,2)=(wf(2)-wf(0))/(x(2)-x(0)) !d w* / d X a(4,3)=(wf(3)-wf(0))/(beta(3)-beta(0)) !d w* / d beta a(4,4)=(wf(4)-wf(0))/(w(4)-w(0)) !d w* / d w ENDIF c DO i=1,ppara c WRITE(*,2000)a(i,:) c ENDDO c PRINT* b(1,1)=r/rsol-1.d0 ; b(1,2)=l/lsol-1.d0 ; b(1,3)=zsx_zc-zsx_sol IF(w_rot /= 0.d0)THEN b(1,4)=ABS(w_rot)-w_sol ELSE b(1,4)=0.d0 ENDIF d=ABS(b(1,:)) d(3)=d(3)/zsx_sol ; d(4)=d(4)/tr_period(p_sol,.TRUE.) PRINT*,'calibration en valeurs absolues' WRITE(*,2)b 2 FORMAT('d_r=',es9.2,', d_l=',es9.2,', d_Zsx=',es9.2,', d_w=',es9.2) PRINT*,'calibration en valeurs relatives' WRITE(*,2)d d(4)=d(4)*1.d-2 cal_max=MAXVAL(d) IF(cal_max < cal)THEN WRITE(*,13)cal_max 13 FORMAT(/,'C BON! modèle calibré à mieux que:',es10.3) un=0 ELSE WRITE(*,14)cal_max 14 FORMAT(/,'C pas encore BON! modèle calibré à:',es10.3,/) CALL gauss_band(a,b,indpc,ppara,ppara,ppara,1,inversible) IF(.NOT.inversible)THEN PRINT*,'jacobien non inversible, Arrêt' ; STOP ENDIF WRITE(*,21)b(1,:) 21 FORMAT('corrections :',5es10.3) B1: DO IF(ABS(b(1,1)) > 0.1D0*alpha)THEN b=b/2.d0 ; CYCLE B1 ELSEIF(ABS(b(1,2)) > 0.1D0*x0)THEN b=b/2.d0 ; CYCLE B1 ELSEIF(ABS(b(1,3)) > 0.1D0*zsx0)THEN b=b/2.d0 ; CYCLE B1 ELSEIF(ABS(b(1,4)) > 0.1D0*ABS(wrot) .AND. wrot /= 0.d0)THEN b=b/2.d0 ; CYCLE B1 ENDIF EXIT B1 ENDDO B1 WRITE(*,22)b(1,:) 22 FORMAT('corrections effectives:',5es10.3) alpha=alpha-b(1,1) x0=x0-b(1,2) zsx0=zsx0-b(1,3) z0=x0*zsx0 y0=1.d0-x0-z0 IF(wrot /= 0.d0)wrot=ABS(wrot)-b(1,4) PRINT* WRITE(*,6)x0,y0,z0,zsx0,alpha,mtot,tr_period(wrot,.FALSE.),wrot 6 FORMAT('nouvelles valeurs pour calibration',/, 1 'X0=',es12.5,', Y0=',es12.5,', Z0=',es12.5,/, 2 'Z0/X0=',es12.5,', alpha=',es12.5,', Mtot=',es12.5,/, 3 'période=',es12.5,' jours =',es12.5', rad/sec',/) IF(alpha < 0.d0)THEN PRINT*,'Erreur l/Hp négatif' ; err=.TRUE. ELSEIF(x0 < 0.d0)THEN PRINT*,'Erreur X < 0' ; err=.TRUE. ELSEIF(x0 > 1.d0)THEN PRINT*,'Erreur X > 1' ; err=.TRUE. ELSEIF(y0 < 0.d0)THEN PRINT*,'Erreur Y < 0' ; err=.TRUE. ELSEIF(y0 > 1.d0)THEN PRINT*,'Erreur Y > 1' ELSEIF(z0 < 0.d0)THEN PRINT*,'Erreur Z < 0' ; err=.TRUE. ELSEIF(z0 > 1.d0)THEN PRINT*,'Erreur Z > 1' ; err=.TRUE. ELSEIF(wrot < 0.d0)THEN PRINT*,'Erreur sur W < 0' ; err=.TRUE. ENDIF w_rot=sign(wrot,w_rot) ; unit='rad/s' CALL write_nl un=1 ENDIF IF(err)THEN x0=0.d0 ; x0=1.d0/x0 ; STOP ENDIF c on dispose 0 ou 1 dans mon_modele.tempo suivant que la calibration c est ou non acquise chaine= TRIM(nom_fich2)//'.tempo' OPEN(unit=25,form='formatted',status='unknown',file=chaine) WRITE(25,*)un ; CLOSE(unit=25) STOP CONTAINS c********************************************************************** REAL (kind=dp) FUNCTION tr_period(p,sens) c transformation de la période jours <==> période rad/sec c sens=true: période jours ==> période rad/sec c sens=false: période rad/sec ==> période jours c Auteur: P.Morel, Département J.D. Cassini, O.C.A. c CESAM2k c-------------------------------------------------------------- USE mod_kind IMPLICIT NONE REAL (kind=dp), INTENT(in) :: p LOGICAL, INTENT(in) :: sens c------------------------------------------------------------- IF(ABS(p) > 0.d0)THEN IF(sens)THEN !période jours ==> période rad/sec tr_period=sign(2.d0*pi/24.d0/3600.d0/ABS(p),p) ELSE !période rad/sec ==> période jours tr_period=sign(2.d0*pi/24.d0/3600.d0/ABS(p),p) ENDIF ELSE tr_period=0.d0 ENDIF RETURN END FUNCTION tr_period END PROGRAM calib2k_pms