c--------------------------------------------------------------------------- subroutine ppxi2_bp(mtot,nchim) c routine du module mod_bp_for_alecian_new ! Version 2.2a (16/02/2007) ! Cette subroutine effectue la lecture des options, des donnees et ! prepare les tableaux necessaires au calcul des ! accelerations (grad2_bp). ! Auteur: ! Georges ALECIAN ! DAEC - UMR 8631, Observatoire de Meudon ! F-92195 MEUDON CEDEX, FRANCE ! Tel: 01 45 07 74 20, + 33 1 45 07 74 20 c Adaptation à CESAM2k B.Pichon, P.Morel c------------------------------------------------------------------------- USE mod_donnees, ONLY : langue, nom_frad USE mod_kind Implicit None Real(DP), Intent(In) :: mtot Integer, Intent(In) :: nchim !** declarations Morel mod. Pichon integer :: i,j,k !** declarations Alecian ! mes variables locales integer, parameter :: nsvp=300,nlist=100 integer :: iunit_svp,izn,ios integer :: icompte,ic,icf,numion integer, dimension(nsvp) :: np,iz,numa integer, dimension(3,nsvp) :: n_cross integer, dimension(nlist) :: mass_svp,numions integer, dimension(nsvp,nlist) :: nps,izs character(1) :: b character(2), dimension(nsvp) :: el character(4), dimension(nsvp,nlist) :: listis character(4), dimension(nsvp) :: listi character(64), dimension(nlist) :: listf real, dimension(nsvp,nlist) :: phis,psis,xis,alphs,acs,bets real, dimension(nsvp) :: phi,psi,xi,alph,ac,bet,absol real(DP) :: mass_st_sol real(DP), dimension(nlist) :: X,F1,F2,F3,F4,F5,F6 ! igrad (si =0 le g_rad des elements non-identifies sera zero) ! (si =1 le g_rad des elements non-identifies sera -gravite) ! tgrad (si =0 le g_rad de tous les elements identifies sera zero) ! (si =1 le g_rad de tous les elements identifies sera -gravite) ! (si =2 le g_rad de tous les elements identifies sera calcule) ! np: nb de protons du noyau ! iz: degre d'ionization (0 pour neutre) ! izn: (non utilisee) ! listi: nom de l'ion (C1 = carbone neutre) ! numion: nb d'ions dans la base !************************** c write(10,*) 'From_Alecian ==> version du 21/03/2007; 15h40' ! unites de lecture iunit_svp=71 ! initialisation des tableaux id_ppxi = 0 isotops = 0 phistar = 0. psistar = 0. xistar = 0. alphstar = -0.5 acstar = 1. betstar = 0. c g_rad pour éléments non identifiés IF(INDEX(nom_frad,'1') /= 0)THEN igrad=1 ELSE igrad=0 ENDIF c écritures SELECT CASE(langue) CASE('english') SELECT CASE(igrad) CASE(0) WRITE(*,1001) ; WRITE(2,1001) 1001 FORMAT('Radiative accelerations according to G.Alécian',/, 1 'The g_rad of unidentified elements will be zero (igrad=0)') CASE(1) WRITE(*,1002) ; WRITE(2,1002) 1002 FORMAT('Radiative accelerations according to G.Alécian',/, 1 'The g_rad of unidentified elements will be -gravity (igrad=1)') END SELECT SELECT CASE(tgrad) CASE(0) WRITE(*,1003) ; WRITE(2,1003) 1003 FORMAT('g_rad of unidentified elements will be zero (tgrad=0)',/) CASE(1) WRITE(*,1004) ; WRITE(2,1004) 1004 FORMAT('g_rad of unidentified elements will be -gravity (tgrad=1)',/) CASE(2) WRITE(*,1005) ; WRITE(2,1005) 1005 FORMAT('g_rad of unidentified elements will be computed (tgrad=2)',/) END SELECT CASE DEFAULT SELECT CASE(igrad) CASE(0) WRITE(*,1) ; WRITE(2,1) 1 FORMAT('Accélérations radiatives selon G.Alécian',/, 1 'Le g_rad des éléments non-identifiés sera zéro (igrad=0)') CASE(1) WRITE(*,2) ; WRITE(2,2) 2 FORMAT('Accélérations radiatives selon G.Alécian',/, 1 'g_rad des éléments non-identifiés sera -gravité (igrad=1)') END SELECT SELECT CASE(tgrad) CASE(0) WRITE(*,3) ; WRITE(2,3) 3 FORMAT('g_rad des éléments identifiés sera zéro (tgrad=0)',/) CASE(1) WRITE(*,4) ; WRITE(2,4) 4 FORMAT('g_rad des éléments identifiés sera -gravité (tgrad=1)',/) CASE(2) WRITE(*,5) ; WRITE(2,5) 5 FORMAT('g_rad des éléments identifiés sera calculé (tgrad=2)',/) END SELECT END SELECT !*** ! Tableau des isotopes de CESAM, matrice (isotops) du genre: ! H1 He3 He4 C12 C13 N14 N15 O16 O17 Fe56 ! H1 1 0 0 0 0 0 0 0 0 0 ! He3 0 1 1 0 0 0 0 0 0 0 ! He4 0 1 1 0 0 0 0 0 0 0 ! C12 0 0 0 1 1 0 0 0 0 0 ! C13 0 0 0 1 1 0 0 0 0 0 ! N14 0 0 0 0 0 1 1 0 0 0 ! N15 0 0 0 0 0 1 1 0 0 0 ! O16 0 0 0 0 0 0 0 1 1 0 ! O17 0 0 0 0 0 0 0 1 1 0 !Fe56 0 0 0 0 0 0 0 0 0 1 j=1 do while (j.le.nchim) icompte=0 do k=1,nchim if(k.lt.j) then isotops(j,k)=0 else if(k.eq.j) then isotops(j,k)=1 else if(zi(j).eq.zi(k)) then isotops(j,k)=1 icompte=icompte+1 else isotops(j,k)=0 end if end if end do ic=1 do while (ic.le.icompte) j=j+1 do k=1,nchim isotops(j,k)=isotops(j-1,k) end do ic=ic+1 end do j=j+1 end do ! lecture de la liste des fichiers contenant les tables phi, psi, etc. OPEN (UNIT = iunit_svp, + FILE = TRIM(nomch)//'datai_gr2p1/liste_svp.data', + STATUS = 'old', + IOStat=ios + ) If ( ios /= 0 ) STOP " Pb- liste_svp.data" icf=0 read(iunit_svp,'(a1)') b read(iunit_svp,'(a1)') b do i=1,100 read(iunit_svp,'(a)',end=90) listf(i) icf=icf+1 end do 90 continue close(iunit_svp) ! determination des masses disponibles dans la base do i=1,icf read(listf(i),'(8x,i3)',iostat=ios) mass_svp(i) if(ios/=0) then print*,'erreur identification de mass_svp => stop' stop end if end do mass_st_sol=NINT(mtot*100.) if((NINT(mass_st_sol) < mass_svp(1)).or. + (NINT(mass_st_sol) > mass_svp(icf))) then print* print*,'Attention mtot hors limites dans ppxi2_bp:' print*,'mtot*100/msol=',mass_st_sol print*,'min*100=',mass_svp(1) print*,'max*100=',mass_svp(icf) print* end if ! lecture de la base des phi et psi do i=1,icf OPEN (UNIT = iunit_svp, + FILE = TRIM(nomch)//'datai_gr2p1/'//listf(i), + STATUS = 'old', + IOStat=ios + ) If ( ios /= 0 ) STOP " Pb- svp_*.data" read(iunit_svp,'(a1)') b numions(i)=0 do j=1,1000 read(iunit_svp,'(3i4,1x,a4,6e12.3)',end=100) + nps(j,i),izs(j,i),izn,listis(j,i), + phis(j,i),psis(j,i),xis(j,i), + alphs(j,i),acs(j,i),bets(j,i) numions(i)=numions(i)+1 end do 100 continue close(unit=iunit_svp) end do 10 format(i4,a1,i4,a1,e10.3,a1,e10.3,a1,i4,a1,a4) ! control des fichiers lus do i=2,icf if(numions(i-1).ne.numions(i)) then print*, 'Anomalie des numions => STOP' stop end if do j=1,numions(i) if(listis(j,i-1).ne.listis(j,i)) then print*, 'Anomalie des listis => STOP' stop end if end do end do ! preparation numion=numions(1) n_cross= 0 !tableau(3,nsvp) np(1:numion) = nps(1:numion,1) ! nb protons du noyau iz(1:numion) = izs(1:numion,1) ! charge ion listi(1:numion) = listis(1:numion,1) ! interpolations do i=1,icf X(i)=float(mass_svp(i)) end do do j=1,numion F1(1:icf) = phis(j,1:icf) F2(1:icf) = psis(j,1:icf) F3(1:icf) = xis(j,1:icf) F4(1:icf) = alphs(j,1:icf) F5(1:icf) = acs(j,1:icf) F6(1:icf) = bets(j,1:icf) phi(j) = FT(mass_st_sol,100,X,F1) psi(j) = FT(mass_st_sol,100,X,F2) xi(j) = FT(mass_st_sol,100,X,F3) alph(j) = FT(mass_st_sol,100,X,F4) ac(j) = FT(mass_st_sol,100,X,F5) bet(j) = FT(mass_st_sol,100,X,F6) end do ! Cross-identification ppxi/cesam. On pose id_ppxi(j)=1 si au moins un ion OK ! On suppose que les proprietes atomiques (pour les g_rad) sont ! les memes entre isotopes d'un meme element. do j=1,nchim do k=1,numion if(zi(j).eq.np(k)) then n_cross(1,k) = iz(k) !correspondance svp/cesam n_cross(2,k) = j !correspondance svp/cesam n_cross(3,k) = n_cross(3,k) +1 ! nb.isotopes id_ppxi(j) = 1 phistar(iz(k),j) = phi(k) psistar(iz(k),j) = psi(k) xistar(iz(k),j) = xi(k) alphstar(iz(k),j) = alph(k) acstar(iz(k),j) = ac(k) betstar(iz(k),j) = bet(k) end if end do end do ! lecture des niveaux atomiques pour le calcul du g_cont (photoionisation) call niv(listi,numion,n_cross,np,iunit_svp) ! lecture de la liste des abondances solaires. OPEN (UNIT = iunit_svp, + FILE = TRIM(nomch)//'datai_gr2p1/solabnum.dat', + STATUS = 'old', + IOStat=ios + ) If ( ios /= 0 ) STOP " Pb- solabnum.dat" ic=0 do k=1,nsvp read(iunit_svp,*,end=101) el(k),numa(k),absol(k) ic=ic+1 end do 101 continue close (iunit_svp) ! concentration par rapport a H pour le calcul de khi dans grad2_bp c print*, 'ic= ',ic do j=1,nchim do i=1,ic if(zi(j).eq.numa(i)) then el_svp(j) = el(i) C_sol(j) = absol(i) exit end if end do end do do j=2,nchim C_sol(j)=C_sol(j)/C_sol(1) end do C_sol(1)=1. end subroutine ppxi2_bp