program read_table implicit none integer :: alpha ! Number of bins integer :: Nvarchimie ! Number of species integer :: nvar ! Equal to Nvarchimie integer :: nion=9 ! Number of non-grain species integer,parameter::dp=kind(1.0d0) ! default !!!! Grains info real(dp), allocatable, dimension(:) :: x_g ! Abundance of grains per bins real(dp), allocatable, dimension(:) :: r_g ! Radius of grains per bins real(dp), allocatable, dimension(:) :: m_g ! Mass of grains per bins !!!! Constants real(dp), parameter :: pi=3.1415927410125732422_dp real(dp), parameter :: mp=1.6726d-24 ! Proton mass in g real(dp), parameter :: me=9.1094d-28 ! Electron mass in g real(dp), parameter :: e=4.803204d-10 ! Electron charge in cgs real(dp), parameter :: Kb = 1.3807d-16 ! Boltzmann constant erg/K real(dp), allocatable, dimension(:) :: q ! Electric charge real(dp), allocatable, dimension(:) :: m ! Mass !!!! bins of grains real(dp), parameter :: rho_s=2.3_dp ! Bulk density g.cc real(dp), parameter :: a_0=0.0375d-4 ! Reference radius cm real(dp), parameter :: a_min=0.0181d-4 ! Minimum radius cm real(dp), parameter :: a_max=0.9049d-4 ! Maximum radius cm real(dp), parameter :: zeta=a_min/a_max ! a_min/a_max real(dp), parameter :: lambda_pow=-3.5d0 ! Coeff power law real(dp) :: rho_gtot ! Grain density real(dp) :: Lp1,Lp3,Lp4 ! resistivites (cf Kunz & Mouschovias 2009, Marchand et al. 2016) real(dp), allocatable, dimension(:) :: sigma real(dp), allocatable, dimension(:) :: zetas real(dp), allocatable, dimension(:) :: phi real(dp), allocatable, dimension(:) :: tau_sn,tau_gp,tau_gm real(dp), allocatable, dimension(:,:) :: tau_inel real(dp), allocatable, dimension(:) :: gamma_zeta,gamma_omega real(dp), allocatable, dimension(:) :: omega,omega_bar real(dp), parameter :: c_l=299792458.d2 ! c in cm/s real(dp) :: sigv, muuu integer :: iB,iH,iT,iX real(dp) :: B,bmaxchimie,nH,T,xi,sigH,sigO,sigP integer :: cmp_X ! For loops integer :: tchimie ! tchmie steps in temperature real(dp) :: dtchimie ! dtchmie step in temperature real(dp) :: tminchimie ! min tchmie integer :: nchimie ! nchmie steps in density real(dp) :: dnchimie ! dnchmie step in density real(dp) :: nminchimie ! min nchmie integer :: bchimie ! bchmie steps in B field real(dp) :: dbchimie ! dbchmie step in B field real(dp) :: bminchimie ! min bchmie integer :: xichimie ! Steps in ionisation rate real(dp) :: dxichimie ! dxichimie step in ion rate real(dp) :: ximinchimie ! min xichimie integer :: nislin,tislin,xiislin ! Linear scale or not real(dp),allocatable,dimension(:,:,:,:,:)::conductivities ! resistivities real(dp),allocatable,dimension(:,:,:,:)::abundances ! abundances integer :: i,j,k,ij real(dp) :: sig_p, sig_perp, sig_h, eta_ohm, eta_hall, eta_ad open(42,file='Table_abundances.dat', status='old') read(42,*) nchimie, tchimie, xichimie, nvar ! Number of steps in density, temperature and ionisation rate read(42,*) nislin, tislin, xiislin ! 1: Constant step in logscale, 0: Not constant step in logscale read(42,*) read(42,*) read(42,*) allocate(abundances(-2:nvar,nchimie,tchimie,xichimie)) ! abundances(-2,:,:,:)=density ! abundances(-1,:,:,:)=temperature ! abundances(0,:,:,:)=ionisation rate ! abundances(n>°,:,:,:)=relative abundance of species n do ij=1,xichimie do i=1,tchimie do j=1,nchimie read(42,*)abundances(-2:nvar,j,i,ij) end do read(42,*) read(42,*) end do end do close(42) nminchimie=(abundances(-2,1,1,1)) !min density tminchimie=(abundances(-1,1,1,1)) !min temperature ximinchimie=(abundances(0,1,1,1)) !min ionisation rate ! Density, temperature and IR steps in log (if nislin, tislin, xiislin = 1) dnchimie=(log10(abundances(-2,nchimie,1,1))-log10(abundances(-2,1,1,1)))/(nchimie-1) dtchimie=(log10(abundances(-1,1,tchimie,1))-log10(abundances(-1,1,1,1)))/(tchimie-1) dxichimie=(log10(abundances(0,1,1,xichimie))-log10(abundances(0,1,1,1)))/(xichimie-1) !Allocating arrays alpha=(nvar-nion)/3 allocate(x_g(alpha)) allocate(r_g(alpha)) allocate(m_g(alpha)) allocate(q(nvar)) allocate(m(nvar)) allocate(sigma(nvar)) allocate(zetas(nvar)) allocate(phi(nvar)) allocate(tau_sn(nvar)) allocate(tau_gp(nvar)) allocate(tau_gm(nvar)) allocate(gamma_zeta(nvar)) allocate(gamma_omega(nvar)) allocate(omega(nvar)) allocate(omega_bar(nvar)) allocate(tau_inel(nvar,nvar)) ! Computation of grains radii and species mass Lp1=dble(lambda_pow+1) Lp3=dble(lambda_pow+3) Lp4=dble(lambda_pow+4) ! Grain radius if (alpha==1) then r_g(1)=a_0 else do i=1,alpha ! cf Marchand et al. 2016 r_g(i)=a_min*zeta**(-dble(i)/dble(alpha)) * & & dsqrt( Lp1/Lp3* (1d0-zeta**(Lp3/dble(alpha)))/(1d0-zeta**(Lp1/dble(alpha)))) end do end if q(:)=1.d0*e ! cations q(1)=-1.d0*e ! electrons do i=nion+1,nvar if (mod(i-nion,3)==1) q(i)=1.d0*e ! g+ if (mod(i-nion,3)==2) q(i)=-1.d0*e ! g- if (mod(i-nion,3)==0) q(i)=0.d0 ! g0 end do m(:)=0.d0 m(1) = me ! e- m(2) = 23.5d0*mp ! ions metalliques m(3) = 29.d0*mp ! ions moleculaires m(4) = 3*mp ! H3+ m(5) = mp ! H+ m(6) = 12.d0*mp ! C+ m(7) = 4.d0*mp ! He+ m(8) = 39.098*mp ! K+ m(9) = 22.99d0*mp ! Na+ do i=1,alpha ! masse des grains m_g(i)=4.d0/3.d0*pi*r_g(i)**3*rho_s m(nion+3*(i-1)+1:nion+3*i)=m_g(i) end do ! values for Btable bminchimie=1.d-10 bmaxchimie=1.d10 ! ok for first core in nimhd. maybe not enough for second core. bchimie=150 dbchimie=(log10(bmaxchimie)-log10(bminchimie))/real((bchimie-1),dp) allocate(conductivities(0:3,1:nchimie,1:tchimie,1:xichimie,1:bchimie)) !conductivities(1,:,:,:)= sigma_// !conductivities(2,:,:,:)= sigma_perp !conductivities(3,:,:,:)= sigma_hall !conductivities(0,:,:,:)= sigma_hall sign ! Computation of resistivities tau_sn = 0.0_dp omega = 0.0_dp sigma = 0.0_dp phi = 0.0_dp zetas = 0.0_dp gamma_zeta = 0.0_dp gamma_omega = 0.0_dp omega_bar = 0.0_dp do iB=1,bchimie do iX=1,xichimie do iT=1,tchimie do iH=1,nchimie nh=abundances(-2,iH,iT,iX) ! density (.cc) of current point B =10.0d0**(log10(bminchimie)+dble(iB-1)*dbchimie) ! Magnetic field T =abundances(-1,iH,iT,iX) ! Temperature xi =abundances(0,iH,iT,iX) ! Ionisation rate do i=1,nion if (i==1) then ! electron sigv=3.16d-11 * (dsqrt(8d0*kB*1d-7*T/(pi*me*1d-3))*1d-3)**1.3d0 tau_sn(i) = 1.d0/1.16d0*(m(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nH/2.d0*sigv) else if (i>=2 .and. i<=nion) then ! ions muuu=m(i)*2d0*mp/(m(i)+2d0*mp) if (i==2 .or. i==3) then sigv=2.4d-9 *(dsqrt(8d0*kB*1d-7*T/(pi*muuu*1d-3))*1d-3)**0.6d0 else if (i==4) then sigv=2d-9 * (dsqrt(8d0*kB*1d-7*T/(pi*muuu*1d-3))*1d-3)**0.15d0 else if (i==5) then sigv=3.89d-9 * (dsqrt(8d0*kB*1d-7*T/(pi*muuu*1d-3))*1d-3)**(-0.02d0) else sigv=1.69d-9 end if tau_sn(i) = 1.d0/1.14d0*(m(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nH/2.d0*sigv) end if omega(i) = q(i)*B/(m(i)*c_l) sigma(i) = abundances(i,iH,iT,iX)*nH*(q(i))**2*tau_sn(i)/m(i) phi(i) = 0.d0 zetas(i) = 0.d0 gamma_zeta(i) = 1.d0 gamma_omega(i) = 1.d0 omega_bar(i) = 0.d0 end do do i=1,alpha ! grains tau_sn(nion+1+3*(i-1))=1.d0/1.28d0*(m_g(i)+2.d0*mp)/(2.d0*mp)*1.d0/(nH/2.d0*(pi*r_g(i)**2*(8.d0*Kb*T/(pi*2.d0*mp))**0.5)) omega(nion+1+3*(i-1)) = q(nion+1+3*(i-1))*B/(m_g(i)*c_l) sigma(nion+1+3*(i-1)) = abundances(nion+1+3*(i-1),iH,iT,iX)*nH*(q(nion+1+3*(i-1)))**2*tau_sn(nion+1+3*(i-1))/m_g(i) tau_sn(nion+2+3*(i-1))=tau_sn(nion+1+3*(i-1)) omega(nion+2+3*(i-1)) = q(nion+2+3*(i-1))*B/(m_g(i)*c_l) sigma(nion+2+3*(i-1)) = abundances(nion+2+3*(i-1),iH,iT,iX)*nH*(q(nion+2+3*(i-1)))**2*tau_sn(nion+2+3*(i-1))/m_g(i) end do sigP=0.d0 sigO=0.d0 sigH=0.d0 do i=1,nvar sigP=sigP+sigma(i) sigO=sigO+sigma(i)/(1.d0+(omega(i)*tau_sn(i))**2) sigH=sigH-sigma(i)*omega(i)*tau_sn(i)/(1.d0+(omega(i)*tau_sn(i))**2) end do conductivities(1,iH,iT,iX,iB)=log10(sigP) conductivities(2,iH,iT,iX,iB)=log10(sigO) conductivities(3,iH,iT,iX,iB)=log10(abs(sigH)) conductivities(0,iH,iT,iX,iB)=sign(1.0_dp,sigH) !conductivities(:,:,:,:) contains the log10(conductivities) eta_ohm = 1d0/sigP eta_hall = sigH/(sigO**2 + sigH**2) eta_ad = sigO/(sigO**2 + sigH**2)-1d0/sigP end do end do end do end do end program