Renaming masse_mol_gaz to mu_mH and masseH to mH

src/ML_prodimo.f90

@@ -195,15 +195,15 @@ subroutine xgb_compute_features()
           feature_Tgas(2,:) = z_grid(:)
        endif
        feature_Tgas(3,:) = Tdust(:)
-       feature_Tgas(4,:) = densite_gaz(:) * masse_mol_gaz / m3_to_cm3 ! g.cm^3
+       feature_Tgas(4,:) = densite_gaz(:) * mu_mH / m3_to_cm3 ! g.cm^3
        feature_Tgas(5:43,:) = J_ML(:,:)
        feature_Tgas(44:47,:) = N_grains(:,:)
        do i=1,n_directions
           feature_Tgas(48+i-1,:) = CD(:,i)  ! CD(n_cells, n_directions)
        enddo
     else if (n_features == 45) then
        feature_Tgas(1,:) = Tdust(:)
-       feature_Tgas(2,:) = densite_gaz(:) * masse_mol_gaz / m3_to_cm3 ! g.cm^3
+       feature_Tgas(2,:) = densite_gaz(:) * mu_mH / m3_to_cm3 ! g.cm^3
        feature_Tgas(3:41,:) = J_ML(:,:)
        feature_Tgas(42:45,:) = N_grains(:,:)
     endif

src/SPH2mcfost.f90

@@ -304,8 +304,8 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g
     call allocate_densities(n_cells_max = n_SPH + n_etoiles) ! we allocate all the SPH particule for libmcfost
     ! Tableau de densite et masse de gaz
     !do icell=1,n_cells
-    !   densite_gaz(icell) = rho(icell) / masse_mol_gaz * m3_to_cm3 ! rho is in g/cm^3 --> part.m^3
-    !   masse_gaz(icell) =  densite_gaz(icell) * masse_mol_gaz * volume(icell)
+    !   densite_gaz(icell) = rho(icell) / mu_mH * m3_to_cm3 ! rho is in g/cm^3 --> part.m^3
+    !   masse_gaz(icell) =  densite_gaz(icell) * mu_mH * volume(icell)
     !enddo
     !masse_gaz(:) = masse_gaz(:) * AU3_to_cm3
 
@@ -314,7 +314,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g
        iSPH = Voronoi(icell)%id
        if (iSPH > 0) then
           masse_gaz(icell)    = massgas(iSPH) * Msun_to_g ! en g
-          densite_gaz(icell)  = masse_gaz(icell) /  (masse_mol_gaz * volume(icell) * AU3_to_m3)
+          densite_gaz(icell)  = masse_gaz(icell) /  (mu_mH * volume(icell) * AU3_to_m3)
        else ! star
           masse_gaz(icell)    = 0.
           densite_gaz(icell)  = 0.
@@ -648,7 +648,7 @@ subroutine Hydro_to_Voronoi_atomic(n_SPH,T_tmp,vt_tmp,mass_gas,mass_ne_on_massga
     ! mask : -1 means skip, 0 means transparent, 1 means compute atomic transfer
     ! ************************************************************************************ !
     use parametres
-    use constantes,   only : masseH
+    use constantes,   only : mH
     use Voronoi_grid, only : Voronoi, volume
     use disk_physics, only : compute_othin_sublimation_radius
     use mem
@@ -673,7 +673,7 @@ subroutine Hydro_to_Voronoi_atomic(n_SPH,T_tmp,vt_tmp,mass_gas,mass_ne_on_massga
     !-> fills element abundances structures for elements
     call alloc_atomrt_grid
     call read_abundance
-    rho_to_nH = 1d3 / masseH / wght_per_H !convert from density to number density nHtot
+    rho_to_nH = 1d3 / mH / wght_per_H !convert from density to number density nHtot
 
     Vxmax = 0
     Vymax = 0

src/benchmarks.f90

@@ -177,8 +177,7 @@ subroutine readMolecule_benchmark2()
 
   read(1,'(a)') mol(1)%name
 
-  read(1,*) molecularWeight
-  masse_mol = masseH * molecularWeight
+  read(1,*) mol(1)%molecularWeight
 
   read(1,*) nLevels, nTrans_tot
 
@@ -304,7 +303,7 @@ subroutine init_benchmark_vanZadelhoff1()
   Tcin = 20._dp
   vfield(:) = 0.0
 
-  masse_mol = 1.0
+  mol(1)%molecularWeight = 1.0
 
   linfall = .true.
   lkeplerian = .false.

src/constants.f90

@@ -31,9 +31,9 @@ module constantes
 
   real(kind=dp), parameter :: Na = 6.022140857e23_dp   ! Nombre d'Avogadro CODATA 2014
   real(kind=dp), parameter :: amu = 1.0_dp/Na          ! atomic mass unit [g]
-  real(kind=dp), parameter :: masseH = 1.007825032231_dp * amu   ! H atom mass [g] CODATA 2014
-  real, parameter :: mu = 2.3                          ! [g]  2.3 following Walker 2004
-  real, parameter :: masse_mol_gaz = mu * masseH
+  real(kind=dp), parameter :: mH = 1.007825032231_dp * amu   ! H atom mass [g] CODATA 2014
+  real, parameter :: mu = 2.3                          ! [mH] 2.3 following Walker 2004
+  real, parameter :: mu_mH = mu * mH                   ! mean molecular weight [g]
   real, parameter :: T_Cmb = 2.7260                    ! K
 
   real(kind=dp), parameter ::  e_ion_hmin = 0.754 * electron_charge !Ionization energy (affinity) Hmin in [J]

src/density.f90

@@ -88,7 +88,7 @@ subroutine define_gas_density()
 
      ! Facteur multiplicatif pour passer en masse de gaz
      ! puis en nH2/AU**3 puis en nH2/m**3
-     cst_gaz(izone) = C * dz%diskmass * dz%gas_to_dust / masse_mol_gaz * Msun_to_g / AU3_to_m3
+     cst_gaz(izone) = C * dz%diskmass * dz%gas_to_dust / mu_mH * Msun_to_g / AU3_to_m3
   enddo
 
   do izone=1, n_zones
@@ -282,7 +282,7 @@ subroutine define_gas_density()
      ! Calcul de la masse de gaz de la zone
      mass = 0.
      do icell = 1, n_cells
-        mass = mass + densite_gaz_tmp(icell) *  masse_mol_gaz * volume(icell)
+        mass = mass + densite_gaz_tmp(icell) *  mu_mH * volume(icell)
      enddo
      mass =  mass * AU3_to_m3 * g_to_Msun
 
@@ -318,7 +318,7 @@ subroutine define_gas_density()
 
   ! Tableau de masse de gaz
   do icell=1,n_cells
-     masse_gaz(icell) =  densite_gaz(icell) * masse_mol_gaz * volume(icell) * AU3_to_m3
+     masse_gaz(icell) =  densite_gaz(icell) * mu_mH * volume(icell) * AU3_to_m3
   enddo
   write(*,*) 'Total  gas mass in model:', real(sum(masse_gaz) * g_to_Msun),' Msun'
 
@@ -442,7 +442,7 @@ subroutine define_dust_density()
 
         do i=1, n_rad
 
-           !write(*,*) "     ", rcyl, rho0*masse_mol_gaz*cm_to_m**2, dust_pop(pop)%rho1g_avg
+           !write(*,*) "     ", rcyl, rho0*mu_mH*cm_to_m**2, dust_pop(pop)%rho1g_avg
            !write(*,*) "s_opt", rcyl, s_opt/1000.
 
            bz : do j=j_start,nz
@@ -678,11 +678,11 @@ subroutine define_dust_density()
               do k=1, n_az
                  rho0 = densite_gaz(cell_map(i,1,k)) ! pour dependance en R : pb en coord sperique
                  !s_opt = rho_g * cs / (rho * Omega)    ! cs = H * Omega ! on doit trouver 1mm vers 50AU
-                 !omega_tau= dust_pop(ipop)%rho1g_avg*(r_grain(l)*mum_to_cm) / (rho * masse_mol_gaz/m_to_cm**3 * H*AU_to_cm)
+                 !omega_tau= dust_pop(ipop)%rho1g_avg*(r_grain(l)*mum_to_cm) / (rho * mu_mH/m_to_cm**3 * H*AU_to_cm)
                  icell = cell_map(i,1,k)
                  rcyl = r_grid(icell)
                  H = dz%sclht * (rcyl/dz%rref)**dz%exp_beta
-                 s_opt = (rho0*masse_mol_gaz*cm_to_m**3  /dust_pop(pop)%rho1g_avg) *  H * AU_to_m * m_to_mum
+                 s_opt = (rho0*mu_mH*cm_to_m**3  /dust_pop(pop)%rho1g_avg) *  H * AU_to_m * m_to_mum
 
                  !write(*,*) "r=", rcyl, "a_migration =", s_opt
 
@@ -1551,7 +1551,7 @@ subroutine read_density_file()
   ! Calcul de la masse de gaz de la zone
   mass = 0.
   do icell=1,n_cells
-     mass = mass + densite_gaz(icell) *  masse_mol_gaz * volume(icell)
+     mass = mass + densite_gaz(icell) *  mu_mH * volume(icell)
   enddo !icell
   mass =  mass * AU3_to_m3 * g_to_Msun
 
@@ -1569,7 +1569,7 @@ subroutine read_density_file()
 
   ! Tableau de masse de gaz
   do icell=1,n_cells
-     masse_gaz(icell) =  densite_gaz(icell) * masse_mol_gaz * volume(icell) * AU3_to_m3
+     masse_gaz(icell) =  densite_gaz(icell) * mu_mH * volume(icell) * AU3_to_m3
   enddo ! icell
   write(*,*) 'Total  gas mass in model:', real(sum(masse_gaz) * g_to_Msun),' Msun'
 
@@ -2010,7 +2010,7 @@ real(kind=dp) function omega_tau(rho,H,l)
   ipop = grain(l)%pop
   !write(*,*) ipop, dust_pop(ipop)%rho1g_avg, rho
   if (rho > tiny_dp) then
-     omega_tau = dust_pop(ipop)%rho1g_avg*(r_grain(l)*mum_to_cm) / (rho * masse_mol_gaz/m_to_cm**3 * H*AU_to_cm)
+     omega_tau = dust_pop(ipop)%rho1g_avg*(r_grain(l)*mum_to_cm) / (rho * mu_mH/m_to_cm**3 * H*AU_to_cm)
   else
      omega_tau = huge_dp
   endif

src/disk_physics.f90

@@ -191,7 +191,7 @@ subroutine equilibre_hydrostatique()
   real, parameter :: gas_dust = 100
 
   M_etoiles = sum(etoile(:)%M) * Msun_to_kg
-  M_mol = masse_mol_gaz * g_to_kg
+  M_mol = mu_mH * g_to_kg
 
   cst = Ggrav * M_etoiles * M_mol / (kb * AU_to_m**2)
 

src/for_astrochem.f90

@@ -7,7 +7,7 @@
   G(:) = compute_UV_field() ! Habing
 
   ! unit conversions for astrochem :
-  gas_density(:) = densite_gaz(1:n_cells) * masse_mol_gaz / m3_to_cm3 ! nH2/m**3 --> g/cm**3
+  gas_density(:) = densite_gaz(1:n_cells) * mu_mH / m3_to_cm3 ! nH2/m**3 --> g/cm**3
 
   do icell=1,n_cells
      dust_mass_density(icell) = sum(densite_pouss(:,icell) * M_grain(:)) ! M_grain en g

src/input.f90

@@ -86,7 +86,7 @@ subroutine readmolecule(imol)
 
   read(1,*) junk
   read(1,*) molecularWeight
-  masse_mol = masseH * molecularWeight
+  mol(imol)%molecularWeight = molecularWeight
 
   read(1,*) junk
   read(1,*) nLevels

src/io_prodimo.f90

@@ -789,7 +789,7 @@ subroutine mcfost2ProDiMo()
     do ri=1, n_rad
        do zj=1,nz
           icell = cell_map(ri,zj,1)
-          dens(ri,zj) =  densite_gaz(icell) * masse_mol_gaz / m3_to_cm3 ! g.cm^-3
+          dens(ri,zj) =  densite_gaz(icell) * mu_mH / m3_to_cm3 ! g.cm^-3
        enddo
     enddo
 

src/mcfost2phantom.f90

@@ -569,7 +569,7 @@ subroutine diffusion_opacity(temp,icell,kappa_diffusion)
           somme2 = somme2 + B*delta_wl
        enddo
        kappa_diffusion = somme2/somme &
-          *cm_to_AU/(densite_gaz(icell)*masse_mol_gaz*(cm_to_m)**3) ! cm^2/g
+          *cm_to_AU/(densite_gaz(icell)*mu_mH*(cm_to_m)**3) ! cm^2/g
        ! check : somme2/somme * cm_to_AU /(masse_gaz(icell)/(volume(icell)*AU_to_cm**3))
     else
        kappa_diffusion = 0.

src/mhd2mcfost.f90

@@ -261,7 +261,7 @@ end subroutine setup_mhd_to_mcfost
    !      close(unit=1)
 
    !      !rho -> nH
-   !      nHtot = nHtot * 1d3 / masseH / wght_per_H
+   !      nHtot = nHtot * 1d3 / mH / wght_per_H
 
    !      write(*,*) "Read ", size(pack(icompute_atomRT,mask=icompute_atomRT>0)), " density zones"
    !      write(*,*) "Read ", size(pack(icompute_atomRT,mask=icompute_atomRT==0)), " transparent zones"

src/mol_transfer.f90

@@ -929,7 +929,7 @@ subroutine init_dust_mol(imol)
            ! AU_to_cm**2 car on veut kappa_abs_LTE en AU-1
            write(*,*) "TODO : the water benchmark 3 needs to be updated for cell pointer in opacity table"
            do icell=1,n_cells
-              kappa_abs_LTE(icell,iTrans) =  kap * (densite_gaz(icell) * cm_to_m**3) * masse_mol_gaz / &
+              kappa_abs_LTE(icell,iTrans) =  kap * (densite_gaz(icell) * cm_to_m**3) * mu_mH / &
                    gas_dust / cm_to_AU
            enddo
 

src/molecular_emission.f90

@@ -33,8 +33,7 @@ module molecular_emission
   real :: correct_Tgas
   logical :: lcorrect_Tgas
 
-  real :: nH2, masse_mol
-  ! masse_mol_gaz sert uniquement pour convertir masse disque en desnite de particule
+  real :: nH2, masse_mol ! masse_mol is the mass of the current molecule
   real(kind=dp), dimension(:,:), allocatable :: kappa_mol_o_freq, kappa_mol_o_freq2 ! n_cells, nTrans
   real(kind=dp), dimension(:,:), allocatable :: emissivite_mol_o_freq,  emissivite_mol_o_freq2 ! n_cells, nTrans
   real, dimension(:,:), allocatable :: tab_nLevel, tab_nLevel2 ! n_cells, nLevels
@@ -73,7 +72,7 @@ module molecular_emission
 
   type molecule
      integer :: n_speed_rt, n_speed_center_rt, n_extraV_rt, nTrans_raytracing, iLevel_max
-     real :: vmin_center_rt, vmax_center_rt, extra_deltaV_rt, abundance
+     real :: vmin_center_rt, vmax_center_rt, extra_deltaV_rt, abundance, molecularWeight
      logical :: lcst_abundance, lline, l_sym_ima
      character(len=512) :: abundance_file, filename
      character(len=32) :: name
@@ -158,9 +157,9 @@ subroutine init_Doppler_profiles(imol)
      ! WARNING : c'est pas un sigma mais un delta, cf Cours de Boisse p47
      ! Consistent avec benchmark
      if (trim(v_turb_unit) == "cs") then
-        v_turb(icell) = sqrt((kb*Tcin(icell) / (masse_mol_gaz * g_to_kg))) * v_turb(icell)
+        v_turb(icell) = sqrt((kb*Tcin(icell) / (mu_mH * g_to_kg))) * v_turb(icell)
      endif
-     sigma2 =  2.0_dp * (kb*Tcin(icell) / (masse_mol * g_to_kg)) + v_turb(icell)**2
+     sigma2 =  2.0_dp * (kb*Tcin(icell) / (mol(imol)%molecularWeight * mH  * g_to_kg)) + v_turb(icell)**2
      v_line(icell) = sqrt(sigma2)
 
      !  write(*,*) "FWHM", sqrt(sigma2 * log(2.)) * 2.  ! Teste OK bench water 1
@@ -421,8 +420,8 @@ subroutine equilibre_LTE_mol(imol)
   !$omp end do
   !$omp  end parallel
 
-  ! write(*,*) real( (sum(masse) * g_to_kg * gas_dust / masse_mol_gaz) / (4.*pi/3. * (rout * AU_to_cm)**3 ) )
-  ! write(*,*) (sum(masse) * g_to_kg * gas_dust / masse_mol_gaz) * abundance * fAul(1) * hp * transfreq(1)
+  ! write(*,*) real( (sum(masse) * g_to_kg * gas_dust / mu_mH) / (4.*pi/3. * (rout * AU_to_cm)**3 ) )
+  ! write(*,*) (sum(masse) * g_to_kg * gas_dust / mu_mH) * abundance * fAul(1) * hp * transfreq(1)
 
   return
 

src/optical_depth.f90

@@ -341,7 +341,7 @@ subroutine compute_column(type, column, lambda)
   real(kind=dp) :: x0,y0,z0, x1,y1,z1, norme, l, u,v,w, l_contrib, l_void_before, CD_units, factor, sum
 
   if (type==1) then
-     CD_units = AU_to_m * masse_mol_gaz / (m_to_cm)**2 ! g/cm^-2 and AU_to_m factor as l_contrib is in AU
+     CD_units = AU_to_m * mu_mH / (m_to_cm)**2 ! g/cm^-2 and AU_to_m factor as l_contrib is in AU
   else if (type==3) then
      CD_units = AU_to_m / (m_to_cm)**2 ! particle/cm^-2 and AU_to_m factor as l_contrib is in AU
   endif

src/output.f90

@@ -1680,7 +1680,7 @@ subroutine write_disk_struct(lparticle_density,lcolumn_density,lvelocity)
 
   dens = 0.0
 
-  dens(:) = densite_gaz(1:n_cells) * masse_mol_gaz / m3_to_cm3 ! nH2/m**3 --> g/cm**3
+  dens(:) = densite_gaz(1:n_cells) * mu_mH / m3_to_cm3 ! nH2/m**3 --> g/cm**3
 
   ! le e signifie real*4
   call ftppre(unit,group,fpixel,nelements,dens,status)

src/read1d_models.f90

@@ -23,13 +23,13 @@ module read1d_models
 	logical, parameter :: lcell_centered = .true.
 
 	contains
-	
-	
+
+
 	subroutine read_model_1d(mod_file)
 		character(len=*), intent(in) :: mod_file
 		integer :: i, nrad_0
 		real(kind=dp) :: Fnorm, dnu, Rmax_c
-		
+
 		grid_type = 2
 		n_rad_in = 1
 		n_az = 1
@@ -40,7 +40,7 @@ subroutine read_model_1d(mod_file)
 		vfield_coord = 3
 		lmagnetized = .false.
 		lcalc_ne = .false.
-		
+
 		open(unit=1,file=mod_file, status="old")
 		read(1,*) atmos_1d%rstar
 		read(1,*) atmos_1d%nr
@@ -91,7 +91,7 @@ subroutine read_model_1d(mod_file)
 
 		n_rad = atmos_1d%nr - 1
 		n_cells = n_rad
-		n_etoiles = 1 !force		
+		n_etoiles = 1 !force
 
 		Rmin = minval(atmos_1d%r) * atmos_1d%rstar * m_to_au
 		Rmax = maxval(atmos_1d%r,mask=(atmos_1d%iz>0)) * atmos_1d%rstar* m_to_au
@@ -127,7 +127,7 @@ subroutine read_model_1d(mod_file)
 		write(*,*) "WARNING distance and map size set to : "
 		distance = Rmax * au_to_pc !pc
 		map_size = 2.005 * Rmax !au
-		write(*,*) distance, ' pc', map_size * 1e3, 'mau'	
+		write(*,*) distance, ' pc', map_size * 1e3, 'mau'
 
 		return
 	end subroutine read_model_1d
@@ -139,7 +139,7 @@ subroutine setup_model1d_to_mcfost()
 
 		call alloc_atomrt_grid()
 		call read_abundance
-		rho_to_nH = 1d3 / masseH / wght_per_H
+		rho_to_nH = 1d3 / mH / wght_per_H
 
 		!- MCFOST is a cell centered radiative transfer code while the 1D spherically symmetric codes
 		! are based on a set of points corresponding to the nodes of multi-dimensional models.
@@ -157,7 +157,7 @@ subroutine setup_model1d_to_mcfost()
 		!allows to properly balance that effect.
 		if (lcell_centered) then
 			!core temperature, read from file as "Tstar"
-			! etoile(1)%T = real(	atmos_1d%T(1) ) ! + something else here 
+			! etoile(1)%T = real(	atmos_1d%T(1) ) ! + something else here
 			icell = n_cells + 1
 			nrad_0 = n_rad
 			icell0 = n_cells
@@ -188,7 +188,7 @@ subroutine setup_model1d_to_mcfost()
 				icompute_atomRT(n_cells) = atmos_1d%iz(n_rad+1)
 			endif
 			!core temperature, read from file as "Tstar"
-			! etoile(1)%T = real(	atmos_1d%T(1) ) ! + irrad ? 
+			! etoile(1)%T = real(	atmos_1d%T(1) ) ! + irrad ?
 			icompute_atomRT(:) = atmos_1d%iz(2:n_cells+1)
 			T(:) = atmos_1d%T(2:n_cells+1)
 			nHtot(:) = atmos_1d%rho(2:n_cells+1) * rho_to_nH
@@ -299,5 +299,5 @@ subroutine check_for_coronal_illumination()
 
 		return
 	end subroutine check_for_coronal_illumination
-	
-end module read1d_models
+
+end module read1d_models

src/read_athena++.f90

@@ -1,6 +1,7 @@
 module read_athena
 
   use parametres
+  use constantes
   use messages
   use mcfost_env
   use grid
@@ -333,7 +334,7 @@ subroutine read_athena_model()
     ! Calcul de la masse de gaz de la zone
     mass = 0.
     do icell=1,n_cells
-       mass = mass + densite_gaz(icell) *  masse_mol_gaz * volume(icell)
+       mass = mass + densite_gaz(icell) *  mu_mH * volume(icell)
     enddo !icell
     mass =  mass * AU3_to_m3 * g_to_Msun
 
@@ -344,7 +345,7 @@ subroutine read_athena_model()
        ! Somme sur les zones pour densite finale
        do icell=1,n_cells
           densite_gaz(icell) = densite_gaz(icell) * facteur
-          masse_gaz(icell) = densite_gaz(icell) * masse_mol_gaz * volume(icell) * AU3_to_m3
+          masse_gaz(icell) = densite_gaz(icell) * mu_mH * volume(icell) * AU3_to_m3
        enddo ! icell
     else
        call error('Gas mass is 0')