From db4a42d15c9f5f9c5da7a4c143f7a65824f3f9e8 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Sat, 28 Sep 2024 21:09:44 +1000 Subject: [PATCH 1/2] Renaming icell_ref to icell1 --- src/SPH2mcfost.f90 | 2 +- src/benchmarks.f90 | 6 ++--- src/cylindrical_grid.f90 | 2 -- src/diffusion.f90 | 10 ++++----- src/disk_physics.f90 | 2 +- src/dust_prop.f90 | 26 ++++++++++++---------- src/dust_ray_tracing.f90 | 27 +++++++++++------------ src/dust_transfer.f90 | 6 ++--- src/gas/atom_transfer.f90 | 22 +++++++++---------- src/gas/escape.f90 | 46 +++++++++++++++++++-------------------- src/mol_transfer.f90 | 2 +- src/optical_depth.f90 | 32 +++++++++++++-------------- src/parameters.f90 | 2 +- src/radiation_field.f90 | 2 +- src/thermal_emission.f90 | 8 +++---- 15 files changed, 97 insertions(+), 98 deletions(-) diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index ab085c764..86b6360e5 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -244,7 +244,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g limit_threshold = (1.0 - SPH_keep_particles) * 0.5 ; - icell_ref = 1 + icell1 = 1 write(*,*) "# Farthest particules :" write(*,*) "x =", minval(x), maxval(x) diff --git a/src/benchmarks.f90 b/src/benchmarks.f90 index bf69b62c1..9efef7162 100644 --- a/src/benchmarks.f90 +++ b/src/benchmarks.f90 @@ -263,7 +263,7 @@ subroutine init_GG_Tau_mol() Tcin(icell) = 30.0 * (r_grid(icell)/100.)**(-0.5) enddo - icell = icell_ref + icell = icell1 write(*,*) "Density @ 100 AU", real(densite_gaz(icell) / 100.**3 * (sqrt(r_grid(icell)**2 + z_grid(icell)) / 100.)**2.75) return @@ -311,8 +311,8 @@ subroutine init_benchmark_vanZadelhoff1() !tab_abundance = abundance - write(*,*) "Density", real(densite_gaz(icell_ref) * & - (r_grid(icell_ref)**2 + z_grid(icell_ref)**2)/rmin**2) + write(*,*) "Density", real(densite_gaz(icell1) * & + (r_grid(icell1)**2 + z_grid(icell1)**2)/rmin**2) return diff --git a/src/cylindrical_grid.f90 b/src/cylindrical_grid.f90 index 76eea7669..6ae8f79ca 100644 --- a/src/cylindrical_grid.f90 +++ b/src/cylindrical_grid.f90 @@ -47,8 +47,6 @@ subroutine build_cylindrical_cell_mapping() ! work also in spherical integer :: i,j,k,icell, ntot, ntot2, alloc_status integer :: istart,iend,jstart,jend,kstart,kend, istart2,iend2,jstart2,jend2,kstart2,kend2 - icell_ref = 1 - istart = 1 iend = n_rad diff --git a/src/diffusion.f90 b/src/diffusion.f90 index d178dd099..3021dde9d 100644 --- a/src/diffusion.f90 +++ b/src/diffusion.f90 @@ -31,7 +31,7 @@ subroutine setDiffusion_coeff(i) cst_Dcoeff = pi/(12.*sigma) - p_icell = icell_ref + p_icell = icell1 do k=1,n_az do j=j_start, nz @@ -90,7 +90,7 @@ subroutine setDiffusion_coeff0(i) cst_Dcoeff = pi/(12.*sigma) - p_icell = icell_ref + p_icell = icell1 do k=1,n_az do j=j_start,nz @@ -662,7 +662,7 @@ subroutine compute_Planck_opacities(icell, Planck_opacity,rec_Planck_opacity) if (lvariable_dust) then p_icell => icell0 else - p_icell => icell_ref + p_icell => icell1 endif somme = 0.0_dp @@ -686,8 +686,8 @@ subroutine compute_Planck_opacities(icell, Planck_opacity,rec_Planck_opacity) somme2 = somme2 + B * kappa(p_icell,lambda) norm = norm + B*delta_wl enddo - rec_Planck_opacity = norm/somme * kappa_factor(icell0) - Planck_opacity = somme2/norm * kappa_factor(icell0) + rec_Planck_opacity = norm/somme * kappa_factor(icell) + Planck_opacity = somme2/norm * kappa_factor(icell) else rec_Planck_opacity = 1e-30 Planck_opacity = 1e-30 diff --git a/src/disk_physics.f90 b/src/disk_physics.f90 index 8aec21a69..2ec9db25c 100644 --- a/src/disk_physics.f90 +++ b/src/disk_physics.f90 @@ -31,7 +31,7 @@ subroutine compute_othin_sublimation_radius() cst=cst_th/dust_pop(1)%T_sub cst=cst_th/1500. - icell = icell_ref + icell = icell1 do lambda=1, n_lambda ! longueur d'onde en metre diff --git a/src/dust_prop.f90 b/src/dust_prop.f90 index 58eac21b4..239ddce74 100644 --- a/src/dust_prop.f90 +++ b/src/dust_prop.f90 @@ -811,10 +811,9 @@ subroutine opacite(lambda, p_lambda, no_scatt) ! c'est pour les prop de diffusion en relatif donc la veleur exacte n'a pas d'importante ldens0 = .false. if (.not.lvariable_dust) then - icell = icell_ref - if (maxval(densite_pouss(:,icell)) < tiny_real) then + if (icell_not_empty > icell1) then ! icell1==1 ldens0 = .true. - densite_pouss(:,icell) = densite_pouss(:,icell_not_empty) + densite_pouss(:,icell1) = densite_pouss(:,icell_not_empty) endif endif @@ -1010,8 +1009,7 @@ subroutine opacite(lambda, p_lambda, no_scatt) ! On remet la densite à zéro si besoin if (ldens0) then - icell = icell_ref - densite_pouss(:,icell) = 0.0_sp + densite_pouss(:,icell1) = 0.0_sp endif if ((ldust_prop).and.(lambda == n_lambda)) then @@ -1055,10 +1053,9 @@ subroutine calc_local_scattering_matrices(lambda, p_lambda) ! c'est pour les prop de diffusion en relatif donc la veleur exacte n'a pas d'importante ldens0 = .false. if (.not.lvariable_dust) then - icell = icell_ref - if (maxval(densite_pouss(:,icell)) < tiny_real) then + if (icell_not_empty > icell1) then ldens0 = .true. - densite_pouss(:,icell) = densite_pouss(:,icell_not_empty) + densite_pouss(:,icell1) = densite_pouss(:,icell_not_empty) endif endif @@ -1216,8 +1213,7 @@ subroutine calc_local_scattering_matrices(lambda, p_lambda) ! see opacite() ! On remet la densite à zéro si besoin if (ldens0) then - icell = icell_ref - densite_pouss(:,icell) = 0.0_sp + densite_pouss(:,icell1) = 0.0_sp endif return @@ -1335,8 +1331,14 @@ subroutine write_dust_prop() allocate(S11_lambda_theta(n_lambda,0:nang_scatt),pol_lambda_theta(n_lambda,0:nang_scatt)) allocate(kappa_grain(n_lambda,n_grains_tot)) - icell = icell_not_empty - p_icell = icell_ref + if (lvariable_dust) then + write(*,*) "Warning: dust is not uniform, picking cell #", icell_not_empty + icell = icell_not_empty + p_icell = icell_not_empty + else + icell = icell_not_empty + p_icell = 1 + endif kappa_lambda=real((kappa(icell,:)*kappa_factor(icell)/AU_to_cm)/(masse(icell)/(volume(icell)*AU_to_cm**3))) ! cm^2/g albedo_lambda=tab_albedo_pos(icell,:) diff --git a/src/dust_ray_tracing.f90 b/src/dust_ray_tracing.f90 index 0489af6e8..b46c2abd0 100644 --- a/src/dust_ray_tracing.f90 +++ b/src/dust_ray_tracing.f90 @@ -503,7 +503,7 @@ subroutine calc_xI_scatt(id,lambda,p_lambda,icell, phik,psup,l,stokes,flag_star) if (lvariable_dust) then p_icell = icell else - p_icell = icell_ref + p_icell = icell1 endif do ibin = 1, RT_n_incl @@ -565,7 +565,7 @@ subroutine calc_xI_scatt_pola(id,lambda,p_lambda,icell,phik,psup,l,stokes,flag_s if (lvariable_dust) then p_icell = icell else - p_icell = icell_ref + p_icell = icell1 endif do ibin = 1, RT_n_incl @@ -670,9 +670,9 @@ subroutine init_dust_source_fct1(lambda,ibin,iaz) !$omp default(none) & !$omp private(icell,p_icell,facteur,kappa_sca,kappa_ext,itype,I_scatt) & !$omp shared(n_cells,energie_photon,volume,n_az_rt,n_theta_rt,kappa,kappa_factor,N_type_flux,lambda,iRT) & - !$omp shared(J_th,xI_scatt,eps_dust1,lsepar_pola,lsepar_contrib,n_Stokes,nb_proc,tab_albedo_pos,lvariable_dust,icell_ref) + !$omp shared(J_th,xI_scatt,eps_dust1,lsepar_pola,lsepar_contrib,n_Stokes,nb_proc,tab_albedo_pos,lvariable_dust,icell1) - p_icell = icell_ref + p_icell = icell1 !$omp do schedule(static,n_cells/nb_proc) do icell=1, n_cells @@ -750,9 +750,9 @@ subroutine init_dust_source_fct2(lambda,p_lambda,ibin) !$omp parallel & !$omp default(none) & !$omp shared(lambda,n_cells,nang_ray_tracing,eps_dust2,I_sca2,J_th,kappa,kappa_factor,eps_dust2_star,lsepar_pola) & - !$omp shared(lsepar_contrib,n_Stokes,nang_ray_tracing_star,nb_proc,tab_albedo_pos,icell_ref,lvariable_dust) & + !$omp shared(lsepar_contrib,n_Stokes,nang_ray_tracing_star,nb_proc,tab_albedo_pos,icell1,lvariable_dust) & !$omp private(icell,p_icell,dir,iscatt,Q,U,kappa_ext) - p_icell = icell_ref + p_icell = icell1 !$omp do schedule(static,n_cells/nb_proc) do icell=1, n_cells if (lvariable_dust) p_icell = icell @@ -829,9 +829,9 @@ subroutine calc_Jth(lambda) !$omp parallel & !$omp default(none) & !$omp shared(n_cells,Tdust,wl,lambda,kappa_abs_LTE,kappa_factor,cst_E,J_th) & - !$omp shared(lvariable_dust,icell_ref) & + !$omp shared(lvariable_dust,icell1) & !$omp private(icell,p_icell,Temp,cst_wl,coeff_exp) - p_icell = icell_ref + p_icell = icell1 !$omp do do icell=1, n_cells if (lvariable_dust) p_icell = icell @@ -1076,7 +1076,7 @@ subroutine calc_Isca_rt2(lambda,p_lambda,ibin) !$omp parallel & !$omp default(none) & - !$omp shared(lvariable_dust,Inu,I_sca2,n_cells,tab_s11_pos,uv0,w0,n_Stokes,icell_ref,energie_photon,volume) & + !$omp shared(lvariable_dust,Inu,I_sca2,n_cells,tab_s11_pos,uv0,w0,n_Stokes,icell1,energie_photon,volume) & !$omp shared(tab_s12_o_s11_pos,tab_s22_o_s11_pos,tab_s33_o_s11_pos,tab_s34_o_s11_pos,tab_s44_o_s11_pos) & !$omp shared(lsepar_pola,tab_k,tab_sin_scatt_norm,lambda,p_lambda,n_phi_I,n_theta_I,nang_ray_tracing,lsepar_contrib) & !$omp shared(s11_save,tab_cosw,tab_sinw,nb_proc,kappa,kappa_factor,tab_albedo_pos) & @@ -1098,8 +1098,7 @@ subroutine calc_Isca_rt2(lambda,p_lambda,ibin) ROP(4,4) = 1.0_dp if (.not.lvariable_dust) then ! we precalculate the s11 as they are all the same - icell = icell_ref - + icell = icell1 do dir=0,1 do iscatt = 1, nang_ray_tracing ! Loop over the specific intensity bins @@ -1122,7 +1121,7 @@ subroutine calc_Isca_rt2(lambda,p_lambda,ibin) enddo ! dir endif ! .not.lvariable_dust - p_icell = icell_ref + p_icell = icell1 !$omp do schedule(static, n_cells/nb_proc) do icell = 1, n_cells !$ id = omp_get_thread_num() + 1 @@ -1269,7 +1268,7 @@ subroutine calc_Isca_rt2_star(lambda,p_lambda,ibin) !$omp parallel & !$omp default(none) & !$omp shared(n_cells,lambda,p_lambda,ibin,I_spec_star,nang_ray_tracing_star,cos_thet_ray_tracing_star) & - !$omp shared(lvariable_dust,r_grid,z_grid,icell_ref,omega_ray_tracing_star,lsepar_pola) & + !$omp shared(lvariable_dust,r_grid,z_grid,icell1,omega_ray_tracing_star,lsepar_pola) & !$omp shared(tab_s11_pos,tab_s12_o_s11_pos,tab_s22_o_s11_pos,tab_s33_o_s11_pos,tab_s34_o_s11_pos,tab_s44_o_s11_pos) & !$omp shared(energie_photon,volume,tab_albedo_pos,kappa,kappa_factor,eps_dust2_star) & !$omp private(id,icell,p_icell,stokes,x,y,z,norme,u,v,w,dir,iscatt,cos_scatt,k,omega,cosw,sinw,RPO,ROP,S) & @@ -1290,7 +1289,7 @@ subroutine calc_Isca_rt2_star(lambda,p_lambda,ibin) stokes(:) = 0.0_dp - p_icell = icell_ref + p_icell = icell1 !$omp do do icell=1, n_cells diff --git a/src/dust_transfer.f90 b/src/dust_transfer.f90 index 529490ffa..bd43af35f 100644 --- a/src/dust_transfer.f90 +++ b/src/dust_transfer.f90 @@ -212,8 +212,8 @@ subroutine transfert_poussiere() if (loptical_depth_to_cell) call write_optical_depth_to_cell(1) write(*,*) "" - write(*,*) "Dust properties in cell #", icell_ref - p_icell = icell_ref + write(*,*) "Dust properties in cell #", icell_not_empty + p_icell = icell_not_empty if (aniso_method==2) write(*,*) "g ", tab_g_pos(p_icell,1) write(*,*) "albedo ", tab_albedo_pos(p_icell,1) if (lsepar_pola.and.(scattering_method == 2)) write(*,*) "polarisability", maxval(-tab_s12_o_s11_pos(:,p_icell,1)) @@ -998,7 +998,7 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif ! Boucle sur les interactions du paquets: diff --git a/src/gas/atom_transfer.f90 b/src/gas/atom_transfer.f90 index 9881052ce..bbc691127 100644 --- a/src/gas/atom_transfer.f90 +++ b/src/gas/atom_transfer.f90 @@ -122,7 +122,7 @@ end subroutine io_write_convergence_maps subroutine nlte_loop_mali() ! ------------------------------------------------------------------------------------ ! - ! Solve the set of statistical equilibrium equations (SEE) with the + ! Solve the set of statistical equilibrium equations (SEE) with the ! Multi-level Accelerated Lambda Iterations method (Rybicki & Hummer 92, Apj 202 209). ! ! By default, the SEE are solved twice. @@ -150,7 +150,7 @@ subroutine nlte_loop_mali() logical :: l_iterate, l_iterate_ne integer :: n_xc real(kind=dp), dimension(:), pointer :: tab_xc - + !Ng's acceleration integer, parameter :: Ng_Ndelay_init = 5 !minimal number of iterations before starting the cycle. !min is 1 (so that conv_speed can be evaluated) @@ -159,7 +159,7 @@ subroutine nlte_loop_mali() real(kind=dp) :: conv_speed_limit integer :: iorder, i0_rest, n_iter_accel, iacc integer :: Ng_Ndelay, ng_index - logical :: lng_turned_on, ng_rest, lconverging, accelerated + logical :: lng_turned_on, ng_rest, lconverging, accelerated logical, parameter :: lextrapolate_electron = .false. !extrapolate electrons with populations real(kind=dp), dimension(:,:), allocatable :: ngtmp @@ -329,7 +329,7 @@ subroutine nlte_loop_mali() ng_index = Neq_ng - mod(n_iter-1,Neq_ng) !overwritten if lng_acceleration. ! NOTE :: the index 1 is hard-coded within the non-LTE loop as it stores the actual (running) value ! of the populations and electronic density. - ! + ! ! goes with maxIter write(*,'(" *** Iteration #"(1I4)"; step #"(1I1)"; threshold: "(1ES11.2E3)"; Nrays: "(1I5))') & n_iter, etape, precision, n_rayons @@ -522,7 +522,7 @@ subroutine nlte_loop_mali() ! reshape(ngpop(1:at%Nlevel,nact,:,:),shape=[n_cells,at%Nlevel,Neq_ng],order=[2,1,3]),&!then flatten ! (/n_cells*at%Nlevel,Neq_ng/)) - call Accelerate(n_cells*at%Nlevel,Ng_Norder,ngtmp) + call Accelerate(n_cells*at%Nlevel,Ng_Norder,ngtmp) !reform in (Nlevel,N_cells,Neq_ng) ngpop(1:at%Nlevel,nact,:,:) = reshape(ngtmp,(/at%Nlevel,n_cells,Neq_ng/)) @@ -1175,7 +1175,7 @@ subroutine atom_line_transfer() if (lany_init4) then call nlte_loop_sobolev() - ! MALI step after Sobolev step ? + ! MALI step after Sobolev step ? ! if (.not.lescape_prob) then ! !to preserve ray-resolution ? ! istep_start = 2 @@ -1220,7 +1220,7 @@ subroutine atom_line_transfer() if (laccretion_shock) then !3 + lg(Facc) = lg(Facc) erg/cm2/s - if (max_Facc>0) write(*,'("max(lgFacc)="(1F7.3)" W/m^2; min(lgFacc)="(1F7.3)" W/m^2")') & + if (max_Facc>0) write(*,'("max(lgFacc)="(1F7.3)" W/m^2; min(lgFacc)="(1F7.3)" W/m^2")') & log10(max_Facc), log10(min_Facc) if (max_Tshock>0) write(*,'("max(Tshock)="(1I8)" K; min(Tshock)="(1I8)" K")') & nint(max_Tshock), nint(min_Tshock) @@ -1245,7 +1245,7 @@ subroutine atom_line_transfer() end subroutine atom_line_transfer subroutine init_dust_temperature() - !lowering too much the treshold might create some convergence issues in the non-LTE pops or in + !lowering too much the treshold might create some convergence issues in the non-LTE pops or in ! the electronic density calculations (0 division mainly for low values). real(kind=dp), parameter :: T_formation = 2000.0 ! [K] logical, dimension(:), allocatable :: ldust @@ -1261,7 +1261,7 @@ subroutine init_dust_temperature() where(ldust) T(:) = Tdust(:) if (any(T < T_formation)) call warning('(atom_transfer) Setting the gas transparent where T < 2000 K.') where (T < T_formation) icompute_atomRT = 0 - !or set the atomic gas temperature to 0 in dust regions ? + !or set the atomic gas temperature to 0 in dust regions ? if (any(icompute_atomRT>0)) then write(*,*) "T:", real(maxval(T,icompute_atomRT>0)), real(minval(T,icompute_atomRT>0)) write(*,*) "T (rho_dust = 0):", real(maxval(T,(icompute_atomRT>0).and..not.ldust)), & @@ -1297,7 +1297,7 @@ subroutine init_dust_atom() if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif call realloc_dust_atom() @@ -1692,7 +1692,7 @@ subroutine emission_line_tau_surface_map(tau,ibin,iaz) do i = 1,npix_x !$ id = omp_get_thread_num() + 1 do j = 1,npix_y - + pixelcenter(:,id) = Icorner(:) + (i-0.5_dp) * dx(:) + (j-0.5_dp) * dy(:) x0(:) = pixelcenter(1,id) diff --git a/src/gas/escape.f90 b/src/gas/escape.f90 index cd7eedae0..e20da9f68 100644 --- a/src/gas/escape.f90 +++ b/src/gas/escape.f90 @@ -126,11 +126,11 @@ subroutine count_neighbours(n) integer, parameter :: order = 1 if ((lcylindrical)) then - call warning("(count_neighbours) Not tested for cylindrical grid yet") + call warning("(count_neighbours) Not tested for cylindrical grid yet") endif if ((lvoronoi)) then - call error("(count_neighbours) Not working for voronoi grid yet") + call error("(count_neighbours) Not working for voronoi grid yet") endif !spherical grid @@ -310,7 +310,7 @@ subroutine mean_velocity_gradient_faster() call progress_bar(ibar) !$omp atomic ibar = ibar+1 - endif + endif enddo !$omp end do !$omp end parallel @@ -347,7 +347,7 @@ subroutine mean_velocity_gradient_faster() ! if (.not.lintersect) cycle ! if (lintersect) then n_rays_star(i_star) = n_rays_star(i_star) + 1 - if (is_inshock(id, 1, i_star, icell, x0, y0, z0, Tchoc, F1, T1)) then + if (is_inshock(id, 1, i_star, icell, x0, y0, z0, Tchoc, F1, T1)) then !fraction actually because all rays touch the star here Tchoc_average(i_star) = Tchoc_average(i_star) + Tchoc * nHtot(icell) n_rays_shock(i_star) = n_rays_shock(i_star) + 1 @@ -503,7 +503,7 @@ subroutine mean_velocity_gradient() call intersect_stars(x0,y0,z0, u,v,w, lintersect_stars, i_star, icell_star)!will intersect - !to do propagate along the ray to get dv_max + !to do propagate along the ray to get dv_max previous_cell = 0 ! unused, just for Voronoi call cross_cell(x0,y0,z0, u,v,w,icell, previous_cell, x1,y1,z1, next_cell,l, l_contrib, l_void_before) @@ -526,7 +526,7 @@ subroutine mean_velocity_gradient() inf : do if (next_cell==icell_star) then !in shock ?? - if (is_inshock(id, iray, i_star, icell_in, xa, xb, xc, Tchoc, F1, T1)) then + if (is_inshock(id, iray, i_star, icell_in, xa, xb, xc, Tchoc, F1, T1)) then dOmega_shock(icell,i_star) = dOmega_shock(icell,i_star) + wei Tchoc_average(i_star) = Tchoc_average(i_star) + Tchoc * nHtot(icell_in) n_rays_shock(i_star) = n_rays_shock(i_star) + 1 @@ -561,7 +561,7 @@ subroutine mean_velocity_gradient() call progress_bar(ibar) !$omp atomic ibar = ibar+1 - endif + endif enddo !$omp end do !$omp end parallel @@ -573,7 +573,7 @@ subroutine mean_velocity_gradient() Tchoc_average = Tchoc_average / rho_shock f_shock = real(n_rays_shock) / real(n_rays_star) endif - + write(*,'("max()="(1ES17.8E3)" s^-1; min()="(1ES17.8E3)" s^-1")') & maxval(mean_grad_v), minval(mean_grad_v,icompute_atomRT>0) write(*,'("max()="(1ES17.8E3)" m; min()="(1ES17.8E3)" m")') & @@ -636,7 +636,7 @@ subroutine nlte_loop_sobolev() logical :: l_iterate, l_iterate_ne integer :: n_xc real(kind=dp), dimension(:), pointer :: tab_xc - + !convergence check type (AtomType), pointer :: at integer(kind=8) :: mem_alloc_local @@ -806,7 +806,7 @@ subroutine nlte_loop_sobolev() stream = 0.0 stream(:) = [(init_sprng(gtype, i-1,nb_proc,seed,SPRNG_DEFAULT),i=1,nb_proc)] - + l_iterate_ne = .false. if( n_iterate_ne > 0 ) then l_iterate_ne = ( mod(n_iter,n_iterate_ne)==0 ) .and. (n_iter>Ndelay_iterate_ne) @@ -942,7 +942,7 @@ subroutine nlte_loop_sobolev() write(*,*) " *** stopping electronic density convergence at iteration ", n_iter n_iterate_ne = 0 endif - endif + endif end if !***********************************************************! @@ -1022,9 +1022,9 @@ subroutine nlte_loop_sobolev() !TO DO: takes time with large grid ! call opacity_atom_bf_loc(icell, n_xc, tab_xc, chi_cont(:,icell), eta_cont(:,icell)) !background continua not needed if no continuum RT - if (limit_mem < 2) then - call calc_contopac_loc(icell,n_xc,tab_xc) - endif + if (limit_mem < 2) then + call calc_contopac_loc(icell,n_xc,tab_xc) + endif end if !if l_iterate end do cell_loop2 !icell @@ -1219,12 +1219,12 @@ subroutine radrates_sobolev_average(id, icell) !pops inversions if (at%n(i,icell)-at%lines(kr)%gij*at%n(j,icell) <= 0.0) cycle atr_loop - + !assumes the underlying radiation is flat across the line. !or take n0 as a function of velocity Icore = Itot(n0,1,id) - + !could be stored in mem. chi_ij = hc_fourPI * at%lines(kr)%Bij * (at%n(i,icell)-at%lines(kr)%gij*at%n(j,icell)) tau_escape = fact_tau * chi_ij * mean_length_scale(icell) / vth @@ -1352,7 +1352,7 @@ function I_sobolev_1ray(id,icell_in,x,y,z,u,v,w,iray,N,lambda) I_sobolev_1ray(:) = 0.0_dp - p_icell => icell_ref + p_icell => icell1 if (lvariable_dust) p_icell => icell ! Will the ray intersect a star @@ -1376,7 +1376,7 @@ function I_sobolev_1ray(id,icell_in,x,y,z,u,v,w,iray,N,lambda) do nat=1, NactiveAtoms do kr=1, ActiveAtoms(nat)%p%nline nl = ActiveAtoms(nat)%p%lines(kr)%Nr - ActiveAtoms(nat)%p%lines(kr)%Nb + 1 - I0_line(1:nl,kr,nat,id) = Icore(ActiveAtoms(nat)%p%lines(kr)%Nb:ActiveAtoms(nat)%p%lines(kr)%Nr) + I0_line(1:nl,kr,nat,id) = Icore(ActiveAtoms(nat)%p%lines(kr)%Nb:ActiveAtoms(nat)%p%lines(kr)%Nr) enddo enddo return @@ -1468,11 +1468,11 @@ subroutine accumulate_radrates_sobolev_1ray(id, icell, iray, dOmega) !pops inversions if (at%n(i,icell)-at%lines(kr)%gij*at%n(j,icell) <= 0.0) cycle atr_loop - + !assumes the underlying radiation is flat across the line. !or take n0 as a function of velocity Icore = maxval(I0_line(:,kr,nact,id)) !shock + star weighted by the solid angle - + !could be stored in mem. chi_ij = hc_fourPI * at%lines(kr)%Bij * (at%n(i,icell)-at%lines(kr)%gij*at%n(j,icell)) tau_escape = fact_tau * chi_ij * ds(iray,id) / vth @@ -1596,7 +1596,7 @@ subroutine xccont(id,icell,iray,domega,rates) if (present(rates)) then if (.not.rates) return else !not present, return, equivalent to .false. - return + return endif !accumulate rates @@ -1700,7 +1700,7 @@ subroutine opt_thick_mali_rates_cont(id,icell,cont) Nb = cont%Nb; Nr = cont%Nr Nl = Nr-Nb + 1 i0 = Nb - 1 - + ! xcc_down = 0.0 ! do l=1, Nl ! if (l==1) then @@ -1749,4 +1749,4 @@ subroutine opt_thick_mali_rates_cont(id,icell,cont) return end subroutine opt_thick_mali_rates_cont -end module escape \ No newline at end of file +end module escape diff --git a/src/mol_transfer.f90 b/src/mol_transfer.f90 index f1d51976d..3bd49ceef 100644 --- a/src/mol_transfer.f90 +++ b/src/mol_transfer.f90 @@ -896,7 +896,7 @@ subroutine init_dust_mol(imol) if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif call realloc_dust_mol(imol) diff --git a/src/optical_depth.f90 b/src/optical_depth.f90 index e38224ada..ae2f53651 100644 --- a/src/optical_depth.f90 +++ b/src/optical_depth.f90 @@ -70,7 +70,7 @@ subroutine physical_length(id,lambda,p_lambda,Stokes,icell,xio,yio,zio,u,v,w,fla if (lvariable_dust) then p_icell => icell0 else - p_icell => icell_ref + p_icell => icell1 endif ! Boucle infinie sur les cellules @@ -215,9 +215,9 @@ subroutine integ_tau(lambda) if (.not.lvariable_dust) then icell = icell_not_empty - if (kappa(icell_ref,lambda) * kappa_factor(icell) > tiny_real) then + if (kappa(icell1,lambda) * kappa_factor(icell) > tiny_real) then write(*,*) " Column density (g/cm^2) = ", real(tau*(masse(icell)/(volume(icell)*AU_to_cm**3))/ & - (kappa(icell_ref,lambda) * kappa_factor(icell)/AU_to_cm)) + (kappa(icell1,lambda) * kappa_factor(icell)/AU_to_cm)) endif endif @@ -233,9 +233,9 @@ subroutine integ_tau(lambda) if (.not.lvariable_dust) then icell = icell_not_empty - if (kappa(icell_ref,lambda) * kappa_factor(icell) > tiny_real) then + if (kappa(icell1,lambda) * kappa_factor(icell) > tiny_real) then write(*,*) " Column density (g/cm^2) = ", real(tau*(masse(icell)/(volume(icell)*AU_to_cm**3))/ & - (kappa(icell_ref,lambda) * kappa_factor(icell)/AU_to_cm)) + (kappa(icell1,lambda) * kappa_factor(icell)/AU_to_cm)) endif endif @@ -283,7 +283,7 @@ subroutine optical_length_tot(id,lambda,Stokes,icell,xi,yi,zi,u,v,w,tau_tot_out, if (lvariable_dust) then p_icell => icell0 else - p_icell => icell_ref + p_icell => icell1 endif ! Boucle infinie sur les cellules @@ -350,9 +350,9 @@ subroutine compute_column(type, column, lambda) do direction = 1, n_directions !$omp parallel default(none) & !$omp shared(densite_gaz,tab_abundance,lVoronoi,Voronoi,direction,column,r_grid,z_grid,phi_grid,n_cells,cross_cell) & - !$omp shared(CD_units,kappa,kappa_factor,lambda,type,test_exit_grid,icell_ref,lvariable_dust) & + !$omp shared(CD_units,kappa,kappa_factor,lambda,type,test_exit_grid,icell1,lvariable_dust) & !$omp private(icell,previous_cell,next_cell,icell0,p_icell,x0,y0,z0,x1,y1,z1,norme,u,v,w,l,l_contrib,l_void_before,factor,sum) - p_icell = icell_ref + p_icell = icell1 !$omp do do icell=1,n_cells if (lVoronoi) then @@ -468,7 +468,7 @@ subroutine integ_ray_mol(id,imol,icell_in,x,y,z,u,v,w,iray,labs, ispeed,tab_spee if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif @@ -637,7 +637,7 @@ subroutine physical_length_mol(imol,iTrans,icell_in,x,y,z,u,v,w, ispeed, tab_spe if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif ! Boucle infinie sur les cellules @@ -766,7 +766,7 @@ subroutine physical_length_mol_Flux(imol,iTrans,icell_in,x,y,z,u,v,w, ispeed, ta if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif ! Boucle infinie sur les cellules @@ -1030,7 +1030,7 @@ subroutine optical_length_tot_mol(imol,icell_in,x,y,z,u,v,w, ispeed, tab_speed, if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif !*** propagation dans la grille @@ -1115,7 +1115,7 @@ subroutine integ_ray_atom(id,icell_in,x,y,z,u,v,w,iray,labs,N,lambda) Itot(:,iray,id) = 0.0_dp - p_icell => icell_ref + p_icell => icell1 if (lvariable_dust) p_icell => icell ! Will the ray intersect a star @@ -1244,7 +1244,7 @@ subroutine physical_length_atom(id,icell_in,x,y,z,u,v,w,N,lambda,tau_threshold,f !flag_sortie = .true. not useful. By default, (x,y,z) are 0 if we don't reach the surface. x = 0.0; y = 0.0; z = 0.0 - p_icell => icell_ref + p_icell => icell1 if (lvariable_dust) p_icell => icell ! Will the ray intersect a star @@ -1358,7 +1358,7 @@ function integ_ray_dust(lambda,icell_in,x,y,z,u,v,w) if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif ! Will the ray intersect a star @@ -1448,7 +1448,7 @@ subroutine define_dark_zone(lambda,p_lambda,tau_max,ldiff_approx) if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif do pk=1, n_az diff --git a/src/parameters.f90 b/src/parameters.f90 index e5839bcd5..4a82b14ce 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -124,7 +124,7 @@ module parametres integer :: n_az, j_start, pj_start ! Nombre de cellules totale integer :: n_cells, nrz, p_n_cells - integer, target :: icell_ref + integer, target :: icell1 = 1 logical :: lregular_theta real :: theta_max, theta_mask_max diff --git a/src/radiation_field.f90 b/src/radiation_field.f90 index ac91394b5..09b8ca70c 100644 --- a/src/radiation_field.f90 +++ b/src/radiation_field.f90 @@ -46,7 +46,7 @@ subroutine save_radiation_field(id,lambda,p_lambda,icell, Stokes, l, x0,y0,z0, if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif if (letape_th) then diff --git a/src/thermal_emission.f90 b/src/thermal_emission.f90 index b7d4bcdf7..9514ee554 100644 --- a/src/thermal_emission.f90 +++ b/src/thermal_emission.f90 @@ -456,7 +456,7 @@ subroutine init_reemission(lheating,dudt) !$omp parallel default(none) & !$omp private(id,icell,T,lambda,integ, Qcool,Qcool0,extra_heating,Qcool_minus_extra_heating,Temp,u_o_dt) & !$omp shared(cst_E,kappa_abs_LTE,kappa_factor,volume,B,lextra_heating,xT_ech,log_Qcool_minus_extra_heating,J0) & - !$omp shared(n_T,n_cells,p_n_cells,n_lambda,tab_Temp,ldudt_implicit,ufac_implicit,dudt,lRE_nLTE,lvariable_dust,icell_ref) + !$omp shared(n_T,n_cells,p_n_cells,n_lambda,tab_Temp,ldudt_implicit,ufac_implicit,dudt,lRE_nLTE,lvariable_dust,icell1) id = 1 ! Pour code sequentiel !$ id = omp_get_thread_num() + 1 @@ -661,7 +661,7 @@ subroutine Temp_LTE(id, icell, Ti, Temp, frac) if (lvariable_dust) then p_icell = icell else - p_icell = icell_ref + p_icell = icell1 endif if (id==0) then @@ -1785,7 +1785,7 @@ subroutine repartition_energie(lambda) if (lvariable_dust) then p_icell => icell else - p_icell => icell_ref + p_icell => icell1 endif cst_wl_max = log(huge_real)-1.0e-4 @@ -1960,7 +1960,7 @@ integer function select_absorbing_grain(lambda,icell, aleat, heating_method) res if (lvariable_dust) then p_icell = icell else - p_icell = icell_ref + p_icell = icell1 endif ! We scale the random number so that it is between 0 and kappa_sca (= last value of CDF) From ce556f22ef56f3aa068ed371daa959615b186fe0 Mon Sep 17 00:00:00 2001 From: Christophe Pinte Date: Sat, 28 Sep 2024 22:24:44 +1000 Subject: [PATCH 2/2] (phantom) Bug fix: using incorrect flag for masking particles --- src/SPH2mcfost.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SPH2mcfost.f90 b/src/SPH2mcfost.f90 index 86b6360e5..54df20aa3 100644 --- a/src/SPH2mcfost.f90 +++ b/src/SPH2mcfost.f90 @@ -575,7 +575,7 @@ subroutine SPH_to_Voronoi(n_SPH, ndusttypes, particle_id, x,y,z,h, vx,vy,vz, T_g if (present(mask)) then if (allocated(mask)) then do icell=1,n_cells - iSPH = Voronoi(icell)%id + iSPH = Voronoi(icell)%original_id if (iSPH > 0) then Voronoi(icell)%masked = mask(iSPH) else @@ -822,7 +822,7 @@ subroutine delete_masked_particles() k=0 do icell=1, n_cells - if (Voronoi(icell)%masked == 1) then + if (Voronoi(icell)%masked == 1) then ! 1 is transparent k=k+1 masse_gaz(icell) = 0. densite_gaz(icell) = 0.