diff --git a/src/benchmarks.f90 b/src/benchmarks.f90 index 49403f31..d710c467 100644 --- a/src/benchmarks.f90 +++ b/src/benchmarks.f90 @@ -284,7 +284,7 @@ subroutine init_HH_30_mol() vfield(icell) = 2.0 * (r_grid(icell)/100.)**(-0.55) enddo - v_turb = 230. + v_turb2 = 230.**2 return @@ -298,7 +298,7 @@ subroutine init_benchmark_vanZadelhoff1() ldust_mol = .false. - v_turb = 150._dp !m.s-1 + v_turb2 = 150._dp**2 !m.s-1 Tdust = 20. Tcin = 20._dp vfield(:) = 0.0 @@ -382,13 +382,13 @@ subroutine init_benchmark_vanzadelhoff2() Tdust(icell) = tmp_T(l-1) + frac * (tmp_T(l) - tmp_T(l-1)) Tcin(icell) = tmp_T(l-1) + frac * (tmp_T(l) - tmp_T(l-1)) vfield(icell) = tmp_v(l-1) + frac * (tmp_v(l) - tmp_v(l-1)) - v_turb(icell) = tmp_vturb(l-1) + frac * (tmp_vturb(l) - tmp_vturb(l-1)) + v_turb2(icell) = (tmp_vturb(l-1) + frac * (tmp_vturb(l) - tmp_vturb(l-1)))**2 enddo enddo ! Conversion vitesses en m.s-1 vfield = vfield * 1.0e3 - v_turb = v_turb * 1.0e3 + v_turb2 = v_turb2 * 1.0e6 ! Conversion part.m-3 densite_gaz(:) = densite_gaz(:) / (cm_to_m)**3 @@ -416,7 +416,7 @@ subroutine init_benchmark_water1() densite_gaz = 1.e4 / (cm_to_m)**3 ! part.m-3 Tcin = 40. vfield = 0.0 - v_turb = 0.0 + v_turb2 = 0.0 linfall = .true. lkeplerian = .false. @@ -445,7 +445,7 @@ subroutine init_benchmark_water2() densite_gaz = 1.e4 / (cm_to_m)**3 ! part.m-3 Tcin = 40. - v_turb = 0.0 + v_turb2 = 0.0 do icell=1, n_cells vfield(icell) = 1e5 * sqrt(r_grid(icell)**2 + z_grid(icell)**2) * AU_to_pc @@ -543,10 +543,10 @@ subroutine init_benchmark_water3() if (rayon < 5.95) then vfield(icell) = 0.0 - v_turb(icell) = 3. + v_turb2(icell) = 9. else vfield(icell) = exp( log_tmp_v(l-1) + frac * (log_tmp_v(l) - log_tmp_v(l-1)) ) - v_turb(icell) = 1. + v_turb2(icell) = 1. endif enddo endif @@ -554,11 +554,11 @@ subroutine init_benchmark_water3() enddo ! Conversion FWHM ---> vitesse - v_turb = v_turb / (2.*sqrt(log(2.))) + v_turb2 = v_turb2 / (2.*sqrt(log(2.)))**2 ! Conversion vitesses en m.s-1 vfield = vfield * 1.0e3 - v_turb = v_turb * 1.0e3 + v_turb2 = v_turb2 * 1.0e6 ! Conversion part.m-3 densite_gaz(:) = densite_gaz(:) / (cm_to_m)**3 diff --git a/src/io_prodimo.f90 b/src/io_prodimo.f90 index c1a33586..4fa7560f 100644 --- a/src/io_prodimo.f90 +++ b/src/io_prodimo.f90 @@ -1962,16 +1962,15 @@ subroutine read_ProDiMo2mcfost(imol) ! Vitesse Doppler do i=1, n_rad do j=1, nz - if (lCII) sigma2 = dvCII(i,j)**2 !/ (2 * log(2.)) - if (lOI) sigma2 = dvOI(i,j)**2 - if (lCO) sigma2 = dvCO(i,j)**2 - if (loH2O) sigma2 = dvoH2O(i,j)**2 - if (lpH2O) sigma2 = dvpH2O(i,j)**2 icell = cell_map(i,j,1) - v_line(icell) = sqrt(sigma2) - ! write(*,*) "FWHM", sqrt(sigma2 * log(2.)) * 2. ! Teste OK bench water 1 - if (sigma2 <=0.) call error("ProDiMo data, dv = 0") + if (lCII) dv_line(icell) = dvCII(i,j)**2 !/ (2 * log(2.)) + if (lOI) dv_line(icell) = dvOI(i,j)**2 + if (lCO) dv_line(icell) = dvCO(i,j)**2 + if (loH2O) dv_line(icell) = dvoH2O(i,j)**2 + if (lpH2O) dv_line(icell) = dvpH2O(i,j)**2 + + sigma2 = dv_line(icell)**2 sigma2_m1 = 1.0_dp / sigma2 sigma2_phiProf_m1(icell) = sigma2_m1 diff --git a/src/mem.f90 b/src/mem.f90 index 2141db20..e72bfcb8 100644 --- a/src/mem.f90 +++ b/src/mem.f90 @@ -296,9 +296,9 @@ subroutine alloc_dynamique(n_cells_max) if (alloc_status > 0) call error('Allocation error vfield') vfield=0.0 !; vx=0.0 ; vy=0.0 - allocate(v_turb(Nc), v_line(Nc), deltaVmax(Nc), stat=alloc_status) + allocate(v_turb2(Nc), dv_line(Nc), deltaVmax(Nc), stat=alloc_status) if (alloc_status > 0) call error('Allocation error sigma2') - v_turb = 0.0 ; v_line = 0.0 ; deltaVmax = 0.0 + v_turb2 = 0.0 ; dv_line = 0.0 ; deltaVmax = 0.0 allocate(tab_dnu_o_freq(Nc), stat=alloc_status) if (alloc_status > 0) call error('Allocation error tab_dnu') diff --git a/src/mol_transfer.f90 b/src/mol_transfer.f90 index 0c10bfc4..0bb1aef0 100644 --- a/src/mol_transfer.f90 +++ b/src/mol_transfer.f90 @@ -997,6 +997,7 @@ subroutine init_molecular_disk(imol) integer, intent(in) :: imol logical, save :: lfirst_time = .true. + real(dp) :: factor integer :: icell ldust_mol = .true. @@ -1024,7 +1025,15 @@ subroutine init_molecular_disk(imol) enddo endif endif - v_turb = vitesse_turb + + if (lvturb_in_cs) then + factor = vitesse_turb**2 + do icell=1, n_cells + v_turb2(icell) = (kb*Tcin(icell) / (mu_mH * g_to_kg)) * factor ! cs**2 * factor + enddo + else + v_turb2(:) = vitesse_turb**2 ! constant vturb + endif endif ! lfirst_time ! Abondance diff --git a/src/molecular_emission.f90 b/src/molecular_emission.f90 index de47a25e..4e4ebf8e 100644 --- a/src/molecular_emission.f90 +++ b/src/molecular_emission.f90 @@ -38,8 +38,9 @@ module molecular_emission 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 - real, dimension(:), allocatable :: v_turb, v_line ! n_cells + real, dimension(:), allocatable :: v_turb2, dv_line ! n_cells, v_turb2 is (v_turb in m/s)**2 + logical :: lvturb_in_cs real :: vitesse_turb, dv, dnu, v_syst character(len=8) :: v_turb_unit integer, parameter :: n_largeur_Doppler = 15 @@ -147,7 +148,7 @@ subroutine init_Doppler_profiles(imol) integer, intent(in) :: imol - real(kind=dp) :: sigma2, sigma2_m1, vmax + real(kind=dp) :: sigma2, sigma2_m1 integer :: icell, n_speed n_speed = mol(imol)%n_speed_rt @@ -156,11 +157,8 @@ subroutine init_Doppler_profiles(imol) ! Utilisation de la temperature LTE de la poussiere comme temperature cinetique ! 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) / (mu_mH * g_to_kg))) * v_turb(icell) - endif - sigma2 = 2.0_dp * (kb*Tcin(icell) / (mol(imol)%molecularWeight * mH * g_to_kg)) + v_turb(icell)**2 - v_line(icell) = sqrt(sigma2) + sigma2 = 2.0_dp * (kb*Tcin(icell) / (mol(imol)%molecularWeight * mH * g_to_kg)) + v_turb2(icell) + dv_line(icell) = sqrt(sigma2) ! write(*,*) "FWHM", sqrt(sigma2 * log(2.)) * 2. ! Teste OK bench water 1 sigma2_m1 = 1.0_dp / sigma2 @@ -172,9 +170,8 @@ subroutine init_Doppler_profiles(imol) ! Echantillonage du profil de vitesse dans la cellule ! 2.15 correspond a l'enfroit ou le profil de la raie faut 1/100 de ! sa valeur au centre : exp(-2.15^2) = 0.01 - vmax = sqrt(sigma2) - tab_dnu_o_freq(icell) = largeur_profile * vmax / (real(n_speed)) - deltaVmax(icell) = largeur_profile * vmax !* 2.0_dp ! facteur 2 pour tirage aleatoire + tab_dnu_o_freq(icell) = largeur_profile * dv_line(icell) / (real(n_speed)) + deltaVmax(icell) = largeur_profile * dv_line(icell) !* 2.0_dp ! facteur 2 pour tirage aleatoire enddo !icell return diff --git a/src/optical_depth.f90 b/src/optical_depth.f90 index aa8ce082..c9a6c75f 100644 --- a/src/optical_depth.f90 +++ b/src/optical_depth.f90 @@ -890,7 +890,7 @@ function local_line_profile(icell,lsubtract_avg, x0,y0,z0,x1,y1,z1,u,v,w,l_void_ ! Nbre de points d'integration en fct du differentiel de vitesse ! compare a la largeur de raie de la cellule de depart - n_vpoints = min(max(2,nint(dv/v_line(icell)*20.)),n_vpoints_max) + n_vpoints = min(max(2,nint(dv/dv_line(icell)*20.)),n_vpoints_max) ! Vitesse projete le long du trajet dans la cellule do ivpoint=2, n_vpoints-1 diff --git a/src/read_param.f90 b/src/read_param.f90 index 458cedec..ea416472 100644 --- a/src/read_param.f90 +++ b/src/read_param.f90 @@ -444,10 +444,17 @@ subroutine read_para(para) read(1,*) lpop, lprecise_pop, lmol_LTE largeur_profile = 15. read(1,*) vitesse_turb, v_turb_unit - if (v_turb_unit(1:2) == "km") then - vitesse_turb = vitesse_turb * km_to_m ! Conversion en m.s-1 - else if (trim(v_turb_unit) /= "cs") then - call error("Turbulence velocity unit not recognised") + if (trim(v_turb_unit) == "cs") then + lvturb_in_cs = .true. + else + lvturb_in_cs = .true. + if (v_turb_unit(1:1) /= "m") then + if (v_turb_unit(1:2) == "km") then + vitesse_turb = vitesse_turb * km_to_m ! Conversion en m.s-1 + else + call error("Turbulence velocity unit not recognised") + endif + endif endif read(1,*) n_molecules