diff --git a/src/ML_prodimo.f90 b/src/ML_prodimo.f90 index de50db64..9f7876ee 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 b69feef6..ec16dac7 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. @@ -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 diff --git a/src/benchmarks.f90 b/src/benchmarks.f90 index 89affda4..49403f31 100644 --- a/src/benchmarks.f90 +++ b/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. diff --git a/src/constants.f90 b/src/constants.f90 index 1a7d3ea2..9c06d417 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 882dcca8..d5d89f54 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 c8bf862e..63a57829 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 7bb7e9bd..856a18bc 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/input.f90 b/src/input.f90 index 0381802a..21662abd 100644 --- a/src/input.f90 +++ b/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 diff --git a/src/io_prodimo.f90 b/src/io_prodimo.f90 index b2d6a064..c1a33586 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 66775cdf..5a504d5f 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/mhd2mcfost.f90 b/src/mhd2mcfost.f90 index 2f8c3cf0..16987a3a 100644 --- a/src/mhd2mcfost.f90 +++ b/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" diff --git a/src/mol_transfer.f90 b/src/mol_transfer.f90 index 73f0fcbc..0c10bfc4 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 1c9063e9..de47a25e 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 @@ -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 diff --git a/src/optical_depth.f90 b/src/optical_depth.f90 index 7ed82208..aa8ce082 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 e125451e..b6417817 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/read1d_models.f90 b/src/read1d_models.f90 index 570f24de..b8b3e38b 100644 --- a/src/read1d_models.f90 +++ b/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 \ No newline at end of file + +end module read1d_models diff --git a/src/read_athena++.f90 b/src/read_athena++.f90 index c4ba400c..d085c308 100644 --- a/src/read_athena++.f90 +++ b/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') diff --git a/src/read_fargo3d.f90 b/src/read_fargo3d.f90 index 15a7f1a7..e2f56091 100644 --- a/src/read_fargo3d.f90 +++ b/src/read_fargo3d.f90 @@ -6,6 +6,7 @@ module read_fargo3d use grid use cylindrical_grid use density + use constantes use stars, only : compute_stellar_parameters implicit none @@ -273,7 +274,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 +285,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 09fc0903..3955c116 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 682c103a..614aee00 100644 --- a/src/read_pluto.f90 +++ b/src/read_pluto.f90 @@ -1,6 +1,7 @@ module read_pluto use parametres + use constantes use messages use mcfost_env use grid @@ -258,7 +259,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 +270,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 a22fefd8..e2a92a82 100644 --- a/src/read_spherical_grid.f90 +++ b/src/read_spherical_grid.f90 @@ -184,15 +184,15 @@ subroutine read_spherical_model(filename) ! do k=1, n_az ! icell = cell_map(i,jj,k) ! T(icell) = T_tmp(i,j,k) - ! nHtot(icell) = rho(i,j,k) * 1d3 / masseH / wght_per_H ! [H/m^3] + ! nHtot(icell) = rho(i,j,k) * 1d3 / mH / wght_per_H ! [H/m^3] ! icompute_atomRT(icell) = dz(i,j,k) ! vfield3d(icell,:) = vtmp(i,j,k,:) ! !-> 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] @@ -219,15 +219,15 @@ subroutine read_spherical_model(filename) icell = cell_map(i,jj,k) endif T(icell) = T_tmp(i,j,k) - nHtot(icell) = rho(i,j,k) * 1d3 / masseH / wght_per_H ! [H/m^3] + nHtot(icell) = rho(i,j,k) * 1d3 / mH / wght_per_H ! [H/m^3] icompute_atomRT(icell) = dz(i,j,k) vfield3d(icell,:) = vtmp(i,j,k,:) !-> 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] @@ -245,9 +245,9 @@ subroutine read_spherical_model(filename) !dust part see read_pluto.f90 ! ********************************** ! mass = sum(masse_gaz) * g_to_Msun - ! mass = sum(densite_gaz * volume) * AU3_to_m3 * masseH * g_to_Msun + ! mass = sum(densite_gaz * volume) * AU3_to_m3 * mH * g_to_Msun if (mass <= 0.0) call error('Gas mass is 0') - ! masse_gaz(:) = masseH * densite_gaz(:) * volume(:) * AU3_to_m3 ! [g] + ! masse_gaz(:) = mH * densite_gaz(:) * volume(:) * AU3_to_m3 ! [g] !--> no normalisation of density. The dust and gas densities are provided in the model. diff --git a/src/stars.f90 b/src/stars.f90 index 50142365..cbe1c60a 100644 --- a/src/stars.f90 +++ b/src/stars.f90 @@ -992,7 +992,7 @@ function is_inshock(id, iray, i_star, icell0, x, y, z, Thp, Tshock, Facc) if (vaccr < 0.0_dp) then !Facc = 1/2 rho vs^3 - Facc = 0.5 * (1d-3 * masseH * rho) * abs(vaccr)**3 + Facc = 0.5 * (1d-3 * mH * rho) * abs(vaccr)**3 Tloc = ( 0.75 * Facc / sigma )**0.25 ! is_inshock = (Tloc > 0.5 * etoile(i_star)%T) is_inshock = (T_hp > 1.0_dp * etoile(i_star)%T) @@ -1003,7 +1003,7 @@ function is_inshock(id, iray, i_star, icell0, x, y, z, Thp, Tshock, Facc) Thp = abs(T_hp) * Tloc endif !assuming mu is 0.5 - Tshock = 0.5 * (3.0/16.0) * (1d-3 * masseH) / kb * vaccr**2 + Tshock = 0.5 * (3.0/16.0) * (1d-3 * mH) / kb * vaccr**2 max_Thp = max(max_Thp, Thp); min_Thp = min(min_Thp, Thp) max_Tshock = max(max_Tshock, Tshock); min_Tshock = min(min_Tshock, Tshock) max_Facc = max(max_Facc,Facc); min_Facc = min(min_Facc, Facc)