Merge pull request #111 from cpinte/UV_field

Option -output_UV_field with Voronoi

src/mem.f90

@@ -81,7 +81,7 @@ subroutine alloc_dust_prop()
      allocate(tab_mueller(4,4,0:nang_scatt,n_grains_tot,n_lambda), stat=alloc_status)
      if (alloc_status > 0) call error('Allocation error tab_mueller')
      tab_mueller = 0
-  else 
+  else
      allocate(tab_s11(0:nang_scatt,n_grains_tot,n_lambda), stat=alloc_status)
      if (alloc_status > 0) call error('Allocation error tab_s11')
      tab_s11 = 0
@@ -121,6 +121,7 @@ subroutine alloc_dynamique(n_cells_max)
   use stars, only : allocate_stellar_spectra
   use thermal_emission, only : allocate_temperature, allocate_thermal_emission, &
        allocate_weight_proba_emission, allocate_thermal_energy
+  use radiation_field, only :  allocate_radiation_field_step1
   use output, only : allocate_origin
 
   integer, intent(in), optional :: n_cells_max
@@ -217,11 +218,11 @@ subroutine alloc_dynamique(n_cells_max)
      if (alloc_status > 0) call error('Allocation error prob_s11_pos')
      prob_s11_pos = 0
 
-     if (lmueller) then 
+     if (lmueller) then
      	allocate(tab_mueller_pos(4,4,0:nang_scatt, p_Nc, p_n_lambda_pos), stat=alloc_status)
      	if (alloc_status > 0) call error('Allocation error tab_mueller_pos')
      	tab_mueller_pos = 0
-     else 
+     else
 	allocate(tab_s11_pos(0:nang_scatt, p_Nc, p_n_lambda_pos), stat=alloc_status)
 	if (alloc_status > 0) call error('Allocation error tab_s11_pos')
 	tab_s11_pos = 0
@@ -269,6 +270,8 @@ subroutine alloc_dynamique(n_cells_max)
   ! **************************************************
   if (lTemp) call allocate_thermal_emission(Nc, p_Nc)
 
+  call allocate_radiation_field_step1(Nc)
+
   ! **************************************************
   ! Tableaux relatifs aux SEDs
   ! **************************************************
@@ -378,7 +381,6 @@ subroutine realloc_dust_mol(imol)
   integer, intent(in) :: imol
 
   integer :: alloc_status, iTrans_min, iTrans_max
-  real :: mem_size
 
   iTrans_min = mol(imol)%iTrans_min
   iTrans_max = mol(imol)%iTrans_max
@@ -444,23 +446,23 @@ subroutine clean_mem_dust_mol()
 
 end subroutine clean_mem_dust_mol
 
-subroutine realloc_dust_atom()
-!this routine should be the mirror of the mol one, except for the tab_lambda which is allocated
-!when the gas atom RT grid is defined (from reading the atomic species).
-!Still, the test on the allocation and the call of clean_mem_dust_mol in init_dust_atom means
-!that theya re not deallocated after temperature calculation. What am I missing ? 
+!******************************************************************************
 
-!   use stars, only : allocate_stellar_spectra !-> not use yet in atom transfer.
+subroutine realloc_dust_atom()
+  ! This routine should be the mirror of the mol one, except for the tab_lambda which is allocated
+  ! when the gas atom RT grid is defined (from reading the atomic species).
+  ! Still, the test on the allocation and the call of clean_mem_dust_mol in init_dust_atom means
+  ! that they are not deallocated after temperature calculation. What am I missing ?
 
+  !   use stars, only : allocate_stellar_spectra !-> not use yet in atom transfer.
 
   integer :: alloc_status
-  real :: mem_size
 
   if (lvariable_dust) then
      write(*,*) " WARNING: sizes of dust transfer could be very big !"
      !TO DO: better storing of quantities / recuction of n_lambda
   endif
-	
+
   !Note: tab_lambda(n_lambda) is already allocated in atomic_transfer
   !	the tab_lambda_* or tab_delta_lambda should be de-allocated when tab_lambda is allocated in atom_rt.
   allocate(tab_lambda_inf(n_lambda), tab_lambda_sup(n_lambda), tab_delta_lambda(n_lambda), &
@@ -521,8 +523,6 @@ end subroutine realloc_dust_atom
 
 !******************************************************************************
 
-!******************************************************************************
-
 subroutine realloc_step2()
 ! Ajout du cas ou les matrices de Mueller sont donnees en entrees
 ! 20/04/2023
@@ -617,7 +617,7 @@ subroutine realloc_step2()
      allocate(tab_mueller(4,4,0:nang_scatt,n_grains_tot,n_lambda2), stat=alloc_status)
      if (alloc_status > 0) call error('Allocation error tab_mueller')
      tab_mueller = 0
-  else 
+  else
      deallocate(tab_s11)
      allocate(tab_s11(0:nang_scatt,n_grains_tot,n_lambda2), stat=alloc_status)
      if (alloc_status > 0) call error('Allocation error tab_s11')
@@ -707,7 +707,7 @@ end subroutine realloc_step2
 subroutine realloc_ray_tracing_scattering_matrix()
 ! Ajout du cas ou les matrices de Mueller sont donnees en entrees
 ! 20/04/2023
- 
+
   integer, parameter :: p_n_lambda2_pos = 1
   integer :: alloc_status
 

src/output.f90

@@ -2392,17 +2392,23 @@ subroutine ecriture_UV_field()
   bitpix=-32
   extend=.true.
 
-  if (l3D) then
-     naxis=3
-     naxes(1)=n_rad
-     naxes(2)=2*nz
-     naxes(3)=n_az
-     nelements=naxes(1)*naxes(2)*naxes(3)
+  if (lVoronoi) then
+     naxis=1
+     naxes(1)=n_cells
+     nelements=naxes(1)
   else
-     naxis=2
-     naxes(1)=n_rad
-     naxes(2)=nz
-     nelements=naxes(1)*naxes(2)
+     if (l3D) then
+        naxis=3
+        naxes(1)=n_rad
+        naxes(2)=2*nz
+        naxes(3)=n_az
+        nelements=naxes(1)*naxes(2)*naxes(3)
+     else
+        naxis=2
+        naxes(1)=n_rad
+        naxes(2)=nz
+        nelements=naxes(1)*naxes(2)
+     endif
   endif
 
   !  Write the required header keywords.
@@ -2434,16 +2440,20 @@ end subroutine ecriture_UV_field
 
 function compute_UV_field() result(G)
 
-  real, dimension(n_cells) :: G
+  real, dimension(:), allocatable :: G
 
   integer :: lambda, l, icell
   integer, parameter :: n=200
 
   real(kind=dp) :: n_photons_envoyes, energie_photon
-  real(kind=dp), dimension(n_lambda,n_cells) :: J
+  real(kind=dp), dimension(:,:), allocatable :: J ! n_lambda,n_cells
   real(kind=dp), dimension(n) :: wl, J_interp
   real(kind=dp), dimension(n_lambda) :: lamb
   real(kind=dp) :: delta_wl
+  integer :: alloc_status
+
+  allocate(J(n_lambda,n_cells), G(n_cells), stat=alloc_status)
+  if (alloc_status /= 0) call error("allocation error in compute_UV_field")
 
   do lambda=1, n_lambda
      n_photons_envoyes = sum(n_phot_envoyes(lambda,:))

src/radiation_field.f90

@@ -149,7 +149,8 @@ subroutine allocate_radiation_field_step1(Nc)
 
   lxJ_abs = lProDiMo.or.lML.or.loutput_UV_field.or.loutput_J
   lxJ_abs_step1 = lRE_nLTE .or. lnRE .or. loutput_J_step1
-  if (lxJ_abs_step1 .or. (lxJ_abs.and.lsed.and.lsed_complete)) then
+
+  if (lxJ_abs_step1 .or. (lxJ_abs.and.lsed)) then
      allocate(xJ_abs(Nc,n_lambda,nb_proc), J0(Nc,n_lambda), stat=alloc_status) ! BIG array
      if (alloc_status > 0) call error('Allocation error xJ_abs')
      xJ_abs=0.0 ; J0 = 0.0

src/read_phantom.f90

@@ -936,8 +936,6 @@ subroutine phantom_2_mcfost(np,nptmass,ntypes,ndusttypes,ldust_moments,n_files,d
                 vz(j) = vzi * uvelocity
              endif
 
-
-
              if (ldust_moments) dust_moments(:,j) = nucleation(1:4,i) ! indexing is different from phantom as I read starting at k0
 
              T_gas(j) = T_gasi

src/thermal_emission.f90

@@ -114,8 +114,6 @@ subroutine allocate_thermal_emission(Nc,p_Nc)
   if (alloc_status > 0) call error('Allocation error kappa_abs_1grain')
   DensE = 0.0 ; DensE_m1 = 0.0
 
-  call allocate_radiation_field_step1(Nc)
-
   allocate(xT_ech(Nc,nb_proc), stat=alloc_status)
   if (alloc_status > 0) call error('Allocation error xT_ech')
   xT_ech = 2

src/wavelengths.f90

@@ -1,7 +1,7 @@
 module wavelengths
 
   use mcfost_env, only : dp, band
-  use parametres, only : n_pop
+  use parametres
   use messages
 
   implicit none
@@ -64,6 +64,10 @@ subroutine init_lambda()
 
   endif
 
+  if (loutput_UV_field .and. lsed_complete .and. (tab_lambda_inf(1) > 0.0912)) then ! 0.2
+     call warning("Shortest wavelengths is larger than 0.0912mum, UV field won't be properly estimated !!!")
+  endif
+
 end subroutine init_lambda
 
 !**********************************************************************
@@ -83,6 +87,10 @@ subroutine init_lambda2()
      tab_delta_lambda(i)= tab_delta_lambda2(i)
   enddo
 
+  if (loutput_UV_field .and. (tab_lambda_inf(1) > 0.0912)) then ! 0.2
+     call warning("Shortest wavelengths is larger than 0.0912mum, UV field won;t be properly estimated !!!")
+  endif
+
   return
 
 end subroutine init_lambda2