From 4de45e7b08b9d6a71a8d1c38b9de2bdd066f4809 Mon Sep 17 00:00:00 2001
From: Bryan Herman <bherman@mit.edu>
Date: Mon, 3 Mar 2014 11:16:50 -0500
Subject: [PATCH] dhats can be reset from cmfd input

---
 src/cmfd_data.F90 | 25 +++++++++++++++++++++----
 1 file changed, 21 insertions(+), 4 deletions(-)

diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90
index df6b38e910..1a35ac09bc 100644
--- a/src/cmfd_data.F90
+++ b/src/cmfd_data.F90
@@ -10,8 +10,6 @@ module cmfd_data
   private
   public :: set_up_cmfd, neutron_balance
 
-  logical :: dhat_reset = .false.
-
 contains
 
 !==============================================================================
@@ -705,7 +703,8 @@ contains
   subroutine compute_dhat()
 
     use constants,  only: CMFD_NOACCEL, ZERO
-    use global,     only: cmfd, cmfd_coremap
+    use global,     only: cmfd, cmfd_coremap, message, dhat_reset
+    use output,     only: write_message
 
     integer :: nx             ! maximum number of cells in x direction
     integer :: ny             ! maximum number of cells in y direction
@@ -825,7 +824,9 @@ contains
               cmfd%dhat(l,g,i,j,k) = dhat
 
               ! check for dhat reset
-              if (dhat_reset) cmfd%dhat(l,g,i,j,k) = ZERO
+              if (dhat_reset) then
+                cmfd%dhat(l,g,i,j,k) = ZERO
+              end if
 
             end do LEAK
 
@@ -837,6 +838,12 @@ contains
 
     end do ZLOOP
 
+    ! write that dhats are zero
+    if (dhat_reset) then
+      message = 'Dhats reset to zero.'
+      call write_message(1)
+    end if
+
   end subroutine compute_dhat
 
 !===============================================================================
@@ -914,6 +921,16 @@ contains
     real(8) :: siga1      ! group 1 abs xs
     real(8) :: siga2      ! group 2 abs xs
     real(8) :: sigs12_eff ! effective downscatter xs
+    real(8) :: a          ! matrix (1,1) element for balance
+    real(8) :: b          ! matrix (1,2) element for balance
+    real(8) :: c          ! matrix (2,1) element for balance
+    real(8) :: d          ! matrix (2,2) element for balance
+    real(8) :: r1         ! right hand side element 1
+    real(8) :: r2         ! right hand side element 2
+    real(8) :: det        ! determinant of balance matrix
+
+    message = 'Correcting neutron balance'
+    call write_message(1)
 
     ! Extract spatial and energy indices from object
     nx = cmfd % indices(1)