cleaned up code

src/cmfd_data.F90

@@ -20,7 +20,8 @@ contains
 
     use cmfd_header,         only: allocate_cmfd
     use constants,           only: CMFD_NOACCEL
-    use global,              only: cmfd, cmfd_coremap, cmfd_downscatter
+    use global,              only: cmfd, cmfd_coremap, cmfd_downscatter, &
+                                   cmfd_fix_balance
 
     ! Check for core map and set it up
     if ((cmfd_coremap) .and. (cmfd%mat_dim == CMFD_NOACCEL)) call set_coremap()
@@ -57,12 +58,16 @@ contains
                             ONE, TINY_BIT
     use error,        only: fatal_error
     use global,       only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes,&
-                            matching_bins, current_batch
+                            matching_bins
     use mesh,         only: mesh_indices_to_bin
     use mesh_header,  only: StructuredMesh
     use string,       only: to_str
     use tally_header, only: TallyObject
 
+#ifdef CMFD_DEBUG
+    use global,       only: current_batch
+#endif
+
     integer :: nx            ! number of mesh cells in x direction
     integer :: ny            ! number of mesh cells in y direction
     integer :: nz            ! number of mesh cells in z direction
@@ -72,7 +77,6 @@ contains
     integer :: k             ! iteration counter for z
     integer :: g             ! iteration counter for g
     integer :: h             ! iteration counter for outgoing groups
-    integer :: l             ! iteration counter for leakage debug
     integer :: ital          ! tally object index
     integer :: ijk(3)        ! indices for mesh cell
     integer :: score_index   ! index to pull from tally object
@@ -82,10 +86,13 @@ contains
     integer :: i_filter_eout ! index for outgoing energy filter
     integer :: i_filter_surf ! index for surface filter
     real(8) :: flux          ! temp variable for flux
-    real(8) :: leak1         ! group 1 leakage 
-    real(8) :: leak2         ! group 2 leakage
     type(TallyObject),    pointer :: t => null() ! pointer for tally object
     type(StructuredMesh), pointer :: m => null() ! pointer for mesh object
+#ifdef CMFD_DEBUG
+    integer :: l             ! iteration counter for leakage debug
+    real(8) :: leak1         ! group 1 leakage 
+    real(8) :: leak2         ! group 2 leakage
+#endif
 
     ! Extract spatial and energy indices from object
     nx = cmfd % indices(1)
@@ -709,9 +716,12 @@ contains
   subroutine compute_dhat()
 
     use constants,  only: CMFD_NOACCEL, ZERO
-    use global,     only: cmfd, cmfd_coremap, message, dhat_reset, current_batch
+    use global,     only: cmfd, cmfd_coremap, message, dhat_reset
     use output,     only: write_message
     use string,     only: to_str
+#ifdef CMFD_DEBUG
+    use global,     only: current_batch
+#endif
 
     integer :: nx             ! maximum number of cells in x direction
     integer :: ny             ! maximum number of cells in y direction
@@ -869,8 +879,8 @@ contains
 
   function get_reflector_albedo(l, g, i, j, k)
 
-    use constants,  only: ALBEDO_REJECT
-    use global,     only: cmfd, cmfd_hold_weights
+    use constants,  only: ONE
+    use global,     only: cmfd
 
     real(8) :: get_reflector_albedo ! reflector albedo
     integer, intent(in) :: i ! iteration counter for x
@@ -892,8 +902,7 @@ contains
     ! Calculate albedo
     if ((shift_idx ==  1 .and. current(2*l  ) < 1.0e-10_8) .or. &
         (shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then
-      albedo = ALBEDO_REJECT 
-      cmfd_hold_weights = .true. 
+      albedo = ONE
     else
       albedo = (current(2*l-1)/current(2*l))**(shift_idx)
     end if
@@ -958,8 +967,6 @@ contains
 
     keff_bal = cmfd % keff_bal
 
-    print *, 'KEFF BALANCE is:', keff_bal
-
     ! Return if not two groups
     if (ng /= 2) return
 
@@ -1007,12 +1014,6 @@ contains
           nsigf12 = cmfd % nfissxs(1,2,i,j,k)
           nsigf22 = cmfd % nfissxs(2,2,i,j,k)
 
-          ! Check for no fission into group 2
-          if (.not.(nsigf12 < 1e-6_8 .and. nsigf22 < 1e-6_8)) then
-            write(OUTPUT_UNIT,'(A,1PE11.4,1X,1PE11.4)') 'Fission in G=2', &
-                  nsigf12,nsigf22
-          end if
-
           ! Compute absorption xs
           siga1 = sigt1 - sigs11 - sigs12
           siga2 = sigt2 - sigs22 - sigs21

src/cmfd_execute.F90

@@ -84,7 +84,7 @@ contains
 
     use global,            only: cmfd_begin, cmfd_on, &
                                  cmfd_reset, cmfd_run,               &
-                                 current_batch, cmfd_hold_weights
+                                 current_batch
 
     ! Check to activate CMFD diffusion and possible feedback
     ! this guarantees that when cmfd begins at least one batch of tallies are
@@ -283,7 +283,6 @@ contains
     integer :: nz       ! maximum number of cells in z direction
     integer :: ng       ! maximum number of energy groups
     integer :: i        ! iteration counter
-    integer :: j        ! iteration counter
     integer :: ijk(3)   ! spatial bin location
     integer :: e_bin    ! energy bin of source particle
     integer :: n_groups ! number of energy groups

src/cmfd_input.F90

@@ -66,6 +66,7 @@ contains
     use xml_interface
     use, intrinsic :: ISO_FORTRAN_ENV
 
+    integer :: i
     integer :: ng
     integer, allocatable :: iarray(:)
     integer, allocatable :: int_array(:)

src/cmfd_solver.F90

@@ -102,9 +102,8 @@ contains
 
     use constants, only: ONE, ZERO
     use error,     only: fatal_error
-    use global,    only: cmfd, message, cmfd_shift, keff, &
-                         cmfd_ktol, cmfd_stol
-    use global,    only: cmfd_write_matrices
+    use global,    only: cmfd, cmfd_shift, keff, cmfd_ktol, cmfd_stol, &
+                         cmfd_write_matrices
 
     logical :: adjoint
 
@@ -180,10 +179,13 @@ contains
   subroutine compute_adjoint()
 
     use error,   only: fatal_error
+#ifndef PETSC
     use global,  only: message
-#ifdef PETSC
+#else
     use global,  only: cmfd_write_matrices
+#endif
 
+#ifdef PETSC
     ! Transpose matrices
     call loss % transpose()
     call prod % transpose()

src/constants.F90

@@ -392,9 +392,6 @@ module constants
   ! constant to represent a zero flux "albedo"
   real(8), parameter :: ZERO_FLUX = 999.0_8
 
-  ! constant to represent albedo rejection
-  real(8), parameter :: ALBEDO_REJECT = 999.0_8
-
   ! constant for writing out no residual
   real(8), parameter :: CMFD_NORES = 99999.0_8
 

src/global.F90

@@ -325,6 +325,9 @@ module global
   ! Flag to reset dhats to zero
   logical :: dhat_reset = .false.
 
+  ! Flag to enforce neutron balance
+  logical :: cmfd_fix_balance = .false.
+
   ! Flag to activate neutronic feedback via source weights
   logical :: cmfd_feedback = .false.
 
@@ -332,9 +335,6 @@ module global
   integer :: n_cmfd_meshes  = 1 ! # of structured meshes
   integer :: n_cmfd_tallies = 3 ! # of user-defined tallies
 
-  ! Flag to hold cmfd weight adjustment factors
-  logical :: cmfd_hold_weights = .false.
-
   ! Eigenvalue solver type
   character(len=10) :: cmfd_solver_type = 'power'