From 9d2a9860143e7ab6ddb0343830fcaed295f1a1c8 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Wed, 11 Dec 2024 20:30:23 +1100 Subject: [PATCH] Renaming masse_mol_gaz to mu_mH --- src/ML_prodimo.f90 | 4 ++-- src/SPH2mcfost.f90 | 6 +++--- src/constants.f90 | 6 +++--- src/density.f90 | 18 +++++++++--------- src/disk_physics.f90 | 2 +- src/for_astrochem.f90 | 2 +- src/io_prodimo.f90 | 2 +- src/mcfost2phantom.f90 | 2 +- src/mol_transfer.f90 | 2 +- src/molecular_emission.f90 | 9 ++++----- src/optical_depth.f90 | 2 +- src/output.f90 | 2 +- src/read_athena++.f90 | 4 ++-- src/read_fargo3d.f90 | 4 ++-- src/read_idefix.f90 | 4 ++-- src/read_pluto.f90 | 4 ++-- src/read_spherical_grid.f90 | 8 ++++---- 17 files changed, 40 insertions(+), 41 deletions(-) diff --git a/src/ML_prodimo.f90 b/src/ML_prodimo.f90 index de50db641..9f7876eef 100644 --- a/src/ML_prodimo.f90 +++ b/src/ML_prodimo.f90 @@ -195,7 +195,7 @@ 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 @@ -203,7 +203,7 @@ subroutine xgb_compute_features() 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 diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index b69feef6c..ec46a38e6 100644 --- a/src/SPH2mcfost.f90 +++ b/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. diff --git a/src/constants.f90 b/src/constants.f90 index 1a7d3ea24..9c06d417d 100644 --- a/src/constants.f90 +++ b/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] diff --git a/src/density.f90 b/src/density.f90 index 882dcca8a..d5d89f547 100644 --- a/src/density.f90 +++ b/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 diff --git a/src/disk_physics.f90 b/src/disk_physics.f90 index c8bf862e0..63a578294 100644 --- a/src/disk_physics.f90 +++ b/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) diff --git a/src/for_astrochem.f90 b/src/for_astrochem.f90 index 7bb7e9bdf..856a18bc9 100644 --- a/src/for_astrochem.f90 +++ b/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 diff --git a/src/io_prodimo.f90 b/src/io_prodimo.f90 index b2d6a0640..c1a335866 100644 --- a/src/io_prodimo.f90 +++ b/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 diff --git a/src/mcfost2phantom.f90 b/src/mcfost2phantom.f90 index 66775cdff..5a504d5f2 100644 --- a/src/mcfost2phantom.f90 +++ b/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. diff --git a/src/mol_transfer.f90 b/src/mol_transfer.f90 index 73f0fcbc8..0c10bfc45 100644 --- a/src/mol_transfer.f90 +++ b/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 diff --git a/src/molecular_emission.f90 b/src/molecular_emission.f90 index 1c9063e9b..82a8435ac 100644 --- a/src/molecular_emission.f90 +++ b/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 @@ -158,7 +157,7 @@ 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 v_line(icell) = sqrt(sigma2) @@ -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 diff --git a/src/optical_depth.f90 b/src/optical_depth.f90 index 7ed822083..aa8ce0825 100644 --- a/src/optical_depth.f90 +++ b/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 diff --git a/src/output.f90 b/src/output.f90 index e125451e2..b6417817c 100644 --- a/src/output.f90 +++ b/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) diff --git a/src/read_athena++.f90 b/src/read_athena++.f90 index c4ba400c0..f55ccc189 100644 --- a/src/read_athena++.f90 +++ b/src/read_athena++.f90 @@ -333,7 +333,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 +344,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') diff --git a/src/read_fargo3d.f90 b/src/read_fargo3d.f90 index 15a7f1a75..5cedbdc3d 100644 --- a/src/read_fargo3d.f90 +++ b/src/read_fargo3d.f90 @@ -273,7 +273,7 @@ subroutine read_fargo3d_files() ! 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 @@ -284,7 +284,7 @@ subroutine read_fargo3d_files() ! 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') diff --git a/src/read_idefix.f90 b/src/read_idefix.f90 index 09fc0903e..3955c1164 100644 --- a/src/read_idefix.f90 +++ b/src/read_idefix.f90 @@ -210,7 +210,7 @@ subroutine read_idefix_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 @@ -221,7 +221,7 @@ subroutine read_idefix_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') diff --git a/src/read_pluto.f90 b/src/read_pluto.f90 index 682c103a1..d6c208c7c 100644 --- a/src/read_pluto.f90 +++ b/src/read_pluto.f90 @@ -258,7 +258,7 @@ subroutine read_pluto_files() ! 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 @@ -269,7 +269,7 @@ subroutine read_pluto_files() ! 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') diff --git a/src/read_spherical_grid.f90 b/src/read_spherical_grid.f90 index a22fefd83..8ed34213d 100644 --- a/src/read_spherical_grid.f90 +++ b/src/read_spherical_grid.f90 @@ -191,8 +191,8 @@ subroutine read_spherical_model(filename) ! !-> wrapper for dust RT. ! !-> taking into account proper weights assuming only molecular gas in dusty regions ! if (rho_dust(i,j,k) > 0.0) then ! dusty region - ! densite_gaz(icell) = rho(i,j,k) * 1d3 / masse_mol_gaz !total molecular gas density in H2/m^3 - ! densite_pouss(:,icell) = rho_dust(i,j,k) * 1d3 / masse_mol_gaz ! [m^-3] + ! densite_gaz(icell) = rho(i,j,k) * 1d3 / mu_mH !total molecular gas density in H2/m^3 + ! densite_pouss(:,icell) = rho_dust(i,j,k) * 1d3 / mu_mH ! [m^-3] ! disk_zone(1)%diskmass = disk_zone(1)%diskmass + rho_dust(i,j,k) * volume(icell) ! else !No dust. ! densite_gaz(icell) = nHtot(icell) * wght_per_H !total atomic gas density in [m^-3] @@ -226,8 +226,8 @@ subroutine read_spherical_model(filename) !-> wrapper for dust RT. !-> taking into account proper weights assuming only molecular gas in dusty regions if (rho_dust(i,j,k) > 0.0) then ! dusty region - densite_gaz(icell) = rho(i,j,k) * 1d3 / masse_mol_gaz !total molecular gas density in H2/m^3 - densite_pouss(:,icell) = rho_dust(i,j,k) * 1d3 / masse_mol_gaz ! [m^-3] + densite_gaz(icell) = rho(i,j,k) * 1d3 / mu_mH !total molecular gas density in H2/m^3 + densite_pouss(:,icell) = rho_dust(i,j,k) * 1d3 / mu_mH ! [m^-3] disk_zone(1)%diskmass = disk_zone(1)%diskmass + rho_dust(i,j,k) * volume(icell) else !No dust. densite_gaz(icell) = nHtot(icell) * wght_per_H !total atomic gas density in [m^-3]