c**************************************************************** c subroutine gr3interpatmas59cgm(teff0,logg0,tauTe,Ttaue,tauPe,Ptaue, c 1 taulaf,rholaf,teffaf) SUBROUTINE piau3(teff0,logg0,tauTe,Ttaue,tauPe,Ptaue, 1 taulaf,rholaf,teffaf) c lois T(tau) et P(tau) interpolees sur les relations atmospheriques Atlas12 c pour les objets de population I avec theorie de la convection CGM/Bernkopf c alpha=0.5. c Composition solaire As05 ou As09. c Auteur : Laurent Piau, SAp-Saclay c entree : c teff : temperature effective c grav : gravite c sortie : c P(tau) : relation P de tau a la valeur la plus proche c de Teff et de logg c ou c T(tau) : relation T de tau a la valeur la plus proche c de Teff et de logg c Et valeurs de tau & rho a la limite atmospherique c Possibilite d'interpoler lineairement en temperature sur les tableaux c de donnees aux quatre plus proches voisins : cf consignes encadrement & decommenter c ICI VERSION (HAL/GR)INTERPATM3 (analogue a (HAL/GR)INTERPATM2): INTERPOLATION AUX QUATRES TABLEAUX VOISINS c MAIS DE PLUS RELATION RECALIBREE AVEC LA LOI EMPIRIQUE INSPIREE DE VANDENBERG 2008 & MUELLER 1974 c Variables intermediaires : c matatm : matrice atmospherique : matrice de vecteurs atmospheriques c vecatm : vecteur atmospherique : structure contenant Teff, log(g) c et les lois P(tau) & T(tau) c Les valeurs Teff et log(g) sont aux c points de grille NextGen. c Domaine couvert Teff / logg c c logg 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 c Teff c c 4000 + + + + + + + + + + + c | + + + + + + + + + + + c | + + + + + + + + + + + c 6400 + + + + + + + + + + + USE mod_donnees, only : nom_chemin USE mod_numerique, only : linf USE mod_variables, ONLY : sortie implicit none integer, parameter :: itauP = 72 ! Nombre de point P/tau NG integer, parameter :: itauT = 72 ! Nombre de point T/tau NG type atmstruct real (kind=dp) :: teff,logg real (kind=dp) :: taula,rhola,teffa real (kind=dp) :: tauP(itauP) real (kind=dp) :: Ptau(itauP) real (kind=dp) :: tauT(itauT) real (kind=dp) :: Ttau(itauT) end type atmstruct c type (atmstruct) atmstruct1 type matatmty type (atmstruct), dimension(25,11) :: vecatm end type matatmty type (matatmty), SAVE :: matatm real (kind=dp) :: teff0,logg0,taulaf,rholaf,teffaf real (kind=dp) :: yy1(25) real (kind=dp) :: yy2(11) real (kind=dp) :: tauPe(itauP) real (kind=dp) :: Ptaue(itauP) real (kind=dp) :: tauTe(itauT) real (kind=dp) :: Ttaue(1,itauT) real (kind=dp) :: teffl1,loggl2,teffl1p1,loggl2p1, 1 taulaint1,taulaint2,rholaint1,rholaint2 real (kind=dp) :: tauPint1(itauP) real (kind=dp) :: tauPint2(itauP) real (kind=dp) :: Ptauint1(itauP) real (kind=dp) :: Ptauint2(itauP) real (kind=dp) :: tauTint1(itauT) real (kind=dp) :: tauTint2(itauT) real (kind=dp) :: Ttauint1(itauT) real (kind=dp) :: Ttauint2(itauT) integer,SAVE :: i,im,jm,l1=1,l2=1 logical :: ok CHARACTER (len=100) :: fich_piau c--------------------------------------------------------------------------- c vérification de l'existence des fichiers de données fich_piau=TRIM(nom_chemin)//'PIAU/AT12PopIAs05Ku08pres.dat' INQUIRE(file=TRIM(fich_piau),exist=ok) IF(.NOT.ok)THEN PRINT*,'sortie, fichier inconnu: ',TRIM(fich_piau) CALL sortie('piau3 1') ENDIF fich_piau=TRIM(nom_chemin)//'PIAU/AT12PopIAs05Ku08temp.dat' INQUIRE(file=TRIM(fich_piau),exist=ok) IF(.NOT.ok)THEN PRINT*,'sortie, fichier inconnu: ',TRIM(fich_piau) CALL sortie('piau3 2') ENDIF fich_piau=TRIM(nom_chemin)//'PIAU/AT12PopIAs05Ku08rhotau.dat' INQUIRE(file=TRIM(fich_piau),exist=ok) IF(.NOT.ok)THEN PRINT*,'sortie, fichier inconnu: ',TRIM(fich_piau) CALL sortie('piau3 3') ENDIF open(unit=41,form='formatted',status='old', 1 file=TRIM(nom_chemin)//'PIAU/AT12PopIAs05Ku08pres.dat') open(unit=51,form='formatted',status='old', 1 file=TRIM(nom_chemin)//'PIAU/AT12PopIAs05Ku08temp.dat') open(unit=71,form='formatted',status='old', 1 file=TRIM(nom_chemin)//'PIAU/AT12PopIAs05Ku08rhotau.dat') c teff0=6000. c logg0=4.00 teff0=teff0/100. c RELATION P(tau) c Lecture des donnees atmospheriques modeles de geantes Teff 4000 -> 6400 K pas de 100 K c et logg = 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 do jm = 1,11 do im = 1,25 read(41,100)matatm%vecatm(im,jm)%teff, 1 matatm%vecatm(im,jm)%logg c write(6,101)matatm%vecatm(im,jm)%teff, c 1 matatm%vecatm(im,jm)%logg do i=1,itauP read(41,200)matatm%vecatm(im,jm)%tauP(i), 1 matatm%vecatm(im,jm)%Ptau(i) enddo c do i=1,itauP c write(6,200)matatm%vecatm(im,jm)%tauP(i), c 1 matatm%vecatm(im,jm)%Ptau(i) c enddo read(41,*) c do i=1,itauT c read(41,200)matatm%vecatm(im,jm)%tauT(i), c 1 matatm%vecatm(im,jm)%Ttau(i) c enddo enddo ! im enddo ! jm c pause c RELATION T(tau) c Lecture des donnees atmospheriques modeles sur les geantes Teff 4000 -> 6400 K pas de 100 K c et logg = 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 : http:// do jm = 1,11 do im = 1,25 read(51,100)matatm%vecatm(im,jm)%teff, 1 matatm%vecatm(im,jm)%logg read(71,102)matatm%vecatm(im,jm)%teff, 1 matatm%vecatm(im,jm)%logg, 1 matatm%vecatm(im,jm)%taula, 1 matatm%vecatm(im,jm)%rhola c write(6,102)matatm%vecatm(im,jm)%teff, c 1 matatm%vecatm(im,jm)%logg, c 1 matatm%vecatm(im,jm)%taula, c 1 matatm%vecatm(im,jm)%rhola do i=1,itauT read(51,200)matatm%vecatm(im,jm)%tauT(i), 1 matatm%vecatm(im,jm)%Ttau(i) c write(6,200)matatm%vecatm(im,jm)%tauT(i), c 1 matatm%vecatm(im,jm)%Ttau(i) enddo read(51,*) enddo ! im enddo ! jm c pause 100 format(f4.0,f5.2) 101 format(f4.0,1x,f4.1) 102 format(f4.0,1x,f5.2,1x,1pd12.5,1x,1pd12.5) 200 format(1x,1pd12.5,1x,1pd12.5) close(unit=41) close(unit=51) close(unit=71) c Choix des structures de Teff et logg les plus proches des donnees do im = 1,25 yy1(im)=matatm%vecatm(im,1)%teff enddo do jm = 1,11 yy2(jm)=matatm%vecatm(1,jm)%logg enddo call linf(teff0,yy1,25,l1) call linf(logg0,yy2,11,l2) c write(6,*)'Teff encadrants',matatm%vecatm(l1,l2)%teff, c 1 teff0*100.,matatm%vecatm(l1+1,l2)%teff c write(6,*)'logg encadrants',matatm%vecatm(l1,l2)%logg, c 1 logg0,matatm%vecatm(l1,l2+1)%logg c Si structure atmospherique au point de grille le plus proche : decommenter c if(abs(teff0-matatm%vecatm(l1,l2)%teff) .gt. c 1 abs(teff0-matatm%vecatm(l1+1,l2)%teff)) l1=l1+1 c if(abs(logg0-matatm%vecatm(l1,l2)%logg) .gt. c 1 abs(logg0-matatm%vecatm(l1,l2+1)%logg)) l2=l2+1 c Si structure atmospherique au point de grille le plus proche : decommenter if((teff0 .ge. 64).or.(teff0 .le. 40) .or. (logg0 .le. 0.5) 1 .or. (logg0 .ge. 5.5)) then write(6,*) 1 'Extrapolation des tables de structure atmospherique NextGen' write(2,*) 1 'Extrapolation des tables de structure atmospherique NextGen' tauPe=matatm%vecatm(l1,l2)%tauP Ptaue=matatm%vecatm(l1,l2)%Ptau tauTe=matatm%vecatm(l1,l2)%tauT Ttaue(1,:)=matatm%vecatm(l1,l2)%Ttau taulaf=matatm%vecatm(l1,l2)%taula rholaf=matatm%vecatm(l1,l2)%rhola teffaf=matatm%vecatm(l1,l2)%teff*100. c write(6,*)'Teff & Teff approchante',teff0*100., c 1 matatm%vecatm(l1,l2)%teff*100,l1 c write(6,*)'logg & logg approchant',logg0, c 1 matatm%vecatm(l1,l2)%logg,l2 c write(6,*)matatm%vecatm(l1,l2)%tauP c write(6,*)matatm%vecatm(l1,l2)%tauT c write(6,*)'Tau & Rho associes', c 1 matatm%vecatm(l1,l2)%taula, c 1 matatm%vecatm(l1,l2)%rhola c pause goto 1000 endif write(6,*)'Teff & Teff approchante',teff0*100., 1 matatm%vecatm(l1,l2)%teff*100,l1 write(6,*)'logg & logg approchant',logg0, 1 matatm%vecatm(l1,l2)%logg,l2 c write(6,*)matatm%vecatm(l1,l2)%tauP c write(6,*)matatm%vecatm(l1,l2)%tauT c write(6,*)'Tau & Rho associes', c 1 matatm%vecatm(l1,l2)%taula, c 2 matatm%vecatm(l1,l2)%rhola c pause c Consigne encadrement : soit : structure atmospherique au point de grille le plus proche c soit : interpolations aux quatre points de grille encadrants c Structure atmospherique au point de grille le plus proche : decommenter c tauPe=matatm%vecatm(l1,l2)%tauP c Ptaue=matatm%vecatm(l1,l2)%Ptau c tauTe=matatm%vecatm(l1,l2)%tauT c Ttaue(1,:)=matatm%vecatm(l1,l2)%Ttau c taulaf=matatm%vecatm(l1,l2)%taula c rholaf=matatm%vecatm(l1,l2)%rhola c teffaf=matatm%vecatm(l1,l2)%teff*100. c Structure atmospherique au point de grille le plus proche : decommenter c Si interpolation aux quatre points de grille encadrants : decommenter teffl1=matatm%vecatm(l1,l2)%teff teffl1p1=matatm%vecatm(l1+1,l2)%teff loggl2=matatm%vecatm(l1,l2)%logg loggl2p1=matatm%vecatm(l1,l2+1)%logg tauPint1=matatm%vecatm(l1,l2)%tauP 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2)%tauP-matatm%vecatm(l1,l2)%tauP) tauPint2=matatm%vecatm(l1,l2+1)%tauP 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2+1)%tauP-matatm%vecatm(l1,l2+1)%tauP) tauPe=tauPint1+((logg0-loggl2)/(loggl2p1-loggl2)) 1 *(tauPint2-tauPint1) Ptauint1=matatm%vecatm(l1,l2)%Ptau 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2)%Ptau-matatm%vecatm(l1,l2)%Ptau) Ptauint2=matatm%vecatm(l1,l2+1)%Ptau 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2+1)%Ptau-matatm%vecatm(l1,l2+1)%Ptau) Ptaue=Ptauint1+((logg0-loggl2)/(loggl2p1-loggl2)) 1 *(Ptauint2-Ptauint1) tauTint1=matatm%vecatm(l1,l2)%tauT 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2)%tauT-matatm%vecatm(l1,l2)%tauT) tauTint2=matatm%vecatm(l1,l2+1)%tauT 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2+1)%tauT-matatm%vecatm(l1,l2+1)%tauT) tauTe=tauTint1+((logg0-loggl2)/(loggl2p1-loggl2)) 1 *(tauTint2-tauTint1) Ttauint1=matatm%vecatm(l1,l2)%Ttau 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2)%Ttau-matatm%vecatm(l1,l2)%Ttau) Ttauint2=matatm%vecatm(l1,l2+1)%Ttau 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2+1)%Ttau-matatm%vecatm(l1,l2+1)%Ttau) Ttaue(1,:)=Ttauint1+((logg0-loggl2)/(loggl2p1-loggl2)) 1 *(Ttauint2-Ttauint1) taulaint1=matatm%vecatm(l1,l2)%taula 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2)%taula-matatm%vecatm(l1,l2)%taula) taulaint2=matatm%vecatm(l1,l2+1)%taula 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2+1)%taula-matatm%vecatm(l1,l2+1)%taula) taulaf=taulaint1+((logg0-loggl2)/(loggl2p1-loggl2)) 1 *(taulaint2-taulaint1) rholaint1=matatm%vecatm(l1,l2)%rhola 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2)%rhola-matatm%vecatm(l1,l2)%rhola) rholaint2=matatm%vecatm(l1,l2+1)%rhola 1 +((teff0-teffl1)/(teffl1p1-teffl1)) 2 *(matatm%vecatm(l1+1,l2+1)%rhola-matatm%vecatm(l1,l2+1)%rhola) rholaf=rholaint1+((logg0-loggl2)/(loggl2p1-loggl2)) 1 *(rholaint2-rholaint1) teffaf=teff0*100. c Si interpolation aux quatre points de grille encadrants 1000 continue teff0=teff0*100. c do i=1,itauT c write(6,*)i c write(6,*)'Avant offset HM74',teff0,tauTe(i),Ttaue(1,i) c call hm74interpatmas59cgm(teff0,tauTe(i),Ttaue(1,i),thm74,tkusol) c write(6,*)'Apres offset HM74',teff0,tauTe(i),Ttaue(1,i) c pause c enddo return c end subroutine gr3interpatmas59cgm END SUBROuTINE piau3