c**************************************************************** c subroutine gr1interpatmas59cgm(teff0,logg0,tauTe,Ttaue,tauPe,Ptaue, c 1 taulaf,rholaf,teffaf) SUBROUTINE piau1(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 la theorie CGM & Bernkopf: 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)INTERPATM1 : TABLEAUX LES PLUS PROCHES 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_kind 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 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) integer, SAVE :: i,im,jm,l1=1,l2=1 logical, SAVE :: init=.true. logical :: ok CHARACTER (len=100) :: fich_piau c------------------------------------------------------------------------- teff0=teff0/100. if(init)then init=.false. 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('piau1 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('piau1 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('piau1 3') ENDIF c PRINT*,ok ; PAUSE'ici' 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 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 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(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,101)matatm%vecatm(im,jm)%teff, c 1 matatm%vecatm(im,jm)%logg 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) endif 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 if(abs(teff0-matatm%vecatm(l1,l2)%teff) .gt. 1 abs(teff0-matatm%vecatm(l1+1,l2)%teff)) l1=l1+1 if(abs(logg0-matatm%vecatm(l1,l2)%logg) .gt. 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' 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 do i=1,itauT c write(6,200)matatm%vecatm(l1,l2)%tauT(i) c enddo c write(6,*)'Tau & Rho associes', c 1 matatm%vecatm(l1,l2)%taula, c 1 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 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 Structure atmospherique au point de grille le plus proche : decommenter c Si interpolation aux quatre points de grille encadrants : decommenter C teffl1=matatm%vecatm(l1,l2)%teff C teffl1p1=matatm%vecatm(l1+1,l2)%teff C loggl2=matatm%vecatm(l1,l2)%logg C loggl2p1=matatm%vecatm(l1,l2+1)%logg c tauPint1=matatm%vecatm(l1,l2)%tauP c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2)%tauP-matatm%vecatm(l1,l2)%tauP) c tauPint2=matatm%vecatm(l1,l2+1)%tauP c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2+1)%tauP-matatm%vecatm(l1,l2+1)%tauP) c tauPe=tauPint1+((logg0-loggl2)/(loggl2p1-loggl2)) c 1 *(tauPint2-tauPint1) c Ptauint1=matatm%vecatm(l1,l2)%Ptau c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2)%Ptau-matatm%vecatm(l1,l2)%Ptau) c Ptauint2=matatm%vecatm(l1,l2+1)%Ptau c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2+1)%Ptau-matatm%vecatm(l1,l2+1)%Ptau) c Ptaue=Ptauint1+((logg0-loggl2)/(loggl2p1-loggl2)) c 1 *(Ptauint2-Ptauint1) c tauTint1=matatm%vecatm(l1,l2)%tauT c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2)%tauT-matatm%vecatm(l1,l2)%tauT) c tauTint2=matatm%vecatm(l1,l2+1)%tauT c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2+1)%tauT-matatm%vecatm(l1,l2+1)%tauT) c tauTe=tauTint1+((logg0-loggl2)/(loggl2p1-loggl2)) c 1 *(tauTint2-tauTint1) c Ttauint1=matatm%vecatm(l1,l2)%Ttau c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2)%Ttau-matatm%vecatm(l1,l2)%Ttau) c Ttauint2=matatm%vecatm(l1,l2+1)%Ttau c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2+1)%Ttau-matatm%vecatm(l1,l2+1)%Ttau) c Ttaue(1,:)=Ttauint1+((logg0-loggl2)/(loggl2p1-loggl2)) c 1 *(Ttauint2-Ttauint1) c taulaint1=matatm%vecatm(l1,l2)%taula c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2)%taula-matatm%vecatm(l1,l2)%taula) c taulaint2=matatm%vecatm(l1,l2+1)%taula c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2+1)%taula-matatm%vecatm(l1,l2+1)%taula) c taulaf=taulaint1+((logg0-loggl2)/(loggl2p1-loggl2)) c 1 *(taulaint2-taulaint1) c rholaint1=matatm%vecatm(l1,l2)%rhola c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2)%rhola-matatm%vecatm(l1,l2)%rhola) c rholaint2=matatm%vecatm(l1,l2+1)%rhola c 1 +((teff0-teffl1)/(teffl1p1-teffl1)) c 2 *(matatm%vecatm(l1+1,l2+1)%rhola-matatm%vecatm(l1,l2+1)%rhola) c rholaf=rholaint1+((logg0-loggl2)/(loggl2p1-loggl2)) c 1 *(rholaint2-rholaint1) c teffaf=teff0*100. c Si interpolation aux quatre points de grille encadrants teff0=teff0*100. return c end subroutine gr1interpatmas59cgm END SUBROUTINE piau1