diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cf52075721..91cd57ae29 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2341,17 +2341,20 @@ contains if (run_CE) then ! By default, isotropic-in-lab is not used - allocate(mat % p0(n)) - mat % p0(:) = .false. + if (list_iso_lab % size() > 0) then + mat % has_isotropic_nuclides = .true. + allocate(mat % p0(n)) + mat % p0(:) = .false. - ! Apply isotropic-in-lab treatment to specified nuclides - do j = 1, list_iso_lab % size() - do k = 1, n - if (names % data(k) == list_iso_lab % data(j)) then - mat % p0(k) = .true. - end if + ! Apply isotropic-in-lab treatment to specified nuclides + do j = 1, list_iso_lab % size() + do k = 1, n + if (names % data(k) == list_iso_lab % data(j)) then + mat % p0(k) = .true. + end if + end do end do - end do + end if end if ! Check to make sure either all atom percents or all weight percents are diff --git a/src/material_header.F90 b/src/material_header.F90 index eda0b42302..a3c14f8cab 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -57,7 +57,8 @@ module material_header logical :: fissionable = .false. logical :: depletable = .false. - ! enforce isotropic scattering in lab + ! enforce isotropic scattering in lab for specific nuclides + logical :: has_isotropic_nuclides = .false. logical, allocatable :: p0(:) contains @@ -302,7 +303,6 @@ contains real(8) :: awr integer, allocatable :: new_nuclide(:) real(8), allocatable :: new_density(:) - logical, allocatable :: new_p0(:) character(:), allocatable :: name_ name_ = to_f_string(name) @@ -339,11 +339,6 @@ contains if (n > 0) new_density(1:n) = m % atom_density call move_alloc(FROM=new_density, TO=m % atom_density) - allocate(new_p0(n + 1)) - if (n > 0) new_p0(1:n) = m % p0 - new_p0(n + 1) = .false. - call move_alloc(FROM=new_p0, TO=m % p0) - ! Append new nuclide/density k = nuclide_dict % get(to_lower(name_)) m % nuclide(n + 1) = k diff --git a/src/physics.F90 b/src/physics.F90 index cb309f93ac..fd6a0087c9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -419,19 +419,22 @@ contains ! Set event component p % event = EVENT_SCATTER - ! Sample new outgoing angle for isotropic in lab scattering - if (materials(p % material) % p0(i_nuc_mat)) then + ! Sample new outgoing angle for isotropic-in-lab scattering + associate (mat => materials(p % material)) + if (mat % has_isotropic_nuclides) then + if (materials(p % material) % p0(i_nuc_mat)) then + ! Sample isotropic-in-lab outgoing direction + uvw_new(1) = TWO * prn() - ONE + phi = TWO * PI * prn() + uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) + uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) + p % mu = dot_product(uvw_old, uvw_new) - ! Sample isotropic-in-lab outgoing direction - uvw_new(1) = TWO * prn() - ONE - phi = TWO * PI * prn() - uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) - uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) - p % mu = dot_product(uvw_old, uvw_new) - - ! Change direction of particle - p % coord(1) % uvw = uvw_new - end if + ! Change direction of particle + p % coord(1) % uvw = uvw_new + end if + end if + end associate end subroutine scatter