Re-writing non uniform v_turb

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,22 +543,22 @@ 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
 
   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

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

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')

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

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

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

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