From 0465aaf58d1e6469e5f36f24ec02362fdc6c5a28 Mon Sep 17 00:00:00 2001
From: Wiarda <wiardada@ornl.gov>
Date: Fri, 5 Nov 2021 12:23:32 -0400
Subject: [PATCH] Make most of XCT independent of SAMMY global parameters

---
 sammy/src/blk/AuxGridHelper_M.f90          |   2 +-
 sammy/src/blk/Fixedi_common.f90            |   3 -
 sammy/src/blk/Varyr_common.f90             |  10 -
 sammy/src/clq/ArtificalCross_M.f90         |   5 +-
 sammy/src/cro/mcro0.f90                    |   5 +-
 sammy/src/inp/minp05.F                     |   2 +
 sammy/src/mlb/mmlb0.f90                    |  10 +-
 sammy/src/rec/mrec0.f                      |   4 +-
 sammy/src/rec/mrec3.f90                    |   2 +-
 sammy/src/sammy/CMakeLists.txt             |   1 -
 sammy/src/the/CrossSectionCalcDriver_M.f90 |  31 +-
 sammy/src/the/CrossSectionCalculator_M.f90 | 103 +++++--
 sammy/src/the/ZeroKCrossCorrections_M.f90  |  42 ++-
 sammy/src/the/mthe0.f90                    |  17 +-
 sammy/src/xct/XctCrossCalcImpl_M.f90       |  79 ++++-
 sammy/src/xct/XctCrossCalc_M.f90           |  13 +-
 sammy/src/xct/mxct0.f90                    | 328 ++-------------------
 sammy/src/xct/mxct02.f90                   | 219 --------------
 sammy/src/xct/mxct03.f90                   |  16 +-
 sammy/src/xct/mxct04.f90                   |  42 +--
 sammy/src/xct/mxct05.f90                   |  34 +--
 sammy/src/xct/mxct11.f90                   |   5 +-
 sammy/src/xct/mxct17.f90                   |  31 +-
 sammy/src/xct/mxct18.f90                   |   2 +-
 sammy/src/xct/mxct21.f90                   |   5 +-
 sammy/src/xct/mxct22.f90                   |   5 +-
 sammy/src/xct/mxct23.f90                   |   2 +-
 sammy/src/xct/mxct26.f90                   |   9 +-
 sammy/src/xct/mxct31.f90                   |   1 -
 sammy/src/xct/mxct32.f90                   |   1 -
 30 files changed, 300 insertions(+), 729 deletions(-)
 delete mode 100644 sammy/src/blk/Varyr_common.f90

diff --git a/sammy/src/blk/AuxGridHelper_M.f90 b/sammy/src/blk/AuxGridHelper_M.f90
index c4c232f4a..acdc7de1d 100644
--- a/sammy/src/blk/AuxGridHelper_M.f90
+++ b/sammy/src/blk/AuxGridHelper_M.f90
@@ -51,7 +51,7 @@ module AuxGridHelper_M
            if (.not.C_ASSOCIATED(grid%instance_ptr)) return
 
            ii  = ioff-1
-           if (ii.lt.0) ii = 0
+           if (ii.lt.0) ii = 0      
            call grid%setRowMax(ii)
        end subroutine setAuxGridRowmax
 
diff --git a/sammy/src/blk/Fixedi_common.f90 b/sammy/src/blk/Fixedi_common.f90
index d0f2a50c0..07e55e10f 100644
--- a/sammy/src/blk/Fixedi_common.f90
+++ b/sammy/src/blk/Fixedi_common.f90
@@ -194,8 +194,6 @@ module fixedi_m
      integer,pointer :: Mxwrec  => lfdim(139)
      integer,pointer :: Jtheta  => lfdim(140)
      integer,pointer :: Kaddcr  => lfdim(141)
-     integer,pointer :: Nrfil3  => lfdim(142)
-     integer,pointer :: Npfil3  => lfdim(143)
      integer,pointer :: Nnnrpi  => lfdim(144)
      integer,pointer :: Krefit  => lfdim(145)
      
@@ -251,7 +249,6 @@ module fixedi_m
      integer,pointer :: Nudeng  => lfdim(187)
      integer,pointer :: Nudtim  => lfdim(188)
      integer,pointer :: Nudmax  => lfdim(189)
-     integer,pointer :: Kkxlmn  => lfdim(190)
 
      ! indexer on lfdim covers up to 190 (next should be 191
      ! old group 9
diff --git a/sammy/src/blk/Varyr_common.f90 b/sammy/src/blk/Varyr_common.f90
deleted file mode 100644
index e76e53645..000000000
--- a/sammy/src/blk/Varyr_common.f90
+++ /dev/null
@@ -1,10 +0,0 @@
-
-module varyr_common_m
-
-      IMPLICIT NONE
-      double precision, save :: Su
-      double precision, save :: Squ     
-      double precision, save :: Etz
-      double precision, save :: Elz
-
-end module varyr_common_m
diff --git a/sammy/src/clq/ArtificalCross_M.f90 b/sammy/src/clq/ArtificalCross_M.f90
index bf9800f4a..087f9263a 100644
--- a/sammy/src/clq/ArtificalCross_M.f90
+++ b/sammy/src/clq/ArtificalCross_M.f90
@@ -27,12 +27,13 @@ module ArtificalCross_M
   end type ArtificalCross
 
 contains
-subroutine ArtificalCross_setRange(this, emin, emax, num)
+subroutine ArtificalCross_setRange(this, emin, emax, num, solve)
     class(ArtificalCross) :: this
     real(kind=8)::emin,emax
     integer::num
+    logical::solve
 
-    call CrossSectionCalculator_setRange(this, emin, emax, num)
+    call CrossSectionCalculator_setRange(this, emin, emax, num, solve)
     this%Jdirac = 0
 
     this%X2 = emax
diff --git a/sammy/src/cro/mcro0.f90 b/sammy/src/cro/mcro0.f90
index 5d789c691..d69ce2751 100644
--- a/sammy/src/cro/mcro0.f90
+++ b/sammy/src/cro/mcro0.f90
@@ -7,8 +7,7 @@
       use ifwrit_m, only : Ks_Res, Ksolve
       use oopsch_common_m, only : Nowwww, Segmen
       use rsl7_m, only : Set_Kws_Xct, Figure_Kws_Cro
-      use AuxGridHelper_M, only : setAuxGridOffset,   &
-                                  setAuxGridRowMax
+
       IMPLICIT None
       integer::I,Kone, Idimen
       external Idimen
@@ -21,8 +20,6 @@
       Segmen(3) = 'O'
       Nowwww = 0
 
-      call setAuxGridOffset(1) ! reset starting point for auxillary grid
-      call setAuxGridRowMax(0)
 !
       Ks_Res = Ksolve
 !
diff --git a/sammy/src/inp/minp05.F b/sammy/src/inp/minp05.F
index 2dc2147b3..0d359ec3a 100644
--- a/sammy/src/inp/minp05.F
+++ b/sammy/src/inp/minp05.F
@@ -299,6 +299,8 @@ C
          Kredwa = 1
       ELSE IF (Number.EQ.185) THEN
          Kgauss = 1
+         write(6,*) "GAUSSIAN RESONANCES ARE ASSUMED not yet supported"
+         STOP '[STOP in Fixdef in inp/minp05.f  # 2a]'
       ELSE IF (Number.EQ.161) THEN
          Leeb = 1
 C
diff --git a/sammy/src/mlb/mmlb0.f90 b/sammy/src/mlb/mmlb0.f90
index 1775eba32..8f8fd924e 100644
--- a/sammy/src/mlb/mmlb0.f90
+++ b/sammy/src/mlb/mmlb0.f90
@@ -16,11 +16,9 @@ module MlbwCalc_m
       use oopsch_common_m, only : Nowwww, Segmen
       use cbro_common_m, only : Segnam  
       use rsl7_m     
-      use AuxGridHelper_M, only : setAuxGridOffset, setAuxGridRowMax
-      IMPLICIT None      
-      integer::I,Kone, Idimen
 
-      external Idimen
+      IMPLICIT None      
+      integer::I,Kone
 
 
       IF (Krmatx.EQ. 1) WRITE (6,99999)
@@ -37,10 +35,6 @@ module MlbwCalc_m
       Segmen(2) = 'L'
       Segmen(3) = 'B'
       Nowwww = 0
-!
-! *** guesstimate size of array needed for SAMMY-MLBW
-      call setAuxGridOffset(1) ! reset starting point for auxillary grid
-      call setAuxGridRowMax(0)
 
       CALL Figure_Kws_Cro (Kone)
 !
diff --git a/sammy/src/rec/mrec0.f b/sammy/src/rec/mrec0.f
index 3bfa94b63..57047ec63 100644
--- a/sammy/src/rec/mrec0.f
+++ b/sammy/src/rec/mrec0.f
@@ -58,8 +58,6 @@ C *** find array sizes for endf/b-vi file 3 cross sections
          end if
          call zeroKCalc%endfReader%crossAddedToStellarAverage(tab1Data)
       End if
-      Nrfil3 = 1
-      Npfil3 = 1
 C
 C *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-RECONSTRUCT
       CALL Estrec (N1, N2, NA, N3, NN, N5, N6, Nesave, Nesub, Nreact,
@@ -234,7 +232,7 @@ C
 C *** the following are not used but must be defined anyway
 C
 C *** zero ***
-      K0 = 2*(Nrfil3+Npfil3)
+      K0 = 4
 C
 C *** one ***
       K1 = Mres
diff --git a/sammy/src/rec/mrec3.f90 b/sammy/src/rec/mrec3.f90
index c789e335c..4cc20d374 100644
--- a/sammy/src/rec/mrec3.f90
+++ b/sammy/src/rec/mrec3.f90
@@ -12,7 +12,6 @@ contains
 ! *** This file is modified from routines in segment XCT, to generate
 ! ***    cross sections one energy at a time, with no derivatives wanted
 !
-      use varyr_common_m, only : Su, Squ
       use xct4_m
       use XctCrossCalc_M
       IMPLICIT None
@@ -20,6 +19,7 @@ contains
       class(XctCrossCalc)::xct
       real(kind=8):: Ssseee(*), Eee
       integer::I
+      real(kind=8)::Su, Squ
 !
       Su = Eee
       Squ = dSQRT(Su)
diff --git a/sammy/src/sammy/CMakeLists.txt b/sammy/src/sammy/CMakeLists.txt
index ef175daaa..efbe7fc24 100644
--- a/sammy/src/sammy/CMakeLists.txt
+++ b/sammy/src/sammy/CMakeLists.txt
@@ -555,7 +555,6 @@ APPEND_SET(SAMMY_SOURCES
             ../blk/Rpijnk_common.f90
             ../blk/Rpires_common.f90
             ../blk/Rpirrr_common.f90
-            ../blk/Varyr_common.f90
             ../blk/Kzznew_common.f90
             ../blk/Aaarfs_common.f90
             ../blk/Z00001_common.f90
diff --git a/sammy/src/the/CrossSectionCalcDriver_M.f90 b/sammy/src/the/CrossSectionCalcDriver_M.f90
index 584cd40d4..da01bc102 100644
--- a/sammy/src/the/CrossSectionCalcDriver_M.f90
+++ b/sammy/src/the/CrossSectionCalcDriver_M.f90
@@ -88,6 +88,7 @@ module CrossSectionCalcDriver_M
      !!  - Itzero is t0 to be adjusted
      !!  - do we need to calculated derivatives
      !!    Note: wantDeriv is false if we don't do a fit. If pup'ed parameters exist, derivatives for them are calculated
+     !!  - Kpolar are the fission width given in polar coordinates
      !!
      subroutine CrossSectionCalcDriver_setUpCalculator(this, Krmatx,Kkclqx, Kkkclq, niso, covData, resParData, &
                                                        radFitFlags,Ilzero, Itzero, wantDeriv, doShiftRes, Kpolar)
@@ -173,12 +174,21 @@ module CrossSectionCalcDriver_M
         call this%calculator%setEnergyIndependent(this%Kcros, Twomhb, this%kwcoul, Etac)
      end subroutine
 
+     !!
+     !! Set additional parameters for the calculator
+     !!
+     !!  - llmax maximum number of Clebsch-Gordon coefficients
+     !!  - Kfinit finite-size corrections flag
+     !!  - wantSelfIndicate  do we need to calculate self-indicated cross section data
+     !!  - Kkfis do we have fission
+     !!  - debug  do we want  extra debug output
+     !!
      subroutine CrossSectionCalcDriver_setAddtionalParams (this,  lllmax, Kfinit, wantSelfIndicate, Kkkfis, debug)
           class(CrossSectionCalcDriver)::this
           logical::debug             ! do we want  extra debug output
           integer::lllmax            ! maximum number of Clebsch-Gordon coefficients
           integer::Kfinit            ! finite-size corrections flag
-          integer::Kkkfis            ! what fission channels for  cro
+          integer::Kkkfis            ! what fission channels for  cro and xct
           logical::wantSelfIndicate  ! do we need to calculate self-indicated cross section data
 
           if (Kkkfis.ne.0.and.associated(this%croCalc)) then
@@ -204,6 +214,7 @@ module CrossSectionCalcDriver_M
         empty = .true.
         call this%calcData%setNotSetReturnsZero(empty)
         call this%calcData%setAccumulate(empty)
+        call this%calcDataSelf%setAccumulate(empty)
 
 
         if (this%calculatorInit) then
@@ -217,11 +228,13 @@ module CrossSectionCalcDriver_M
         empty = .false.
         call this%calcData%setNotSetReturnsZero(empty)
         call this%calcData%setAccumulate(empty)
+        call this%calcDataSelf%setAccumulate(empty)
      end subroutine
 
 
      !!
      !! Set up the desired cross section data if using xct algorithm
+     !!  - calc the object on which to set the information
      !!
      subroutine setup_desiredCross(this, calc)
         class(CrossSectionCalcDriver)::this
@@ -339,6 +352,7 @@ module CrossSectionCalcDriver_M
      !!
      !! Do we want to get the abundance from the
      !! spin group or from the isotope
+     !!  - setAbund if true the abundance is retrieved from the spin group instead of the isotope
      !!
      subroutine CrossSectionCalcDriver_getAbundanceFromSpinGroup(this, setAbund)
          class(CrossSectionCalcDriver)::this
@@ -349,13 +363,24 @@ module CrossSectionCalcDriver_M
          end if
      end subroutine
 
-     subroutine CrossSectionCalcDriver_setRange(this, emin, emax, num)
+
+     !!
+     !! Set the range and number of energy points
+     !! This gives the calculator a change to reserve data
+     !! in the DerivativeHandler object or other needed data.
+     !!  - emin the lowest energy point to be calculated
+     !!  - emax the upper energy point to be calculated
+     !!  - num the number of  energy points
+     !!  - solve do we do a solve step, i.e. all derivatives are needed
+     !!
+     subroutine CrossSectionCalcDriver_setRange(this, emin, emax, num, solve)
         class(CrossSectionCalcDriver) :: this
         real(kind=8)::emin,emax
         integer::num
+        logical::solve
 
         if (this%calculatorInit) then
-           call this%calculator%setRange(emin, emax, num)
+           call this%calculator%setRange(emin, emax, num, solve)
         end if
   end subroutine
 
diff --git a/sammy/src/the/CrossSectionCalculator_M.f90 b/sammy/src/the/CrossSectionCalculator_M.f90
index aeccd8aaf..16c4d784c 100644
--- a/sammy/src/the/CrossSectionCalculator_M.f90
+++ b/sammy/src/the/CrossSectionCalculator_M.f90
@@ -70,8 +70,8 @@ module CrossSectionCalculator_M
      procedure, pass(this) :: setRange => CrossSectionCalculator_setRange   ! set the range and reserve enough  space in crossData (the latter for efficiency only)
      procedure, pass(this) :: initialize => CrossSectionCalculator_initialize
      procedure, pass(this) :: destroy => CrossSectionCalculator_destroy
-     procedure, pass(this) :: calcCross => CrossSectionCalculator_calcCross
-     procedure, pass(this) :: setDerivFlag  => CrossSectionCalculator_setDerivFlag
+     procedure, pass(this) :: calcCross => CrossSectionCalculator_calcCross   ! calculate the cross section
+     procedure, pass(this) :: setDerivFlag  => CrossSectionCalculator_setDerivFlag  ! check whether a certain derivative neeeds to be calculated
   end type CrossSectionCalculator
 
 contains
@@ -116,6 +116,7 @@ subroutine CrossSectionCalculator_setEchan(this)
      end do
   end do
 end subroutine
+
 !!
 !! Set the boundary condition for each group and channel
 !! - Twomhb sqrt(2m)/hba
@@ -140,6 +141,7 @@ subroutine CrossSectionCalculator_setBound(this, Twomhb, kwcoul, Etac)
      end do
   end do
 end subroutine
+
 !!
 !! Get the abundance for the indicated spin group
 !! - igr the group for which to get the abundance
@@ -207,17 +209,20 @@ subroutine CrossSectionCalculator_calcCross(this, ener, Ipoten)
     if (this%enerSq.gt.0.0d0) this%enerSq = dSQRT(this%enerSq)
 end subroutine
 
-!
-!  For a given spin group, determine the number of varied parameters
-!  (including their resonance and channel.
-!
-!  Please note that the index of the varied parameter in Inum does
-!  not at all be the same order as the position in the covariance, i.e. Inum(index, 1)
-!
-!  Thus, if for example gamma width data  for a given spin group are fitted together,
-!  there will be an Inum entry for all gamma width. However, for all of the gamma width
-!  the position in the covariance, i.e. Inum(index,1) will be identical
-!
+!!
+!!  For a given spin group, determine the number of varied parameters
+!!  (including their resonance and channel.
+!!
+!!  Please note that the index of the varied parameter in Inum does
+!!  not at all be the same order as the position in the covariance, i.e. Inum(index, 1)
+!!
+!!  Thus, if for example gamma width data  for a given spin group are fitted together,
+!!  there will be an Inum entry for all gamma width. However, for all of the gamma width
+!!  the position in the covariance, i.e. Inum(index,1) will be identical
+!!
+!!  - ires the starting resonance for the spin group (it is assumed resonance are ordered by spingroup)
+!!  - igr the spin group for which to get the data
+!!  - fillIt if present and false don't fill the this%Inum array (true is assumed if not present)
 subroutine CrossSectionCalculator_getParamPerSpinGroup(this, ires, igr,  fillIt)
     class(CrossSectionCalculator) :: this
     integer,intent(inout)::ires  ! starting resonance for this group   
@@ -268,17 +273,43 @@ subroutine CrossSectionCalculator_getParamPerSpinGroup(this, ires, igr,  fillIt)
     end do    
 end subroutine
 
-logical function CrossSectionCalculator_setDerivFlag(this, iflag, doSolve) result(set)
+!!
+!! Determine if a derivativate needs to be calculated:
+!!  - If we are solving, all flagged and pup'ed parameter derivatives are needed
+!!  - If we are not solving, only pup'ed parmeters need  to be considered
+!!
+!! parameters are:
+!!   - iflag, the fit flag, aka the position in the covariance matrix
+!!   - doSolve do we want to solve. If not, we only need derivatives for pup'ed parameters
+!!   _ inside optional parameter (assumed false) that is used if calling the function from inside this class
+logical function CrossSectionCalculator_setDerivFlag(this, iflag, doSolve, inside) result(set)
     class(CrossSectionCalculator) :: this
     integer::iflag
     logical::doSolve
 
-    set = .false.
-    if (iflag.le.0) return
+    logical::inside
+    logical::insideReal
+
+    optional inside
+
+    insideReal = .false.
+    if (present(inside)) insideReal = inside
+
+    set = .false.   ! assume false
+    if (iflag.le.0) return   ! not flagged
+
+    if (this%covariance%isPupedParameter(iflag)) then  ! always if puped
+        set = .true.
+        return
+     end if
+
+    if (.not.doSolve) return   ! don't want to solve and it's not a pup'ed parameter
+
+    if (insideReal) then
+        if  (this%onlyPupDerivs)  return ! inside doSolve is true, if we have pup'ed resonance parameters
+    end if
+
     set = .true.
-    if (doSolve) return
-    set = .false.
-    if (this%covariance%isPupedParameter(iflag)) set = .true.
 end  function CrossSectionCalculator_setDerivFlag
 
 !!
@@ -305,21 +336,21 @@ subroutine CrossSectionCalculator_Which_Derivs(this)
    this%Ifzzz = .false.
    if (.not.this%wantDerivs) return
 
-   set = this%setDerivFlag(this%Itzero, this%wantDerivs)
+   set = this%setDerivFlag(this%Itzero, this%wantDerivs, .true.)
    if (set) this%Ifzzz = .true.
-   set = this%setDerivFlag(this%Ilzero, this%wantDerivs)
+   set = this%setDerivFlag(this%Ilzero, this%wantDerivs, .true.)
    if (set) this%Ifzzz = .true.
 
    DO ig=1,this%resData%getNumSpinGroups()
       call this%resData%getSpinGroupInfo(spinInfo, ig)
       DO Ich=1,spinInfo%getNumChannels()
-         set = this%setDerivFlag(this%radiusData%getTrueFitFlag(Ig, Ich), this%wantDerivs)
+         set = this%setDerivFlag(this%radiusData%getTrueFitFlag(Ig, Ich), this%wantDerivs, .true.)
          if (set) then
             this%Ifrad = .true.
             this%Ifradt = .true.
          end if
          if (this%Ifrad.and.this%Ifradt) exit
-         set = this%setDerivFlag(this%radiusData%getEffFitFlag(Ig, Ich), this%wantDerivs)
+         set = this%setDerivFlag(this%radiusData%getEffFitFlag(Ig, Ich), this%wantDerivs, .true.)
          if(set) then
             this%Ifrad = .true.
          end if
@@ -330,7 +361,7 @@ subroutine CrossSectionCalculator_Which_Derivs(this)
 
    do is = 1, this%resData%getNumIso()
       call this%resData%getIsoInfo(isoInfo, Is)
-      set = this%setDerivFlag(isoInfo%getFitOption(), this%wantDerivs)
+      set = this%setDerivFlag(isoInfo%getFitOption(), this%wantDerivs, .true.)
       if (set) then
          this%Ifiso = .true.
          exit
@@ -343,7 +374,7 @@ subroutine CrossSectionCalculator_Which_Derivs(this)
         if (.not.this%resData%hasRexInfo(ig, Ich)) cycle
         call this%resData%getRextInfoByGroup(rextInfo, ig, Ich)
         DO Is = 1, rextInfo%getNrext()
-           set = this%setDerivFlag(rextInfo%getIflSammyIndex(Is), this%wantDerivs)
+           set = this%setDerivFlag(rextInfo%getIflSammyIndex(Is), this%wantDerivs, .true.)
            if( set) then
               this%Ifext = .true.
               exit
@@ -393,19 +424,32 @@ end subroutine
 !!  - emin the lowest energy point to be calculated
 !!  - emax the upper energy point to be calculated
 !!  - num the number of  energy points
+!!  - solve do we do a solve step, i.e. all derivatives are needed
 !!
-subroutine CrossSectionCalculator_setRange(this, emin, emax, num)
+subroutine CrossSectionCalculator_setRange(this, emin, emax, num, solve)
    class(CrossSectionCalculator) :: this
    real(kind=8)::emin,emax
    integer::num
+   logical::solve
+
+   integer::itot, irow, nsig
 
-   integer::itot, irow
+   nsig = this%crossData%getNnnsig()
+   if (this%reactType.eq.6) then
+       if( nsig.ne.1) then
+           write(6,*)" Expected Nnsig to be 1 if calculating eta"
+           write(21,*)" Expected Nnsig to be 1 if calculating eta"
+           stop
+       end if
+       call this%crossData%setNnsig(2)
+       nsig = 2
+   end if
 
     itot = 1
-    if (this%wantDerivs) then
+    if (solve.or.this%covariance%getPupedParam().gt.0) then
        itot = this%covariance%getNumTotalParam() + 1
      end if
-     irow = this%crossData%getNnnsig() * num
+     irow =  nsig * num
      call this%crossData%reserve(irow, itot)
 
      if (this%wantCrossPerSpin) then
@@ -422,6 +466,7 @@ end subroutine
 !! - needAngular is angular distribution desired
 !! - Itzero fit flag for t0
 !! - Ilzero fit flag for flight path length
+!! - doShiftRes do we want to SHIFT RESONANCE ENERGIES VIA SHIFT FACTOR
 subroutine CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
     class(CrossSectionCalculator) :: this
     type(SammyRMatrixParameters)::pars
diff --git a/sammy/src/the/ZeroKCrossCorrections_M.f90 b/sammy/src/the/ZeroKCrossCorrections_M.f90
index ccd575a71..51e58b891 100644
--- a/sammy/src/the/ZeroKCrossCorrections_M.f90
+++ b/sammy/src/the/ZeroKCrossCorrections_M.f90
@@ -47,16 +47,18 @@ contains
   !!              In practice this is only used for the fake  cross section data, to determine the
   !!              position of the delta-function peak)
   !!  - expData needed as an argument for grid and userGrid to get the enery
-  subroutine ZeroKCrossCorrections_resonanceRecon(this, grid, userGrid, expData)
+  !!  - solve do we do a solve step, i.e. all derivatives are needed
+  subroutine ZeroKCrossCorrections_resonanceRecon(this, grid, userGrid, expData, solve)
     use SammyGridAccess_M
     use GridData_M
     class(ZeroKCrossCorrections)::this
     type(SammyGridAccess)::grid, userGrid
     type(GridDataList)::expData
+    logical::solve
 
-    real(kind=8)::ener, emin, emax
+    real(kind=8)::ener, emin, emax, val
     integer::numEl, numUser
-    integer::ipos, iel, ig, ngr, is
+    integer::ipos, iel, ig, ngr, is, iso
 
     call this%driver%calcData%nullify()
     call this%driver%calcDataSelf%nullify()
@@ -65,7 +67,7 @@ contains
     numUser = userGrid%getNumEnergies(expData)
     emin =  userGrid%getEnergy(1, expData)
     emax =  userGrid%getEnergy(numUser, expData)
-    call this%driver%setRange(emin, emax, numEl)
+    call this%driver%setRange(emin, emax, numEl, solve)
 
     ipos = 0
     this%dataStart = 0
@@ -81,10 +83,10 @@ contains
           WRITE (21,10000) ener
        END IF
 
-       if (ener.le.0.0d0.and..not.this%wantNeg) then
+       if (ener.le.0.0d0.and..not.this%wantNeg) then           
            if (this%moreCorrections) then
               ipos = ipos + 1
-              call this%driver%calcCross(ipos, 0.0d0, this%Ipoten)   ! set data to zero
+              call this%driver%calcCross(ipos, 0.0d0, this%Ipoten)   ! set data to zero             
            end if
            cycle
         end if
@@ -92,10 +94,24 @@ contains
        ipos = ipos + 1
        if (this%dataStart.eq.0) this%dataStart = ipos
 
-       if (.not.(this%artificalGrid.and.ener.gt.this%artificalGridMax)) then
+       if (this%artificalGrid.and.ener.gt.this%artificalGridMax) then
            ! do nothing  all data are already set to null for
            ! artificalGrid and ener.gt.this%artificalGridMax
-          call this%driver%calcCross(ipos, ener, this%Ipoten)
+
+              ! SAMMY did not zero the cross section for isotopes 2, 3, ...
+              ! but reused the previous value
+              ! todo: delete this to make the code correct
+              if (ipos.gt.1) then
+                 do Iso=2,this%driver%calculator%numIso
+                    do is = 1, this%driver%calcData%getNnnsig()
+                       val = this%driver%calcData%getDataNs(ipos-1, is, 0, iso)
+                       if (val.eq.0.0d0) cycle
+                       call this%driver%calcData%addDataNs(ipos, is, 0, iso, val)
+                    end do
+                 end do
+              end if
+       else
+          call this%driver%calcCross(ipos, ener, this%Ipoten)       
        end if
     end do
 
@@ -103,6 +119,7 @@ contains
       this%summedOverIsotopes = .true.
     end if
 
+
   end subroutine
 
   !!
@@ -133,7 +150,7 @@ contains
   end subroutine
 
   !!
-  !! Save the cross section data only int Theory
+  !! Save the cross section data only into Theory
   !!
   !! - Theory object to save the data to
   !! - nnsig the number of sections in Theory
@@ -514,6 +531,13 @@ contains
   END function
 
 
+  !!
+  !! Add paramagnetic cross section
+  !!  - grid   energy grid on which to calculate the cross section
+  !!  - userGrid  grid from which to calculate the emin and emax value
+  !!  - expData needed as an argument for grid and userGrid to get the enery
+  !!  - covariance covariance data, needed to find the position for the derivatives
+  !!  - wantDeriv do we need derivatives (i.e. are we solving).
   subroutine ZeroKCrossCorrections_AddParam(this,  grid, expData, covariance, wantDeriv)
      use SammyGridAccess_M
      use paramagnetic_cross_m
diff --git a/sammy/src/the/mthe0.f90 b/sammy/src/the/mthe0.f90
index bd9ef229f..7fbba54b6 100644
--- a/sammy/src/the/mthe0.f90
+++ b/sammy/src/the/mthe0.f90
@@ -1,5 +1,6 @@
 module mthe0_M
   contains
+
 !
 !
 ! --------------------------------------------------------------
@@ -110,13 +111,9 @@ module mthe0_M
       call setAuxGridOffset(1) ! reset starting point for auxillary grid
       call setAuxGridRowMax(0)
 
-      ! reconstruct the cross section at 0K on the desired grid
-      !
-      ! for now this will only do something for MLBW or SLBW    
-      call zeroKCalc%resonanceRecon(grid, usergrid, expData)
 
       ! pick up some additional information
-      ! and sizing
+      ! and sizing, needed in downstream modules
       IF (Kkkclq.NE.0) THEN
          CALL Samclq_0
       ELSE IF (Krmatx.EQ.0) THEN
@@ -131,11 +128,19 @@ module mthe0_M
          STOP '[STOP in Samthe_0 in the/mthe0.f]'
       END If
 
+      ! reconstruct the cross section at 0K on the desired grid
+      !
+      ! move this after the call to Samxct_0
+      ! as they can update Ksolve
       if (Ksolve.eq.2) then
           wantDeriv = .false.
       else
          wantDeriv = .true.
       end if
+      call zeroKCalc%resonanceRecon(grid, usergrid, expData, wantDeriv)
+
+
+
 
       if (kcros.eq.6) then
           call zeroKCalc%Fix_Eta(grid, expData, covData, wantDeriv) ! calculate eta if needed
@@ -192,7 +197,7 @@ module mthe0_M
       end if
       if (.not.calcDataSelfInit) then
           call calcDataSelf%initialize()
-          calcDataSelfInit = .true.
+          calcDataSelfInit = .true.         
       end if
 
       ! initialize the global zeroZCalc object
diff --git a/sammy/src/xct/XctCrossCalcImpl_M.f90 b/sammy/src/xct/XctCrossCalcImpl_M.f90
index 22dd3878d..484b254ff 100644
--- a/sammy/src/xct/XctCrossCalcImpl_M.f90
+++ b/sammy/src/xct/XctCrossCalcImpl_M.f90
@@ -17,27 +17,93 @@ module XctCrossCalcImpl_M
 
   type, extends(XctCrossCalc) :: XctCrossCalcImpl      
      contains
+     procedure, pass(this) :: setRange => XctCrossCalcImpl_setRange   ! set the range and reserve enough  space in crossData (the latter for efficiency only)
      procedure, pass(this) :: calcCross => XctCrossCalcImpl_calcCross
      procedure, pass(this) :: setEnergyIndependent => XctCrossCalcImpl_setEnergyIndependent  ! set energy independent values using current parameter values
   end type XctCrossCalcImpl
 
 contains
 
+  subroutine XctCrossCalcImpl_setRange(this, emin, emax, num, solve)
+     use exploc_common_m, only : A_Iccoul, A_Idcoul
+     use fixedi_m, only : Nnnsig, Ndasig
+     class(XctCrossCalcImpl) :: this
+     real(kind=8)::emin,emax
+     integer::num
+     logical::solve
+
+     integer::itot, irow, ii
+
+      call this%crossData%setNnsig(Nnnsig)      
+      call CrossSectionCalculator_setRange(this, emin, emax, num, solve)
+      itot = 1
+      if (solve.or.this%covariance%getPupedParam().gt.0) then
+         itot = this%covariance%getNumTotalParam()
+       end if
+
+       ! an error in sammy that does not set these right
+       this%crossSelfWhy  = .false.
+       do ii = Ndasig + 1, itot
+          this%crossSelfWhy(ii) = .true.  ! these should be parameter that are not shared
+       end do
+
+      if (this%IfCoul.GT.0) THEN
+         ii  = 2*this%Ntotc*this%resData%getNumSpinGroups()*num
+         call allocate_real_data(A_Iccoul, ii)
+         ii = ii * this%covariance%getNumTotalParam()
+         call allocate_real_data(A_Idcoul, ii)
+      end if
+
+      if (this%wantSelfIndicate) then
+         call this%crossDataSelf%nullify()
+         call this%crossDataSelf%setNnsig(1)
+         call this%crossDataSelf%reserve(num, itot+1)  
+      end if
+  end subroutine
+
   subroutine XctCrossCalcImpl_calcCross(this, ener, Ipoten)
      use xct2_m
-     use varyr_common_m, only : Su, Squ
+     use xct3_m
+     use xct4_m
+     use mxct26_m
+     use mxct18_m
+     use fixedr_m, only : Elzero, Tttzzz
+     use ifwrit_m, only : Nd_Xct
+     use exploc_common_m, only : A_Iccoul, A_Idcoul
      class(XctCrossCalcImpl) ::this
      real(kind=8)::ener
      integer::ipoten
 
+     integer::itot, ngroup
 
+     call CrossSectionCalculator_calcCross(this, ener, Ipoten)
 
-     real(kind=8)::val
-     integer::i, itot, ns, is
-     logical(C_BOOL)::accu
+     this%Etz = 0.0d0
+     IF (this%Itzero.GT.0) this%Etz = dabs(ener)/Tttzzz/Elzero
 
-     call CrossSectionCalculator_calcCross(this, ener, Ipoten)
+     this%Elz = 0.0d0
+     IF (this%ILzero.GT.0) this%Elz = this%enerSq/Elzero
+
+     IF (this%hasDirectCapture) CALL Find_Drcpt (this, dAbs(ener))
+
+     ! make cross section data
+     IF (Nd_Xct.NE.0) THEN
+        CALL N_D_Zcross (this)
+     ELSE
+        CALL Zcross (this, Ipoten)
+     END IF
 
+     ! ************ Store Coul if needed
+     IF (this%IfCoul.GT.0) THEN
+         itot = 0
+         if (this%wantDerivs) itot = this%covariance%getNumTotalParam()
+         ngroup = this%resData%getNumSpinGroups()
+         CALL Store_Coul (A_Iccoul , A_Idcoul ,   &
+           this%angInternal, this%Ccoulx , this%row , itot, this%ntotc, ngroup )
+     END IF
+
+     ! ************    Set the particular cross sections needed for this run
+     CALL Zwhich (this)
   end subroutine
 
   subroutine XctCrossCalcImpl_setEnergyIndependent(this, reactType, Twomhb, kwcoul, Etac)
@@ -47,7 +113,7 @@ contains
       use exploc_common_m, only : A_Iprmsc, I_Iflmsc, I_Ijkmsc
       use par_parameter_names_common_m, only : Nammsc
       use ifwrit_m, only : Kadddc
-      use fixedi_m, only : Nummsc
+      use fixedi_m, only : Nummsc      
       class(XctCrossCalcImpl) :: this
       integer::kwcoul, reactType
       real(kind=8)::Twomhb, Etac
@@ -86,7 +152,6 @@ contains
          call Babb (this, .true.)
          CALL Babbga (this, kwcoul )
       end if
-
   end subroutine
 
 
diff --git a/sammy/src/xct/XctCrossCalc_M.f90 b/sammy/src/xct/XctCrossCalc_M.f90
index b720be945..b2d1d8204 100644
--- a/sammy/src/xct/XctCrossCalc_M.f90
+++ b/sammy/src/xct/XctCrossCalc_M.f90
@@ -18,9 +18,8 @@ module XctCrossCalc_M
       logical,allocatable,dimension(:)::crossSelfWhy  ! reproduce a SAMMY bug for self-indication experiments. To Do fix the bug instead
 
       integer::lllmax = 0                  ! maximum number of Clebsch-Gordon coefficients
-      integer::Kfinit = 0                  ! finite-size corrections flag
       logical::wantSelfIndicate = .false.  ! do we need to calculate self-indicated cross section data
-      integer::Kssmsc                      ! is multiple scattering desired latter in the code (affects which 0K data are needed)
+
       integer::Kpolar = 0                  ! are fission width given in polar form
 
       integer::ntotc    ! the maximum number of channels (excluding gamma width) in any spin group
@@ -87,7 +86,8 @@ module XctCrossCalc_M
 
      real(kind=8),allocatable,dimension(:,:)::Ccoulx
 
-
+     real(kind=8):: Etz   ! partial sqrt(E) wrt  tzero, energy dependent
+     real(kind=8):: Elz   ! partial sqrt(E) wrt eLzero, energy dependent
 
      ! direct capture arrays. Todo update the direct capture to read and store data in C++ class
      logical::hasDirectCapture
@@ -118,8 +118,8 @@ subroutine XctCrossCalc_setAddtionalParams(this,  lllmax, Kfinit, wantSelfIndica
      class(XctCrossCalc)::this
      integer::lllmax            ! maximum number of Clebsch-Gordon coefficients
      integer::Kfinit            ! finite-size corrections flag
+     integer::Kssmsc            ! flag for scattering correction
      logical::wantSelfIndicate  ! do we need to calculate self-indicated cross section data
-     integer::Kssmsc
      integer::numIso
      logical::debug
 
@@ -136,15 +136,14 @@ subroutine XctCrossCalc_setAddtionalParams(this,  lllmax, Kfinit, wantSelfIndica
          end if
      end if
 
-     this%Kfinit  = Kfinit
+
      if (Kfinit.ne.0) this%Ifcros = .true. ! need all cross sections
      this%lllmax = lllmax
-     this%Kssmsc = Kssmsc
 
      this%Kslow = 1
      this%C_G_Kxlmn = 1
      if (this%lllmax.gt.0) then
-        IF (this%reactType.EQ.7 .OR. this%reactType.EQ.11 .OR. this%Kssmsc.GT.0) THEN
+        IF (this%reactType.EQ.7 .OR. this%reactType.EQ.11 .OR. Kssmsc.GT.0) THEN
            if (.not.allocated(this%Alj)) then
                ng = this%resData%getNumSpinGroups()
                ntotTrig = (this%ntotc*(this%ntotc+1))/2
diff --git a/sammy/src/xct/mxct0.f90 b/sammy/src/xct/mxct0.f90
index 0eeff0472..44093febf 100644
--- a/sammy/src/xct/mxct0.f90
+++ b/sammy/src/xct/mxct0.f90
@@ -5,34 +5,19 @@ module xct_m
 !
       Subroutine Samxct_0(xct)
 !
-      use fixedi_m, only : Ifdif, Iq_Val, Kiniso, Kkkiso, Kkxlmn, Lllmax, Ncrsss, Ndasig, Ndbsig, Niniso, Nnniso, Nrext, Nrfil3, &
-                           Ntotc, Ntriag, numcro, Npfil3
-      use ifwrit_m, only : Ifcoul, Kadddc, Ks_Res, Ksolve, ktzero, Nd_Xct, Nnpar, Nnparx
+      use fixedi_m, only : Iq_Iso, Iq_Val, Kiniso, Kkkiso, Nnniso, Niniso
+      use ifwrit_m, only : Ifcoul, Ks_Res, Ksolve, Nd_Xct, Nnpar
       use exploc_common_m
       use array_sizes_common_m
       use oopsch_common_m, only : Nowwww, Segmen
-      use EndfData_common_m
-      use AuxGridHelper_M
-      use SammyGridAccess_M
-      use AllocateFunctions_m
-      use rsl7_m
-      use xct1_m
-      use xct2_m
-      use mxct27_m
-      use mcro8_m
-      use mxct25_m, only : Kount_Clbsch, Clbsch
-      use mxct30_m, only : Clbsch_Slow
+      use EndfData_common_m, only : covData
+      use rsl7_m, only : Figure_Kws_Xct, Set_Kws_Xct
       use SammyLptPrinting_m     
-      use XctCrossCalc_M
-      use array_sizes_common_m, only : calcData
-      IMPLICIT none
-      type(SammyGridAccess)::grid, auxGrid
+      use XctCrossCalc_M      
+      IMPLICIT none      
       character(len=80)::line
-      integer::iflagMatch
       integer::Idimen
-      integer::K_Coul_N, Lllmmm, Mxany, N, Neight, Nfive1, Nfive1x, Nfive2, Nfive3, Nfive3x, Nfive4, Nfour
-      integer::Ifinal, Krext, Nfour1, Nfour2, Ng, Nnine, Nsix, Nthr1, Nthr2, Nthr3, ntwo1, numElAux, Nw1, NcrssxO
-      integer::Nfprrr
+      integer:: Low, I
       class(XctCrossCalc)::xct   ! temporarily here so that energy indepdent code can move in steps
       external Idimen
 !
@@ -44,21 +29,7 @@ module xct_m
       Segmen(2) = 'C'
       Segmen(3) = 'T'
       Nowwww = 0
-      NcrssxO = 0
-      if (any(xct%Ifcros)) NcrssxO = 1
-
-      call grid%initialize()
-      call grid%setParameters(numcro, ktzero)
-      call grid%setToExpGrid(expData)
-      call auxGrid%initialize()
-      call auxGrid%setParameters(numcro, ktzero)
-      call auxGrid%setToAuxGrid(expData)
-      call setAuxGridOffset(1)
-      call setAuxGridRowMax(0)
-      numElAux = auxGrid%getNumEnergies(expData)
 
-      iflagMatch = radFitFlags%matchFitFlag()
-      CALL Initil
 !
       IF (Nd_Xct.NE.0) THEN
          WRITE (6,99998)
@@ -71,304 +42,41 @@ module xct_m
       END IF
 !
 !
-! *** organize for broadening
-      Nnparx = Nnpar
-      Nnpar  = Ndasig + Ndbsig
-!
-!
-! *** find max number of res for any one spin group
-      Mxany = resParData%getNumResonances()
-!
-! *** find Nrfil3, Npfil3 for added cross sections (endf/b-vi file 3)
-!     Now done in C++
-      Npfil3 = 1
-      Nrfil3 = 1
-!
-! *** Count how many non-zero elements are in Xlmn
-      Kkxlmn = xct%C_G_Kxlmn
-!
-!
 ! *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-XCT
       Ifcoul = xct%IfCoul
-      CALL Estxct (xct, Ntwo1, Nthr1, Nthr2, Nthr3, Nfour, Nfour1, Nfour2, &
-         Nfive1, Nfive2, Nfive3, Nfive4, Nfive1x, Nfive3x, Nsix, Neight, &
-         Nnine, Nw1, Mxany, Nfprrr, K_Coul_N, numElAux, NcrssxO)
-!
-! *** Zero ***
-      N = Ndasig
-!
-      IF (Ifcoul.NE.0) THEN
-         call make_A_Iccoul(K_Coul_N)
-         IF (Nnpar.GT.0) THEN
-            call make_A_Idcoul(K_Coul_N*Nnpar)
-         END IF
-      END IF
-!
+
+      CALL Figure_Kws_Xct (Low)
+
+      Iq_Iso = xct%numIso
+
+      I = Idimen (0  ,  0, '          0     ,  0')
 !
 ! *** one ***
       CALL Set_Kws_Xct
 !
-!     Derivative list is set up Work to
-!     allow for combined data values that are not
-!     in the resonance data yet and therefore
-!     are not yet added by setUpList
-!
-
-      Ifinal = Idimen (1, 1, 'Ifinal    1, 1')
-!     Setting Ifinal here allows us to use the "Vsigxx" space, but then
-!        come back to it later via "I = Idimen (Ifinal, -1)"
-!
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
       IF (Ksolve.NE.2 .OR. covData%getPupedParam().GT.0) THEN
-         Ks_Res = 0
-         Krext = Nrext
-         IF (Nrext.EQ.0) Krext = 1
+         Ks_Res = 0        
+         Nnpar = covData%getNumTotalParam()
       ELSE
          Ks_Res = 2
+         Nnpar = 0
       END IF
-!
-! *** Two ***
-!
-! - - - - - - - - - - - - - - - - <      
-! - - - - - - - - - - - - - - - - >
-!
-! *** four ***
-!
-! *** five ***
-!
-!
-! - - - - - - - - - - - - - - - - - - - - - - <
-! *** six ***
-!                                Nsix = Nfpres
-! - - - - - - - - - - - - - - - - <
-! CALL Abpart, Abpga
-! - - - - - - - - - - - - - - - - >
-! *** eight ***
-!
-! *** nine ***
-      IF (Ifdif.NE.0.and.2*Ntriag.gt.N) then
-         N = 2*Ntriag
-      ENd if
-      N  = Ntotc
-! CALL Sinsix
-! - - - - - - - - - - - - - - - - <
-! *** eleven ***
-! CALL Setr
-! CALL Yinvrs
-! CALL Setxqx
-! - - - - - - - - - - - - - - - - >
-!
-! *** twelve ***
-      N = (Ncrsss-2)
-      IF (N.EQ.0) N = 1
-      IF (NcrssxO.EQ.0) N = 1
-      if (n.lt.ntotc) n = ntotc
-! CALL Sectio
-! CALL Setqri
-! CALL Settri
-! CALL Derres
-! CALL Dercap
-! CALL Dereee
-! CALL Derext
-! CALL Setpxr
-! CALL Derrho
-! CALL Derrad
-! CALL Derwid
-! CALL Dertze
-! CALL Derzzz
-! CALL Deriso
-! - - - - - - - - - - - - - - - - - - - - - - >
-!
-!
-      Lllmmm = Lllmax
-      IF (Lllmax.EQ.0) Lllmmm = 1
-      CALL Work (    xct, calcData , calcDataSelf,  Lllmmm)
-! *** SBROUTINE Work generates theory and derivatives
-! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
-!
-!
+
       Ks_Res = Ksolve
       Ksolve = 0
       IF (Ks_Res.EQ.2) THEN
          IF (covData%getPupedParam().eq.0) Ksolve = 2
       END IF
 !
-      Ifinal = Idimen (Ifinal, -1, '          Ifinal, -1')
-      Nnpar = Nnparx
       Kiniso = Kkkiso
       Niniso = Nnniso
       IF (Iq_Val.NE.0) Niniso = Iq_Val
 
-      call grid%destroy()
-      call auxgrid%destroy()
       RETURN
 !
       END
-!
-!
-! --------------------------------------------------------------
-!
-      SUBROUTINE Estxct (xct, Ntwo1, Nthr1, Nthr2, Nthr3, Nfour, Nfour1, &
-         Nfour2, Nfive1, Nfive2, Nfive3, Nfive4, Nfive1x, Nfive3x, &
-         Nsix, Neight, Nnine, Nw1, Mxany, Nfprrr, K_Coul_N, &
-         numElAux, Ncrssx)
-!
-! *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-XCT
-!
-      use fixedi_m, only : Ifdif, Iq_Iso, Iq_Val, Kkxlmn, Lllmax, Nres, Mres, Ncrsss,  &
-                           Ndasig, Ndbsig, Nfpres, Nnniso, Npfil3, Nrfil3, Ntotc, Ntriag,      &
-                           Kshift
-      use ifwrit_m, only : Ifcoul, Kcros, Kpiece, Ksolve, Kssmsc, Nd_Xct,Nnpar
-      use XctCrossCalc_M
-      use EndfData_common_m, only : resParData, radFitFlags
-      use rsl7_m, only : Figure_Kws_Xct
-      IMPLICIT none
-      class(XctCrossCalc)::xct   ! temporarily here so that energy indepdent code can move in steps
-      integer::Ntwo1, Nthr1, Nthr2, Nthr3, Nfour, Nfour1, &
-               Nfour2, Nfive1, Nfive2, Nfive3, Nfive4, Nfive1x, Nfive3x, &
-               Nsix, Neight, Nnine, Nw1, Mxany, Nfprrr, K_Coul_N, &
-               numElAux, Ncrssx
-       integer::Idimen
-       integer::I,K, K1, K2, K3,K4, K6, K7, K8, Ke, Kw, low, N,Nthr4
-       external Idimen
-!
-! *** Andy
-      K = 2*(Npfil3+Nrfil3)
-!
-! *** one
-      CALL Figure_Kws_Xct (Low)
-      K = resParData%getNumSpinGroups()*2 + K
-!
-! *** Two
-      K1 = Mres
-      Ntwo1 = 1
-      IF (Kshift.NE.0) Ntwo1 = Mres
-      K2 = Ntwo1
-!
-! *** three
-      Iq_Iso = Nnniso
-      IF ( (Kcros.EQ.7 .OR. Kcros.EQ.11) .OR. Kssmsc.GT.0) THEN
-         IF (Iq_Val.GT.Iq_Iso) Iq_Iso = Iq_Val         
-         Nthr1 = Lllmax*Iq_Iso
-         Nthr2 = Nthr1*Ndasig
-         Nthr3 = Kkxlmn
-         IF (Nthr1.EQ.0) Nthr1 = 1
-         IF (Nthr2.EQ.0) Nthr2 = 1
-         IF (Nthr3.EQ.0) Nthr3 = 1
-         N = Nthr1 + Nthr2 + 2*Nthr3
-         K3 = N
-      ELSE
-         Nthr1 = 1
-         Nthr2 = 1
-         Nthr3 = 1
-         Nthr4 = 1
-         K3  = 4
-      END IF
-!
-! *** four
-      Nfour = 1
-      IF (xct%wantDerivs) Nfour = Ntriag*Nfpres
-      Nfour1 = Ntriag*Ntotc*resParData%getNumResonances()
-      IF (Nfour1.EQ.0) Nfour1 = 1
-      Nfprrr = radFitFlags%getNumTotalVaried()
-      IF (Nfprrr.EQ.0) Nfprrr = 1
-      Nfour2 = Ntriag*Nfprrr*resParData%getNumSpinGroups()
-      IF (Nfour2.EQ.0) Nfour2 = 1
-      K4 = Nfour*4 + Nfour1 + 2*Nfour2
-!
-! *** five
-      Nfive1 = Ncrsss*resParData%getNumSpinGroups()
-      IF (Ncrssx.EQ.0) Nfive1 = 1
-      IF (Ksolve.NE.2 .AND. Ndasig+Ndbsig.NE.Nnpar) THEN
-         STOP '[STOP in Estxct in xct/mxct0.f]'
-      END IF
-      Nfive2 = Nfive1*(Ndasig+Ndbsig)
-      IF (Nfive2.EQ.0) Nfive2 = 1
-      IF (Ncrssx.EQ.0) Nfive2 = 1
-      IF (Kcros.NE.7 .AND. Kcros.NE.11 .AND. Kssmsc.EQ.0) THEN
-         Nfive3 = 1
-         IfCoul = 0
-      ELSE
-         Nfive3 = 2*Ntotc*Ntotc*resParData%getNumSpinGroups()         
-      END IF
-      Nfive4 = Nfive3*Nnpar
-      IF (Nfive1.EQ.0) Nfive1 = 1
-      IF (Nfive2.EQ.0) Nfive2 = 1
-      IF (Nfive3.EQ.0) Nfive3 = 1
-      IF (Nfive4.EQ.0) Nfive4 = 1
-      Nfive1x = 1
-      Nfive3x = 1
-      IF (Nd_Xct.EQ.1) THEN
-         Nfive1x = Nfive1
-         Nfive3x = Nfive3
-      END IF
-      K4 = K4 + Nfive1 + Nfive2 + Nfive3 + Nfive4 + Nfive1x + Nfive3x
-!
-! *** six
-      nsix = 0
-      IF (xct%wantDerivs) Nsix = Nfpres
-      IF (Nsix.EQ.0) Nsix = 1
-      K6 = 3*Mres + Nsix
-!
-! *** seven
-      K7 = 2*Mres + 2*Nsix
-!
-! *** eight
-      Neight = 1
-      IF (Kpiece.EQ.1) Neight = numElAux*resParData%getNumSpinGroups()
-      IF (Neight.EQ.0) Neight = 1
-      K8 = Neight
-!
-! *** nine
-      IF (Ifdif.EQ.0) THEN
-         Nnine = 1
-      ELSE
-         Nnine = 2*Ntriag
-      END IF
-      K8 = K8 + 5*Ntotc + Nnine
-!
-! *** ten
-      K8 = K8 +  3*Ntotc + 2*Ntriag + 2*Ntotc*Ntotc + 2*Ntriag
-!
-! *** eleven
-      Ke = 4*Ntriag
-!
-! *** twelve
-      N = Ncrsss - 2 + 2*Ncrsss*Ntriag
-      IF (Ncrssx.EQ.0) N = 3
-      Kw = N + 2*Ntriag + 4*Ntriag**2
-      Kw = Kw + 3*Ntotc + Ntotc**2 + 4*Ntriag
-      Kw = Kw + 2*Ntotc + 2*Ntotc**2 + &
-           6*Ntriag*resParData%getNumSpinGroups() + Ncrsss
-      Nw1 = 1
-      IF (Kpiece.EQ.1) Nw1 = numElAux
-!
-      K8 = K8 + MAX0(Ke,Kw)
-      K6 = K6 + MAX0(K7,K8)
-      K4 = K4 + MAX0(K6,Nw1)
-      K3 = K3 + K4
-      K1 = K1 + MAX0(K2,K3)
-      K = K + K1
-!
-! *** ZERO
-      IF (Ifcoul.NE.0) THEN
-         K_Coul_N = 2*Ntotc*resParData%getNumSpinGroups()*numElAux
-         K = K + K_Coul_N
-         N = Ndasig
-         IF (N.GT.0) K = K + K_Coul_N *N
-      ELSE
-         K_Coul_N = 1
-      END IF
-!
-! *** Alltogether
-      Low = Low + K
-! Previously this was a check that we have enough array space
-! But now many of the arrays have moved to allocatable array
-! and the test is no longer meaningful
-      I = Idimen (0  ,  0, '          0     ,  0')
-!
-      RETURN
-      END
+
 !
 !
 ! --------------------------------------------------------------
diff --git a/sammy/src/xct/mxct02.f90 b/sammy/src/xct/mxct02.f90
index 5d958f92b..5f909089d 100644
--- a/sammy/src/xct/mxct02.f90
+++ b/sammy/src/xct/mxct02.f90
@@ -1,226 +1,7 @@
 !
 module xct2_m
   contains
-!
-! --------------------------------------------------------------
-!
-      SUBROUTINE Work (calc, derivs,   derivsSelf,  Lllmmm)
-!
-! *** PURPOSE -- Generate theoretical cross sections "theory" and partial
-! ***            Derivatives "dasig"
-!
-      use fixedi_m
-      use ifwrit_m
-      use exploc_common_m
-      use samxxx_common_m
-      use oopsch_common_m
-      use fixedr_m
-      use varyr_common_m, only : Elz, Etz, Su, Squ
-      use cbro_common_m
-      use lbro_common_m
-      use EndfData_common_m
-      use SammyGridAccess_M
-      use DerivativeHandler_M
-      use xct3_m
-      use xct4_m
-      use convert_to_transmission_m
-      ! use xct6_m
-      use mxct27_m
-      use mxct26_m
-      use mxct06_m
-      use mxct18_m
-      use Zgauss_m
-      use XctCrossCalc_M
-      IMPLICIT none
-
-      integer(4), intent(in):: Lllmmm
-      real(8):: Zero, A, Gbx, Theoryx
-      integer(4):: Jdat, Idrcp, Ipoten, Iw, irow, istart, &
-                   ng, numEl, TotalNdasig
-      integer(4) :: Iipar, iso, Jcount, Jsig,isoReal
-      real(8) :: val
-
-      LOGICAL Ywhich
-      type(SammyGridAccess)::grid
-      type(DerivativeHandler)::derivs, derivsSelf
-      integer::iflagMatch, isoOur, is
-      logical::wantNegative, wantDeriv
-      class(XctCrossCalc)::calc
-      logical(C_BOOL)::accu
-
-!
-!      DIMENSION W...(...Ndatb)
-!
-!
-      DATA Zero /0.0d0/
-!
-!
-      call grid%initialize()
-      call grid%setParameters(numcro, ktzero)
-      call grid%setToAuxGrid(expData)
-      numEl = grid%getNumEnergies(expData)
-
-      ! if using Leal-Hwang doppler broadening (Kkkdop.eq.1)
-      ! don't need to extend to negative energies
-      wantNegative = .true.
-      if (Kkkdop.eq.1) wantNegative = .false.
-
-      accu = .true.
-      call derivs%setAccumulate(accu)
-      call derivsSelf%setAccumulate(accu)
-      calc%crossSelfWhy  = .false.
-      do is = Ndasig + 1, Ndbsig
-         calc%crossSelfWhy(is) = .true.  ! these should be parameter that are not shared
-      end do
-
-      Ywhich = Ydoppr.OR.Yresol.OR.Yangle.OR.Yssmsc.OR.Yaverg.OR. Maxwel.EQ.1 .OR. Knocor.EQ.1
-
-      wantDeriv = Ksolve.NE.2
-      if ( covData%getPupedParam().gt.0) wantDeriv = .true.
 
-!
-      Iw = 0
-      IF (Ksindi.GT.0 .AND. Kcros.EQ.8) THEN
-         Iw = 1
-      END IF
-      call derivs%nullify()
-      call derivs%setNnsig(Nnnsig)
-      if (calc%reactType.eq.6) then
-          if( Nnnsig.ne.1) then
-              write(6,*)" Expected Nnsig to be 1 if calculating eta"
-              write(21,*)" Expected Nnsig to be 1 if calculating eta"
-              stop
-          end if
-          call calc%crossData%setNnsig(2)
-      end if
-
-      call derivs%reserve(numEl * calc%crossData%getNnnsig(), Ndasig + Ndbsig + 1)
-
-      IF (Iw.EQ.1.or.Ksitmp.gt.0) THEN
-         call derivsSelf%nullify()
-         call derivsSelf%setNnsig(1)
-         call derivsSelf%reserve(numEl, Ndasig + Ndbsig + 1)
-      end if
-!      
-      irow = 0
-      istart = 0
-!
-! *** (ordering is as variable # is assigned)
-!
-!
-! *** Start Coulomb if needed
-      IF (IfCoul.GT.0) THEN
-         ng = resParData%getNumSpinGroups()
-         CALL Zero_Array (A_Iccoul , 2*Ntotc*Ng*numEl)
-         IF (Nnpar.GT.0) CALL Zero_Array (A_Idcoul , &
-                             2*Ntotc*Ng*numEl*Nnpar)
-      END IF
-!
-      Idrcp = 1
-! *** Calculate cross sections and derivatives for all
-! ***       positive energies
-      DO Jdat=1,numEl
-!
-
-         IF ((Jdat/100000)*100000.EQ.Jdat) WRITE (6,10000) Jdat
-10000    FORMAT ('     *** on data point number', I10)
-         Su = grid%getEnergy(Jdat, expData)
-
-         if (Su.le.0) then
-            if (.not.wantNegative) then
-                if (Ywhich) then
-                  irow = irow + 1
-                  do iso = 1, Nnnsig
-                    call derivs%setToZero((irow-1)*Nnnsig+iso, Ndasig+ndbsig+1)
-                    if(Iw.eq.1.or.Ksitmp.gt.0) THEN
-                        call derivsSelf%setToZero((irow-1)*Nnnsig+iso, Ndasig+ndbsig+1)
-                    end if
-                  end do
-                end if
-
-                cycle
-            end if
-         end if
-         irow = irow + 1
-         calc%row = irow
-         calc%ener = grid%getEnergy(Jdat, expData)         
-
-!
-            IF (Su.GT.Emax .AND. Kartgd.EQ.1) THEN
-               ! if Kartgd.EQ.1 and more than one isotope
-               ! SAMMY did not zero the cross section for isotopes 2, 3, ...
-               ! but reused the previous value
-               ! todo: delete this to make the code correct
-               do Iso=2,Iq_Iso
-                  do is = 1, Nnnsig
-                     val = derivs%getDataNs(irow-1, is, 0, iso)
-                     call derivs%addDataNs(irow, is, 0, iso, val)
-                  end do
-               end do
-               GO TO 20
-            END IF
-            IF (Su.LT.Zero) Su = - Su
-
-            IF (Kadddc.NE.0) CALL Find_Drcpt (calc, Su)
-!
-! ********* Start regular calculation
-            Squ = dSQRT(Su)
-            calc%enerSq = Squ
-! ********* Su  = E = m Dist^2 Elzero^2 / [2 (t-Tzero)^2]
-! ********* Squ = Tttzzz * Dist * Elzero / (t-Tzero)
-! *********       Tttzzz = sqrt(m/2) * conversion factors = 72.3
-!
-            IF (Itzero.GT.0) Etz = Su/Tttzzz/Elzero
-! ********* Etz = sqrt(E)/(t-Tzero) for use in Dereee & elsewhere
-! *********     = (partial of sqrt(E)) wrt (Tzero)
-!
-            IF (ILzero.GT.0) ELz = Squ/Elzero
-! ********* eLz = sqrt(E)/Elzero = (partial of sqrt(E)) wrt (Elzero)
-!
-! ********* Write potential scattering if desired
-            Ipoten = 0
-            IF ((Jdat.EQ.1 .OR. Jdat.EQ.numEl) .AND. Kpoten.NE.0) &
-               CALL Write_Potential (Su, Ipoten)
-!
-            IF (Kgauss.EQ.1) THEN
-! ************ Want dummy Gaussian resonances
-               CALL Zgauss (resparData, val, grid%getEnergy(Jdat, expData))
-               call calc%crossData%addDataNs(calc%row, 1, 0, 1, val)
-            ELSE
-!
-! ************ Generate cross sections and derivatives
-               IF (Nd_Xct.NE.0 .AND. Ksolve.NE.2) THEN
-                  CALL N_D_Zcross (calc)
-               ELSE
-                  CALL Zcross (calc, Ipoten)
-               END IF
-!
-! ************ Store Coul if needed
-               IF (IfCoul.GT.0) THEN
-                  CALL Store_Coul (A_Iccoul , A_Idcoul ,  &
-                     calc%angInternal, calc%Ccoulx , Jdat)
-               END IF
-!
-! ************    Set the particular cross sections needed for this run
-               CALL Zwhich (calc)
-            END IF
-!
-   20       CONTINUE
-
-!
-      END DO
-
-! *** end of do-loop on energies (Jdat)
-!
-!
-      call grid%destroy()
-      accu = .false.
-      call derivs%setAccumulate(accu)
-      call derivsSelf%setAccumulate(accu)
-
-!
-      RETURN
-      END
 !
 !
 ! ______________________________________________________________________
diff --git a/sammy/src/xct/mxct03.f90 b/sammy/src/xct/mxct03.f90
index 00cddd376..b8eaa292c 100644
--- a/sammy/src/xct/mxct03.f90
+++ b/sammy/src/xct/mxct03.f90
@@ -13,15 +13,11 @@ module xct3_m
 !
 ! *** PURPOSE -- Calculate numerically the partial derivatives 
 ! ***               of the cross section wrt R-matrix parameters
-!
-      use fixedi_m, only : Kpolar
-      use exploc_common_m
-      use SammyResonanceInfo_M
+!   
       use RMatResonanceParam_M
       use xct1_m
       use xct5_m
-      use mxct06_m      
-      use AllocateFunctions_m
+      use mxct06_m
       IMPLICIT none
 !
 !
@@ -123,7 +119,7 @@ module xct3_m
 ! ***                GAMMA-SUB-N (PARTICLE WIDTHS)
                      M2 = M - 2
                      ichan = spinInfo%getWidthForChannel(M)
-                     IF (M2.EQ.1 .OR. Kpolar.NE.1) THEN
+                     IF (M2.EQ.1 .OR. calc%Kpolar.NE.1) THEN
                         B = resonance%getWidth(ichan)
                         call resonance%setWidth(ichan, &
                                      B*(One+U_Increment))
@@ -180,8 +176,6 @@ module xct3_m
 !
 ! *** PURPOSE -- Calculate numerically the partial derivatives 
 !
-      use fixedi_m
-      use ifwrit_m
       IMPLICIT none
 !
       class(XctCrossCalc)::calc
@@ -190,8 +184,8 @@ module xct3_m
       real(kind=8)::X, Two_X
       integer::K, Nchan, Nchanx
 !
-      IF (Ncrssx.GT.0) THEN
-         DO K=1,Ncrsss
+      IF (any(calc%Ifcros)) THEN
+         DO K=1,calc%ntotc+1
             un = unpertCross(k, Igrp)
             per =  calc%crossInternal(K, igrp, 0)
             val = - (un-per)/X
diff --git a/sammy/src/xct/mxct04.f90 b/sammy/src/xct/mxct04.f90
index 5041b993c..2812fe918 100644
--- a/sammy/src/xct/mxct04.f90
+++ b/sammy/src/xct/mxct04.f90
@@ -10,21 +10,14 @@ module xct4_m
 ! ***            AND THE ( PARTIAL DERIVATIVES OF THE CROSS SECTION
 ! ***            WITH RESPECT TO THE VARIED PARAMETERS ) = Deriv
 !
-      use over_common_m
-      use oops_common_m
-      use fixedi_m
-      use ifwrit_m
-      use exploc_common_m
-      use EndfData_common_m
       use xct5_m
       use mxct06_m
       use SammySpinGroupInfo_M
       use XctCrossCalc_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      logical::doRadDeriv
-      integer::igr, ich
-      type(SammySpinGroupInfo)::spinInfo
+      IMPLICIT None
+
       class(XctCrossCalc)::calc
+      integer::ipoten
 !
 !
 ! *** Generate Pr and Pi = Partial of R wrt U-parameters
@@ -35,34 +28,7 @@ module xct4_m
 ! *** Generate Pgar & Pgai = partial of R wrt (Gamma-x) *
 ! ***                        partial (Gamma_x) wrt radius
 !
-      doRadDeriv = .false.
-      do igr = 1, resparData%getNumSpinGroups()
-         call resParData%getSpinGroupInfo(spinInfo, igr)
-         DO Ich=1,spinInfo%getNumChannels()
-            ii = radFitFlags%getTrueFitFlag(Igr, Ich)
-            IF (ii.NE.0) THEN
-               if (Ksolve.ne.2) then
-                   doRadDeriv = .true.
-               else if (covData%isPupedParameter(ii)) then
-                  doRadDeriv = .true.
-               end if
-            end if
-            if (doRadDeriv) exit
-            ii = radFitFlags%getEffFitFlag(Igr, Ich)
-            IF (ii.NE.0) THEN
-               if (Ksolve.ne.2) then
-                   doRadDeriv = .true.
-               else if (covData%isPupedParameter(ii)) then
-                  doRadDeriv = .true.
-               end if
-            end if
-         end do
-      end do
-      if (Ks_Res.eq.2.or.resParData%getNumResonances().eq.0) then
-         doRadDeriv = .false.
-      end if
-
-      IF (doRadDeriv) THEN
+      IF (calc%Ifrad) THEN
           CALL Abpga (calc)
       END IF
 !
diff --git a/sammy/src/xct/mxct05.f90 b/sammy/src/xct/mxct05.f90
index 23c178df2..8a204e85c 100644
--- a/sammy/src/xct/mxct05.f90
+++ b/sammy/src/xct/mxct05.f90
@@ -11,14 +11,9 @@ module xct5_m
 ! ***            Also generate Pr and Pi = partial of R wrt U-parameters
 ! ***            And Prer & Prei = partial of R wrt E
 !
-      use fixedi_m
-      use ifwrit_m
-      use varyr_common_m, only : Su
-      use EndfData_common_m
       use SammyResonanceInfo_M
       use RMatResonanceParam_M
       use SammySpinGroupInfo_M
-      !use templc_common_m
       use XctCrossCalc_M
       IMPLICIT None
 !
@@ -39,19 +34,19 @@ module xct5_m
 ! ***      AND Alphai = Gamgam/2 / ( Ditto )
 !
       calc%needAlphai = .false.
-      DO N=1,resParData%getNumResonances()
+      DO N=1,calc%resData%getNumResonances()
          calc%Xden(N) = Zero
          calc%Alphar(N) = Zero
          calc%Alphai(N) = Zero
-         call resParData%getResonanceInfo(resInfo, N)
+         call calc%resData%getResonanceInfo(resInfo, N)
          IF (.not.resInfo%getIncludeInCalc()) THEN
             calc%Difen(N) = Zero
          ELSE
             igrp = resInfo%getSpinGroupIndex()
-            call resParData%getSpinGroupInfo(spinInfo, igrp)
+            call calc%resData%getSpinGroupInfo(spinInfo, igrp)
             IF (spinInfo%getIncludeInCalc()) THEN
-               call resParData%getRedResonance(resonance, resInfo)
-               calc%Difen(N) = resonance%getEres() - Su
+               call calc%resData%getRedResonance(resonance, resInfo)
+               calc%Difen(N) = resonance%getEres() - dAbs(calc%ener)
                if (calc%doResShift) then
                   calc%Difen(N) =  calc%Difen(N) + calc%Xx(N)
                end if
@@ -120,18 +115,18 @@ module xct5_m
          END DO
       END IF
 
-      IF (Itzero.NE.0 .OR. Ilzero.NE.0) THEN
+      IF (calc%Itzero.NE.0 .OR. calc%Ilzero.NE.0) THEN
 !         
          calc%Prer = 0.0d0
          calc%Prei = 0.0d0
 
-         DO Ig=1,resParData%getNumSpinGroups()  
-            call resParData%getSpinGroupInfo(spinInfo, ig)
-            DO N=1,resParData%getNumResonances()
-               call resParData%getResonanceInfo(resInfo, N)
+         DO Ig=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinInfo, ig)
+            DO N=1,calc%resData%getNumResonances()
+               call calc%resData%getResonanceInfo(resInfo, N)
                igrp = resInfo%getSpinGroupIndex()
                IF (igrp.EQ.Ig) THEN
-                  call resParData%getRedResonance(resonance, resInfo)
+                  call calc%resData%getRedResonance(resonance, resInfo)
                   Aa = calc%Alphar(N)**2 - calc%Alphai(N)**2
                   Bb = Two*calc%Alphar(N)*calc%Alphai(N)
                   Ij = 0
@@ -170,9 +165,6 @@ module xct5_m
 ! *** PURPOSE -- generate Pgar and Pgai = partial of R wrt unvaried
 ! ***               width-parameters
 !
-      use fixedi_m
-      use ifwrit_m
-      use EndfData_common_m
       use SammySpinGroupInfo_M
       use XctCrossCalc_M
       IMPLICIT none
@@ -204,7 +196,7 @@ module xct5_m
          call calc%resData%getResonanceInfo(resInfo, ires)
          IF (resInfo%getIncludeInCalc()) THEN             
             ign = resInfo%getSpinGroupIndex()
-            call resParData%getSpinGroupInfo(spinInfo, ign)           
+            call calc%resData%getSpinGroupInfo(spinInfo, ign)
             DO Mchan=1,spinInfo%getNumChannels()
                Ifl = calc%radiusData%getTrueFitFlag(Ign, Mchan)
                IF (Ifl.GT.0) THEN
@@ -218,7 +210,7 @@ module xct5_m
                       end if
                   end do
 
-                  DO Ij=1,Ntriag
+                  DO Ij=1,calc%ntriag
                      IF (calc%Bga(IJ,Mchan,Ires).NE.Zero) THEN
                         calc%Pgar(IJ,Ix,Ign) = calc%Pgar(Ij,Ix,Ign) + &
                            calc%Bga(IJ,Mchan,Ires)*calc%Alphar(Ires)
diff --git a/sammy/src/xct/mxct11.f90 b/sammy/src/xct/mxct11.f90
index 507a3ad6f..838d9c701 100644
--- a/sammy/src/xct/mxct11.f90
+++ b/sammy/src/xct/mxct11.f90
@@ -136,7 +136,6 @@ contains
 ! *** generate Deriv(k,Itzero) & Deriv(k,Ilzero); ditto Derivx
 ! *** ie Deriv(k,i) = partial Crss(k) wrt either Tzero or LZero
 !
-      use varyr_common_m, only : Etz, Elz
       use mxct09_m
       use SammySpinGroupInfo_M
       IMPLICIT None
@@ -207,12 +206,12 @@ contains
                            end select
                            IF (Itz.GT.0) THEN
                              val = calc%angInternal(Ix,Ichanx,Ichan, Igr, Itz)
-                             val = val + A1 + A2*Etz
+                             val = val + A1 + A2*calc%Etz
                              calc%angInternal(Ix,Ichanx,Ichan, Igr, Itz) = val
                            end if
                            IF (iLz.GT.0) THEN
                               val = calc%angInternal(Ix,Ichanx,Ichan, Igr, iLz)
-                              val = val + A1 + A2*Elz
+                              val = val + A1 + A2*calc%Elz
                               calc%angInternal(Ix,Ichanx,Ichan, Igr, iLz) = val
                            end if
                         end do
diff --git a/sammy/src/xct/mxct17.f90 b/sammy/src/xct/mxct17.f90
index ebb271d4d..d96c9a0a0 100644
--- a/sammy/src/xct/mxct17.f90
+++ b/sammy/src/xct/mxct17.f90
@@ -11,7 +11,6 @@ contains
 ! *** Purpose -- Generate derivatives of cross section wrt sqrt(E) via Rho
 ! *** Note that [partial Rho wrt sqrt(E)] = effective radius (Zke is k without sqrt(E) term)
 !
-      use varyr_common_m, only : Elz, Etz
       use SammyResonanceInfo_M
       use RMatResonanceParam_M
 !
@@ -52,7 +51,7 @@ contains
             call spinInfo%getChannelInfo(channelInfo, Ichan)
             call calc%resData%getChannel(channelI, channelInfo)
             IF (Itz.GT.0) THEN               
-               Zz = Etz*channelI%getApe()*calc%Zke(Ichan, Igr)
+               Zz = calc%Etz*channelI%getApe()*calc%Zke(Ichan, Igr)
                do ix = 1, 2
                    val = calc%angInternal(Ix,Ichan,Ichan, Igr, Itz)
                    val = val + Zz*calc%Dsfx(Ix,Ichan,Ichan)
@@ -60,7 +59,7 @@ contains
                end do
             END IF
             IF (ILz.GT.0) THEN
-               Zz = ELz*channelI%getApe()*calc%Zke(Ichan, Igr)
+               Zz = calc%ELz*channelI%getApe()*calc%Zke(Ichan, Igr)
                do ix = 1, 2
                   val = calc%angInternal(Ix,Ichan,Ichan, Igr, ILz)
                   val = val + Zz*calc%Dsfx(Ix,Ichan,Ichan)
@@ -85,14 +84,14 @@ contains
                      do ix = 1, 2
                         IF (Itz.GT.0) THEN
                            val = calc%angInternal(Ix,Jchan,Ichan,Igr, Itz)
-                           val = val +  Etz*R1 * calc%Dstx(Ix,Jchan,Ichan)
-                           val = val +  Etz*R2 * calc%Dstx(Ix,Jchan,Ichan)
+                           val = val +  calc%Etz*R1 * calc%Dstx(Ix,Jchan,Ichan)
+                           val = val +  calc%Etz*R2 * calc%Dstx(Ix,Jchan,Ichan)
                            calc%angInternal(Ix,Jchan,Ichan, Igr, Itz) = val
                         end if
                         IF (ILz.GT.0) THEN
                             val = calc%angInternal(Ix,Jchan,Ichan, Igr, ILz)
-                            val = val + Elz*R1*calc%Dstx(Ix,Jchan,Ichan)
-                            val = val + Elz*R2*calc%Dstx(Ix,Jchan,Ichan)
+                            val = val + calc%Elz*R1*calc%Dstx(Ix,Jchan,Ichan)
+                            val = val + calc%Elz*R2*calc%Dstx(Ix,Jchan,Ichan)
                             calc%angInternal(Ix,Jchan,Ichan, Igr, ILz)  = val
                         end if
                      end do
@@ -114,7 +113,6 @@ contains
 ! *** Purpose -- Generate derivatives of cross section wrt sqrt(E) via Rho
 ! *** Note that [partial Rho wrt sqrt(E)] = true radius (Zke is k without sqrt(E) term)
 !
-      use varyr_common_m, only : Elz, Etz
       use SammyResonanceInfo_M
       use RMatResonanceParam_M
       use SammySpinGroupInfo_M
@@ -190,12 +188,12 @@ contains
                do ix = 1, 2
                   if (Itz.gt.0) then
                       val = calc%angInternal(Ix,Ichan,Ichan, Igr, Itz)
-                      val = val + Zz * Etz * calc%Dstx(Ix,Ichan,Ichan)
+                      val = val + Zz * calc%Etz * calc%Dstx(Ix,Ichan,Ichan)
                       calc%angInternal(Ix,Ichan,Ichan, Igr, Itz) = val
                   end if
                   if (Ilz.gt.0) then
                       val = calc%angInternal(Ix,Ichan,Ichan, Igr, Ilz)
-                      val = val + Zz * Elz * calc%Dstx(Ix,Ichan,Ichan)
+                      val = val + Zz * calc%Elz * calc%Dstx(Ix,Ichan,Ichan)
                       calc%angInternal(Ix,Ichan,Ichan, Igr, Ilz) = val
                   end if
                end do
@@ -210,15 +208,15 @@ contains
                 do ix = 1, 2
                   IF (Itz.GT.0) THEN
                      val = calc%angInternal(Ix,Jchan,Ichan, Igr, Itz)
-                     val = val + Etz * R1 * calc%Dstx(ix,Jchan,Ichan)
-                     val = val + Etz * R2 * calc%Dstx(ix,Jchan,Ichan)
+                     val = val + calc%Etz * R1 * calc%Dstx(ix,Jchan,Ichan)
+                     val = val + calc%Etz * R2 * calc%Dstx(ix,Jchan,Ichan)
                      calc%angInternal(Ix,Jchan,Ichan,Igr, Itz) = val
                      ! ??? double-counting here?
                   end if
                   IF (Ilz.GT.0) THEN
                      val = calc%angInternal(Ix,Jchan,Ichan,Igr, Ilz)
-                     val = val + Elz * R1 * calc%Dstx(ix,Jchan,Ichan)
-                     val = val + Elz * R2 * calc%Dstx(ix,Jchan,Ichan)
+                     val = val + calc%Elz * R1 * calc%Dstx(ix,Jchan,Ichan)
+                     val = val + calc%Elz * R2 * calc%Dstx(ix,Jchan,Ichan)
                      calc%angInternal(Ix,Jchan,Ichan, Igr, Ilz) = val
                   end if
                 end do
@@ -239,7 +237,6 @@ contains
 ! ***    & Elzero) * (E/4pi) that comes directly from the 1/E term (1/k**2)
 ! ***    in the formula for cross section
 !
-      use varyr_common_m, only : Elz, Etz
       use SammySpinGroupInfo_M
 !
       class(XctCrossCalc)::calc
@@ -253,11 +250,11 @@ contains
          val = calc%crossInternal(I, Igr, 0)
          A = calc%Ddddtl(I)*Dgoj - val*Two/calc%enerSq
          IF (calc%Itzero.GT.0)  then
-            val = A*Etz + calc%crossInternal(I,Igr, calc%Itzero)
+            val = A*calc%Etz + calc%crossInternal(I,Igr, calc%Itzero)
             calc%crossInternal(I,Igr, calc%Itzero ) = val
          end if
          IF (calc%Ilzero.GT.0) then
-            val = A*Elz + calc%crossInternal(I,Igr, calc%Ilzero)
+            val = A*calc%Elz + calc%crossInternal(I,Igr, calc%Ilzero)
             calc%crossInternal(I,Igr, calc%Ilzero) = val
          end if
       END DO
diff --git a/sammy/src/xct/mxct18.f90 b/sammy/src/xct/mxct18.f90
index af30e1cd1..45da6f66f 100644
--- a/sammy/src/xct/mxct18.f90
+++ b/sammy/src/xct/mxct18.f90
@@ -24,7 +24,7 @@ contains
       IF (calc%ener.EQ.Zero) RETURN
 !
       IF ((calc%reactType.EQ. 9 .OR. (calc%reactType.LT.7 .AND. (Kssmsc.EQ.0 .OR.   &
-         Kssmsc.EQ.-1))) .AND. Kaverg.NE.2) THEN
+         Kssmsc.EQ.-1))) .AND. Kaverg.NE.2) THEN        
 !
 ! ***    Here for no scattering of any kind, only one type of cross section   
          CALL Prtclr (calc)
diff --git a/sammy/src/xct/mxct21.f90 b/sammy/src/xct/mxct21.f90
index 0db8e1fc6..e2fd49046 100644
--- a/sammy/src/xct/mxct21.f90
+++ b/sammy/src/xct/mxct21.f90
@@ -4,7 +4,6 @@ IMPLICIT NONE
 
 private Jxnnn, Jxmmm
 
-! Todo: Don't use global parameters  Iq_Iso, Iq_Val
 
 contains
     integer function Jxnnn (Nb,Na,N, Ntotc, Ngroup)
@@ -26,7 +25,7 @@ contains
 ! *** Purpose -- set  coefficient of Legendre polynomial
 ! ***                         P-sub-(L-1) for Isotope Iso
 !
-      use fixedi_m, only : Iq_Iso, Iq_Val
+      use fixedi_m, only :  Iq_Val
       use SammySpinGroupInfo_M
       use mdat9_m
 
@@ -46,7 +45,7 @@ contains
 
       CALL Findpr (calc%C_G_Kxlmn, Klmn)
 !
-      DO Iq=1,Iq_Iso
+      DO Iq=1,calc%numIso
          IF (Iq_Val.NE.0) THEN
             Iso = Iso_Qv(Iq)
             C2 = Cmlab(2,Iq)
diff --git a/sammy/src/xct/mxct22.f90 b/sammy/src/xct/mxct22.f90
index 5320df1bd..e468f3fef 100644
--- a/sammy/src/xct/mxct22.f90
+++ b/sammy/src/xct/mxct22.f90
@@ -4,7 +4,6 @@ module mxct22_m
 
   private Jxnnn, Jxmmm
 
-! Todo: Don't use global parameters  Iq_Iso, Iq_Val
 contains
 
   integer function Jxnnn (Nb,Na,N, Ntotc, Ngroup)
@@ -36,7 +35,7 @@ contains
     ! ***
     ! *** Note    -- Enter this routine only if we want derivatives
     !
-    use fixedi_m, only : Iq_Iso, Iq_Val
+    use fixedi_m, only : Iq_Val
     use SammySpinGroupInfo_M
     use SammyIsoInfo_M
     use mdat9_m
@@ -58,7 +57,7 @@ contains
     ngroup = calc%resData%getNumSpinGroups()
     CALL Findpr (calc%C_G_Kxlmn, Klmn)
     !
-    DO Iq=1,Iq_Iso
+    DO Iq=1,calc%numIso
        IF (Iq_Val.NE.0) THEN
           Iso = Iso_Qv(Iq)
           C2 = Cmlab(2,Iq)
diff --git a/sammy/src/xct/mxct23.f90 b/sammy/src/xct/mxct23.f90
index c7c1480b5..e57ed9326 100644
--- a/sammy/src/xct/mxct23.f90
+++ b/sammy/src/xct/mxct23.f90
@@ -1,7 +1,7 @@
 module mxct23_m
 use XctCrossCalc_M
 implicit none
-! Note: To do: don't use the global parameter Ksindi, Numdet or Nfissl
+
 contains
 !
 !
diff --git a/sammy/src/xct/mxct26.f90 b/sammy/src/xct/mxct26.f90
index 1dfc4e858..cb9b59116 100644
--- a/sammy/src/xct/mxct26.f90
+++ b/sammy/src/xct/mxct26.f90
@@ -62,18 +62,15 @@ contains
 !
 ! ______________________________________________________________________
 !
-      SUBROUTINE Store_Coul (Ccoul, Dcoul, angInternal, Ccoulx, Jdat)
-      use fixedi_m, only : Ntotc, Ngroup
-      use ifwrit_m, only : Nnpar
-      use  EndfData_common_m     
+      SUBROUTINE Store_Coul (Ccoul, Dcoul, angInternal, Ccoulx, Jdat, Nnpar, Ntotc, Ngroup)
       real(kind=8),dimension(:,:,:,:,0:)::angInternal
-      integer::Jdat
+      integer::Jdat, Nnpar, Ntotc, Ngroup
       real(kind=8)::Ccoul(2,Ntotc,Ngroup,*),   &
              Dcoul(2,Ntotc,Nnpar,Ngroup,*)
       real(kind=8)::Ccoulx(:,:)
       integer::igroup, Nn, Ix, Iipar
 
-      DO Igroup=1,resParData%getNumSpinGroups()
+      DO Igroup=1,Ngroup
          DO Nn=1,Ntotc
             DO Ix=1,2
                Ccoul(Ix,Nn,Igroup,Jdat) =    &
diff --git a/sammy/src/xct/mxct31.f90 b/sammy/src/xct/mxct31.f90
index 0c70d0dbc..43d524228 100644
--- a/sammy/src/xct/mxct31.f90
+++ b/sammy/src/xct/mxct31.f90
@@ -3,7 +3,6 @@ use XctCrossCalc_M
 implicit none
 public Setleg_Slow, Find_Kountr_Jx_Slow
 
-! Todo: Delete dependence on Iq_val, number of distinct Q values
 
 contains
 
diff --git a/sammy/src/xct/mxct32.f90 b/sammy/src/xct/mxct32.f90
index 0f23b213e..e39ed654f 100644
--- a/sammy/src/xct/mxct32.f90
+++ b/sammy/src/xct/mxct32.f90
@@ -4,7 +4,6 @@ module mxct32_m
 
   public Derleg_Slow
 
-  ! Todo: Delete dependence on Iq_val, number of distinct Q values
 contains
 
   integer function Jxlmn(Mb,Ma,M, Ntotc)
-- 
GitLab