From bbbc46a7736db280923a2ec3498fca5227cccf7a Mon Sep 17 00:00:00 2001
From: Wiarda <wiardada@ornl.gov>
Date: Fri, 17 Sep 2021 11:39:30 -0400
Subject: [PATCH] Add an implementation layer for cro and xct so the file
 structure does not need to change too much Take out I_Iflpol in favor of
 using the fit flag for the fission width directly Move A_Ibga to
 XctCrossCalc_M

---
 sammy/src/amr/mamr.f90                     |   6 +-
 sammy/src/amr/mamr2.f90                    |   1 -
 sammy/src/amx/mamx.f90                     |   6 +-
 sammy/src/blk/Exploc_common.f90            |   5 -
 sammy/src/blk/Templc_common.f90            |   4 -
 sammy/src/clq/ArtificalCross_M.f90         |   6 +-
 sammy/src/cro/CroCrossCalcImpl_M.f90       | 199 +++++++++++++++++++
 sammy/src/cro/CroCrossCalc_M.f90           |  27 ++-
 sammy/src/cro/mcro0.f90                    | 137 ++++++-------
 sammy/src/cro/mcro1.f90                    | 146 ++------------
 sammy/src/cro/mcro2.f90                    |  28 ---
 sammy/src/cro/mcro2a.f90                   |  55 ++++-
 sammy/src/cro/mcro4.f90                    |   2 +-
 sammy/src/cro/mcro8.f90                    | 126 +++++-------
 sammy/src/end/mout.f                       |  14 +-
 sammy/src/fin/mfin0.f90                    |   3 +-
 sammy/src/fin/mfin3.f90                    |  22 +-
 sammy/src/fin/mfin4.f90                    |   4 +-
 sammy/src/fin/mfin5.f90                    |  15 +-
 sammy/src/inp/minp01.f                     |   1 -
 sammy/src/mlb/MlbwCrossCalc_M.f90          |   6 +-
 sammy/src/new/mnew0.f                      |   1 -
 sammy/src/old/ReadCovarianceParts.f90      |  17 +-
 sammy/src/old/mold0.f                      |   4 +-
 sammy/src/old/mold2.f                      |   9 +-
 sammy/src/old/mold5.f                      |   7 +-
 sammy/src/par/mpar0.f90                    |   2 -
 sammy/src/par/mpar11.f90                   |  75 -------
 sammy/src/rec/mrec0.f                      |  14 +-
 sammy/src/rec/mrec2.f90                    |   2 +-
 sammy/src/rec/mrec3.f                      |   3 +-
 sammy/src/ref/mrfs0.f                      |   1 -
 sammy/src/ref/mrfs2.f                      |   1 -
 sammy/src/ref/mrfs4.f                      |   2 -
 sammy/src/ref/mwrt0.f                      |   1 -
 sammy/src/sammy/CMakeLists.txt             |   2 +
 sammy/src/the/CrossSectionCalcDriver_M.f90 |  70 +++++--
 sammy/src/the/CrossSectionCalculator_M.f90 |  11 +-
 sammy/src/the/DerivativeHandler.cpp        |   3 +-
 sammy/src/the/ZeroKCrossCorrections_M.f90  |   8 +-
 sammy/src/the/mthe0.f90                    |  28 ++-
 sammy/src/the/mthe1.f90                    |   7 +-
 sammy/src/xct/XctCrossCalcImpl_M.f90       |  68 +++++++
 sammy/src/xct/XctCrossCalc_M.f90           | 172 +++++++++++++++-
 sammy/src/xct/mxct0.f90                    |  99 ++++-----
 sammy/src/xct/mxct01.f90                   | 221 +++++++++++----------
 sammy/src/xct/mxct02.f90                   |  44 +---
 sammy/src/xct/mxct03.f90                   |  30 +--
 sammy/src/xct/mxct04.f90                   |  16 +-
 sammy/src/xct/mxct05.f90                   |  27 ++-
 sammy/src/xct/mxct06.f90                   |   3 +
 sammy/src/xct/mxct09.f90                   |   1 +
 sammy/src/xct/mxct10.f90                   |   1 +
 sammy/src/xct/mxct18.f90                   |  10 +-
 sammy/src/xct/mxct20.f90                   |  24 ++-
 sammy/src/xct/mxct21.f90                   | 115 +++++++----
 sammy/src/xct/mxct22.f90                   |  70 +++++--
 sammy/src/xct/mxct25.f90                   | 124 +++++++-----
 sammy/src/xct/mxct26.f90                   |  65 +++---
 sammy/src/xct/mxct30.f90                   |  56 ++++--
 sammy/src/xct/mxct31.f90                   |  67 +++++--
 sammy/src/xct/mxct32.f90                   |  59 ++++--
 62 files changed, 1378 insertions(+), 975 deletions(-)
 create mode 100644 sammy/src/cro/CroCrossCalcImpl_M.f90
 create mode 100644 sammy/src/xct/XctCrossCalcImpl_M.f90

diff --git a/sammy/src/amr/mamr.f90 b/sammy/src/amr/mamr.f90
index d4d594fb7..9aebc03b8 100644
--- a/sammy/src/amr/mamr.f90
+++ b/sammy/src/amr/mamr.f90
@@ -112,7 +112,7 @@
 ! ***         VALUES, ETC.
          CALL Rdcov (  A_Iprbrd , I_Iflbrd ,                       &
             A_Iprdet , I_Ifldet , I_Iigrde ,                       &
-            A_Ipolar , I_Iflpol ,                                  &
+            A_Ipolar ,                                             &
             A_Iprmsc , I_Iflmsc , I_Irdmsc , I_Ijkmsc , A_Ietaee , &
             A_Iprpmc , I_Iflpmc , I_Isopmc ,                       &
             A_Iprorr , I_Iflorr , A_Icrnch , A_Iedets ,            &
@@ -128,7 +128,7 @@
          Isingl = Idimen (Nsingl, 1, 'Nsingl, 1')
          Nsingl = Nsingl*2
          CALL Rdcovx ( A_Iprbrd , I_Iflbrd ,                       &
-            A_Ipolar , I_Iflpol , A_Iprmsc , I_Iflmsc , I_Irdmsc , &
+            A_Ipolar , A_Iprmsc , I_Iflmsc , I_Irdmsc ,            &
             A_Iprpmc , I_Iflpmc , I_Isopmc ,                       &
             A_Iprorr , I_Iflorr , A_Icrnch , A_Iedets ,            &
                        A_Iseses , A_Iesese , A_Isgdts ,            &
@@ -170,7 +170,7 @@
 ! *** Write THE "COVARIANCE" FILE
       CALL Wrcov (  A_Iprbrd , I_Iflbrd , &
          A_Iprdet , I_Ifldet , I_Iigrde , &
-         A_Ipolar , I_Iflpol , &
+         A_Ipolar ,  &
          A_Iprmsc , I_Iflmsc , I_Irdmsc , I_Ijkmsc , A_Ietaee , &
          A_Iprpmc , I_Iflpmc , I_Isopmc , &
          A_Iprorr , I_Iflorr , A_Icrnch , A_Iedets , &
diff --git a/sammy/src/amr/mamr2.f90 b/sammy/src/amr/mamr2.f90
index b79b570ab..73b31a1eb 100644
--- a/sammy/src/amr/mamr2.f90
+++ b/sammy/src/amr/mamr2.f90
@@ -39,7 +39,6 @@ module amr2
 !
          IF (Kpolar.EQ.1) THEN
             call make_A_Ipolar(2*Nres)
-            call make_I_Iflpol(2*Nres)
          END IF
       END IF
 !
diff --git a/sammy/src/amx/mamx.f90 b/sammy/src/amx/mamx.f90
index 410c1f4fb..a6c7bd1c5 100644
--- a/sammy/src/amx/mamx.f90
+++ b/sammy/src/amx/mamx.f90
@@ -114,7 +114,7 @@
 ! &&&    from old/mold2.f
          CALL Rdcov (  A_Iprbrd , I_Iflbrd ,  &
             A_Iprdet , I_Ifldet , I_Iigrde ,  &
-            A_Ipolar , I_Iflpol , &
+            A_Ipolar ,                        &
             A_Iprmsc , I_Iflmsc , I_Irdmsc , I_Ijkmsc , A_Ietaee , &
             A_Iprpmc , I_Iflpmc , I_Isopmc , &
             A_Iprorr , I_Iflorr , A_Icrnch , A_Iedets , &
@@ -132,7 +132,7 @@
          Nsingl = Nsingl*2
 ! &&&    from old/mold5.f
          CALL Rdcovx ( A_Iprbrd , I_Iflbrd ,  &
-            A_Ipolar , I_Iflpol , A_Iprmsc , I_Iflmsc , I_Irdmsc , &
+            A_Ipolar , A_Iprmsc , I_Iflmsc , I_Irdmsc , &
             A_Iprpmc , I_Iflpmc , I_Isopmc , &
             A_Iprorr , I_Iflorr , A_Icrnch , A_Iedets , &
                        A_Iseses , A_Iesese , A_Isgdts , &
@@ -176,7 +176,7 @@
 ! *** Write THE "COVARIANCE" FILE
       CALL Wrcov (  A_Iprbrd , I_Iflbrd , &
          A_Iprdet , I_Ifldet , I_Iigrde , &
-         A_Ipolar , I_Iflpol , &
+         A_Ipolar ,  &
          A_Iprmsc , I_Iflmsc , I_Irdmsc , I_Ijkmsc , A_Ietaee , &
          A_Iprpmc , I_Iflpmc , I_Isopmc , &
          A_Iprorr , I_Iflorr , A_Icrnch , A_Iedets , &
diff --git a/sammy/src/blk/Exploc_common.f90 b/sammy/src/blk/Exploc_common.f90
index 5df315dba..3c9fa4663 100644
--- a/sammy/src/blk/Exploc_common.f90
+++ b/sammy/src/blk/Exploc_common.f90
@@ -41,7 +41,6 @@ module exploc_common_m
      real(kind=8),allocatable,dimension(:)::A_Ipolar 
      
      ! old group '*'
-     integer,allocatable,dimension(:)::I_Iflpol
      real(kind=8),allocatable,dimension(:),target::A_Iprmsc 
      integer,allocatable,dimension(:)::I_Iflmsc 
      real(kind=8),allocatable,dimension(:)::A_Idemsc
@@ -256,10 +255,6 @@ module exploc_common_m
        call allocate_real_data(A_Ipolar,want)
      end subroutine make_A_Ipolar
 
-     subroutine make_I_Iflpol(want)
-       integer::want
-       call allocate_integer_data(I_Iflpol,want)
-     end subroutine make_I_Iflpol
 
      subroutine make_A_Iprmsc(want)
        integer::want
diff --git a/sammy/src/blk/Templc_common.f90 b/sammy/src/blk/Templc_common.f90
index c60ee4ee0..96abf345c 100644
--- a/sammy/src/blk/Templc_common.f90
+++ b/sammy/src/blk/Templc_common.f90
@@ -7,11 +7,8 @@ module templc_common_m
 !
       IMPLICIT NONE
       
-      real(kind=8),allocatable,dimension(:,:)::A_Ibr
-      real(kind=8),allocatable,dimension(:,:)::A_Ibi
       real(kind=8),allocatable,dimension(:,:)::A_Ipr
       real(kind=8),allocatable,dimension(:,:)::A_Ipi
-      real(kind=8),allocatable,dimension(:)::A_Ibga
       real(kind=8),allocatable,dimension(:)::A_Ipgar
       real(kind=8),allocatable,dimension(:)::A_Ipgai
       real(kind=8),allocatable,dimension(:)::A_Ialphr
@@ -50,7 +47,6 @@ module templc_common_m
       real(kind=8),allocatable,dimension(:)::A_Ipxri
       real(kind=8),allocatable,dimension(:)::A_Ixxxxr
       real(kind=8),allocatable,dimension(:)::A_Ixxxxi
-      real(kind=8),allocatable,dimension(:)::A_Ixx
       real(kind=8),allocatable,dimension(:)::A_Icccll
       real(kind=8),allocatable,dimension(:)::A_Idddll
       real(kind=8),allocatable,dimension(:)::A_Ixlmn
diff --git a/sammy/src/clq/ArtificalCross_M.f90 b/sammy/src/clq/ArtificalCross_M.f90
index 4de7e5290..bf9800f4a 100644
--- a/sammy/src/clq/ArtificalCross_M.f90
+++ b/sammy/src/clq/ArtificalCross_M.f90
@@ -203,15 +203,15 @@ subroutine ArtificalCross_calcCross(this, ener, Ipoten)
     call this%crossData%setSharedValNs(this%row, 1, 0, Sigxxx)
 end subroutine
 
-subroutine ArtificalCross_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero)
+subroutine ArtificalCross_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
     class(ArtificalCross) :: this
     type(SammyRMatrixParameters)::pars
     type(CovarianceData)::cov
     type(AdjustedRadiusData)::rad
     integer,intent(in)::niso, Itzero, Ilzero
-    logical,intent(in)::needAngular
+    logical,intent(in)::needAngular, doShiftRes
 
-    call CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero)
+    call CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
 end subroutine
 subroutine ArtificalCross_destroy(this)
     class(ArtificalCross) :: this
diff --git a/sammy/src/cro/CroCrossCalcImpl_M.f90 b/sammy/src/cro/CroCrossCalcImpl_M.f90
new file mode 100644
index 000000000..29a6cbc88
--- /dev/null
+++ b/sammy/src/cro/CroCrossCalcImpl_M.f90
@@ -0,0 +1,199 @@
+!
+!  Add an extra layer between CroCrossCalc and the
+!  existing cro methods, so that we can pass
+!  a CroCrossCalc object to them
+!  without having to change too many
+!  implementations
+!
+module CroCrossCalcImpl_M
+  use, intrinsic :: ISO_C_BINDING
+  use CroCrossCalc_M
+  use DerivativeHandler_M
+  use SammyRMatrixParameters_M
+  use CovarianceData_M
+  use AdjustedRadiusData_M
+  use XctCrossCalc_M
+  implicit none
+
+  type, extends(CroCrossCalc) :: CroCrossCalcImpl
+     contains
+     procedure, pass(this) :: calcCross => CroCrossCalcImpl_calcCross
+     procedure, pass(this) :: setEnergyIndependent => CroCrossCalcImpl_setEnergyIndependent  ! set energy independent values using current parameter values
+  end type CroCrossCalcImpl
+
+contains
+
+   subroutine CroCrossCalcImpl_calcCross(this, ener, Ipoten)
+      use workCro_M
+      use varyr_common_m, only : Su, Squ
+      class(CroCrossCalcImpl) ::this
+      real(kind=8)::ener
+      integer::ipoten
+
+
+      real(kind=8)::val
+      integer::i, itot, ns, is
+      logical(C_BOOL)::accu
+
+      call CrossSectionCalculator_calcCross(this, ener, Ipoten)
+      if (ener.eq.0.0d0) return
+
+      IF (this%ener.LT.0.0d0) this%ener = -this%ener
+
+      Su = this%ener
+      Squ = this%enerSq
+
+      this%A_Isigxx = 0.0d0
+      this%A_Idasig = 0.0d0
+      this%A_Idbsig = 0.0d0
+
+
+      call Work_Cro (this,  Ipoten, this%A_Isigxx, this%A_Idasig, this%A_Idbsig)    
+
+      if (ener.lt.0.0d0) then
+         itot = this%covariance%getNumTotalParam() + 1
+         if ( .not.this%wantDerivs) itot = 1
+         ns = this%crossData%getNnnsig()
+
+         accu = .false.
+         call this%crossData%setAccumulate(accu)
+         do is = 1, ns
+            do i = 0, itot
+              val = this%crossData%getDataNs(this%row, is, i, 1)
+              if (val.ne.0.0d0) then
+                  val = -val
+                  call this%crossData%addDataNs(this%row, is, i, 1, val)
+              end if
+            end do
+         end do
+         accu = .true.
+         call this%crossData%setAccumulate(accu)
+      end if
+
+      this%ener = ener
+   end subroutine
+
+   subroutine CroCrossCalcImpl_setEnergyIndependent(this, reactType, Twomhb, kwcoul, Etac)
+       use fixedi_m, only : Kiniso, Kkkiso, Niniso, Nnniso, Nnnsig, Nres, Nrext, Ntotc, Kshift, Ndasig, Ndbsig, Ngbout, Nfpres, Kpolar
+       use ifwrit_m, only : Ks_Res, Ksolve
+       use exploc_common_m, only : A_Ibound, A_Iechan, A_Ipolar, A_Izke
+       use oopsch_common_m, only : Nowwww, Segmen
+       use EndfData_common_m, only : covData, resParData, radFitFlags
+       use templc_common_m
+       use AllocateFunctions_m
+       use xct1_m
+       use mcro8_m
+
+       class(CroCrossCalcImpl) :: this
+       integer::kwcoul, reactType
+       real(kind=8)::Twomhb, Etac
+
+       integer::n_ix1
+       integer::Krext, Mxany, N, N2, N3, N4, N6, N8, NA, Nb, NN
+       integer::ns
+ !
+
+       call CrossSectionCalculator_setEnergyIndependent(this, reactType, Twomhb, kwcoul, Etac)
+       if (this%wantDerivs) then
+          call Babb (this,.false.)
+       end if
+
+       ! all the reset, except for the call to fixx and babb will
+       ! go away after we did the proper allocating in CroCrossCalc
+       !
+       ! in order to avoid circular dependencies in the code that will be deleted, some
+       ! is just copied here so that it compiles
+       !
+
+       ns = this%crossData%getNnnsig()
+       call allocate_real_data(this%A_Isigxx, ns)
+       call allocate_real_data(this%A_Idasig, Ndasig * ns)
+       call allocate_real_data(this%A_Idbsig, Ndbsig * ns)
+
+
+
+       !
+       ! *** find array sizes for xx
+       Mxany = this%resData%getNumResonances()
+       !
+       ! *** find array sizes for added cross section from endf/b-vi
+       !     No done in C++
+
+
+       Na = 1
+       Nb = 1
+       IF (Kshift.NE.0) Na = Nres
+       IF (Kshift.NE.0) Nb = Mxany*Mxany*this%resData%getNumSpinGroups()
+       NN = (Ntotc*(Ntotc+1))/2
+       IF (this%wantDerivs) THEN
+          N2 = NN*Nfpres
+          N3 = Ngbout
+          N4 = Nfpres
+       ELSE
+          N2 = 1
+          N3 = 1
+          N4 = 1
+       END IF
+       N6 = 1
+       N8 = 1
+
+
+
+
+       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
+       Krext = Nrext
+       IF (Nrext.EQ.0) Krext = 1
+       !
+       ! *** two ***
+       N = N2
+
+       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - <
+       ! *** three ***
+       n_ix1 = NA
+       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - >
+       !
+       ! *** five ***
+       call allocate_real_data(A_Ialphr, Nres)
+       call allocate_real_data(A_Ialphi, Nres)
+       call allocate_integer_data(I_Inot, Nres)
+       !
+       ! - - - - - - - - - - - - - - - - <
+       ! *** six ***
+       call allocate_real_data(A_Idifen, Nres)
+       call allocate_real_data(A_Ixden, Nres)
+       ! CALL Abpart
+       ! - - - - - - - - - - - - - - - - >
+       !
+       ! *** seven ***
+       N = Ntotc
+       call allocate_real_data(A_Ics, N)
+       call allocate_real_data(A_Isi, N)
+       call allocate_real_data(A_Idphi, N)
+       call allocate_real_data(A_Iz, N)
+       N = Nn
+       call allocate_real_data(A_Iwr, N)
+       call allocate_real_data(A_Iwi, N)
+       call allocate_real_data(A_Ixxxxr, N)
+       call allocate_real_data(A_Ixxxxi, N)
+       N = Nn*2
+       call allocate_real_data(A_Irmat, N)
+       call allocate_real_data(A_Irinv, N)
+       N = Ntotc*Ntotc
+       call allocate_real_data(A_Ipwrr, N)
+       call allocate_real_data(A_Ipwri, N)
+       call allocate_real_data(A_Ixqr, N)
+       call allocate_real_data(A_Ixqi, N)
+       call allocate_real_data(A_Itr, N)
+       call allocate_real_data(A_Iti, N)
+       N = Ntotc
+       call allocate_real_data(A_Isphr, N)
+       call allocate_real_data(A_Isphi, N)
+       call allocate_real_data(A_Iphr, N)
+       call allocate_real_data(A_Iphi, N)
+       N = Nn*Nn
+       call allocate_real_data(A_Iqr, N)
+       call allocate_real_data(A_Iqi, N)
+       N = N6
+   end subroutine
+
+end module CroCrossCalcImpl_M
diff --git a/sammy/src/cro/CroCrossCalc_M.f90 b/sammy/src/cro/CroCrossCalc_M.f90
index b6683f79a..2d7aa1cc6 100644
--- a/sammy/src/cro/CroCrossCalc_M.f90
+++ b/sammy/src/cro/CroCrossCalc_M.f90
@@ -5,30 +5,47 @@ module CroCrossCalc_M
   use SammyRMatrixParameters_M
   use CovarianceData_M
   use AdjustedRadiusData_M
+  use XctCrossCalc_M
+  use AllocateFunctions_m
   implicit none
 
-  type, extends(CrossSectionCalculator) :: CroCrossCalc
+  type, extends(XctCrossCalc) :: CroCrossCalc
      contains
+     procedure, pass(this) :: setUpDerivativeList => CroCrossCalc_setUpDerivativeList    ! set up  crossData, depending on number of isotopes
      procedure, pass(this) :: initialize => CroCrossCalc_initialize
      procedure, pass(this) :: destroy => CroCrossCalc_destroy
   end type CroCrossCalc
 
 contains
+subroutine CroCrossCalc_setUpDerivativeList(this)
+   class(CroCrossCalc) :: this
 
-subroutine CroCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero)
+   call CrossSectionCalculator_setUpDerivativeList(this)
+
+   if (this%reactType.eq.6) then   !  This is for eta-data only
+      call this%crossData%setNnsig(2)
+   end if
+
+   if (this%wantCrossPerSpin) then
+       call this%crossPerSpin%nullify()
+   end if
+end subroutine
+
+subroutine CroCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
     class(CroCrossCalc) :: this
     type(SammyRMatrixParameters)::pars
     type(CovarianceData)::cov
     type(AdjustedRadiusData)::rad
     integer,intent(in)::niso, Itzero, Ilzero
-    logical,intent(in)::needAngular
+    logical,intent(in)::needAngular, doShiftRes
 
-    call CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero)
+    call XctCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
 end subroutine
+
 subroutine CroCrossCalc_destroy(this)
     class(CroCrossCalc) :: this
 
-    call CrossSectionCalculator_destroy(this)
+    call XctCrossCalc_destroy(this)
 end subroutine
 
 end module CroCrossCalc_M
diff --git a/sammy/src/cro/mcro0.f90 b/sammy/src/cro/mcro0.f90
index 96db93dc4..9a0ca819c 100644
--- a/sammy/src/cro/mcro0.f90
+++ b/sammy/src/cro/mcro0.f90
@@ -3,26 +3,21 @@
 !
       SUBROUTINE Samcro_0
 !
-      use oops_common_m
-      use fixedi_m
-      use ifwrit_m
-      use exploc_common_m
-      use array_sizes_common_m
-      use oopsch_common_m
-      use cbro_common_m
-      use lbro_common_m
-      use EndfData_common_m
+      use fixedi_m, only : Iq_Iso, Jcros, Kiniso, Kkkiso, Niniso, Nnniso, Nnnsig, Npfil3, Nres, Nrext, Nrfil3, Ntotc, needResDerivs
+      use ifwrit_m, only : Ks_Res, Ksolve
+      use exploc_common_m, only : A_Ibound, A_Iechan, A_Ipolar, A_Izke
+      use oopsch_common_m, only : Nowwww, Segmen
+      use EndfData_common_m, only : covData, resParData, radFitFlags
       use templc_common_m
       use AllocateFunctions_m
-      use rsl7_m
-      use xct_m
-      use xct1_m
-      use mxct27_m
+      use rsl7_m, only : Set_Kws_Xct
+      use mcro8_m
       use AuxGridHelper_M, only : setAuxGridOffset,   &
                                   setAuxGridRowMax
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      IMPLICIT None
       real(kind=8),allocatable,dimension(:)::A_Ix1
       integer::n_ix1
+      integer::Krext, Mxany, N, N2, N3, N4, N6, N8, NA, Nb, NN
 !
 !
       WRITE (6,99999)
@@ -54,7 +49,7 @@
       IF (Nnniso.NE.1) STOP '[STOP in Samcro_0 in cro/mcro0.f]'
 !
 ! *** find array sizes for xx
-      IF (resParData%getNumResonances().GT.0) CALL Findmx (Mxany)
+      MxAny = resParData%getNumResonances()
 !
 ! *** find array sizes for added cross section from endf/b-vi
 !     Now done in C++
@@ -67,7 +62,7 @@
 !
 ! *** one ***
       CALL Set_Kws_Xct
-      call calcData%setUpList(resParData, radFitFlags, Iq_Iso)
+      !call calcData%setUpList(resParData, radFitFlags, Iq_Iso)
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
 
       Krext = Nrext
@@ -76,86 +71,72 @@
 ! *** two ***
       N = N2
 !
-      call allocate_real_data(A_Ixx, Nres)
+      !call allocate_real_data(A_Ixx, Nres)
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - <
 ! *** three ***
       n_ix1 = NA
       if (Mxany.gt.NA) n_ix1 = Mxany
-      ! needs to be dimensioned at least Mxany
-      ! as we zero the array to size Mxanny
-      ! in subroutine Fixx
-      call allocate_real_data(A_Ix1, n_ix1)
-      IF (Nres.GT.0) CALL Fixx (   &
-         A_Ibound , A_Iechan ,   &
-         A_Izke  , A_Ixx    , A_Ix1   ,  Mxany)
-! *** Sbroutine Fixx sets Xx = energy shift
-      deallocate(A_Ix1)
+
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 !
 !
-!   Fals is passed to babb since it is used to set parameters NOT for numerical differentiation
-      IF (Ksolve.NE.2 .AND. needResDerivs) CALL Babb(   &
-         A_Ipolar , I_Iflpol , A_Ibr , A_Ibi , A_Ixx ,.false.)
-! *** Sbroutine Babb_cro generates energy-independent portion of
-! ***    partial derivatives
-!
 !
 ! *** four ***
 ! - - - - - - - - - - - - - - - - - - - - - - <
 !
 ! *** five ***
-      call allocate_real_data(A_Ialphr, Nres)
-      call allocate_real_data(A_Ialphi, Nres)
-      call allocate_integer_data(I_Inot, Nres)
+      !call allocate_real_data(A_Ialphr, Nres)
+      !call allocate_real_data(A_Ialphi, Nres)
+      !call allocate_integer_data(I_Inot, Nres)
 !
 ! - - - - - - - - - - - - - - - - <
 ! *** six ***
-      call allocate_real_data(A_Idifen, Nres)
-      call allocate_real_data(A_Ixden, Nres)
+      !call allocate_real_data(A_Idifen, Nres)
+      !call allocate_real_data(A_Ixden, Nres)
 ! CALL Abpart
 ! - - - - - - - - - - - - - - - - >
 !
 ! *** seven ***
       N = Ntotc
-      call allocate_real_data(A_Ics, N)
-      call allocate_real_data(A_Isi, N)
-      call allocate_real_data(A_Idphi, N)
-      call allocate_real_data(A_Iz, N)
+      !call allocate_real_data(A_Ics, N)
+      !call allocate_real_data(A_Isi, N)
+      !call allocate_real_data(A_Idphi, N)
+      !call allocate_real_data(A_Iz, N)
       N = Nn
-      call allocate_real_data(A_Iwr, N)
-      call allocate_real_data(A_Iwi, N)
-      call allocate_real_data(A_Ixxxxr, N)
-      call allocate_real_data(A_Ixxxxi, N)
+      !call allocate_real_data(A_Iwr, N)
+      !call allocate_real_data(A_Iwi, N)
+      !call allocate_real_data(A_Ixxxxr, N)
+      !call allocate_real_data(A_Ixxxxi, N)
       N = Nn*2
-      call allocate_real_data(A_Irmat, N)
-      call allocate_real_data(A_Irinv, N)
+      !call allocate_real_data(A_Irmat, N)
+      !call allocate_real_data(A_Irinv, N)
       N = Ntotc*Ntotc
-      call allocate_real_data(A_Ipwrr, N)
-      call allocate_real_data(A_Ipwri, N)
-      call allocate_real_data(A_Ixqr, N)
-      call allocate_real_data(A_Ixqi, N)
-      call allocate_real_data(A_Itr, N)
-      call allocate_real_data(A_Iti, N)
+      !call allocate_real_data(A_Ipwrr, N)
+      !call allocate_real_data(A_Ipwri, N)
+      !call allocate_real_data(A_Ixqr, N)
+      !call allocate_real_data(A_Ixqi, N)
+      !call allocate_real_data(A_Itr, N)
+      !call allocate_real_data(A_Iti, N)
       N = Ntotc
-      call allocate_real_data(A_Isphr, N)
-      call allocate_real_data(A_Isphi, N)
-      call allocate_real_data(A_Iphr, N)
-      call allocate_real_data(A_Iphi, N)
+      !call allocate_real_data(A_Isphr, N)
+      !call allocate_real_data(A_Isphi, N)
+      !call allocate_real_data(A_Iphr, N)
+      !call allocate_real_data(A_Iphi, N)
       N = Nn*Nn
-      call allocate_real_data(A_Iqr, N)
-      call allocate_real_data(A_Iqi, N)
+      !call allocate_real_data(A_Iqr, N)
+      !call allocate_real_data(A_Iqi, N)
       N = N6
 ! CALL Parsh
 ! - - - - - - - - - - - - - - - - - - - - - - >
 !
 !
 ! *** eight ***
-      CALL Work_Cro (I_Iflmsc ,   &
-             A_Isigxx , A_Idasig , A_Idbsig)
-! *** Sbroutine Work_Cro generates theory and derivatives
-      deallocate(A_Ics)
-      deallocate(A_Isi)
-      deallocate(A_Idphi)
+!      actual routine that does the reconstruction is called from
+!      ZeroKCrossCorrections
+!
+      !deallocate(A_Ics)
+      !deallocate(A_Isi)
+      !deallocate(A_Idphi)
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 !
 !
@@ -173,22 +154,18 @@
 !
 ! *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-Cross section
 !
-      use fixedi_m
-      use ifwrit_m
-      use broad_common_m
-      use lbro_common_m
-      use EndfData_common_m
-      use SammyGridAccess_M
+      use fixedi_m, only : Jcros, Kshift, Nfpres, Ngbout, Ngroup, Npfil3, Nres, Ntotc, Nrfil3, numcro
+      use ifwrit_m, only : Ksolve
       use rsl7_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      IMPLICIT None
 
-      type(SammyGridAccess)::grid
+      integer::Na, Nb, N2, N3, N4, NN, N6, N8, Mxany
+      integer::I, Idimen
+      integer::K, K1, K2, K3, K4, K5, K6, K7, K8, Ka, Kone,numE
+      external Idimen
 !
 !
-      call grid%initialize()
-      call grid%setParameters(numcro, ktzero)
-      call grid%setToAuxGrid(expData)
-      numE = grid%getNumEnergies(expData)
+
       Na = 1
       Nb = 1
       IF (Kshift.NE.0) Na = Nres
@@ -204,9 +181,7 @@
          N4 = 1
       END IF
       N6 = 1
-      IF (Kpiece.EQ.1) N6 = numE*Ngroup
-      N8 = 1
-      IF (Kpiece.EQ.1) N8 = numE
+      N8 = 1     
 !
 ! *** one ***
       CALL Figure_Kws_Cro (Kone)
@@ -247,6 +222,6 @@
       K = Idimen (K,  1, 'K,  1')
       I = Idimen (K, -1, 'K, -1')
       I = Idimen (0, 0, '0, 0')
-      call grid%destroy()
+
       RETURN
       END
diff --git a/sammy/src/cro/mcro1.f90 b/sammy/src/cro/mcro1.f90
index f8b16614d..f3ca37546 100644
--- a/sammy/src/cro/mcro1.f90
+++ b/sammy/src/cro/mcro1.f90
@@ -1,137 +1,35 @@
+module workCro_M
+contains
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Work_Cro (Iflmsc,   &
-         Sigxxx, Dasigx, Dbsigx)
+      SUBROUTINE Work_Cro (calc, Ipoten, Sigxxx, Dasigx, Dbsigx)
 !
 ! *** PURPOSE -- GENERATE THEORETICAL DATA "theory" AND PARTIAL
 ! ***            DERIVATIVES "GB"
 !
-      use oops_common_m
-      use fixedi_m
-      use ifwrit_m
-      use oopsch_common_m
-      use fixedr_m
-      use varyr_common_m
-      use lbro_common_m
-      use mxct27_m
-      use EndfData_common_m
-      use SammyGridAccess_M
-      use array_sizes_common_m, only : calcData
-      use convert_to_transmission_m
-      use normalize_and_background
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      LOGICAL Ywhich
-!
-      type(SammyGridAccess)::grid
-      logical::wantDeriv
-      DIMENSION Iflmsc(*)
-      DIMENSION Sigxxx(Nnnsig,*),   &
-         Dasigx(Nnnsig,*), Dbsigx(Nnnsig,Ndbxxx,*)
-!
-!      DIMENSION Iflmsc(Nummsc)
-!
-      DATA Zero /0.0d0/
-!
-      call grid%initialize()
-      call grid%setParameters(numcro, ktzero)
-      call grid%setToAuxGrid(expData)
-      numE = grid%getNumEnergies(expData)
-!
-      Ywhich = Ydoppr.OR.Yresol
-
-      wantDeriv = Ksolve.NE.2
-      if ( covData%getPupedParam().gt.0) wantDeriv = .true.
-!
-      Kkkkkk = 0
-      Kkkmin = 0
-      call calcData%nullify()
-      call calcData%setNnsig(Nnnsig)
-!
-! *** Calculate cross sections and derivatives for all
-! ***      positive energies
-         DO J=1,numE
-            Jdat = J
-            do Jsig = 1,Nnnsig
-               Jcount = (Jdat-1)*Nnnsig + Jsig
-               call calcData%setToZero(Jcount, Ndasig + Ndbsig)            
-            end do
-
-            CALL Zero_Array (Sigxxx, Nnnsig*Nnniso)
-            IF (Ndasig.GT.0) CALL Zero_Array (Dasigx, Nnnsig*Ndasig)
-            IF (Ndbsig.GT.0) CALL Zero_Array (Dbsigx,   &
-                                               Nnnsig*Ndbsig*Nnniso)
-            Su = grid%getEnergy(J, expData)
-            IF (Kkkdop.EQ.1 .AND. Su.LE.Zero) cycle
-!
-            IF (Su.NE.Zero) THEN
-!
-               IF (Su.LT.Zero) Su = - Su
-               Squ = dSQRT(Su)
-               Ipoten = 0
-               IF (Kpoten.NE.0 .AND. (J.EQ.1 .OR. J.EQ.numE)) THEN
-                  Ipoten = 1
-               END IF
-!
-! *********    Generate cross sections and derivatives
-               IF (grid%getEnergy(J, expData).NE.Zero) THEN
-                  CALL Zparsh (Ipoten)
-                  IF (grid%getEnergy(J, expData).LT.Zero) THEN
-                     CALL Fix_Negative_E (Sigxxx, Dasigx, Dbsigx)
-                  END IF
-               ENDIF
-!
-            END IF
+      use fixedi_m, only : Nnnsig, Ndasig, Ndbsig
+      use CroCrossCalc_M
+      use mcro2a_m
+      IMPLICIT None
 
+      class(CroCrossCalc)::calc
+      integer::Ipoten
+      real(kind=8)::Sigxxx(Nnnsig),  Dasigx(Nnnsig,*), Dbsigx(Nnnsig,Ndbsig,*)
+      integer::Jsig, Iipar
 
-            IF (Maxwel.EQ.1 .OR. Knocor.EQ.1.or.   &
-                Ydoppr.or.Yresol .OR. Yangle.or.   &
-                Kfake.eq.1)  then
-                Kkkkkk = Jdat
-            else
-                Kkkkkk = Kkkkkk + 1
-            end if
+      CALL Zparsh (calc, Ipoten)
 
-            call calcData%addCalculatedData(Kkkkkk, Nnnsig, ndasig,   &
-                         ndbsig, -1, Sigxxx, Dasigx, Dbsigx)
-!
+      do Jsig = 1, calc%crossData%getNnnsig()
+          call calc%crossData%setSharedValNs(calc%row, Jsig, 0, Sigxxx(Jsig))
+          do Iipar = 1, Ndasig
+              call calc%crossData%setSharedValNs(calc%row, Jsig, Iipar, Dasigx(Jsig, Iipar))
+          end do
+          DO Iipar=1,Ndbsig
+              call calc%crossData%setSharedValNs(calc%row,Jsig, Iipar + Ndasig, Dbsigx(Jsig, Iipar, 1))
+          end do
       end do
-!
-      call grid%destroy()
-      RETURN
-      END
-!
-!
-! --------------------------------------------------------------
-!
-      SUBROUTINE Fix_Negative_E (Sigxxx, Dasigx, Dbsigx)
-      use fixedi_m
-      use ifwrit_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      DIMENSION Sigxxx(Nnnsig,*), Dasigx(Nnnsig,*),   &
-        Dbsigx(Nnnsig,Ndbxxx,*)
-      DO Iso=1,Nnniso
-         DO N=1,Nnnsig
-            Sigxxx(N,Iso) = - Sigxxx(N,Iso)
-         END DO
-      END DO
-      IF (Ksolve.NE.2) THEN
-         IF (Ndasig.GT.0) THEN
-            DO Ipar=1,Ndasig
-               DO N=1,Nnnsig
-                  Dasigx(N,Ipar) = - Dasigx(N,Ipar)
-               END DO
-            END DO
-         END IF
-         IF (Ndbsig.GT.0) THEN
-            DO Iso=1,Nnniso
-               DO Ipar=1,Ndbsig
-                  DO N=1,Nnnsig
-                     Dbsigx(N,Ipar,Iso) = - Dbsigx(N,Ipar,Iso)
-                  END DO
-               END DO
-            END DO
-         END IF
-      END IF
+
       RETURN
       END
+end module workCro_M
diff --git a/sammy/src/cro/mcro2.f90 b/sammy/src/cro/mcro2.f90
index 58fcb06ca..e45240af1 100644
--- a/sammy/src/cro/mcro2.f90
+++ b/sammy/src/cro/mcro2.f90
@@ -2,34 +2,6 @@
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Zparsh (Ipoten)
-!
-! *** PURPOSE -- FORM THE CROSS SECTION Sigma
-! ***            AND THE ( PARTIAL DERIVATIVES OF THE CROSS SECTION
-! ***            WITH RESPECT TO THE VARIED PARAMETERS ) = SJkL
-!
-      use over_common_m
-      use oops_common_m
-      use fixedi_m
-      use ifwrit_m
-      use exploc_common_m
-      use array_sizes_common_m
-      use templc_common_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-!
-!
-      CALL Abpart_Cro (   &
-         A_Ialphr , A_Ialphi ,   &
-         A_Idifen , A_Ixden ,   &
-         A_Idifma , I_Inot  , I_Inotu , A_Ixx )
-!
-      CALL Parsh (   &
-         A_Izke  ,   &
-         Ipoten,   &
-         A_Isigxx , A_Idasig , A_Idbsig )
-!
-      RETURN
-      END
 !
 !
 ! --------------------------------------------------------------
diff --git a/sammy/src/cro/mcro2a.f90 b/sammy/src/cro/mcro2a.f90
index 28e597d6f..97a4a503e 100644
--- a/sammy/src/cro/mcro2a.f90
+++ b/sammy/src/cro/mcro2a.f90
@@ -1,10 +1,41 @@
+module mcro2a_m
+contains
 !
+     SUBROUTINE Zparsh (calc, Ipoten)
+!
+! *** PURPOSE -- FORM THE CROSS SECTION Sigma
+! ***            AND THE ( PARTIAL DERIVATIVES OF THE CROSS SECTION
+! ***            WITH RESPECT TO THE VARIED PARAMETERS ) = SJkL
+!
+     use over_common_m
+     use oops_common_m
+     use fixedi_m
+     use ifwrit_m
+     use exploc_common_m
+     use templc_common_m
+     use CroCrossCalc_M
+     IMPLICIT DOUBLE PRECISION (a-h,o-z)
+     class(CroCrossCalc)::calc
+!
+!
+     CALL Abpart_Cro (  calc,  &
+        A_Ialphr , A_Ialphi ,   &
+        A_Idifen , A_Ixden ,   &
+        A_Idifma , I_Inot  , I_Inotu)
+!
+     CALL Parsh (   &
+        A_Izke  ,   &
+        Ipoten,   &
+        calc%A_Isigxx , calc%A_Idasig , calc%A_Idbsig )
+!
+     RETURN
+     END
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Abpart_Cro (Alphar,   &
+      SUBROUTINE Abpart_Cro (calc, Alphar,   &
          Alphai, Difen, Xden, Difmax,   &
-         Not, Notu, Xx)
+         Not, Notu)
 !
 ! *** PURPOSE -- GENERATE Upr AND Upi = ENERGY-DEPENDENT Pieces OF
 ! ***            PR AND PI = PARTIAL OF R WRT U-PARAMETERS
@@ -18,15 +49,17 @@
       use EndfData_common_m
       use RMatResonanceParam_M
       use templc_common_m
+      use CroCrossCalc_M
       IMPLICIT None
 !
+      type(CroCrossCalc)::calc
       type(SammyResonanceInfo)::resInfo
       type(SammySpinGroupInfo)::spinInfo
       type(RMatResonance)::resonance
       real(kind=8)::   &
          Alphar(*),Alphai(*),   &
          Difen(*), Xden(*),   &
-         Difmax(*), Xx(*)
+         Difmax(*)
       integer:: Not(*), Notu(*)
       real(kind=8)::Zero, Two
       real(kind=8)::Aa, G2, G3
@@ -51,7 +84,9 @@
          Alphai(I)= Zero
          Not(I)   = 1
          Difen(I) = resonance%getEres() - Su
-         Difen(I) = Difen(I) + Xx(I)
+         if (calc%doResShift) then
+            Difen(I) = Difen(I) + calc%Xx(I)
+         end if
          IF (dABS(Difen(I)).LT.100.0D0*Difmax(I)) Not(I) = 0
 
          call resParData%getResonanceInfo(resInfo, I)
@@ -137,11 +172,11 @@
                DO I=1,Ntotc
                   DO J=1,I
                      Ij = Ij + 1
-                     IF (A_Ibr(Ij,K).NE.Zero) then
-                         A_Ipr(Ij,K) = A_Ibr(Ij,K)*Upr(K)
+                     IF (calc%Br(Ij,K).NE.Zero) then
+                         A_Ipr(Ij,K) = calc%Br(Ij,K)*Upr(K)
                      END IF
-                     IF (A_Ibi(Ij,K).NE.Zero) then
-                         A_Ipi(Ij,K) = A_Ibi(Ij,K)*Upi(K)
+                     IF (calc%Bi(Ij,K).NE.Zero) then
+                         A_Ipi(Ij,K) = calc%Bi(Ij,K)*Upi(K)
                      END IF
                   END DO
                END DO
@@ -289,7 +324,7 @@
             Agoj = VarAbn*spinInfo%getGFactor()
 ! ***       TOTAL CROSS SECTIONS
             nent = spinInfo%getNumEntryChannels()
-            next = spinInfo%getNumExitChannels()
+            next = spinInfo%getNumExitChannels()          
             IF (Kcros.EQ.1) CALL Total (Agoj, nent, Ntotnn,   &
                A_Ics, A_Isi,   &
                A_Idphi , A_Iwr , A_Iwi , A_Ipwrr , A_Ipwri , A_Itr ,   &
@@ -331,3 +366,5 @@
 !
       RETURN
       END
+
+end module mcro2a_m
diff --git a/sammy/src/cro/mcro4.f90 b/sammy/src/cro/mcro4.f90
index cf8d97dbc..78ff910cd 100644
--- a/sammy/src/cro/mcro4.f90
+++ b/sammy/src/cro/mcro4.f90
@@ -919,7 +919,7 @@
                END DO
                if (Notu(M).gt.0) then
                   Dasigx(1,Notu(M)) = Dasigx(1,Notu(M)) +   &
-                                      Fourpi*Agoj*S/Su
+                                      Fourpi*Agoj*S/Su               
                else
                   Ifl = -1 * Notu(M) - Ndasig
                   Dbsigx(1,Ifl,iso) = Dbsigx(1,Ifl,Iso) +   &
diff --git a/sammy/src/cro/mcro8.f90 b/sammy/src/cro/mcro8.f90
index 2ee3437e7..bc35b69a1 100755
--- a/sammy/src/cro/mcro8.f90
+++ b/sammy/src/cro/mcro8.f90
@@ -1,59 +1,22 @@
+module mcro8_m
+contains
 !
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Findmx (Mxany)
-!
-!     *** PURPOSE -- find the maximim number of resonances for any spin group
-!
-      use fixedi_m,only:Kshift
-      use EndfData_common_m
-      use SammySpinGroupInfo_M
-      use RMatResonanceParam_M
-      IMPLICIT NONE
-      integer::Mxany
-      type(RMatSpinGroup)::spinGroup
-      type(SammySpinGroupInfo)::spinInfo
-      integer::n,nmax
-      
-      Mxany = 0
-      if( resParData%getNumResonances().eq.0) return
-      
-      if (Kshift.eq.0) then
-         Mxany = resParData%getNumResonances()
-         return
-      end if   
-
-      do n = 1, resParData%getNumSpinGroups()
-         call resParData%getSpinGroupInfo(spinInfo, n)
-         call resParData%getSpinGroup(spinGroup, spinInfo)
-         nmax = spinGroup%getNres()
-         IF (nmax.GT.Mxany) Mxany = nmax
-      end do
- 
-      RETURN
-      END
-!
-!
-! --------------------------------------------------------------
-!
-      SUBROUTINE Fixx (Bound, Echan,   &
-         Zke, Xx, X1, Mxany)
+      SUBROUTINE Fixx (calc, XX,  X1)
 !
 ! *** PURPOSE -- set Xx(I) = shift for resonance I when shift factors are
 ! ***              to be used in this manner
 !
-      use fixedi_m
-      use ifwrit_m
-      use EndfData_common_m
       use SammyResonanceInfo_M
       use RMatResonanceParam_M
       use xxx5
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      use CrossSectionCalculator_M
+      IMPLICIT None
 !
-      DIMENSION   &
-         Bound(Ntotc,*), Echan(Ntotc,*),   &
-         Zke(Ntotc,*), Xx(Mres), X1(Mres)
+      real(kind=8) :: X1(:), XX(:)
+      class(CrossSectionCalculator)::calc
       type(SammyResonanceInfo)::resInfo
       type(RMatResonance)::resonance
       type(SammySpinGroupInfo)::spinInfo
@@ -62,35 +25,34 @@
       type(RMatChannelParams)::channel
       type(SammyChannelInfo)::channelInfo
       real(kind=8),allocatable::resTmp(:)
-!
-!      DIMENSION Ntot(Ngroup),
-!     *   Bound(Ntotc,Ngroup), Echan,(Ntotc,Ngroup),
-!     *   Zke(Ntotc,Ngroup), Xx(Nres), X1(Nres)
-!
-      DATA Small /0.00001D0/, Zero /0.0d0/
+
+      real(kind=8),parameter:: Small = 0.00001D0, Zero = 0.0d0
+      integer::number, nres
+
+      real(kind=8)::beta, del, eres, Rho, X
+      integer::I, Ich, Ichan, Jsmall, K, KI, L, maxr, minr, Nnnn, Ii
       DATA Number /40/
 !
-      CALL Zero_Array (Xx, Mxany)
-      CALL Zero_Array (X1, Mxany)
-!
-      IF (Kshift.EQ.0) RETURN
+
+
+      Xx = 0.0d0
+      if( calc%resData%getNumResonances().eq.0) return
 !
       WRITE (21,20000)
-20000 FORMAT (' *** CAUTION *** Sub routine Fixx is not working properly   &
-       yet for', /,   &
-              '                 other-than-isolated resonances.  Proceed   &
-       with care.', /,   &
-              '                 Contact the author (N M Larson) if help   &
-      is needed.')
+20000 FORMAT (' *** CAUTION *** Sub routine Fixx is not working properly yet for', /,   &
+              '                 other-than-isolated resonances.  Proceed  with care.', /,   &
+              '                 Contact the author (N M Larson) if help  is needed.')
 !
+       X1 = 0.0d0
+
       DO Number=1,5
          Jsmall = 0
          maxr = 0
-         DO K=1,Ngroup
-            call resParData%getSpinGroupInfo(spinInfo, K)
+         DO K=1, calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinInfo, K)
             minr = maxr + 1
-            do i = minr, resParData%getNumResonances()
-               call resParData%getResonanceInfo(resInfo, i)
+            do i = minr, calc%resData%getNumResonances()
+               call calc%resData%getResonanceInfo(resInfo, i)
                if( resInfo%getSpinGroupIndex() .ne.k) exit
                maxr = i
             end do
@@ -99,8 +61,8 @@
                KI = 0
                Nnnn = K
                DO I=minr, maxr
-                  call resParData%getResonanceInfo(resInfo, i)
-                  call  resParData%getRedResonance(resonance, resInfo)
+                  call calc%resData%getResonanceInfo(resInfo, i)
+                  call  calc%resData%getRedResonance(resonance, resInfo)
                   eres = resonance%getEres()
                   II = 0
                   KI = KI + 1
@@ -109,17 +71,17 @@
                      II = II + Ich
                      ichan = spinInfo%getWidthForChannel(ich)
                      call spinInfo%getChannelInfo(channelInfo, Ich)
-                     call resParData%getParticlePairInfo(   &
+                     call calc%resData%getParticlePairInfo(   &
                          pairInfo,   &
                          channelInfo%getParticlePairIndex())
-                     call resParData%getParticlePair(pair, pairInfo)
-                     call resParData%getChannel(channel, channelInfo)
+                     call calc%resData%getParticlePair(pair, pairInfo)
+                     call calc%resData%getChannel(channel, channelInfo)
                      IF (pair%getCalcShift()) THEN
-                        IF (dABS(eres).GE.Echan(Ich,Nnnn)) THEN
+                        IF (dABS(eres).GE.calc%Echan(Ich,Nnnn)) THEN
                            L = channel%getL()
-                           Rho = channel%getApt()*Zke(Ich,Nnnn)*dSQRT(   &
-                              dABS(eres)+xx(I)-Echan(Ich,Nnnn) )
-                           X = Sf (Rho, L, Bound(Ich,Nnnn))
+                           Rho = channel%getApt()*calc%Zke(Ich,Nnnn)*dSQRT(   &
+                              dABS(eres)+Xx(I)-calc%Echan(Ich,Nnnn) )
+                           X = Sf (Rho, L, calc%Bound(Ich,Nnnn))
                            beta = resonance%getWidth(ichan)
                            beta = beta * beta
                            X1(I) = X1(I) + X*Beta
@@ -139,10 +101,11 @@
       WRITE (6,10300) Number
 10300 FORMAT (' Search for energy-shift did not converge in', I3,   &
            ' iterations.  Values for E-shift are: ')
-      allocate(resTmp(resParData%getNumResonances()))
-      do i = 1, resParData%getNumResonances()
-         call resParData%getResonanceInfo(resInfo, i)
-         call  resParData%getResonance(resonance, resInfo)
+      nres = calc%resData%getNumResonances()
+      allocate(resTmp(nres))
+      do i = 1, calc%resData%getNumResonances()
+         call calc%resData%getResonanceInfo(resInfo, i)
+         call  calc%resData%getResonance(resonance, resInfo)
          resTmp(i) = resonance%getEres()
       end do
       WRITE (6,10600) (resTmp(I),Xx(I),i=1,Nres)
@@ -153,10 +116,10 @@
       WRITE (6,10500) Number
 10500 FORMAT (' Search for energy-shift converged in', I3,   &
          ' iterations.  Values for E-shift are:')
-      allocate(resTmp(resParData%getNumResonances()))
-      do i = 1, resParData%getNumResonances()
-         call resParData%getResonanceInfo(resInfo, i)
-         call  resParData%getResonance(resonance, resInfo)
+      allocate(resTmp(calc%resData%getNumResonances()))
+      do i = 1, calc%resData%getNumResonances()
+         call calc%resData%getResonanceInfo(resInfo, i)
+         call  calc%resData%getResonance(resonance, resInfo)
          resTmp(i) = resonance%getEres()
       end do
       WRITE (6,10600) (resTmp(I),Xx(I),I=1,Nres)
@@ -165,3 +128,4 @@
 !
       RETURN
       END
+end module mcro8_m
diff --git a/sammy/src/end/mout.f b/sammy/src/end/mout.f
index 45b7582be..5a351b5e9 100644
--- a/sammy/src/end/mout.f
+++ b/sammy/src/end/mout.f
@@ -4,7 +4,7 @@ C --------------------------------------------------------------
 C
       SUBROUTINE Outpar (Iftit , Nunit , Parbrd, Iflbrd,
      *   Pardet, Ifldet, Igrdet,
-     *   Polar , Iflpol, Parmsc, Iflmsc, Iradms, Ijkmsc,
+     *   Polar , Parmsc, Iflmsc, Iradms, Ijkmsc,
      *   Etaeee, Parpmc, Iflpmc, Isopmc,
      *   Parorr, Iflorr, Ecrnch, Endets, Sesese, Eseses, Sigdts,
      *   Parrpi, Iflrpi, Parnbk, Iflnbk,
@@ -39,7 +39,7 @@ C
       type(ResonanceCovariance)::physCov
       DIMENSION Parbrd(*), Iflbrd(*),
      *   Pardet(*), Ifldet(*), Igrdet(*),
-     *   Polar(2,*), Iflpol(2,*),  Parmsc(*),
+     *   Polar(2,*), Parmsc(*),
      *   Iflmsc(*), Iradms(*), Ijkmsc(*), Etaeee(*),
      *   Parpmc(4,*), Iflpmc(4,*), Isopmc(*), Parorr(*), Iflorr(*),
      *   Ecrnch(*), Endets(*), Sigdts(*), Sesese(*), Eseses(*),
@@ -54,7 +54,7 @@ C
 C
 C      DIMENSION Parbrd(Numbrd), Iflbrd(Numbrd),
 c     *   Pardet(Numdet), Ifldet(Numdet), Igrdet(Ngroup),
-C     *   Polar(2,Nres), Iflpol(2,Nres),
+C     *   Polar(2,Nres),
 C     *   Parmsc(Nummsc), Iflmsc(Nummsc), Iradms(Ngroup), Ijkmsc(Nummsc),
 C     *   Etaeee(Mjetan),
 C     *   Parpmc(4,Numpmc), Iflpmc(4,Numpmc), Isopmc(Numpmc),
@@ -330,8 +330,8 @@ C ***       Here if want fission gamma widths in polar coordinates
                            Iif(J) = resInfo%getChannelFitOption(j-1)
                         end if
                      END DO
-                     Iif(6) = Iflpol(1,N)
-                     Iif(7) = Iflpol(2,N)                     
+                     Iif(6) = resInfo%getChannelFitOption(3)
+                     Iif(7) = resInfo%getChannelFitOption(4)
                      CALL Setflg (Iif, Mmax+4)
                      CALL Pwrite (pkenPr, resonance,
      *                  Polar(1,N), Nunit, Mmax)
@@ -381,8 +381,8 @@ Ccccc          FORMAT (/' EXCLUDED RESONANCES .....')
                            Iif(J) = resInfo%getChannelFitOption(j-1)
                         end if
                      END DO
-                     Iif(6) = Iflpol(1,N)
-                     Iif(7) = Iflpol(2,N)
+                     Iif(6) = resInfo%getChannelFitOption(3)
+                     Iif(7) = resInfo%getChannelFitOption(4)
                      pkenPr = getPkenPr(resonance%getEres(), reduced) 
                      CALL Setflg (Iif, Mmax+4)
                      CALL Pwrite (pkenPr, resonance,
diff --git a/sammy/src/fin/mfin0.f90 b/sammy/src/fin/mfin0.f90
index 006c2e75b..148eb6169 100644
--- a/sammy/src/fin/mfin0.f90
+++ b/sammy/src/fin/mfin0.f90
@@ -133,7 +133,7 @@ module fin
             Nowrt = 1
             CALL Wrcov   (A_Iprbrd , I_Iflbrd , &
                A_Iprdet , I_Ifldet , I_Iigrde , &
-               A_Ipolar , I_Iflpol , &
+               A_Ipolar ,   &
                A_Iprmsc , I_Iflmsc , I_Irdmsc , I_Ijkmsc , A_Ietaee , &
                A_Iprpmc , I_Iflpmc , I_Isopmc , &
                A_Iprorr , I_Iflorr , A_Icrnch , A_Iedets , &
@@ -214,7 +214,6 @@ module fin
                CALL Oldord ( A_Iprbrd , I_Iflbrd , A_Idebrd ,  &
                   I_Ifexcl , &                  
                   A_Iprdet , I_Ifldet , A_Idedet , I_Iigrde , &
-                  I_Iflpol , &
                   A_Iprmsc , I_Iflmsc , A_Idemsc , I_Irdmsc , I_Ijkmsc , &
                              A_Ietaee , &
                   A_Iprpmc , I_Iflpmc , A_Idepmc , I_Isopmc , &
diff --git a/sammy/src/fin/mfin3.f90 b/sammy/src/fin/mfin3.f90
index f9ecce0b9..77babde1c 100644
--- a/sammy/src/fin/mfin3.f90
+++ b/sammy/src/fin/mfin3.f90
@@ -7,7 +7,6 @@ module fin3
       SUBROUTINE Oldord (Parbrd, Iflbrd, Delbrd, &
          If_Excl, &
          Pardet, Ifldet, Deldet, Igrdet,                 &
-         Iflpol,  &
          Parmsc, Iflmsc, Delmsc, Iradms, Ijkmsc, Etaeee, &
          Parpmc, Iflpmc, Delpmc, Isopmc, &
          Parorr, Iflorr, Delorr, Ecrnch, Endets, Sigdts, &
@@ -37,7 +36,6 @@ module fin3
       DIMENSION Parbrd(*), Iflbrd(*), Delbrd(*),  &
          If_Excl(Ntotc,*), &
          Pardet(*), Ifldet(*), Deldet(*), Igrdet(*), &
-         Iflpol(2,*), &
          Parmsc(*), Iflmsc(*), Delmsc(*), Iradms(*), Ijkmsc(*), &
                     Etaeee(*), &
          Parpmc(4,*), Iflpmc(4,*), Delpmc(4,*), Isopmc(*), &
@@ -83,7 +81,7 @@ module fin3
       END IF
 !
       IF (resParData%getNumResonances().GT.0) THEN
-         CALL Ordres (Iflpol,      &
+         CALL Ordres (             &
             Uncxxx, Ratio, If_Pub, &
             relevantData,  resetFitFlags)
       ELSE
@@ -207,7 +205,7 @@ module fin3
 !
 ! ----------------------------------------------------------------
 !
-      SUBROUTINE Ordres (  Iflpol, &
+      SUBROUTINE Ordres ( &
          Uncxxx, Ratio, &
          If_Pub, relevantData, resetFitFlags)
 !
@@ -236,8 +234,7 @@ module fin3
       type(SammyRMatrixParameters)::relevantData
       type(ResonanceCovariance)::physCov
       integer,allocatable::ifits(:)
-      DIMENSION  Iflpol(2,*), &
-                 Uncxxx(Ntotc2,*), Ratio(*)
+      DIMENSION  Uncxxx(Ntotc2,*), Ratio(*)
       Tab = CHAR(9)
       Ifk = If_Pub*Kpblsh
       IF (Ifk.GT.0) then
@@ -263,12 +260,7 @@ module fin3
 
          
 !     rewrite fit flag data to parameter flags, i.e. 0, 1, 3
-!     insteat of number of varied parameter         
-         IF (Kpolar.NE.0.and.resetFitFlags) THEN
-             DO I=1,2
-                IF (Iflpol(I,N).GT.0) Iflpol(I,N) = 1
-             END DO
-         END IF
+!     insteat of number of varied parameter               
          
          Ig = resInfo%getSpinGroupIndex()
          call resParData%getSpinGroupInfo(spinInfo, Ig)
@@ -316,12 +308,6 @@ module fin3
             IgroupW = -1 * IgroupW
          end if
          IF (IgroupW.LT.-9) IgroupW = Kkkgrp - IgroupW
-         IF (Kpolar.NE.0.and.n2.ge.5) THEN
-            Ifits(4) = Iflpol(1,N)
-            Ifits(5) = Iflpol(2,N)
-            if (Ifits(4).gt.0) Ifits(4) = 1
-            if (Ifits(5).gt.0) Ifits(5) = 1
-         END IF
 
          if(resetFitFlags) then
             call resInfo%setEnergyFitOption(ifits(1))
diff --git a/sammy/src/fin/mfin4.f90 b/sammy/src/fin/mfin4.f90
index caece6356..eefb7e41e 100644
--- a/sammy/src/fin/mfin4.f90
+++ b/sammy/src/fin/mfin4.f90
@@ -60,7 +60,7 @@ module fin4
 ! write partial width data (indicated by last .false. in the argument list)         
       CALL Outpar (   Iftit, Nunit    , A_Iprbrd , I_Iflbrd , &
                              A_Iprdet , I_Ifldet , I_Iigrde , &
-       A_Ipolar , I_Iflpol , A_Iprmsc , I_Iflmsc , &
+       A_Ipolar , A_Iprmsc , I_Iflmsc , &
        I_Irdmsc , I_Ijkmsc , A_Ietaee , &
        A_Iprpmc , I_Iflpmc , I_Isopmc , A_Iprorr , &
        I_Iflorr , A_Icrnch , A_Iedets , A_Iseses , A_Iesese , A_Isgdts , &
@@ -76,7 +76,7 @@ module fin4
 ! write reduced width data (indicated by last .true. in the argument list)         
       CALL Outpar (   Iftit, Nunit    , A_Iprbrd , I_Iflbrd , &
                              A_Iprdet , I_Ifldet , I_Iigrde , &
-       A_Ipolar , I_Iflpol , A_Iprmsc , I_Iflmsc , &
+       A_Ipolar , A_Iprmsc , I_Iflmsc , &
        I_Irdmsc , I_Ijkmsc , A_Ietaee , &
        A_Iprpmc , I_Iflpmc , I_Isopmc , A_Iprorr , &
        I_Iflorr , A_Icrnch , A_Iedets , A_Iseses , A_Iesese , A_Isgdts , &
diff --git a/sammy/src/fin/mfin5.f90 b/sammy/src/fin/mfin5.f90
index 917642433..556ca074c 100644
--- a/sammy/src/fin/mfin5.f90
+++ b/sammy/src/fin/mfin5.f90
@@ -439,7 +439,7 @@ module fin5
 !
       SUBROUTINE Wrcov ( Parbrd, Iflbrd, &
          Pardet, Ifldet, Igrdet,         &
-         Polar , Iflpol, Parmsc, Iflmsc, Iradms, Ijkmsc, Etaeee, &
+         Polar , Parmsc, Iflmsc, Iradms, Ijkmsc, Etaeee, &
          Parpmc, Iflpmc, Isopmc, &
          Parorr, Iflorr, Ecrnch, Endets, Sesese, Eseses, Sigdts, &
          Parrpi, Iflrpi, Parudr, Ifludr, Parnbk, Iflnbk, &
@@ -470,7 +470,7 @@ module fin5
       type(ResonanceCovariance)::physCov
       DIMENSION Parbrd(*), Iflbrd(*), &
          Pardet(*), Ifldet(*), Igrdet(*), &
-         Polar(2,*), Iflpol(2,*), &
+         Polar(2,*),     &
          Parmsc(*), Iflmsc(*), Iradms(*), Ijkmsc(*), Etaeee(*), &
          Parpmc(4,*), Iflpmc(4,*), Isopmc(*), &
          Parorr(*), Iflorr(*), Ecrnch(*), Endets(*), &
@@ -484,7 +484,7 @@ module fin5
 !
 !      DIMENSION Parbrd(Numbrd), Iflbrd(Numbrd),
 !     *   Pardet(Numdet), Ifldet(Numdet), Igrdet(Ngroup),
-!     *   Polar(2,Nres), Iflpol(2,Nres),
+!     *   Polar(2,Nres),
 !     *   Parmsc(Nummsc), Iflmsc(Nummsc), Iradms(Ngroup), Ijkmsc(Nummsc),
 !     *   Parpmc(4,Numpmc), Iflpmc(4,Numpmc), Isopmc(Numpmc),
 !     *   Parorr(Numorr), Iflorr(Numorr), Ecrnch(Numorr-11),
@@ -561,7 +561,7 @@ module fin5
                     i = 1, covData%getNumParam())
 !
       IF (resParData%getNumResonances().GT.0) THEN
-         allocate(fitFlags(ntotc2, resParData%getNumResonances()))
+         allocate(fitFlags(max(ntotc2,2), resParData%getNumResonances()))
          fitFlags(:,:) = 0
          iresTmpSize = (Ntotc*(Ntotc+1))/2
          allocate(resTmpData(iresTmpSize,resParData%getNumResonances()))
@@ -678,7 +678,12 @@ module fin5
             WRITE (Iu64) (( Polar(J,I),J=1,2),I=1, &
                resParData%getNumResonances())
 !#            WRITE (Iu64) 'Iflpol', Nres, 2
-            WRITE (Iu64) ((Iflpol(J,I),J=1,2),I=1, &
+            do Ires=1,resParData%getNumResonances()
+                call resParData%getResonanceInfo(resInfo,  Ires)
+                fitFlags(1,Ires) = resInfo%getChannelFitOption(3)
+                fitFlags(2,Ires) = resInfo%getChannelFitOption(4)
+            end do
+            WRITE (Iu64) ((fitFlags(J,I),J=1,2),I=1, &
                 resParData%getNumResonances())
          END IF
          deallocate(fitFlags)
diff --git a/sammy/src/inp/minp01.f b/sammy/src/inp/minp01.f
index 93a157243..d6d0cd32f 100644
--- a/sammy/src/inp/minp01.f
+++ b/sammy/src/inp/minp01.f
@@ -383,7 +383,6 @@ C *** one PAR
          N = 2*Nres
          IF (Kpolar.EQ.0) N = 1
          call make_A_Ipolar(N)
-         call make_I_Iflpol(N)
       END IF
 C
 C *** two PAR
diff --git a/sammy/src/mlb/MlbwCrossCalc_M.f90 b/sammy/src/mlb/MlbwCrossCalc_M.f90
index 097e0105f..79e3c91e4 100644
--- a/sammy/src/mlb/MlbwCrossCalc_M.f90
+++ b/sammy/src/mlb/MlbwCrossCalc_M.f90
@@ -25,20 +25,20 @@ module MlbwCrossCalc_M
   end type MlbwCrossCalc
 
 contains
-subroutine MlbwCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero)
+subroutine MlbwCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
     class(MlbwCrossCalc) :: this
     type(SammyRMatrixParameters)::pars
     type(CovarianceData)::cov
     type(AdjustedRadiusData)::rad
     integer,intent(in)::niso, Itzero, Ilzero
-    logical,intent(in)::needAngular
+    logical,intent(in)::needAngular, doShiftRes
 
     type(VariedParameterInfo)::info
     integer::itot, ipup, ntot, ntotTrig
     logical(C_BOOL)::countCombined
 
 
-    call CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero)
+    call CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
 
     call info%initialize(pars, cov, rad)
 
diff --git a/sammy/src/new/mnew0.f b/sammy/src/new/mnew0.f
index c95877b62..be2e49e2c 100644
--- a/sammy/src/new/mnew0.f
+++ b/sammy/src/new/mnew0.f
@@ -384,7 +384,6 @@ C     write the par file but with the new fit values
       CALL Oldord (    A_Iprbrd , I_Jflbrd , A_Idebrd ,
      *      I_Ifexcl ,
      *      A_Iprdet , I_Jfldet , A_Idedet , I_Iigrde ,
-     *                 I_Iflpol ,
      *      A_Iprmsc , I_Jflmsc , A_Idemsc , I_Irdmsc , I_Ijkmsc ,
      *                 A_Ietaee ,
      *      A_Iprpmc , I_Jflpmc , A_Idepmc , I_Isopmc ,
diff --git a/sammy/src/old/ReadCovarianceParts.f90 b/sammy/src/old/ReadCovarianceParts.f90
index 40e7c0a21..1d3885a8f 100644
--- a/sammy/src/old/ReadCovarianceParts.f90
+++ b/sammy/src/old/ReadCovarianceParts.f90
@@ -463,13 +463,11 @@ contains
      !! @param old  are the data in old format 
      !! @param resParData the resonance data object to fill
      !! @param Polar  the polar values
-     !! @param iflpol the flags for the polar values
      !!
-     subroutine readResonanceData(old, resParData, Polar, iflpol)
+     subroutine readResonanceData(old, resParData, Polar)
        use fixedi_m
        type(SammyRMatrixParameters)::resParData
        real(kind=8):: Polar(2,*)
-       integer::Iflpol(2,*)
        logical::old
 
        integer,allocatable:: fitFlags(:,:), igroup(:)
@@ -487,7 +485,7 @@ contains
        call readArrayData(old, ggaTmp, nres, iu62)
        call readArrayData2d(old, gsiTmp, Ntotc, Nres, iu62)
 
-       allocate(fitFlags(Ntotc2, Nres))
+       allocate(fitFlags(max(2,Ntotc2), Nres))   ! make space for polar fi
        allocate(igroup(Nres))
        READ (Iu62) ((fitFlags(M,I),M=1,Ntotc2),I=1,Nres)
        READ (Iu62) (Igroup(I),I=1,Nres)
@@ -556,7 +554,6 @@ contains
 
        deallocate(gbetprTmp)
        deallocate(betaprTmp)
-       deallocate(fitFlags)
        deallocate(igroup)
        deallocate(eres)
        deallocate(ggaTmp)
@@ -565,8 +562,14 @@ contains
        ! Polar versions of fission widths, if they exist
        IF (Kpolar.EQ.1) THEN
           call readArrayData2d(old, Polar,   2, Nres, iu62)
-          READ (Iu62) ((Iflpol(M,I),M=1,2),I=1,Nres)
-        END IF       
+          READ (Iu62) ((fitFlags(M,I),M=1,2),I=1,Nres)
+          do i = 1, nres
+             call resParData%getResonanceInfo(resInfo,  I)
+             call resInfo%setChannelFitOption(3, fitFlags(1,I))
+            call resInfo%setChannelFitOption(4, fitFlags(2,I))
+          end do
+        END IF
+         deallocate(fitFlags)
      end subroutine readResonanceData
      
      !!
diff --git a/sammy/src/old/mold0.f b/sammy/src/old/mold0.f
index 146c12124..4494af6ea 100644
--- a/sammy/src/old/mold0.f
+++ b/sammy/src/old/mold0.f
@@ -58,7 +58,7 @@ C ***       previous run, including covariance matrices & parameter values
 C ***    Read the COVariance file, extract information
          CALL Rdcov (  A_Iprbrd , I_Iflbrd , 
      *      A_Iprdet , I_Ifldet , I_Iigrde ,
-     *      A_Ipolar , I_Iflpol ,
+     *      A_Ipolar ,
      *      A_Iprmsc , I_Iflmsc , I_Irdmsc , I_Ijkmsc , A_Ietaee ,
      *      A_Iprpmc , I_Iflpmc , I_Isopmc ,
      *      A_Iprorr , I_Iflorr , A_Icrnch , A_Iedets ,
@@ -78,7 +78,7 @@ C ***    Note that call to Rdcovx should not need updating because obsolete
          Nsingl = Nsingl*2
 C ***    Read the very old COVariance file, extract information
          CALL Rdcovx ( A_Iprbrd , I_Iflbrd ,
-     *      A_Ipolar , I_Iflpol , A_Iprmsc , I_Iflmsc , I_Irdmsc ,
+     *      A_Ipolar , A_Iprmsc , I_Iflmsc , I_Irdmsc ,
      *      A_Iprpmc , I_Iflpmc , I_Isopmc ,
      *      A_Iprorr , I_Iflorr , A_Icrnch , A_Iedets ,
      *                 A_Iseses , A_Iesese , A_Isgdts ,
diff --git a/sammy/src/old/mold2.f b/sammy/src/old/mold2.f
index 545069ff4..99d43368a 100644
--- a/sammy/src/old/mold2.f
+++ b/sammy/src/old/mold2.f
@@ -99,7 +99,7 @@ C --------------------------------------------------------------
 C
       SUBROUTINE Rdcov ( Parbrd, Iflbrd, 
      *   Pardet, Ifldet, Igrdet,
-     *   Polar , Iflpol, Parmsc, Iflmsc, Iradms, Ijkmsc, Etaeee,
+     *   Polar , Parmsc, Iflmsc, Iradms, Ijkmsc, Etaeee,
      *   Parpmc, Iflpmc, Isopmc,
      *   Parorr, Iflorr, Ecrnch, Endets,
      *                   Sesese, Eseses, Sigdts,
@@ -155,7 +155,7 @@ C
 C#      CHARACTER*6 What
       DIMENSION Parbrd(*), Iflbrd(*),
      *   Pardet(*), Ifldet(*), Igrdet(*),
-     *   Polar(2,*), Iflpol(2,*),
+     *   Polar(2,*),
      *   Parmsc(*), Iflmsc(*), Iradms(*), Ijkmsc(*), Etaeee(*),
      *   Parpmc(4,*), Iflpmc(4,*), Isopmc(*),
      *   Parorr(*), Iflorr(*), Ecrnch(*), Endets(*), Sesese(*),
@@ -168,7 +168,7 @@ C#      CHARACTER*6 What
 C
 C      DIMENSION Parbrd(Numbrd), Iflbrd(Numbrd), 
 c     *   Pardet(Numdet), Ifldet(Numdet), Igrdet(Ngroup),
-C     *   Polar(2,Nres), Iflpol(2,Nres),
+C     *   Polar(2,Nres),
 C     *   Parmsc(Nummsc), Iflmsc(Nummsc), Iradms(Ngroup), Ijkmsc(Nummsc),
 C     *   Parpmc(Numpmc), Iflpmc(Numpmc), Isopmc(Numpmc),
 C     *   Parorr(Numorr), Iflorr(Numorr), Ecrnch(Numorr-11),
@@ -205,8 +205,7 @@ C *** RESONANCE PARAMETERS
       IF (Nres.GT.0) THEN
          call readResonanceData(.false.,
      *        resParData,
-     *        Polar,
-     *        iflpol)
+     *        Polar)
       END IF
 C
 C *** R-EXTERNAL PARAMETERS
diff --git a/sammy/src/old/mold5.f b/sammy/src/old/mold5.f
index c2830063a..d03f439f5 100644
--- a/sammy/src/old/mold5.f
+++ b/sammy/src/old/mold5.f
@@ -4,7 +4,7 @@ C --------------------------------------------------------------
 C
       SUBROUTINE  Rdcovx ( Parbrd, Iflbrd, 
      *     Iflext,
-     *     Polar , Iflpol, Parmsc, Iflmsc, Iradms,
+     *     Polar , Parmsc, Iflmsc, Iradms,
      *     Parpmc, Iflpmc, Isopmc,
      *     Parorr, Iflorr, Ecrnch, Endets,
      *                     Sesese, Eseses, Sigdts,
@@ -44,7 +44,7 @@ C this is not a common block, this contains functions
       real(kind=8)  :: dum1
 C
       DIMENSION Parbrd(*), Iflbrd(*),
-     *   Polar(2,*), Iflpol(2,*), Parmsc(*), Iflmsc(*), Iradms(*),
+     *   Polar(2,*), Parmsc(*), Iflmsc(*), Iradms(*),
      *   Parpmc(4,*), Iflpmc(4,*), Isopmc(*),
      *   Parorr(*), Iflorr(*), Ecrnch(*), Endets(*), Sesese(*),
      *   Eseses(*), Sigdts(*), Parrpi(*), Iflrpi(*),
@@ -95,8 +95,7 @@ C *** isotopic abundances etc.
       IF (Nres.GT.0) THEN
            call readResonanceData(.true.,
      *        resParData,
-     *        Polar,
-     *        iflpol)
+     *        Polar)
       END IF
 C
 C *** R-EXTERNAL PARAMETERS
diff --git a/sammy/src/par/mpar0.f90 b/sammy/src/par/mpar0.f90
index 63610c673..21a8d097b 100644
--- a/sammy/src/par/mpar0.f90
+++ b/sammy/src/par/mpar0.f90
@@ -141,14 +141,12 @@ module par_m
 ! *** Sub routine Order reorders resonances according to J-Pi groups
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >
 !
-      IF (Kpolar.EQ.1.AND. Nres.NE.0) CALL Fixpol (I_Iflpol )
       CALL Fix (I_Iflbrd ,                                      &
          I_Ifldet , I_Iflmsc ,                                  &
          I_Irdmsc , I_Iflpmc , I_Iflorr , I_Iflrpi , I_Iflnbk , &
          I_Iflbgf , I_Ifldtp , I_Iflusd)
 ! *** Routine Fix sets Flag = (Parameter Number)
 !
-      IF (Kpolar.EQ.1.AND. Nres.NE.0) CALL Flgpol (I_Iflpol )
 !
       IF (Kdecpl.NE.0.AND. Nres.NE.0) CALL Dddcov
 ! *** Routine Dddcov  updates zero uncertainties in resInfo%getXVal
diff --git a/sammy/src/par/mpar11.f90 b/sammy/src/par/mpar11.f90
index 53eddb915..64dffc8fe 100644
--- a/sammy/src/par/mpar11.f90
+++ b/sammy/src/par/mpar11.f90
@@ -282,81 +282,6 @@ module par11_m
       END
 !
 !
-! --------------------------------------------------------------
-!
-      SUBROUTINE Fixpol (Iflpol)
-! *** Purpose -- Generate flags for polar version of fission widths
-!
-      use fixedi_m
-      use SammyResonanceInfo_M
-      use RMatResonanceParam_M
-      use EndfData_common_m
-      IMPLICIT NONE
-!
-      INTEGER::Iflpol(2,*)
-
-      type(SammyResonanceInfo)::resInfo
-      type(RMatResonance)::resonance
-      integer::n, numChan
-
-      DO N=1,resParData%getNumResonances()
-         call resParData%getResonanceInfo(resInfo, N)
-         call resParData%getResonance(resonance, resInfo)
-         numChan = resonance%getNumChan()
-
-         Iflpol(1,N) = 0
-         Iflpol(2,N) = 0
-
-! if lrf=3, channel 3 and 4 are fission width 1 and 2
-!     Sammy here assumes the same is true if setting iflpol
-!        resInfo%getChannelFitOption(1) -> gamma
-!        resInfo%getChannelFitOption(2) -> capture
-!        resInfo%getChannelFitOption(3) -> fission 1
-!        resInfo%getChannelFitOption(4) -> fission 2
-         if (numChan.ge.3) Iflpol(1,N) = resInfo%getChannelFitOption(3)
-         if (numChan.ge.4) Iflpol(2,N) = resInfo%getChannelFitOption(4)
-         
-         IF (Iflpol(1,N).NE.0 .OR. Iflpol(2,N).NE.0) then
-            if (numChan.ge.3)  call resInfo%setChannelFitOption(3, 1)
-            if (numChan.ge.4)  call resInfo%setChannelFitOption(4, 1)
-         end if
-      END DO
-      RETURN
-      END
-!
-!
-! --------------------------------------------------------------
-!
-      SUBROUTINE Flgpol (Iflpol)
-!
-! *** Purpose -- fix flags for polar version of fission widths
-!
-      use SammyResonanceInfo_M
-      use RMatResonanceParam_M
-      use EndfData_common_m
-      IMPLICIT none
-!
-      integer:: Iflpol(2,*)
-
-      type(SammyResonanceInfo)::resInfo
-      type(RMatResonance)::resonance
-      integer::n, numChan
-
-      DO N=1,resParData%getNumResonances()
-         call resParData%getResonanceInfo(resInfo, N)
-         call resParData%getResonance(resonance, resInfo)
-         numChan = resonance%getNumChan()
-         
-         IF (Iflpol(1,N).NE.0.and.numChan.ge.3) then
-            Iflpol(1,N) = resInfo%getChannelFitOption(3)
-         end if
-         IF (Iflpol(2,N).NE.0.and.numChan.ge.3) then
-            Iflpol(2,N) = resInfo%getChannelFitOption(4)
-         end if
-      END DO
-      RETURN
-      END
-!
 !
 ! --------------------------------------------------------------
 !
diff --git a/sammy/src/rec/mrec0.f b/sammy/src/rec/mrec0.f
index 8b6f5a9f9..e5a507516 100644
--- a/sammy/src/rec/mrec0.f
+++ b/sammy/src/rec/mrec0.f
@@ -17,11 +17,13 @@ C
       use lbro_common_m
       use AllocateFunctions_m
       use GridData_M
+      use EndfData_common_m
       use mrec8_m
       use mdata_M
       use mrec2_m
       use mthe0_M
       use Tab1_M
+      use mcro8_m
       use array_sizes_common_m, only : zeroKCalc, zeroKCalcInit
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
       type(GridData)::energb
@@ -29,7 +31,7 @@ C
       real(kind=8),allocatable,dimension(:)::A_Isig2(:,:)
       type(Tab1)::tab1Data
       real(kind=8)::val
-      real(kind=8),allocatable,dimension(:)::A_ix1, A_Ienode, A_Iwnode
+      real(kind=8),allocatable,dimension(:)::A_Ienode, A_Iwnode
 C
 C
       WRITE (6,99999)
@@ -42,7 +44,7 @@ C
       CALL Initil
 C
 C *** find array sizes for xx
-      CALL Findmx (Mxany)
+      Mxany = resparData%getNumResonances()
       call setup_zeroK()
 
 C
@@ -72,15 +74,7 @@ C
 C *** two ***
       Maxnpu = Nemax
 C
-      call allocate_real_data(A_Ixx, Mres)
 C - - - - - - - - - - - - - - - - - - - - - - - - - - - <
-C *** three ***
-
-      call allocate_real_data(A_Ix1, max(Na, Mxany))
-      CALL Fixx (A_Ibound , A_Iechan ,
-     *   A_Izke , A_Ixx, A_Ix1 , Mxany)
-C *** SBROUTNE Fixx sets Xx = energy shift
-      deallocate(A_Ix1)
 C - - - - - - - - - - - - - - - - - - - - - - - - - - - <
 C
 C *** five ***
diff --git a/sammy/src/rec/mrec2.f90 b/sammy/src/rec/mrec2.f90
index 8d2eba94f..79a0b1a45 100644
--- a/sammy/src/rec/mrec2.f90
+++ b/sammy/src/rec/mrec2.f90
@@ -18,7 +18,7 @@ contains
       use fixedi_m
       use fixedr_m
       use AllocateFunctions_m
-      use mrec6_m
+      use mrec6_m      
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 !
       DIMENSION Crss(*), Enode(*), Widnod(*)
diff --git a/sammy/src/rec/mrec3.f b/sammy/src/rec/mrec3.f
index 981894b9f..2e27395e3 100644
--- a/sammy/src/rec/mrec3.f
+++ b/sammy/src/rec/mrec3.f
@@ -15,6 +15,7 @@ C
       use varyr_common_m
       use ifsubs_common
       use xct4_m
+      use array_sizes_common_m, only : zeroKCalc
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 C
       DIMENSION Ssseee(*), Crss(*)
@@ -31,7 +32,7 @@ C
       Ifradt = 1
       I = 0
 C *** generate cross sections pieces
-      CALL Zcross (Nnndrc, I, Kount_Helmut)
+      CALL Zcross (zeroKCalc%driver%xctCalc, Nnndrc, I, Kount_Helmut)
 C
 C *** set the individual cross sections
       CALL Indivi (Crss, Ssseee, Su)
diff --git a/sammy/src/ref/mrfs0.f b/sammy/src/ref/mrfs0.f
index 608804d7b..314ee38c5 100644
--- a/sammy/src/ref/mrfs0.f
+++ b/sammy/src/ref/mrfs0.f
@@ -55,7 +55,6 @@ C
 C *** collect values needed for parameter file
       Korder = Idimen (1, 1, '1, 1')
       CALL Cnvpar ( A_Iprbrd , I_Iflbrd ,
-     *   I_Iflpol ,
      *   A_Iprmsc , I_Iflmsc , A_Idemsc , I_Irdmsc ,
      *   A_Iprpmc , I_Iflpmc , A_Idepmc , I_Isopmc ,
      *   A_Iprorr , I_Iflorr , A_Ideorr , A_Icrnch , A_Iedets ,
diff --git a/sammy/src/ref/mrfs2.f b/sammy/src/ref/mrfs2.f
index 112053148..9e17e7d13 100644
--- a/sammy/src/ref/mrfs2.f
+++ b/sammy/src/ref/mrfs2.f
@@ -166,7 +166,6 @@ C *** one
          N = 2*Nres
          IF (Kpolar.EQ.0) N = 1
          call make_A_Ipolar(N)
-         call make_I_Iflpol(N)
       END IF
 C
 C *** two
diff --git a/sammy/src/ref/mrfs4.f b/sammy/src/ref/mrfs4.f
index 0ec9e41bf..cedf71538 100644
--- a/sammy/src/ref/mrfs4.f
+++ b/sammy/src/ref/mrfs4.f
@@ -3,7 +3,6 @@ C
 C --------------------------------------------------------------
 C
       SUBROUTINE Cnvpar (  Parbrd, Iflbrd,
-     *     Iflpol,
      *     Parmsc, Iflmsc, Delmsc, Iradms,
      *     Parpmc, Iflpmc, Delpmc, Isopmc,
      *     Parorr, Iflorr, Delorr, Ecrnch, Endets, Sigdts,
@@ -35,7 +34,6 @@ C
       type(SammyIsoInfo)::isoInfo
       logical(C_BOOL)::reduced, combined
       DIMENSION                  Parbrd(*), Iflbrd(*),
-     *              Iflpol(2,*),
      *     Parmsc(*), Iflmsc(*), Delmsc(*), Iradms(*),
      *     Parpmc(4,*), Iflpmc(4,*), Delpmc(4,*), Isopmc(*),
      *     Parorr(*), Iflorr(*), Delorr(*), Ecrnch(*),
diff --git a/sammy/src/ref/mwrt0.f b/sammy/src/ref/mwrt0.f
index 185f45a75..b2238e76e 100644
--- a/sammy/src/ref/mwrt0.f
+++ b/sammy/src/ref/mwrt0.f
@@ -50,7 +50,6 @@ C *** Write Parameter file
                CALL Oldord ( A_Iprbrd , I_Iflbrd , A_Idebrd , 
      *            I_Ifexcl ,
      *            A_Iprdet , I_Ifldet , A_Idedet , I_Iigrde ,
-     *            I_Iflpol ,
      *            A_Iprmsc , I_Iflmsc , A_Idemsc , I_Irdmsc , I_Ijkmsc ,
      *                       A_Ietaee ,
      *            A_Iprpmc , I_Iflpmc , A_Idepmc , I_Isopmc ,
diff --git a/sammy/src/sammy/CMakeLists.txt b/sammy/src/sammy/CMakeLists.txt
index 3c6ff618d..5e20453d1 100644
--- a/sammy/src/sammy/CMakeLists.txt
+++ b/sammy/src/sammy/CMakeLists.txt
@@ -78,6 +78,7 @@ APPEND_SET(SAMMY_SOURCES
             ../cro/mnrm1.f90
             ../cro/mnrm2.f90
             ../cro/CroCrossCalc_M.f90
+            ../cro/CroCrossCalcImpl_M.f90
 
             ../dat/mdat0.f90
             ../dat/mdat1.f90
@@ -505,6 +506,7 @@ APPEND_SET(SAMMY_SOURCES
             ../xct/mxct31.f90
             ../xct/mxct32.f90
             ../xct/XctCrossCalc_M.f90
+            ../xct/XctCrossCalcImpl_M.f90
 
             ../xxx/mxxxa.f90
             ../xxx/mxxxb.f90
diff --git a/sammy/src/the/CrossSectionCalcDriver_M.f90 b/sammy/src/the/CrossSectionCalcDriver_M.f90
index e4961fa02..0819414fe 100644
--- a/sammy/src/the/CrossSectionCalcDriver_M.f90
+++ b/sammy/src/the/CrossSectionCalcDriver_M.f90
@@ -1,8 +1,8 @@
 module CrossSectionCalcDriver_M
    use CrossSectionCalculator_M
    use MlbwCrossCalcImpl_m
-   use CroCrossCalc_M
-   use XctCrossCalc_M
+   use CroCrossCalcImpl_M
+   use XctCrossCalcImpl_M
    use ArtificalCross_M
    use DerivativeHandler_M
    use GridData_M
@@ -21,8 +21,8 @@ module CrossSectionCalcDriver_M
 
      type(ArtificalCross),pointer::artificalCalc => null()
      type(MlbwCrossCalcImpl),pointer,private::mlbwCalc => null()
-     type(CroCrossCalc),pointer,private::croCalc => null()
-     type(XctCrossCalc),pointer,private::xctCalc => null()
+     type(CroCrossCalcImpl),pointer,private::croCalc => null()
+     type(XctCrossCalcImpl),pointer::xctCalc => null()  ! note: private will come back after conversion of xct  codes
 
      logical::wantCrossPerSpin = .false.  ! do we want cross section per spin group
      type(GridData)::crossPerSpin         ! if used: energy, group 1, group 1 + group 2, group 1 + group 2 + group 3,  ... repeated for each Nnnsig value )
@@ -34,6 +34,8 @@ module CrossSectionCalcDriver_M
      integer::Kssmsc ! is multiple scattering desired later in the code (affects which 0K data are needed)
      integer::Itzero, Ilzero  ! fit flag for t0
 
+     integer::Kaptur,  Kaverg, Krecon
+
      contains
      procedure, pass(this) :: initialize => CrossSectionCalcDriver_initialize
      procedure, pass(this) :: setUpCalculator => CrossSectionCalcDriver_setUpCalculator  ! select the type of cross section reconstructor
@@ -44,6 +46,9 @@ module CrossSectionCalcDriver_M
      procedure, pass(this) :: getWantDerivs => CrossSectionCalcDriver_getWantDerivs  ! do we want derivatives
      procedure, pass(this) :: getNumVaried => CrossSectionCalcDriver_getNumVaried    ! get the number of total varied parameters
      procedure, pass(this) :: setWantCrossPerSpin => CrossSectionCalcDriver_setWantCrossPerSpin  ! indicate whether we want cross section per spin
+
+     procedure, pass(this) :: setAddtionalParams => CrossSectionCalcDriver_setAddtionalParams   ! set additional parameters important to ony some formalisms
+                                                                                                ! call after first call to setUpCalculator
   end type CrossSectionCalcDriver
 
   private setup_desiredCross
@@ -84,15 +89,16 @@ module CrossSectionCalcDriver_M
      !!  - 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
      !!
-     subroutine CrossSectionCalcDriver_setUpCalculator(this, Krmatx,Kkclqx, Kkkclq, niso, covData, resParData, radFitFlags,Ilzero, Itzero, wantDeriv)
+     subroutine CrossSectionCalcDriver_setUpCalculator(this, Krmatx,Kkclqx, Kkkclq, niso, covData, resParData, &
+                                                       radFitFlags,Ilzero, Itzero, wantDeriv, doShiftRes, Kpolar)
         class(CrossSectionCalcDriver)::this
         type(SammyRMatrixParameters)::resParData
         type(CovarianceData)::covData
         type(AdjustedRadiusData)::radFitFlags
-        integer::Ilzero, Itzero, Krmatx, Kkclqx, Kkkclq, niso
-        logical::wantDeriv
-        type(VariedParameterInfo)::info
 
+        integer::Ilzero, Itzero, Krmatx, Kkclqx, Kkkclq, niso, Kpolar
+        logical::wantDeriv, doShiftRes
+        type(VariedParameterInfo)::info
 
         integer::nn
         logical::needAngular
@@ -108,13 +114,13 @@ module CrossSectionCalcDriver_M
         if (.not.this%calculatorInit) then
            if (Kkkclq.NE.0) THEN
                 allocate(this%artificalCalc)
-                call this%artificalCalc%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero)
+                call this%artificalCalc%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero, doShiftRes)
                 this%artificalCalc%Kkclqx = Kkclqx
                 this%artificalCalc%Kkkclq = Kkkclq
                 this%calculator => this%artificalCalc
            ELSE IF (Krmatx.EQ.1 .OR. Krmatx.EQ.-1) THEN
                 allocate(this%mlbwCalc)
-                call this%mlbwCalc%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero)
+                call this%mlbwCalc%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero, doShiftRes)
                 if (Krmatx.eq.1)  then
                    this%mlbwCalc%wantMlbw = .true.
                 else if  (Krmatx.eq.-1) then
@@ -125,8 +131,9 @@ module CrossSectionCalcDriver_M
                this%calculator => this%mlbwCalc
            ELSE IF (Krmatx.EQ.0) THEN
               allocate(this%croCalc)
-              call this%croCalc%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero)
+              call this%croCalc%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero, doShiftRes)
               this%calculator => this%croCalc
+              this%croCalc%Kpolar = Kpolar
            ELSE IF (Krmatx.EQ.2) THEN
                 nisoOur = niso
                 if (needAngular) then
@@ -136,8 +143,10 @@ module CrossSectionCalcDriver_M
                     call info%destroy()
                 end if
                 allocate(this%xctCalc)
-                call this%xctCalc%initialize(resParData, covData, radFitFlags, nisoOur, needAngular, Itzero, Ilzero)
+                call this%xctCalc%initialize(resParData, covData, radFitFlags, nisoOur, needAngular, Itzero, Ilzero, doShiftRes)
                 this%calculator => this%xctCalc
+                this%xctCalc%crossDataSelf%instance_ptr = this%calcDataSelf%instance_ptr
+                this%xctCalc%Kpolar = Kpolar
            ELSE
                write(0,*)" Invalid selection for cross section calculator"
                stop
@@ -163,6 +172,17 @@ module CrossSectionCalcDriver_M
         call this%calculator%setEnergyIndependent(this%Kcros, Twomhb, this%kwcoul, Etac)
      end subroutine
 
+     subroutine CrossSectionCalcDriver_setAddtionalParams (this,  lllmax, Kfinit, wantSelfIndicate, 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
+          logical::wantSelfIndicate  ! do we need to calculate self-indicated cross section data
+
+          if (.not.associated(this%xctCalc)) return
+          call this%xctCalc%setAddtionalParams (lllmax, Kfinit, wantSelfIndicate, this%Kssmsc, debug)
+     end subroutine
+
      !!
      !! Calculate the cross section (of the type given by Kcros in the initialization routine)
      !!
@@ -199,16 +219,17 @@ module CrossSectionCalcDriver_M
      !! Set up the desired cross section data if using xct algorithm
      !!
      subroutine setup_desiredCross(this, calc)
-        use fixedi_m, only : Kaptur
-        use ifwrit_m, only : Kaverg, Kfinit, Krecon
-
         class(CrossSectionCalcDriver)::this
         class(XctCrossCalc)::calc
 
-        integer::Ifc3,nreact, I
+        integer::Ifc3,nreact, I, Ifdif
+        logical::doSelfIndicate
+
+        doSelfIndicate = .false.
+        if (associated(this%xctCalc)) doSelfIndicate = this%xctCalc%wantSelfIndicate
 
         calc%Ifcros = .true.
-        if (Krecon.GT.0) return ! want all cross sections
+        if (this%Krecon.GT.0) return ! want all cross sections
 
         nreact = size(calc%Ifcros)
 
@@ -221,8 +242,8 @@ module CrossSectionCalcDriver_M
         ! multiple-scattering corrections  .or. self-shielding
         ! .or. Bondarenko averaging)
         ! for which we need total
-         IF (this%Kcros.EQ.1 .OR. Kfinit.NE.0 .OR. this%Kssmsc.GT.0 &
-            .OR. this%Kssmsc.EQ.-2 .OR. Kaverg.EQ.2) THEN
+         IF (this%Kcros.EQ.1 .OR. doSelfIndicate.OR. this%Kssmsc.GT.0 &
+            .OR. this%Kssmsc.EQ.-2 .OR. this%Kaverg.EQ.2) THEN
              calc%Ifcros(1) = .true.
              calc%Ifcros(2) = .true.
              Ifc3 = 1
@@ -234,7 +255,7 @@ module CrossSectionCalcDriver_M
         ! reaction or inelastic or fission cross section
         IF (this%Kcros.EQ.3) THEN
             Ifc3 = 1
-            IF (Kaptur.EQ.1) calc%Ifcros(2) = .true.
+            IF (this%Kaptur.EQ.1) calc%Ifcros(2) = .true.
                 ! Kaptur=1 means add Reich Moore eliminated capture
                 ! channel information, so really calculate
                 ! [absorp] minus [non-gamma channels]
@@ -278,6 +299,15 @@ module CrossSectionCalcDriver_M
               calc%Ifcros(I) = .true.
            END DO
         END IF
+
+
+        !  differential elastic cross section .or. multiple-scattering
+        ! corrections for which we need differential elastic
+        IF (this%Kcros.EQ.7 .OR. this%Kssmsc.GT.0) Ifdif = 1
+        IF (this%Kcros.EQ.11) Ifdif = 2   !  angular distribution for reaction cross section
+
+        if (associated(this%xctCalc))  this%xctCalc%Ifdif = Ifdif
+        if (associated(this%croCalc))  this%croCalc%Ifdif = Ifdif
      end subroutine
 
 
diff --git a/sammy/src/the/CrossSectionCalculator_M.f90 b/sammy/src/the/CrossSectionCalculator_M.f90
index c3bfb2515..6769fc53f 100644
--- a/sammy/src/the/CrossSectionCalculator_M.f90
+++ b/sammy/src/the/CrossSectionCalculator_M.f90
@@ -22,6 +22,8 @@ module CrossSectionCalculator_M
 
      logical::needAngular  ! do we need angular distributions
 
+     logical::doResShift  = .false.  ! do we want to SHIFT RESONANCE ENERGIES VIA SHIFT FACTOR
+
      integer::numIso   ! the number of isotopes for which cross section are calculated, or 1 if sum over isotopes was already performed
 
      logical::Ifext  ! do we have any external r-matrix parameters to fit
@@ -66,7 +68,7 @@ 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) :: calcCross => CrossSectionCalculator_calcCross     
   end type CrossSectionCalculator
 
 contains
@@ -193,6 +195,7 @@ subroutine CrossSectionCalculator_calcCross(this, ener, Ipoten)
     if ( .not.this%wantDerivs) itot = 1
     call this%crossData%reserveColumns(this%row, itot)
 
+
     this%ener = ener
 
     this%enerSq  = ener
@@ -386,13 +389,13 @@ end subroutine
 !! - needAngular is angular distribution desired
 !! - Itzero fit flag for t0
 !! - Ilzero fit flag for flight path length
-subroutine CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero)
+subroutine CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
     class(CrossSectionCalculator) :: this
     type(SammyRMatrixParameters)::pars
     type(CovarianceData)::cov
     type(AdjustedRadiusData)::rad
     integer,intent(in)::niso, Itzero, Ilzero
-    logical,intent(in)::needAngular
+    logical,intent(in)::needAngular, doShiftRes
 
     type(VariedParameterInfo)::info
     integer::itot, ipup, ntot, itmp
@@ -409,6 +412,8 @@ subroutine CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAng
     this%numIso = niso
     this%row = 1
 
+    this%doResShift = doShiftRes
+
     call info%initialize(pars, cov, rad)
     ! get the number of independently fitted resonance parameters
     countCombined = .false.
diff --git a/sammy/src/the/DerivativeHandler.cpp b/sammy/src/the/DerivativeHandler.cpp
index 13d1c9592..c2e5e6239 100644
--- a/sammy/src/the/DerivativeHandler.cpp
+++ b/sammy/src/the/DerivativeHandler.cpp
@@ -6,6 +6,7 @@ namespace sammy{
        int numIso = params.getNumIso();       
        if (numIso < num) numIso = num;  // this can happen if we want to keep track of channels included in the calculation with unique Echan values
        if (num == 1) numIso = 1;        // indicates that we can already sum and don't need to keep track of isotopes
+
        if (numIso <= 1) {
            if (getGridNumber()  == 0) addGrid();
            return;
@@ -143,7 +144,7 @@ namespace sammy{
       addData(row, col, posOur, value, current);
   }
 
-  int DerivativeHandler::getNumberIsotopes() const{
+  int DerivativeHandler::getNumberIsotopes() const{      
       return getGridNumber();
   }
 
diff --git a/sammy/src/the/ZeroKCrossCorrections_M.f90 b/sammy/src/the/ZeroKCrossCorrections_M.f90
index 27b80ec61..89f4898a0 100644
--- a/sammy/src/the/ZeroKCrossCorrections_M.f90
+++ b/sammy/src/the/ZeroKCrossCorrections_M.f90
@@ -79,7 +79,13 @@ contains
           WRITE (21,10000) ener
        END IF
 
-       if (ener.lt.0.0d0.and..not.this%wantNeg) cycle
+       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
+           end if
+           cycle
+        end if
 
        ipos = ipos + 1
        if (this%dataStart.eq.0) this%dataStart = ipos
diff --git a/sammy/src/the/mthe0.f90 b/sammy/src/the/mthe0.f90
index aaf7ae13c..9c90f37ad 100644
--- a/sammy/src/the/mthe0.f90
+++ b/sammy/src/the/mthe0.f90
@@ -20,7 +20,7 @@ module mthe0_M
       use AuxGridHelper_M
       use mthe1_m
       use Samclq_m
-      use Odfpcs_M
+      use Odfpcs_M      
       use EndfData_common_m, only : covData, resParData, radFitFlags
       use array_sizes_common_m, only : zeroKCalc, zeroKCalcInit,      &
                                        calcData, calcDataSelf,        &
@@ -39,6 +39,7 @@ module mthe0_M
       logical::wantDeriv
       type(SammyGridAccess)::grid, userGrid
       integer::ngroup
+
 !
           WRITE (line,99999)
 99999     FORMAT (' *** SAMMY-THEOR 15 Nov 07 ***')
@@ -121,7 +122,7 @@ module mthe0_M
       ELSE IF (Krmatx.EQ.1 .OR. Krmatx.EQ.-1) THEN
          CALL  Sammlb_0
       ELSE IF (Krmatx.EQ.2) THEN
-         CALL Samxct_0
+         CALL Samxct_0(zeroKCalc%driver%xctCalc)
       ELSE
          WRITE (6,10100)
 10100    FORMAT (' ??? No call here! ???')
@@ -165,15 +166,15 @@ module mthe0_M
       use array_sizes_common_m, only :  zeroKCalc, zeroKCalcInit,       &
                                         calcData, calcDataSelf,         &
                                         calcDataInit, calcDataSelfInit
-      use fixedi_m, only : Numiso, Ktruet, Nnniso, Kaddcr
+      use fixedi_m, only : Numiso, Ktruet, Nnniso, Kaddcr,lllmax, Kshift, Kaptur, Kpolar
       use ifwrit_m, only : Ilzero, Itzero, Kcros, Knocor, Kpoten, Krmatx, Ksolve,  &
                            Kssmsc, kwcoul, Maxwel, Knocor, Kkkdop, Kkclqx, Kkkclq, &
-                           Kartgd, Kpiece
+                           Kartgd, Kpiece, Kfinit, Ksitmp, Ksindi, Kaverg, Krecon
       use fixedr_m, only : emax
-      use lbro_common_m, only : Yangle, Yaverg, Ydoppr, Yresol, Yssmsc, Ytrans
+      use lbro_common_m, only : Yangle, Yaverg, Ydoppr, Yresol, Yssmsc, Ytrans, debug
       use  namfil_common_m, only : Faddcr
       implicit none
-      logical::wantDeriv
+      logical::wantDeriv, wantSelfIndicate, doResShift
 
       ! make sure the object to take the data exists
       if (.not.calcDataInit) then
@@ -193,6 +194,9 @@ module mthe0_M
                       kwcoul, Kssmsc, Kcros)
            zeroKCalc%driverInit = .true.
 
+           zeroKCalc%driver%Kaptur = Kaptur
+           zeroKCalc%driver%Kaverg = Kaverg
+           zeroKCalc%driver%Krecon = Krecon
            zeroKCalcInit = .true.
        end if
 
@@ -215,9 +219,17 @@ module mthe0_M
          else
             wantDeriv = .true.
          end if
+
+         doResShift = .false.
+         if( Kshift.ne.0)  doResShift = .true.
+
          call zeroKCalc%driver%setUpCalculator(Krmatx, Kkclqx, Kkkclq,  Nnniso,  &
-                     covData, resParData,                      &
-                     radFitFlags,Ilzero, Itzero, wantDeriv)
+                     covData, resParData,                                        &
+                     radFitFlags,Ilzero, Itzero, wantDeriv, doResShift, Kpolar)
+         wantSelfIndicate = .false.
+         if (Ksindi.GT.0 .AND. Kcros.EQ.8) wantSelfIndicate = .true.
+         if (Ksitmp.GT.0) wantSelfIndicate = .true.
+         call zeroKCalc%driver%setAddtionalParams(lllmax, Kfinit, wantSelfIndicate, debug)
 
          ! by default we get the abundance (as expected) from the
          ! isotope
diff --git a/sammy/src/the/mthe1.f90 b/sammy/src/the/mthe1.f90
index c62518eb4..9a937be1d 100644
--- a/sammy/src/the/mthe1.f90
+++ b/sammy/src/the/mthe1.f90
@@ -13,8 +13,7 @@ module mthe1_m
       use fixedr_m, only : Emax, Emin
       use ifwrit_m, only : Kdecpl, Kscut, Ksolve, ktzero, Ndat
       use broad_common_m, only : Dopple, Iesopr
-      use templc_common_m, only : I_Inotu, A_Ibr, A_Ibi,   &
-                                  A_Ipr, A_Ipi, Upi, Upr
+      use templc_common_m, only : I_Inotu,   A_Ipr, A_Ipi, Upi, Upr
       use AllocateFunctions_m
       use EndfData_common_m
       use SammyResonanceInfo_M
@@ -203,15 +202,13 @@ module mthe1_m
          end if
 
          ! todo: Sometimes ntotc is set wrong
-         ! and A_Ibr, A_Ibi, A_Ipr, and A_Ipi
+         ! and  A_Ipr, and A_Ipi
          ! are queried up to that number
          if( ntotres.lt.Ntotc) then
             ntotres = Ntotc
          end if
          ntotres = (ntotres*(ntotres+1))/2
          if (Napres.eq.0) ntotres = 0
-         call reallocate_real_data_2d(A_Ibr, ntotres, 0, Napres, 0)
-         call reallocate_real_data_2d(A_Ibi, ntotres, 0, Napres, 0)
          call reallocate_real_data_2d(A_Ipr, ntotres, 0, Napres, 0)
          call reallocate_real_data_2d(A_Ipi, ntotres, 0, Napres, 0)
          call allocate_real_data(Upi, Napres)
diff --git a/sammy/src/xct/XctCrossCalcImpl_M.f90 b/sammy/src/xct/XctCrossCalcImpl_M.f90
new file mode 100644
index 000000000..b573602f4
--- /dev/null
+++ b/sammy/src/xct/XctCrossCalcImpl_M.f90
@@ -0,0 +1,68 @@
+!
+!  Add an extra layer between XctCrossCalc and the
+!  existing xct methods, so that we can pass
+!  a XctCrossCalc object to them
+!  without having to change too many
+!  implementations
+!
+module XctCrossCalcImpl_M
+  use, intrinsic :: ISO_C_BINDING
+  use CrossSectionCalculator_M
+  use DerivativeHandler_M
+  use XctCrossCalc_M
+  use SammyRMatrixParameters_M
+  use CovarianceData_M
+  use AdjustedRadiusData_M
+  implicit none
+
+  type, extends(XctCrossCalc) :: XctCrossCalcImpl
+     contains
+     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_calcCross(this, ener, Ipoten)
+     use xct2_m
+     use varyr_common_m, only : Su, Squ
+     class(XctCrossCalcImpl) ::this
+     real(kind=8)::ener
+     integer::ipoten
+
+
+
+     real(kind=8)::val
+     integer::i, itot, ns, is
+     logical(C_BOOL)::accu
+
+     call CrossSectionCalculator_calcCross(this, ener, Ipoten)
+!     if (ener.eq.0.0d0) then
+
+!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
+
+
+  end subroutine
+
+  subroutine XctCrossCalcImpl_setEnergyIndependent(this, reactType, Twomhb, kwcoul, Etac)
+      use xct1_m
+      class(XctCrossCalcImpl) :: this
+      integer::kwcoul, reactType
+      real(kind=8)::Twomhb, Etac
+
+
+      call XctCrossCalc_setEnergyIndependent(this, reactType, Twomhb, kwcoul, Etac)
+
+      if (this%wantDerivs) then
+         call Babb (this, .true.)
+         CALL Babbga (this, kwcoul )
+      end if
+
+  end subroutine
+
+
+end module XctCrossCalcImpl_M
diff --git a/sammy/src/xct/XctCrossCalc_M.f90 b/sammy/src/xct/XctCrossCalc_M.f90
index 20a0ee143..d392f8e41 100644
--- a/sammy/src/xct/XctCrossCalc_M.f90
+++ b/sammy/src/xct/XctCrossCalc_M.f90
@@ -5,6 +5,7 @@ module XctCrossCalc_M
   use SammyRMatrixParameters_M
   use CovarianceData_M
   use AdjustedRadiusData_M
+  use AllocateFunctions_m
   implicit none
 
   type, extends(CrossSectionCalculator) :: XctCrossCalc
@@ -12,24 +13,159 @@ module XctCrossCalc_M
      ! this is assuming first width is gamma and is not counted, and all other width are as in normal sammy definition
      logical,allocatable,dimension(:)::Ifcros(:)
 
+      type(DerivativeHandler)::crossDataSelf     ! the self indicated cross section and derivative calculated for a given energy
+
+      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
+      integer::ntriag   ! size of lower triangular matrix over width, i.e. (ntotc*(ntotc+1))/2
+
+      logical::debug   ! do we want debug print out
+
+      integer::Kslow = 0     ! should we use the slow or fast version of Clebsch-Gordon calculation
+      integer::C_G_Kxlmn = 1
+
+      integer::Ifdif = 0    ! Ifdif = 1 if need differential elastic, Ifdif = 2 if need differential reaction
+      integer::Ifcoul = 0   ! do we need to calculate coulomb data
+
+      real(kind=8),allocatable,dimension(:,:)::Alj   ! used to count the number of Clebsch-Gordon coefficients
+      real(kind=8),allocatable,dimension(:)::Xx      ! 0.0 or if SHIFT RESONANCE ENERGIES VIA SHIFT FACTOR, the factor
+      real(kind=8),allocatable,dimension(:)::XxHelper    ! helper array to calculate Xx
+
+     real(kind=8),allocatable,dimension(:,:)::Br, Bi    ! energy independent part of derivatives for resonance parameters
+     real(kind=8),allocatable,dimension(:,:,:)::Bga !E ENERGY-INDEPENDENT PORTION OF PARTIAL of R wrt reduced-width-amplitudes (for those which are not varied)
+
+     ! note these are temporary also here so that we can
+     ! refactor without dependencies on array_sizes_common
+     real(kind=8),allocatable,dimension(:)::A_Isigxx, A_Idasig, A_Idbsig
+     real(kind=8),allocatable,dimension(:)::A_Isigsi, A_Idasis, A_Idbsis
      contains
+     procedure, pass(this) :: setUpDerivativeList => XctCrossCalc_setUpDerivativeList    ! set up  crossData, depending on number of isotopes
+     procedure, pass(this) :: setAddtionalParams => XctCrossCalc_setAddtionalParams
+     procedure, pass(this) :: setEnergyIndependent => XctCrossCalc_setEnergyIndependent  ! set energy independent values using current parameter values
+     procedure, pass(this) :: calcCross => XctCrossCalc_calcCross
      procedure, pass(this) :: initialize => XctCrossCalc_initialize
      procedure, pass(this) :: destroy => XctCrossCalc_destroy
   end type XctCrossCalc
 
 contains
+subroutine XctCrossCalc_setAddtionalParams(this,  lllmax, Kfinit, wantSelfIndicate, Kssmsc, debug)
+     use mxct25_m
+     class(XctCrossCalc)::this
+     integer::lllmax            ! maximum number of Clebsch-Gordon coefficients
+     integer::Kfinit            ! finite-size corrections flag
+     logical::wantSelfIndicate  ! do we need to calculate self-indicated cross section data
+     integer::Kssmsc
+     integer::numIso
+     logical::debug
+
+     integer::ntotTrig, ng
 
-subroutine XctCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero)
+    this%debug = debug
+
+
+     this%wantSelfIndicate = wantSelfIndicate
+     if ( wantSelfIndicate) then
+         numIso = this%crossData%getNumberIsotopes()
+         if (this%crossDataSelf%getNumberIsotopes().ne.numIso) then
+            call this%crossDataSelf%setUpList(this%resData, this%radiusData, numIso)
+         end if
+     end if
+
+     this%Kfinit  = Kfinit
+     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 (.not.allocated(this%Alj)) then
+               ng = this%resData%getNumSpinGroups()
+               ntotTrig = (this%ntotc*(this%ntotc+1))/2
+               allocate(this%Alj(ntotTrig, ng))
+           end if
+           CALL Kount_Clbsch (this%Zke, this%Alj, this%Kslow, this%resData, this%C_G_Kxlmn, this%Lllmax, this%ntotc)
+        END IF
+    end if
+end subroutine
+subroutine XctCrossCalc_setEnergyIndependent(this, reactType, Twomhb, kwcoul, Etac)
+    use mcro8_m
+    use mxct26_m
+    class(XctCrossCalc) :: this
+    integer::kwcoul, reactType
+    real(kind=8)::Twomhb, Etac
+
+    call CrossSectionCalculator_setEnergyIndependent(this, reactType, Twomhb, kwcoul, Etac)
+
+    ! calculate energy shift if needed
+    if (this%doResShift) then
+         call Fixx (this, this%XX,  this%XxHelper)
+    end if
+
+    ! determine whether we need coulomb interaction
+    call Find_If_Coulomb (this, this%IfCoul, this%Ifdif)
+end subroutine
+subroutine XctCrossCalc_calcCross(this, ener, Ipoten)
+    class(XctCrossCalc):: this
+    real(kind=8)::ener
+    integer::ipoten
+
+    integer::itot
+
+    call CrossSectionCalculator_calcCross(this, ener, Ipoten)
+    if (this%wantSelfIndicate) then
+       itot = this%covariance%getNumTotalParam() + 1
+       if ( .not.this%wantDerivs) itot = 1
+       call this%crossDataSelf%reserveColumns(this%row, itot)
+    end if
+end subroutine
+subroutine XctCrossCalc_setUpDerivativeList(this)
+   class(XctCrossCalc) :: this
+   integer::itot, i, ii, Nnnsig,numIso
+   type(SammySpinGroupInfo)::spinInfo
+
+   call CrossSectionCalculator_setUpDerivativeList(this)
+
+
+   ! initially set all shared isotope IDs to 1  
+   numIso = this%crossData%getNumberIsotopes()
+   if (numIso.gt.1) then
+      itot = this%covariance%getNumTotalParam() + 1
+      do i = 1, itot
+          ii = this%crossData%getIsotopeForShared(i)
+          if (ii.gt.0) then
+              call this%crossData%addSharedColumn(i, 1)
+              if (this%wantSelfIndicate) then
+                 call this%crossDataSelf%addSharedColumn(i, 1)
+              end if
+          end if
+      end do
+    end if   
+
+    if (this%wantSelfIndicate) then
+        call this%crossDataSelf%nullify()
+        call this%crossDataSelf%setNnsig(1)
+    end if
+end subroutine
+subroutine XctCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
     class(XctCrossCalc) :: this
     type(SammyRMatrixParameters)::pars
     type(CovarianceData)::cov
     type(AdjustedRadiusData)::rad
     integer,intent(in)::niso, Itzero, Ilzero
-    logical,intent(in)::needAngular
     type(VariedParameterInfo)::info
     integer::ntot
+    logical,intent(in)::needAngular, doShiftRes
 
-    call CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero)
+    integer::itot,  ipup, irad, nres
+    logical(C_BOOL)::countCombined
+
+    call CrossSectionCalculator_initialize(this, pars, cov, rad, niso, needAngular, Itzero, Ilzero, doShiftRes)
     call info%initialize(pars, cov, rad)
 
      ntot = info%getMaxChannel()
@@ -38,12 +174,42 @@ subroutine XctCrossCalc_initialize(this, pars, cov, rad, niso, needAngular, Itze
      this%Ifcros = .true.  ! assume cross sections are to be calculated for all channels
 
     call info%destroy()
+
+    call info%initialize(pars, cov, rad)
+    this%ntotc = info%getMaxChannel() -1  ! for most of the internal arrays the gamma channel is not used in the internal arrays
+    this%ntriag = (this%ntotc*(this%ntotc+1))/2
+    countCombined = .true.
+    call info%getNumAll(itot, ipup, countCombined)
+    irad = this%radiusData%getNumTotalVaried()
+    nres = this%resData%getNumResonances()
+    call info%destroy()
+
+    if (doShiftRes) then
+       call allocate_real_data(this%Xx, pars%getNumResonances())
+       call allocate_real_data(this%XxHelper, pars%getNumResonances())
+    end if
+
+    if (itot.gt.0) then
+       call reallocate_real_data_2d(this%Br, this%ntriag, 0, itot, 0)
+       call reallocate_real_data_2d(this%Bi, this%ntriag, 0, itot, 0)
+    end if
+
+    if (irad.gt.0) then
+       if( nres.gt.0) allocate(this%Bga(this%ntriag, this%ntotc, nres))
+    end if
 end subroutine
 subroutine XctCrossCalc_destroy(this)
     class(XctCrossCalc) :: this
 
     call CrossSectionCalculator_destroy(this)
     deallocate(this%Ifcros)
+
+    if (allocated(this%Alj))      deallocate(this%Alj)
+    if (allocated(this%Xx))       deallocate(this%Xx)
+    if (allocated(this%XxHelper)) deallocate(this%XxHelper)
+    if (allocated(this%Br))       deallocate(this%Br)
+    if (allocated(this%Bi))       deallocate(this%Bi)
+    if(allocated(this%Bga))       deallocate(this%Bga)
 end subroutine
 
 end module XctCrossCalc_M
diff --git a/sammy/src/xct/mxct0.f90 b/sammy/src/xct/mxct0.f90
index 40b762ef7..3f2e184e6 100644
--- a/sammy/src/xct/mxct0.f90
+++ b/sammy/src/xct/mxct0.f90
@@ -3,7 +3,7 @@ module xct_m
    contains
 ! --------------------------------------------------------------
 !
-      Subroutine Samxct_0
+      Subroutine Samxct_0(xct)
 !
       use oops_common_m
       use fixedi_m
@@ -22,16 +22,22 @@ module xct_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 SammyLptPrinting_m     
+      use XctCrossCalc_M
       use array_sizes_common_m, only : calcData
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      IMPLICIT none
       type(SammyGridAccess)::grid, auxGrid
-      real(kind=8)::XxTmp(1)
-      real(kind=8),allocatable,dimension(:)::A_Ialj
-      real(kind=8),allocatable,dimension(:)::A_Ix1
       character(len=80)::line
       integer::iflagMatch
-      logical::calcRadDeriv
+      integer::Idimen
+      integer::K_Coul_N, Kslow, 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
+      integer::Nnndrc
+      class(XctCrossCalc)::xct   ! temporarily here so that energy indepdent code can move in steps
+      external Idimen
 !
 !
           WRITE (line,99999)
@@ -72,8 +78,7 @@ module xct_m
 !
 !
 ! *** find max number of res for any one spin group
-      Mxany = 0
-      IF (resParData%getNumResonances().GT.0) CALL Findmx (Mxany)
+      Mxany = resParData%getNumResonances()
 !
 ! *** find Nrfil3, Npfil3 for added cross sections (endf/b-vi file 3)
 !     Now done in C++
@@ -81,13 +86,8 @@ module xct_m
       Nrfil3 = 1
 !
 ! *** Count how many non-zero elements are in Xlmn
-      IF (Kcros.EQ.7 .OR. Kcros.EQ.11 .OR. Kssmsc.GT.0) THEN
-         ng = resParData%getNumSpinGroups()
-         call allocate_real_data(A_Ialj, Ntriag*ng)
-         CALL Kount_Clbsch (A_Izke , A_Ialj, Kslow)
-      ELSE
-         Kkxlmn = 1
-      END IF
+      Kkxlmn = xct%C_G_Kxlmn
+      Kslow = xct%Kslow
 !
       IF (Kadddc.NE.0) THEN
 ! ***    Scan direct-capture file, figure dimensions et al
@@ -97,7 +97,7 @@ module xct_m
       END IF
 !
 ! *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-XCT
-      CALL Find_If_Coulomb (A_Izeta , IfCoul, Ifdif)
+      Ifcoul = xct%IfCoul
       CALL Estxct (Ntwo1, Nthr1, Nthr2, Nthr3, Nfour, Nfour1, Nfour2, &
          Nfive1, Nfive2, Nfive3, Nfive4, Nfive1x, Nfive3x, Nsix, Neight, &
          Nnine, Nw1, Mxany, Nfprrr, K_Coul_N, numElAux)
@@ -140,16 +140,7 @@ module xct_m
 !
 ! *** Two ***
 !
-      call allocate_real_data(A_Ixx, resParData%getNumResonances())
 ! - - - - - - - - - - - - - - - - <      
-      IF (resParData%getNumResonances().GT.0)  then
-         call allocate_real_data(A_Ix1, Mxany) ! was ntwo1, but Fixx assumes Mxany
-         CALL Fixx ( A_Ibound , A_Iechan , &
-                    A_Izke , A_Ixx, A_Ix1, Mxany)
-         deallocate(A_Ix1)
-      end if
-! *** SBROUTINE Fixx sets Xx = energy shift
-
 ! - - - - - - - - - - - - - - - - >
 !
 !
@@ -162,47 +153,22 @@ module xct_m
       call allocate_real_data(A_Ixlmn, Nthr3)
       call allocate_integer_data(I_Ikxlmn, Nthr3*3)
       IF (Kcros.EQ.7 .OR. Kcros.EQ.11 .OR. Kssmsc.GT.0) THEN
-         IF (Kslow.EQ.0) THEN
-            CALL Clbsch  ( A_Ixlmn, A_Ialj   , I_Ikxlmn)
+         IF (xct%Kslow.EQ.0) THEN
+            CALL Clbsch  ( resparData, A_Ixlmn, xct%Alj   , I_Ikxlmn)
          ELSE
-            CALL Clbsch_Slow (A_Ixlmn, A_Ialj   , I_Ikxlmn)
+            CALL Clbsch_Slow (resparData, A_Ixlmn, xct%Alj   , I_Ikxlmn)
          END IF
       END IF
-      if (allocated(A_Ialj)) deallocate(A_Ialj)
 !
 ! *** four ***
       N = Nfour
 !                     Nfour = (Ntriag) * (Nfpres)
       N   = Nfour1
-      call allocate_real_data(A_Ibga, N)
       N   = Nfour2
       call allocate_real_data(A_Ipgar, N)
       call allocate_real_data(A_Ipgai, N)
 !
-      IF (needResDerivs) THEN
-         CALL Babb (   A_Ipolar , I_Iflpol , &
-                       A_Ibr    , A_Ibi    , XxTmp, .true.)
-! ***    SBROUTINE Babb GENERATES ENERGY-INDEPENDENT PORTION OF
-! ***       PARTIAL DERIVATIVES
-      END IF
-!
-      calcRadDeriv = .false.
-      if (Ksolve.eq.2) then  ! only if we have pup'ed rad paramters
-         if (radFitFlags%getNumPupedVaried(covData).gt.0) then
-           calcRadDeriv = .true.
-         end if
-      else  ! otherwise if we have any varied or pup'ed rad parameters
-         if( radFitFlags%getNumTotalVaried().gt.0) then
-            calcRadDeriv = .true.
-         end if
-      end if
-      IF (calcRadDeriv) THEN
-         CALL Babbga ( A_Ibound , A_Iechan , &
-           A_Izke  ,    A_Izeta  , A_Ibga  )
-! ***    Sbroutine Babbga generates energy-independent portion of partial
-! ***      derivatives wrt unvaried widths (for use in figuring partial
-! ***      derivatives wrt varied channel radii)
-      END IF
+
 !
 ! *** five ***
       call allocate_real_data(A_Icrss, Nfive1)
@@ -328,7 +294,7 @@ module xct_m
 !
       Lllmmm = Lllmax
       IF (Lllmax.EQ.0) Lllmmm = 1
-      CALL Work (    calcData , calcDataSelf,                    &
+      CALL Work (    xct, calcData , calcDataSelf,               &
           A_Isigxx , A_Idasig , A_Idbsig , A_Isigsi , A_Idasis , &
           A_Idbsis , I_Iisopa ,                                  &
           A_Iedrcp , A_Icdrcp ,                                  &
@@ -365,14 +331,21 @@ module xct_m
 !
 ! *** GUESSTIMATE SIZE OF ARRAY NEEDED FOR SAMMY-XCT
 !
-      use fixedi_m
-      use ifwrit_m
-      use array_sizes_common_m
-      use broad_common_m
-      use lbro_common_m
-      use EndfData_common_m
-      use rsl7_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      use fixedi_m, only : Ifdif, Iq_Iso, Iq_Val, Kkxlmn, Lllmax, Nres, Mres, Ncrssx, Ncrsss,  &
+                           Ndasig, Ndbsig, Nfpres, Nnniso, Npfil3, Nrfil3, Ntotc, Ntriag,      &
+                           Kshift, needResDerivs
+      use ifwrit_m, only : Ifcoul, Kcros, Kpiece, Ksolve, Kssmsc, Nd_Xct,Nnpar
+      !use lbro_common_m
+      use EndfData_common_m, only : resParData, radFitFlags
+      use rsl7_m, only : Figure_Kws_Xct
+      IMPLICIT none
+      integer::Ntwo1, Nthr1, Nthr2, Nthr3, Nfour, Nfour1, &
+               Nfour2, Nfive1, Nfive2, Nfive3, Nfive4, Nfive1x, Nfive3x, &
+               Nsix, Neight, Nnine, Nw1, Mxany, Nfprrr, K_Coul_N, &
+               numElAux
+       integer::Idimen
+       integer::I,K, K1, K2, K3,K4, K6, K7, K8, Ke, Kw, low, N,Nthr4
+       external Idimen
 !
 ! *** Andy
       K = 2*(Npfil3+Nrfil3)
diff --git a/sammy/src/xct/mxct01.f90 b/sammy/src/xct/mxct01.f90
index 07d9d7300..363f3305e 100755
--- a/sammy/src/xct/mxct01.f90
+++ b/sammy/src/xct/mxct01.f90
@@ -5,72 +5,82 @@ module xct1_m
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Babb (Polar, Iflpol, Br, Bi, Xx, xct)
+      SUBROUTINE Babb (calc, xct)
 !
 ! *** Purpose -- Generate energy-independent portion of partial
 ! ***            derivatives of R with respect to U-parameters
-!
-      use fixedi_m
-      use ifwrit_m
-      use templc_common_m, only : I_Inotu
-      use EndfData_common_m
+!     
       use SammyResonanceInfo_M
       use RMatResonanceParam_M
+      use XctCrossCalc_M
       IMPLICIT None
 !
       type(SammyResonanceInfo)::resInfo
       type(RMatResonance)::resonance
       type(SammySpinGroupInfo)::spinInfo
+      class(XctCrossCalc)::calc
       logical::xct
       integer::Iuse
-      real(kind=8)::Br(:,:), Bi(:,:)
-      real(kind=8):: Polar(2,*), Xx(*)
-      integer::  Iflpol(2,*)
-      real(kind=8)::Zero, Two
       real(kind=8)::Aa, arg, Bb, D, D1, f1, f2
       integer::ichan, iFis1, iFis2, Iflr, Igam, Igr, Ipar, Index
       integer::Ires, J, Jk, K, Kj, Kk, M, Mmax, Mmax2, Mmm
-      real(kind=8)::W2, W3
-!
-!      DIMENSION Polar(2,Nres), Iflpol(2,Nres)
+      real(kind=8)::W2, W3,  polar2
+      integer::iresStart,ii, ntotc, iflpol2
 !
-      DATA Zero /0.0d0/, Two /2.0d0/
+      real(kind=8)::Zero=0.0d0, Two=2.0d0
 !
 !
-      Ipar = 0
-      Br = 0.0d0
-      Bi = 0.0d0
-      DO Ires=1,resParData%getNumResonances()
-         call resParData%getResonanceInfo(resInfo,  Ires)
-         Igr = resInfo%getSpinGroupIndex() 
-         IF (resInfo%getIncludeInCalc()) THEN
-            call resParData%getRedResonance(resonance, resInfo)
-            call resparData%getSpinGroupInfo(spinInfo, igr)
-            IF (.not.spinInfo%getIncludeInCalc()) cycle
+      Ipar = 0      
+
+      ! make loop over spin groups, so that we
+      ! can call getParamPerSpinGroup to ensure parameters
+      ! are counted consistently
+      iresStart = 0
+      ntotc = 0
+      DO igr = 1, calc%resData%getNumSpinGroups()
+         call calc%getParamPerSpinGroup(iresStart, igr)
+         call calc%resData%getSpinGroupInfo(spinInfo, igr)
+
+
+         do ii = 1, calc%inumSize
+            Ipar = Ipar + 1  ! count all
+
+            ! set it to zero per parameter instead of once in the beginning
+            ! is because we might have fitted parameter
+            ! but none of them are resonance parameters
+            ! In that case Br and Bi are not allocated
+            !
+            calc%Br(:,Ipar) = 0.0d0
+            calc%Bi(:,Ipar) = 0.0d0
+            Iflr = calc%Inum(ii,1)
+            IF (.not.calc%covariance%contributes(Iflr)) cycle ! but don't calculate all
+            M = calc%Inum(ii,2)
+            Ires = calc%Inum(ii,3)
+
+            call calc%resData%getResonanceInfo(resInfo,  Ires)
+            call calc%resData%getRedResonance(resonance, resInfo)
 
             Mmax = spinInfo%getNumChannels()
             Mmax2 = spinInfo%getNumResPar()
+            if (Mmax.gt.Ntotc) Ntotc = Mmax
             iGam = spinInfo%getGammaWidthIndex()
+            iflPol2 = 0
+            if (calc%Kpolar.NE.0) then
+               iflPol2 = resInfo%getChannelFitOption(4)
+               W2 = resonance%getWidth(3)
+               W3 = resonance%getWidth(4)
+               Polar2 = dSQRT(W2**2 + W3**2)
+            end if
+
 !
-            DO M=1,Mmax2
-               if(m.eq.1) then
-                  iflr = resInfo%getEnergyFitOption()
-               else
-                  iflr = resInfo%getChannelFitOption(m-1)
-               end if
-               IF (Iflr.GT.0) THEN
-                  Ipar = Ipar + 1  ! count all                  
-                  IF (covData%contributes(Iflr)) THEN  ! but don't calculate all
-                     if (Iflr.ne.abs(I_Inotu(Ipar))) then
-                         STOP 'Count of varied resonance parameter inconsistent in mxct01'
-                     end if
+
+
 !
                      IF (M.EQ.1) THEN
 ! ***                   HERE U-PARAMETER IS RESONANCE ENERGY
-                        if (xct) then
-                           arg = resonance%getEres()
-                        else
-                           arg = resonance%getEres() + Xx(Ires)
+                        arg = resonance%getEres()
+                        if (.not.xct.and.calc%doResShift) then
+                           arg =  arg + calc%Xx(Ires)
                         end if
                         D1 = Two*DSQRT(dABS(arg))
                         Kj = 0
@@ -82,8 +92,8 @@ module xct1_m
                               f2 = resonance%getWidth(ichan)
                               Kj = Kj + 1
                               D = f1*f2*D1
-                              Br(Kj,Ipar) =  D
-                              Bi(Kj,Ipar) = -D - D
+                              calc%Br(Kj,Ipar) =  D
+                              calc%Bi(Kj,Ipar) = -D - D
                            END DO
                         END DO
 !
@@ -100,8 +110,8 @@ module xct1_m
                               Kj = Kj + 1
                               D = f1*f2*D1
                               D = D + D
-                              Br(Kj,Ipar) = - D - D
-                              Bi(Kj,Ipar) =   D
+                              calc%Br(Kj,Ipar) = - D - D
+                              calc%Bi(Kj,Ipar) =   D
                            END DO
                         END DO
 !
@@ -116,8 +126,8 @@ module xct1_m
                                  ichan = spinInfo%getWidthForChannel(k)
                                  D = resonance%getWidth(ichan)
                                  IF (K.EQ.Mmm) D = D + D
-                                 Br(Kj,Ipar) = D
-                                 Bi(Kj,Ipar) = D
+                                 calc%Br(Kj,Ipar) = D
+                                 calc%Bi(Kj,Ipar) = D
                               END IF
                            END DO
                            IF (K.NE.Mmax) THEN
@@ -129,14 +139,14 @@ module xct1_m
                                       spinInfo%getWidthForChannel(k)
                                     D = resonance%getWidth(ichan)
                                     IF (K.EQ.Mmm) D = D + D
-                                    Br(Jk,Ipar) = D
-                                    Bi(Jk,Ipar) = D
+                                    calc%Br(Jk,Ipar) = D
+                                    calc%Bi(Jk,Ipar) = D
                                  END IF
                               END DO
                            END IF
                         END DO
 !
-                        IF (Kpolar.NE.0 .AND. Mmm.EQ.3) THEN   ! >>> mmm.eq.3 (xct) mmm.gt.1 otherwise
+                        IF (calc%Kpolar.NE.0 .AND. Mmm.EQ.3) THEN   ! >>> mmm.eq.3 (xct) mmm.gt.1 otherwise
 ! ***                      Here we use polar coordinants for fission
 ! ***
                            Kj = 0
@@ -147,52 +157,52 @@ module xct1_m
                            DO K=1,Mmax
                               DO J=1,K
                                  Kj = Kj + 1
-                                 Aa = Br(Kj,Ipar-1)
-                                 Bb = Br(Kj,Ipar)
+                                 Aa = calc%Br(Kj,Ipar-1)
+                                 Bb = calc%Br(Kj,Ipar)
                                  IF (Aa.NE.Zero .OR. Bb.NE.Zero) THEN
                                     W2 = resonance%getWidth(iFis1)
                                     W3 = resonance%getWidth(iFis2)
-                                    Br(KJ,Ipar-1) = -W3* Aa +W2*Bb
-                                    IF (Iflpol(2,Ires).EQ.0) THEN
-                                       Br(Kj,Ipar) = Zero
+                                    calc%Br(KJ,Ipar-1) = -W3* Aa +W2*Bb
+                                    IF (Iflpol2.EQ.0) THEN
+                                       calc%Br(Kj,Ipar) = Zero
                                     ELSE
-                                       Br(Kj,Ipar) = ( W2*Aa + w3*Bb ) &
-                                          / Polar(2,Ires)
+                                       calc%Br(Kj,Ipar) = ( W2*Aa + w3*Bb ) &
+                                          / Polar2
                                     END IF
                                  END IF
-                                 Aa = Bi(Kj,Ipar-1)
-                                 Bb = Bi(Kj,Ipar)
+                                 Aa = calc%Bi(Kj,Ipar-1)
+                                 Bb = calc%Bi(Kj,Ipar)
                                  IF (Aa.NE.Zero .OR. Bb.NE.Zero) THEN
                                     W2 = resonance%getWidth(iFis1)
                                     W3 = resonance%getWidth(iFis2)
-                                    Bi(Kj,Ipar-1) = -W3*Aa + W2*Bb
-                                    IF (Iflpol(2,Ires).EQ.0) THEN
-                                       Bi(Kj,Ipar) = Zero
+                                    calc%Bi(Kj,Ipar-1) = -W3*Aa + W2*Bb
+                                    IF (Iflpol2.EQ.0) THEN
+                                       calc%Bi(Kj,Ipar) = Zero
                                     ELSE
-                                       Bi(Kj,Ipar) = ( W2* Aa + W3*Bb ) &
-                                          / Polar(2,Ires)
+                                       calc%Bi(Kj,Ipar) = ( W2* Aa + W3*Bb ) &
+                                          / Polar2
                                     END IF
                                  END IF
                               END DO
                            END DO
                         END IF
                      END IF
-                  END IF
-               END IF
-            END DO
-         END IF
-      END DO
-!
+
+
+
+         END do  ! end of loop over  parameters
+      END DO   ! end of loop over spin groups
+
       if (xct) then
       DO m=1,Ipar
          Kj = 0
          DO K=1,Ntotc
             DO J=1,K
                Kj = Kj + 1
-               IF (Br(Kj,M).NE.Zero) Br(Kj,M) = &
-                    Two*Br(Kj,m)
-               IF (Bi(Kj,M).NE.Zero) Bi(Kj,M) = &
-                    Two*Bi(Kj,M)
+               IF (calc%Br(Kj,M).NE.Zero) calc%Br(Kj,M) = &
+                    Two*calc%Br(Kj,m)
+               IF (calc%Bi(Kj,M).NE.Zero) calc%Bi(Kj,M) = &
+                    Two*calc%Bi(Kj,M)
             END DO
          END DO
       END DO
@@ -204,19 +214,19 @@ module xct1_m
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Babbga (Bound, Echan, Zke, Zeta, Bga)
+      SUBROUTINE Babbga (calc, Kwcoul)
 !
 ! *** PURPOSE -- GENERATE ENERGY-INDEPENDENT PORTION OF PARTIAL of R wrt
 ! ***            reduced-width-amplitudes (for those which are not varied)
 !
       use sammy_CoulombSelector_I
-      use fixedi_m
-      use ifwrit_m
-      use EndfData_common_m
       use SammyResonanceInfo_M
       use RMatResonanceParam_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      use XctCrossCalc_M
+      IMPLICIT none
 !
+      class(XctCrossCalc)::calc
+      integer::Kwcoul
       type(SammyResonanceInfo)::resInfo
       type(RMatResonance)::resonance
       type(SammySpinGroupInfo)::spinInfo
@@ -224,20 +234,27 @@ module xct1_m
       type(RMatChannelParams)::channel
       type(SammyParticlePairInfo)::pairInfo
       type(RMatParticlePair)::pair
-      DIMENSION &
-         Bound(Ntotc,*), Echan(Ntotc,*), Zke(Ntotc,*),  &
-         Zeta(Ntotc,*), Bga(Ntriag,Ntotc,*)
+
 !
-      DATA Zero /0.0d0/
+      real(kind=8),parameter:: Zero = 0.0d0
+      real(kind=8)::Cosphi, D, Dgda, Dp, Dphi, Ds, eres, Rhox
+      real(kind=8)::Eta, Hi, Hr, P, Q, Rho, Sinphi, su, width
+      integer::ichan, Iffy, igr, Ires, Ishift, J, Jk, Kj
+      integer::Iflr, K, Mmm, Kk, Lsp, M, Mmax, nres
 !
-      CALL Zero_Array (Bga, Ntriag*Ntotc*Nres)
+      if (calc%radiusData%getNumTotalVaried().eq.0) return
+
+      nres = calc%resData%getNumResonances()
+      if(nres.eq.0) return
+
+      calc%Bga = Zero
 !
-      DO Ires=1,resParData%getNumResonances()
-         call resParData%getResonanceInfo(resInfo,  Ires)
+      DO Ires=1,calc%resData%getNumResonances()
+         call calc%resData%getResonanceInfo(resInfo,  Ires)
          IF (resInfo%getIncludeInCalc()) THEN
-            call  resParData%getRedResonance(resonance, resInfo)
+            call  calc%resData%getRedResonance(resonance, resInfo)
             igr = resInfo%getSpinGroupIndex()
-            call resParData%getSpinGroupInfo(spinInfo, igr)
+            call calc%resData%getSpinGroupInfo(spinInfo, igr)
             Mmax = spinInfo%getNumChannels()
             eres = resonance%getEres()
             Su = dABS(eres)
@@ -250,11 +267,11 @@ module xct1_m
                end if
                IF (Iflr.LE.0) THEN
                   call spinInfo%getChannelInfo(channelInfo, Mmm)
-                  call resParData%getChannel(channel, channelInfo)
-                  call resParData%getParticlePairInfo( &
+                  call calc%resData%getChannel(channel, channelInfo)
+                  call calc%resData%getParticlePairInfo( &
                      pairInfo, &
                      channelInfo%getParticlePairIndex())
-                 call resParData%getParticlePair(pair, pairInfo)
+                 call calc%resData%getParticlePair(pair, pairInfo)
                  if (pair%getCalcShift()) then
                      Ishift = 1
                  else
@@ -265,16 +282,16 @@ module xct1_m
                   width = resonance%getWidth(ichan)
                   IF (pair%getPnt().EQ.1 .AND. width.NE.Zero) THEN
                      Lsp = channel%getL()
-                     Q = dSQRT( dABS(eres-Echan(Mmm,Igr)) )
-                     Rho  = channel%getApt()*Zke (Mmm,Igr)*Q
-                     Rhox = Zke (Mmm,Igr)*Q
-                     IF (Zeta(Mmm,Igr).EQ.Zero) THEN
-                        CALL Pgh (Rho, Lsp, Bound(Mmm,Igr), Hr, Hi, P, &
+                     Q = dSQRT( dABS(eres-calc%Echan(Mmm,Igr)) )
+                     Rho  = channel%getApt()*calc%Zke (Mmm,Igr)*Q
+                     Rhox = calc%Zke (Mmm,Igr)*Q
+                     IF (calc%Zeta(Mmm,Igr).EQ.Zero) THEN
+                        CALL Pgh (Rho, Lsp, calc%Bound(Mmm,Igr), Hr, Hi, P, &
                            Dp, Ds, Ishift, Iffy)
                      ELSE
-                        Eta  = Zeta(Mmm,Igr)/Q
+                        Eta  = calc%Zeta(Mmm,Igr)/Q
                         CALL f_sammy_columb_Pghcou(Kwcoul, &
-                           Rho, Lsp, Bound(Mmm,Igr), &
+                           Rho, Lsp, calc%Bound(Mmm,Igr), &
                            Hr, Hi, &
                            P, Dp, Ds, Ishift, Iffy, Eta, &
                            Sinphi, Cosphi, Dphi, 0)
@@ -294,7 +311,7 @@ module xct1_m
                               ichan = spinInfo%getWidthForChannel(k)
                               D = resonance%getWidth(ichan)
                               IF (K.EQ.Mmm) D = D + D
-                              Bga(Kj,Mmm,Ires) = D
+                              calc%Bga(Kj,Mmm,Ires) = D
                            END IF
                         END DO
                         IF (K.LE.Mmax) THEN
@@ -305,14 +322,14 @@ module xct1_m
                               ichan = spinInfo%getWidthForChannel(k)
                               D = resonance%getWidth(ichan)
                               IF (K.EQ.Mmm) D = D + D
-                              Bga(Jk,Mmm,Ires) = D
+                              calc%Bga(Jk,Mmm,Ires) = D
                               END IF
                            END DO
                         END IF
                      END DO
-                     DO Jk=1,Ntriag
-                        IF (Bga(Jk,Mmm,Ires).NE.Zero) Bga(Jk,Mmm,Ires) = &
-                                                   Bga(Jk,Mmm,Ires)*Dgda
+                     DO Jk=1,calc%ntriag
+                        IF (calc%Bga(Jk,Mmm,Ires).NE.Zero) calc%Bga(Jk,Mmm,Ires) = &
+                                                   calc%Bga(Jk,Mmm,Ires)*Dgda
                      END DO
                   END IF
                END IF
diff --git a/sammy/src/xct/mxct02.f90 b/sammy/src/xct/mxct02.f90
index 3bd0e9c25..8d94997bb 100644
--- a/sammy/src/xct/mxct02.f90
+++ b/sammy/src/xct/mxct02.f90
@@ -4,7 +4,7 @@ module xct2_m
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Work (derivs,   derivsSelf,                      &
+      SUBROUTINE Work (calc, derivs,   derivsSelf,                &
          Sigxxx, Dasigx, Dbsigx, Sigsin, Dasigs,                  &
          Dbsigs, Isopar,                                          &
          Edrcpt, Cdrcpt, Xdrcpt, Ndrcpt, Nnndrc, Lllmmm, Kslow)
@@ -12,7 +12,6 @@ module xct2_m
 ! *** PURPOSE -- Generate theoretical cross sections "theory" and partial
 ! ***            Derivatives "dasig"
 !
-      use oops_common_m
       use fixedi_m
       use ifwrit_m
       use exploc_common_m
@@ -32,6 +31,10 @@ module xct2_m
       use convert_to_transmission_m
       ! use xct6_m
       use mxct27_m
+      use mxct26_m
+      use mxct06_m
+      use mxct18_m
+      use XctCrossCalc_M
       IMPLICIT none
 
       real(8), intent(in)::                       &
@@ -51,7 +54,8 @@ module xct2_m
       type(SammyGridAccess)::grid
       type(DerivativeHandler)::derivs, derivsSelf
       integer::iflagMatch
-      logical::wantNegative, wantTrans, wantDeriv
+      logical::wantNegative, wantDeriv
+      class(XctCrossCalc)::calc
 
       DIMENSION   &
          Edrcpt(Nnndrc,*), Cdrcpt(Nnndrc,*), Xdrcpt(*), Ndrcpt(*)
@@ -76,18 +80,6 @@ module xct2_m
       wantNegative = .true.
       if (Kkkdop.eq.1) wantNegative = .false.
 
-      ! do we need to convert to transmission
-      ! will only be done if no further doppler broadening is needed
-      wantTrans = .false.
-      IF (Ytrans .AND. Kcros.LT.7 .AND. Ktruet.EQ.0) wantTrans = .true.
-      if (Ydoppr) then
-          ! don't do transmission conversion
-          ! except if Maxwellian averages or energy-averages with no
-          ! correction terms
-          if (.not.(Maxwel.EQ.1 .OR. Knocor.EQ.1)) then
-             wantTrans = .false.
-          end if
-     end if
 
 
       Ywhich = Ydoppr.OR.Yresol.OR.Yangle.OR.Yssmsc.OR.Yaverg.OR. Maxwel.EQ.1 .OR. Knocor.EQ.1
@@ -104,27 +96,11 @@ module xct2_m
       call derivs%nullify()
       call derivs%setNnsig(Nnnsig)
 
-      if (Iq_Iso.gt.1) then
-         TotalNdasig = Nfpres + Nfpext + radFitFlags%getNumTotalVaried() + Nfpiso
-         iflagMatch = radFitFlags%matchFitFlag()
-         if(iflagMatch.gt.0) TotalNdasig = TotalNdasig + 1
-
-          DO Iipar=1,TotalNdasig
-            call derivs%addSharedColumn(Iipar, 1)
-          end do
-      end if
-      call  derivs%setUpList(resParData, radFitFlags, Iq_Iso)
       call derivs%reserve(numEl * Nnnsig, Ndasig + Ndbsig + 1)
 
       IF (Iw.EQ.1.or.Ksitmp.gt.0) THEN
          call derivsSelf%nullify()
          call derivsSelf%setNnsig(1)
-         if (Iq_Iso.gt.1) then
-            DO Iipar=1,TotalNdasig
-                call derivsSelf%addSharedColumn(Iipar, 1)
-            end do
-         end if
-         call  derivsSelf%setUpList(resParData, radFitFlags, Iq_Iso)
          call derivsSelf%reserve(numEl, Ndasig + Ndbsig + 1)
       end if
 
@@ -220,9 +196,9 @@ module xct2_m
 !
 ! ************ Generate cross sections and derivatives
                IF (Nd_Xct.NE.0 .AND. Ksolve.NE.2) THEN
-                  CALL N_D_Zcross (Kount_Helmut)
+                  CALL N_D_Zcross (calc, Kount_Helmut)
                ELSE
-                  CALL Zcross (Nnndrc, Ipoten, Kount_Helmut)
+                  CALL Zcross (calc, Nnndrc, Ipoten, Kount_Helmut)
                END IF
 !
 ! ************ Store Coul if needed
@@ -232,7 +208,7 @@ module xct2_m
                END IF
 !
 ! ************    Set the particular cross sections needed for this run
-               CALL Zwhich (Sigxxx, Dasigx, Dbsigx, Sigsin, &
+               CALL Zwhich (calc, Sigxxx, Dasigx, Dbsigx, Sigsin, &
                   Dasigs, Dbsigs, Theoryx, Su,              &
                   grid%getEnergy(Jdat, expData), Lllmmm, Kslow)
                IF (Kfake.EQ.1) THEN
diff --git a/sammy/src/xct/mxct03.f90 b/sammy/src/xct/mxct03.f90
index c712c626d..58f49d869 100644
--- a/sammy/src/xct/mxct03.f90
+++ b/sammy/src/xct/mxct03.f90
@@ -4,7 +4,7 @@ module xct3_m
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE N_D_Zcross (Kount_Helmut)
+      SUBROUTINE N_D_Zcross (calc, Kount_Helmut)
 !
 ! *** PURPOSE -- Calculate numerically the partial derivatives 
 ! ***               of the cross section wrt R-matrix parameters
@@ -19,13 +19,17 @@ module xct3_m
       use EndfData_common_m
       use SammyResonanceInfo_M
       use RMatResonanceParam_M
+      use xct1_m
+      use xct5_m
+      use mxct06_m
+      use XctCrossCalc_M
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 !
 !
       type(SammyResonanceInfo)::resInfo
       type(RMatResonance)::resonance
       type(SammySpinGroupInfo)::spinInfo
-      real(kind=8)::XxTmp(1)
+      class(XctCrossCalc)::calc
 !
       DATA Zero /0.0d0/, One /1.0d0/
       DATA U_Increment /0.0001d0/
@@ -44,13 +48,11 @@ module xct3_m
 ! *** First, for the original parameter values --
 !     *** Generate energy-independent pieces
 !   True is passed to babb since it is used to set parameters for numerical differentiation      
-      CALL Babb (   &
-        A_Ipolar , I_Iflpol , A_Ibr, &
-        A_Ibi , Xxtmp, .true.)
-       CALL Abpart (  &
-       A_Ialphr , A_Ialphi , A_Ibr    , &
-       A_Ibi    , A_Ipr    , A_Ipi    , A_Idifen , A_Ixden  , &
-       A_Idifma , I_Inot   , I_Inotu  , A_Ixx    , &
+     CALL Babb (   calc, .true.)
+     CALL Abpart (  calc, &
+       A_Ialphr , A_Ialphi ,  &
+       A_Ipr    , A_Ipi    , A_Idifen , A_Ixden  , &
+       A_Idifma , I_Inot   , I_Inotu  ,                       &
        A_Iprer  , A_Iprei )
 !
 ! *** Form the cross section Crss
@@ -134,14 +136,12 @@ module xct3_m
                       END IF
 !
 ! ***              Generate energy-independent pieces with new parameter
-                   CALL Babb ( A_Ipolar , &
-                               I_Iflpol , A_Ibr    , A_Ibi, &
-                               Xxtmp, .true.)
-                   CALL Abpart (A_Ialphr , &
-                     A_Ialphi , A_Ibr    , A_Ibi    , A_Ipr    , &
+                   CALL Babb ( calc,  .true.)
+                   CALL Abpart (calc, A_Ialphr , &
+                     A_Ialphi , A_Ipr    , &
                      A_Ipi    , A_Idifen , A_Ixden  , &
                      A_Idifma , I_Inot   , &
-                     I_Inotu  , A_Ixx    , A_Iprer  , A_Iprei  )
+                     I_Inotu  , A_Iprer  , A_Iprei  )
 !
 ! ***              Form the cross section Crss with new parameter value
                    CALL Crosss (   &
diff --git a/sammy/src/xct/mxct04.f90 b/sammy/src/xct/mxct04.f90
index 79388415c..5bdbf142d 100644
--- a/sammy/src/xct/mxct04.f90
+++ b/sammy/src/xct/mxct04.f90
@@ -4,7 +4,7 @@ module xct4_m
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Zcross (Nnndrc, Ipoten, Kount_Helmut)
+      SUBROUTINE Zcross (calc, Nnndrc, Ipoten, Kount_Helmut)
 !
 ! *** PURPOSE -- FORM THE CROSS SECTION Crss
 ! ***            AND THE ( PARTIAL DERIVATIVES OF THE CROSS SECTION
@@ -19,21 +19,23 @@ module xct4_m
       use templc_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
+      class(XctCrossCalc)::calc
 !
 !
 ! *** Generate Pr and Pi = Partial of R wrt U-parameters
 ! ***    from Upr and Upi = energy-dependent pieces of those derivs    
       IF (resParData%getNumResonances().NE.0)  then
-         CALL Abpart (   &
-         A_Ialphr , A_Ialphi , A_Ibr    , &
-         A_Ibi    , A_Ipr    , A_Ipi    , A_Idifen , A_Ixden  , &
-         A_Idifma , I_Inot   , I_Inotu  , A_Ixx  , &
-         A_Iprer  , A_Iprei )
+         CALL Abpart (  calc,  &
+         A_Ialphr , A_Ialphi , &
+         A_Ipr    , A_Ipi    , A_Idifen , A_Ixden  , &
+         A_Idifma , I_Inot   , I_Inotu  , A_Iprer  , A_Iprei )
          end if
 !
 ! *** Generate Pgar & Pgai = partial of R wrt (Gamma-x) *
@@ -68,7 +70,7 @@ module xct4_m
 
       IF (doRadDeriv) THEN
           CALL Abpga (   &
-             A_Ialphr, A_Ialphi, A_Ibga , A_Ipgar, A_Ipgai, Nfprrr)
+             A_Ialphr, A_Ialphi, calc%Bga , A_Ipgar, A_Ipgai, Nfprrr)
       END IF
 !
 ! *** FORM THE CROSS SECTION Crss AND THE ( PARTIAL DERIVATIVES OF THE
diff --git a/sammy/src/xct/mxct05.f90 b/sammy/src/xct/mxct05.f90
index d19361a28..e23ea46bc 100644
--- a/sammy/src/xct/mxct05.f90
+++ b/sammy/src/xct/mxct05.f90
@@ -4,9 +4,9 @@ module xct5_m
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Abpart ( &
-         Alphar, Alphai, Br, Bi, Pr, Pi, Difen, Xden, &
-         Difmax, Not, Notu, Xx, Prer, Prei)
+      SUBROUTINE Abpart (calc, &
+         Alphar, Alphai, Pr, Pi, Difen, Xden, &
+         Difmax, Not, Notu, Prer, Prei)
 !
 ! *** Purpose -- Generate Alphar & Alphai = energy-independent bits
 ! ***            and Upr and Upi = Energy-dependent pieces of Pr & Pi
@@ -21,14 +21,15 @@ module xct5_m
       use RMatResonanceParam_M
       use SammySpinGroupInfo_M
       use templc_common_m
+      use XctCrossCalc_M
       IMPLICIT None
 !
 !
-      real(kind=8)::Bi(:,:), Br(:,:), Pr(:, :), Pi(:,:)
+      real(kind=8):: Pr(:, :), Pi(:,:)
       real(kind=8):: &
          Alphar(*), Alphai(*), &
          Difen(*), Xden(*), Difmax(*), &
-         Xx(*), Prer(Ntriag,*), Prei(Ntriag,*)
+         Prer(Ntriag,*), Prei(Ntriag,*)
       integer:: Not(*), Notu(*)
       type(SammyResonanceInfo)::resInfo
       type(SammySpinGroupInfo)::spinInfo
@@ -36,12 +37,13 @@ module xct5_m
       real(kind=8)::Zero, One, Two
       real(kind=8)::Aa, Bb, channelWidthC, channelWidthCPrime, G2, G3
       integer::I, ichan, Ig, igrp, Ij, Ijk, Ijl, Iflr, Ipar, J, M, N, N2, K
+      class(XctCrossCalc)::calc
 !     
 !      DIMENSION 
 !     *     Alphar(Nres), Alphai(Nres),
 !     *     Difen(Nres), Xden(Nres),
 !     *     Difmax(Nres),
-!     *     Not(Nres) , Xx(Nres), Prer(Ntriag,Ngroup),
+!     *     Not(Nres) , Prer(Ntriag,Ngroup),
 !     *     Prei(Ntriag,Ngroup)
 !
       DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
@@ -62,7 +64,10 @@ module xct5_m
             call resParData%getSpinGroupInfo(spinInfo, igrp)
             IF (spinInfo%getIncludeInCalc()) THEN
                call resParData%getRedResonance(resonance, resInfo)
-               Difen(N) = resonance%getEres() - Su + Xx(N)
+               Difen(N) = resonance%getEres() - Su
+               if (calc%doResShift) then
+                  Difen(N) =  Difen(N) + calc%Xx(N)
+               end if
                IF (Dabs(Difen(N)).LT.1.0D8*Difmax(N)) Not(N) = 0
 !cccccccccccccccccccccccccccccccccc
 !   not(n)=0
@@ -153,11 +158,11 @@ module xct5_m
                DO I=1,Ntotc
                   DO J=1,I
                      IJ = IJ + 1
-                     IF (Br(IJ,K).NE.Zero) THEN
-                        Pr(IJ,K) = Br(IJ,K)*Upr(K)
+                     IF (calc%Br(IJ,K).NE.Zero) THEN
+                        Pr(IJ,K) = calc%Br(IJ,K)*Upr(K)
                      END IF
-                     IF (Bi(IJ,K).NE.Zero) THEN
-                        Pi(IJ,K) = Bi(IJ,K)*Upi(K)
+                     IF (calc%Bi(IJ,K).NE.Zero) THEN
+                        Pi(IJ,K) = calc%Bi(IJ,K)*Upi(K)
                      END IF                     
                   END DO
                END DO
diff --git a/sammy/src/xct/mxct06.f90 b/sammy/src/xct/mxct06.f90
index 7755f19fa..7bba9bbd7 100644
--- a/sammy/src/xct/mxct06.f90
+++ b/sammy/src/xct/mxct06.f90
@@ -1,3 +1,5 @@
+module mxct06_m
+contains
 !
 !
 ! --------------------------------------------------------------
@@ -431,3 +433,4 @@
 !
       RETURN
       END
+end module mxct06_m
diff --git a/sammy/src/xct/mxct09.f90 b/sammy/src/xct/mxct09.f90
index e86bdf447..b070dc007 100644
--- a/sammy/src/xct/mxct09.f90
+++ b/sammy/src/xct/mxct09.f90
@@ -75,6 +75,7 @@
       use EndfData_common_m
       use SammyResonanceInfo_M
       use RMatResonanceParam_M
+      use mxct26_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 !
       type(SammySpinGroupInfo)::spinInfo
diff --git a/sammy/src/xct/mxct10.f90 b/sammy/src/xct/mxct10.f90
index 4839eb45d..77e12f92b 100644
--- a/sammy/src/xct/mxct10.f90
+++ b/sammy/src/xct/mxct10.f90
@@ -91,6 +91,7 @@
       use EndfData_common_m
       use SammyResonanceInfo_M
       use RMatResonanceParam_M
+      use mxct26_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 !
       type(SammySpinGroupInfo)::spinInfo
diff --git a/sammy/src/xct/mxct18.f90 b/sammy/src/xct/mxct18.f90
index aee7ee7c4..a6b45a274 100644
--- a/sammy/src/xct/mxct18.f90
+++ b/sammy/src/xct/mxct18.f90
@@ -1,8 +1,11 @@
+module mxct18_m
+use XctCrossCalc_M
+contains
 !
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Zwhich (Sigxxx, Dasigx, Dbsigx, Sigsin, Dasigs,   &
+      SUBROUTINE Zwhich (calc, Sigxxx, Dasigx, Dbsigx, Sigsin, Dasigs,   &
          Dbsigs, Theory, Su, Eb, Lllmmm, Kslow)
 !
 ! *** Purpose -- Set the particular cross sections needed for this run
@@ -12,7 +15,9 @@
       use ifwrit_m
       use exploc_common_m
       use templc_common_m
+      use mxct20_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      class(XctCrossCalc)::calc
       real(kind=8),pointer,dimension(:)::A_Ietax
 !
       DIMENSION Sigxxx(Nnnsig,*), Dasigx(Nnnsig,*),   &
@@ -50,7 +55,7 @@
 ! ***    Here where there is scattering, self-shielding, angular
 ! ***       distributions, or something  involving more than one
 ! ***       type of cross section
-         CALL Diffel ( A_Isiabn ,    &
+         CALL Diffel ( calc, A_Isiabn ,    &
             I_Ifexcl ,   &
             A_Iprdet , I_Ifldet , I_Iigrde , I_Iflmsc , A_Iprpmc ,   &
             I_Iflpmc , I_Isopmc ,   &
@@ -63,3 +68,4 @@
       END IF
       RETURN
       END
+end module mxct18_m
diff --git a/sammy/src/xct/mxct20.f90 b/sammy/src/xct/mxct20.f90
index 9d471141e..83c8bf319 100644
--- a/sammy/src/xct/mxct20.f90
+++ b/sammy/src/xct/mxct20.f90
@@ -1,8 +1,11 @@
+module mxct20_m
+use XctCrossCalc_M
+contains
 !
 !
 ! ______________________________________________________________________
 !
-      SUBROUTINE Diffel (Siabnd, Jfexcl,    &
+      SUBROUTINE Diffel (calc, Siabnd, Jfexcl,    &
                  Pardet, Ifldet, Igrdet, Iflmsc, Parpmc,   &
          Iflpmc, Isopmc, Sigxxx, Dasigx, Dbsigx, Sigsin, Dasigs, Dbsigs,   &
          Isopar, Ccclll, Dddlll, Xlmn  , Kxlmn , Crss  , Deriv , Crssx ,   &
@@ -17,8 +20,13 @@
       use fixedi_m
       use ifwrit_m
       use lbro_common_m
+      use mxct32_m
+      use mxct31_m
+      use mxct22_m
+      use mxct21_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 !
+      class(XctCrossCalc)::calc
       DIMENSION Siabnd(*), Jfexcl(Ntotc,*),   &
          Pardet(*), Ifldet(*), Igrdet(*), Iflmsc(*),   &
          Parpmc(*), Iflpmc(*), Isopmc(*),   &
@@ -42,10 +50,10 @@
 !
 ! ****** first, set the Legendre coefficients:
          IF (Kslow.EQ.0) THEN
-            CALL Setleg ( Sigxxx, Ccclll,   &
+            CALL Setleg ( calc, Sigxxx, Ccclll,   &
                Xlmn, Kxlmn, Crssx, Echan, Cmlab, Iso_Qv, Lllmmm, Eb)
          ELSE
-           CALL Setleg_Slow (Sigxxx, Ccclll,   &
+           CALL Setleg_Slow (calc, Sigxxx, Ccclll,   &
                Xlmn, Kxlmn, Crssx, Echan, Cmlab, Iso_Qv, Lllmmm, Eb)
          END IF
 !
@@ -64,12 +72,12 @@
 !
 ! ****** set the cross sections (non-angle-differential):
 ! ****** beginning with self-indication transmission
-         IF (Ksindi.GT.0 .AND. Kcros.EQ.8) CALL Setsel (Siabnd,    &
+         IF (Ksindi.GT.0 .AND. Kcros.EQ.8) CALL Setsel (calc, Siabnd,    &
             Parpmc, Isopmc, Sigsin, Crss, Su, Eb)
 !
 !
 ! ****** now do other terms
-         CALL Setoth (Pardet, Igrdet,   &
+         CALL Setoth (calc, Pardet, Igrdet,   &
             Parpmc, Isopmc, Sigxxx, Crss, Termf, Termfx, Su, Eb)
       END IF
 !
@@ -83,11 +91,11 @@
 ! ****** first, derivatives of Legendre coefficients
          IF (Ndasig.GT.0) THEN
             IF (Kslow.EQ.0) THEN
-               CALL Derleg (   &
+               CALL Derleg (  calc,  &
                   Sigxxx, Dasigx, Ccclll, Dddlll, Xlmn, Kxlmn, Crssx,   &
                   Derivx, Isopar, Echan, Cmlab, Iso_Qv, Lllmmm, Eb)
             ELSE
-               CALL Derleg_Slow (   &
+               CALL Derleg_Slow (  calc,  &
                   Sigxxx, Dasigx, Ccclll, Dddlll, Xlmn, Kxlmn,   &
                   Crssx, Derivx, Isopar, Echan, Cmlab, Iso_Qv, Lllmmm,   &
                   Eb)
@@ -115,3 +123,5 @@
    40 CONTINUE
       RETURN
       END
+
+end module mxct20_m
diff --git a/sammy/src/xct/mxct21.f90 b/sammy/src/xct/mxct21.f90
index 93b345e20..4f4bcd4db 100644
--- a/sammy/src/xct/mxct21.f90
+++ b/sammy/src/xct/mxct21.f90
@@ -1,30 +1,52 @@
+module mxct21_m
+use XctCrossCalc_M
+IMPLICIT NONE
+
+private Jxnnn, Jxmmm
+
+contains
+    integer function Jxnnn (Nb,Na,N, Ntotc, Ngroup)
+       integer Nb,Na,N, Ntotc, Ngroup
+
+       Jxnnn = (((N-1)*Ntotc+Na-1)*Ntotc+Nb-1)*Ngroup - 1
+    end function
+    integer function Jxmmm (Mb,Ma,M,J, Ntotc, Lllmax)
+       integer Mb,Ma,M,J, Ntotc, Lllmax
+
+       Jxmmm   = (((J+M)*Ntotc+Ma-1)*Ntotc+Mb-1)*Lllmax
+    end function
 !
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Setleg (Sigxxx, Ccclll,   &
+      SUBROUTINE Setleg (calc, Sigxxx, Ccclll,   &
          Xlmn, Kxlmn, Crssx, Echan, Cmlab, Iso_Qv, Lllmmm, Eb)
 !
 ! *** Purpose -- set Ccclll(L,Iso) = coefficient of Legendre polynomial
 ! ***                         P-sub-(L-1) for Isotope Iso
 !
-      use fixedi_m
-      use ifwrit_m
-      use lbro_common_m
-      use EndfData_common_m
+      use fixedi_m, only : Nnnsig, Kkxlmn, Ntotc, Iq_Iso, Iq_Val, Lllmax, Nnniso, Numiso
       use SammySpinGroupInfo_M
       use mdat9_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
 !
-      DIMENSION    &
-         Sigxxx(Nnnsig,*), Ccclll(Lllmmm,*), Xlmn(*), Kxlmn(Kkxlmn),   &
-         Crssx(2,Ntotc,Ntotc,*), Echan(Ntotc,*), Cmlab(3,*), Iso_Qv(*)
+      type(XctCrossCalc)::calc
+      real(kind=8)::    &
+         Sigxxx(Nnnsig,*), Ccclll(Lllmmm,*), Xlmn(*),   &
+         Crssx(2,Ntotc,Ntotc,*), Echan(Ntotc,*), Cmlab(3,*)
+      integer::   Kxlmn(Kkxlmn),  Iso_Qv(*)
       type(SammySpinGroupInfo)::spinMgr, spinNgr
-      DATA Zero /0.0d0/
-      Jxnnn (Nb,Na,N)   = (((N-1)*Ntotc+Na-1)*Ntotc+Nb-1)*Ngroup - 1
-      Jxmmm (Mb,Ma,M,J) = (((J+M)*Ntotc+Ma-1)*Ntotc+Mb-1)*Lllmax
+      real(kind=8),parameter::Zero = 0.0d0
+      integer::Lllmmm
+      real(kind=8)::Eb
+      real(kind=8)::Ai, Ar, Br, C2
+      integer::Iq, Iso, isoMgr, isoNgr, Jx, Jxm, Jxn, Klmn, Kountr
+      integer::L, Mchan, Mchanx, Mgr, Nchan, Nchanx, Ngr
+      integer::ngroup
 !
 !
+      ngroup = calc%resData%getNumSpinGroups()
+
       CALL Zero_Array (Ccclll, Iq_Iso*Lllmax)
       CALL Findpr (Kkxlmn, Klmn)
 !
@@ -36,8 +58,8 @@
             Iso = Iq
             C2 = Zero
          END IF
-         DO Ngr=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo(spinNgr, Ngr) 
+         DO Ngr=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinNgr, Ngr)
             IF (spinNgr%getIncludeInCalc()) THEN               
                isoNgr = spinNgr%getIsotopeIndex() 
                IF (IsoNgr.EQ.Iso .OR. Nnniso.NE.Numiso) THEN
@@ -54,9 +76,9 @@
 ! ***                         reaction under consideration
                            IF (Ar.NE.Zero .OR. Ai.NE.Zero) THEN
 ! --------------------------------------------------------------------
-         Jxn = Jxnnn (Nchanx,Nchan,Ngr)
-         DO Mgr=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo(spinMgr, Mgr)  
+         Jxn = Jxnnn (Nchanx,Nchan,Ngr, Ntotc, Ngroup)
+         DO Mgr=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinMgr, Mgr)
             IF (spinMgr%getIncludeInCalc()) THEN   
                isoMgr = spinMgr%getIsotopeIndex() 
                IF (IsoMgr.EQ.Iso .OR. Nnniso.NE.Numiso) THEN
@@ -69,7 +91,7 @@
                            Br = Ar*Crssx(1,Mchanx,Mchan,Mgr) +   &
                                 Ai*Crssx(2,Mchanx,Mchan,Mgr)
                            IF (Br.NE.Zero) THEN
-                              Jxm = Jxmmm (Mchanx,Mchan,Mgr,Jxn)
+                              Jxm = Jxmmm (Mchanx,Mchan,Mgr,Jxn, Ntotc, Lllmax)
                               DO L=1,Lllmax
                                  Jx = Jxm + L
                                  CALL Find_Kountr_Jx (Kxlmn, Kkxlmn,   &
@@ -109,32 +131,33 @@
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Setsel (Siabnd, Parpmc, Isopmc,   &
+      SUBROUTINE Setsel (calc, Siabnd, Parpmc, Isopmc,   &
            Sigsin, Crss, Su, Eb)
 !
 ! *** Purpose -- set self-indication transmission if needed
 !
-      use fixedi_m
-      use ifwrit_m
-      use lbro_common_m
-      use constn_common_m
-      use EndfData_common_m
+      use fixedi_m, only : Ncrsss, Numiso, Numpmc
+      use constn_common_m, only : Fourpi
       use SammySpinGroupInfo_M
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 !
-      DIMENSION Siabnd(*)
-      DIMENSION Parpmc(*), Isopmc(*), Sigsin(*)
-      DIMENSION Crss(Ncrsss,*)
+      class(XctCrossCalc)::calc
+      real(kind=8)::Siabnd(*)
+      real(kind=8)::Parpmc(*), Sigsin(*)
+      integer::Isopmc(*)
+      real(kind=8)::Crss(Ncrsss,*)
+      real(kind=8)::Su, Eb
       type(SammySpinGroupInfo)::spinInfo
-      DATA Zero /0.0d0/
+      real(kind=8), parameter::Zero=0.0d0
+      integer::Iso, isoN, N, Nnnnis
 !
 !
       Nnnnis = Numiso
       IF (Nnnnis.EQ.0) Nnnnis = 1
       DO Iso=1,Nnnnis
          Sitota = Zero
-         DO N=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo(spinInfo, N) 
+         DO N=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinInfo, N)
             IF (spinInfo%getIncludeInCalc()) THEN
                isoN = spinInfo%getIsotopeIndex()
                IF (Numiso.GT.0 .AND. IsoN.EQ.Iso) THEN
@@ -153,24 +176,26 @@
 !
 ! ______________________________________________________________________
 !
-      SUBROUTINE Setoth ( Pardet,   &
+      SUBROUTINE Setoth ( calc, Pardet,   &
          Igrdet, Parpmc, Isopmc, Sigxxx, Crss, Termf, Termfx, Su, Eb)
 !
 ! *** Purpose -- Set "other" cross sections as needed
 !
-      use fixedi_m
-      use ifwrit_m
-      use lbro_common_m
-      use constn_common_m
-      use EndfData_common_m
+      use fixedi_m, only : Nnnsig, Ncrsss, Lllmax, Nnniso, Numdet, Numiso, Numpmc
+      use ifwrit_m, only : Kaverg, Kssmsc, Nfissl, Ntgrlq, Kcros, Kfinit
+      use lbro_common_m, only : Yangle
+      use constn_common_m, only : Fourpi
       use SammySpinGroupInfo_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
 !
-      DIMENSION Pardet(*),   &
-           Igrdet(*), Parpmc(*), Isopmc(*),   &
+      class(XctCrossCalc)::calc
+      real(kind=8):: Pardet(*),  Parpmc(*),   &
            Sigxxx(Nnnsig,*), Crss(Ncrsss,*), Termf(*), Termfx(*)
+      integer:: Igrdet(*),  Isopmc(*)
       type(SammySpinGroupInfo)::spinInfo
-      DATA Zero /0.0d0/
+      real(kind=8),parameter::Zero = 0.0d0
+      real(kind=8)::Su, Eb, Val
+      real(kind=8)::abnSpin, Terma, Termax, Termff, Termn, Termxx
+      integer::I, Iso, Isox, N, Nnn, isoN
 !
 !
       Nnn = Lllmax + 1
@@ -187,8 +212,8 @@
                Termfx(I) = Zero
             END DO
          END IF
-         DO N=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo(spinInfo, N) 
+         DO N=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinInfo, N)
             IF (spinInfo%getIncludeInCalc()) THEN
                isoN= spinInfo%getIsotopeIndex()
                abnSpin = spinInfo%getAbundance()
@@ -277,7 +302,9 @@
 ! ***               but the Jx do not necessarily occur in order.
 ! ***            Ergo, must search through all of Kxlmn each time.
 !
-      DIMENSION Kxlmn(Kkxlmn)
+      integer::Kxlmn(Kkxlmn)
+      integer::Kkxlmn, Jx, Kountr, K
+      integer::N,M,I, Klmn
       N = 1
       Kountr = 0
       IF (Jx.LT.Kxlmn(1)) THEN
@@ -312,3 +339,5 @@
       END IF
       RETURN
       END
+
+end module mxct21_m
diff --git a/sammy/src/xct/mxct22.f90 b/sammy/src/xct/mxct22.f90
index aa39fb4c5..5899b039a 100644
--- a/sammy/src/xct/mxct22.f90
+++ b/sammy/src/xct/mxct22.f90
@@ -1,8 +1,26 @@
+module mxct22_m
+use XctCrossCalc_M
+implicit none
+
+private Jxnnn, Jxmmm
+contains
+
+     integer function Jxnnn (Nb,Na,N, Ntotc, Ngroup)
+        integer Nb,Na,N, Ntotc, Ngroup
+
+        Jxnnn = (((N-1)*Ntotc+Na-1)*Ntotc+Nb-1)*Ngroup - 1
+     end function
+     integer function Jxmmm (Mb,Ma,M,J, Ntotc, Lllmax)
+        integer Mb,Ma,M,J, Ntotc, Lllmax
+
+        Jxmmm   = (((J+M)*Ntotc+Ma-1)*Ntotc+Mb-1)*Lllmax
+     end function
+
 !
 !
 ! ______________________________________________________________________
 !
-      SUBROUTINE Derleg (   &
+      SUBROUTINE Derleg (calc,   &
          Sigxxx, Dasigx, Ccclll, Dddlll, Xlmn, Kxlmn, Crssx,   &
          Derivx, Isopar, Echan, Cmlab, Iso_Qv, Lllmmm, Eb)
 !
@@ -11,27 +29,35 @@
 ! ***
 ! *** Note    -- Enter this routine only if Ndasig > 0
 !
-      use fixedi_m
-      use ifwrit_m
-      use lbro_common_m
-      use EndfData_common_m
+      use fixedi_m, only : Nnnsig, Kkxlmn, Ntotc, Iq_Iso, Iq_Val, Lllmax, Ndasig, &
+                           Nfpiso, Nnniso, Numiso
+      use ifwrit_m, only : Nnpar
       use SammySpinGroupInfo_M
       use SammyIsoInfo_M
       use mdat9_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      use mxct21_m
 !
-      DIMENSION    &
+      class(XctCrossCalc)::calc
+      integer::Lllmmm
+      real(kind=8)::Eb
+      real(kind=8):: &
          Sigxxx(Nnnsig,*),   &
          Dasigx(Nnnsig,*), Ccclll(Lllmmm,*), Dddlll(Lllmmm,*),   &
-         Xlmn(*), Kxlmn(Kkxlmn), Crssx(2,Ntotc,Ntotc,*),   &
-         Derivx(2,Ntotc,Ntotc,Nnpar,*), Isopar(*),   &
-         Echan(Ntotc,*), Cmlab(3,*), Iso_Qv(*)
+         Xlmn(*), Crssx(2,Ntotc,Ntotc,*),   &
+         Derivx(2,Ntotc,Ntotc,Nnpar,*),   &
+         Echan(Ntotc,*), Cmlab(3,*)
+      integer::  Kxlmn(Kkxlmn),Isopar(*),  Iso_Qv(*)
       type(SammySpinGroupInfo)::spinMgr, spinNgr
       type(SammyIsoInfo)::isoInfo
-      DATA Zero /0.0d0/
-      Jxnnn (Nb,Na,N)   = (((N-1)*Ntotc+Na-1)*Ntotc+Nb-1)*Ngroup - 1
-      Jxmmm (Mb,Ma,M,J) = (((J+M)*Ntotc+Ma-1)*Ntotc+Mb-1)*Lllmax
+      real(kind=8),parameter:: Zero = 0.0d0
+      real(kind=8)::val
+      real(kind=8)::Ai, Ar, Bi, Br, C2, Dai, Dar, Dbi, Dbr, Dr
+      integer::Ifl, Iipar, Iq, Iso, IsoMgr, isoNgr, Jx, Jxm, Jxn
+      integer::Klmn, Kountr, Mchan, Mchanx,  Mgr, Nchan, Nchanx
+      integer::Ngr, L
+      integer::ngroup
 !
+      ngroup = calc%resData%getNumSpinGroups()
       CALL Zero_Array (Dddlll, Lllmmm*Ndasig)
       CALL Findpr (Kkxlmn, Klmn)
 !
@@ -43,8 +69,8 @@
             Iso = Iq
             C2 = Zero
          END IF
-         DO Ngr=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo( spinNgr, Ngr)              
+         DO Ngr=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo( spinNgr, Ngr)
             IF (spinNgr%getIncludeInCalc()) THEN
                isoNgr =  spinNgr%getIsotopeIndex()
                IF (Numiso.LE.0 .OR. isoNgr.EQ.Iso .OR.   &
@@ -59,9 +85,9 @@
                            Ai = Crssx(2,Nchanx,Nchan,Ngr)
                            IF (Ar.NE.Zero .OR. Ai.NE.Zero) THEN
 ! ----------------------------------------------------------
-         Jxn = Jxnnn (Nchanx,Nchan,Ngr)
-         DO Mgr=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo( spinMgr, Mgr) 
+         Jxn = Jxnnn (Nchanx,Nchan,Ngr, Ntotc, Ngroup)
+         DO Mgr=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo( spinMgr, Mgr)
             IF (spinMgr%getIncludeInCalc()) THEN               
                isoMgr = spinMgr%getIsotopeIndex()
                IF (Numiso.LE.0 .OR. IsoMgr.EQ.Iso .OR.   &
@@ -74,7 +100,7 @@
                         DO Mchanx=1,spinMgr%getNumEntryChannels()
                            Br = Crssx(1,Mchanx,Mchan,Mgr)
                            Bi = Crssx(2,Mchanx,Mchan,Mgr)
-                           Jxm = Jxmmm (Mchanx,Mchan,Mgr,Jxn)
+                           Jxm = Jxmmm (Mchanx,Mchan,Mgr,Jxn, Ntotc, Lllmax)
 ! ----------------------------------------------------------
          DO Iipar=1,Ndasig
             Dar = Derivx(1,Nchanx,Nchan,Iipar,Ngr)
@@ -126,14 +152,14 @@
 ! ##################### maybe NOT CORRECT YET FOR Iq_Val>0
       IF (Nfpiso.GT.0) THEN
          DO Iso=1,Numiso
-            call resParData%getIsoInfo(isoInfo, Iso)
+            call calc%resData%getIsoInfo(isoInfo, Iso)
             Ifl = isoInfo%getFitOption()
             IF (Ifl.GT.0) THEN
                Ifl = Ifl
                Isopar(Ifl) = Iso
                DO L=1,Lllmax
                   Dddlll(L,Ifl) = Ccclll(L,Iso)/   &
-                   resParData%getAbundanceByIsotope(Iso)
+                   calc%resData%getAbundanceByIsotope(Iso)
                END DO
             END IF
          END DO
@@ -148,3 +174,5 @@
       END IF
       RETURN
       END
+
+end module mxct22_m
diff --git a/sammy/src/xct/mxct25.f90 b/sammy/src/xct/mxct25.f90
index c78ec8df7..b611e314a 100644
--- a/sammy/src/xct/mxct25.f90
+++ b/sammy/src/xct/mxct25.f90
@@ -1,31 +1,41 @@
+module mxct25_m
+implicit none
+public Kount_Clbsch, Clbsch
+contains
 !
 !
 ! ------------------------------------------------------------------
 !
-      SUBROUTINE Kount_Clbsch (Zke, Alj, Kslow)
+      SUBROUTINE Kount_Clbsch (Zke, Alj, Kslow, resParData, Kkxlmn, Lllmax, Ntotc)
 !
 ! *** Purpose -- Generate Alj, count non-zero elements of Xlmn
 !
-      use fixedi_m
-      use fixedr_m
-      use lbro_common_m
-      use EndfData_common_m
       use SammySpinGroupInfo_M
       use RMatResonanceParam_M
       use SammyChannelInfo_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      use SammyRMatrixParameters_M
+      IMPLICIT None
 !
-      EXTERNAL Del2, Ww, Gll
-      DIMENSION  Zke(Ntotc,*), Alj(Ntriag,*)
+      integer::Kkxlmn, Lllmax, Ntotc
+      real(kind=8)::Zke(:,:), Alj(:,:)
       type(SammySpinGroupInfo)::spinMgr, spinNgr
       type(RMatSpinGroup)::spinGroupNgr, spinGroupMgr
       type(RMatChannelParams)::chanNc, chanNcx
       type(RMatChannelParams)::chanMc, chanMcx
       type(SammyChannelInfo)::chanInfoNc, chanInfoNcx
-      type(SammyChannelInfo)::chanInfoMc, chanInfoMcx
-      DATA Zero /0.0d0/, Tenth /0.10D0/, One /1.0d0/, Two /2.0d0/
-!
-      CALL Zero_Array (Alj, Ntriag*Ngroup)
+      type(SammyChannelInfo)::chanInfoMc, chanInfoMcx     
+      real(kind=8)::A, Aaa, Acn, Acnm, Aj, B, BX, C, Cjn, Cjnx, Cx
+      real(kind=8)::D, Dx, EL, Elm, Elmx, Eln, Elnx
+      integer::isoMgr, IsoNgr, Jn, Jnx, K, Km, Kn, Kountr, L
+      integer::Ll, Llm, Llmx, Lln, Llnx, Lmax, Lmaxx, Lmin
+      integer::Lminx, Mc, Mcx, Mgr, N, Nc, Ncx, Ngr, Nmlmax, Nmlmin
+      real(kind=8)::Spinm, Spinn
+      integer::kslow
+      type(SammyRMatrixParameters)::resParData
+      real(kind=8),parameter::Zero=0.0d0, Tenth=0.10D0, One=1.0d0, Two=2.0d0
+
+!
+      Alj = 0.0d0
 !
       DO Ngr=1,resParData%getNumSpinGroups()
          call resParData%getSpinGroupInfo(spinNgr, Ngr)
@@ -112,8 +122,7 @@
                      call resParData%getSpinGroupInfo(spinMgr, Mgr)
                      call resParData%getSpinGroup(spinGroupMgr, spinMgr) 
                      isoMgr = spinMgr%getIsotopeIndex()
-                     IF (Numiso.LE.0 .OR. IsoNgr.EQ.IsoMgr)   &
-                        THEN
+                     IF (IsoNgr.EQ.IsoMgr) THEN
                         Spinm = dABS(spinGroupMgr%getJ())
                         DO Mc=1,spinMgr%getNumChannels()
                            call spinMgr%getChannelInfo(chanInfoMc, Mc)
@@ -224,9 +233,20 @@
       END
 !
 !
+      integer function Jxnnn (Nb,Na,N, Ntotc, Ngroup)
+         integer Nb,Na,N, Ntotc, Ngroup
+
+         Jxnnn = (((N-1)*Ntotc+Na-1)*Ntotc+Nb-1)*Ngroup - 1
+      end function
+      integer function Jxmmm (Mb,Ma,M,J, Ntotc, Lllmax)
+         integer Mb,Ma,M,J, Ntotc, Lllmax
+
+         Jxmmm   = (((J+M)*Ntotc+Ma-1)*Ntotc+Mb-1)*Lllmax
+      end function
+
 ! ------------------------------------------------------------------
 !
-      SUBROUTINE Clbsch (Xlmn, Alj, Kxlmn)
+      SUBROUTINE Clbsch (resparData, Xlmn, Alj, Kxlmn)
 !
 ! *** Purpose -- Generate portions of Clebsch-Gordan-ry for
 ! ***            combining two spin groups to give differential
@@ -234,32 +254,39 @@
 ! ***            and groups to use the regular version
 ! *** Output  -- Array Xlmn
 !
-      use fixedi_m
-      use fixedr_m
-      use lbro_common_m
-      use EndfData_common_m
+      use fixedi_m, only : Kkxlmn, Lllmax, Ntotc
+      use fixedr_m, only : Spninc
+      use lbro_common_m, only : Debug
       use SammySpinGroupInfo_M
       use RMatResonanceParam_M  
       use SammyChannelInfo_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      use SammyRMatrixParameters_M
 !
-      EXTERNAL Del2, Ww, Gll    
-      DIMENSION Alj(Ntriag,*)
-      DIMENSION Xlmn(*), Kxlmn(*)
+      real(kind=8):: Alj(:,:)
+      real(kind=8):: Xlmn(:)
+      integer::Kxlmn(:)
+      type(SammyRMatrixParameters)::resParData
       type(SammySpinGroupInfo)::spinNgr, spinMgr
       type(RMatSpinGroup)::spinGroupMgr, spinGroupNgr
       type(RMatChannelParams)::chanNc, chanNcx
       type(RMatChannelParams)::chanMc, chanMcx
       type(SammyChannelInfo)::chanInfoNc, chanInfoNcx
       type(SammyChannelInfo)::chanInfoMc, chanInfoMcx
-      DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
-      Jxnnn (Nb,Na,N)   = (((N-1)*Ntotc+Na-1)*Ntotc+Nb-1)*Ngroup - 1
-      Jxmmm (Mb,Ma,M,J) = (((J+M)*Ntotc+Ma-1)*Ntotc+Mb-1)*Lllmax
-!
-      CALL Zero_Array (Xlmn, Kkxlmn)
-      CALL Zero_Integer (Kxlmn, Kkxlmn)
+      integer :: ngrps
+      real(kind=8),parameter::Zero=0.0d0, One=1.0d0, Two=2.0d0
+      real(kind=8)::A, Aaa, Acn, Acn2, Acnm, B, Bx, C, Cjn
+      real(kind=8)::Cjnx, Cx, D, Dx, EL, Elm, Eln, Elnx
+      real(kind=8)::Ggginc, Elmx, Spinm, Spinn
+      integer::isoMgr, isoNgr, Jn, Jnx, Jx, Jxm, Jxn, K, Kk, Km, Kn
+      integer::Kountr, L, Ll, Llm, Lln, Llnx, Lmax, Lmaxx, Llmax, Llmx
+      integer::Lmin, Lminx, Mc, Mcx, Mgr, Nc, Ncx, Ngr, Nmlmax, Nmlmin
+
+!
+      Xlmn = 0.0d0
+      Kxlmn = 0
 !
       Kountr = 0
+      ngrps = resparData%getNumSpinGroups()
       Ggginc = Two*dABS(Spninc) + One
       Ggginc = One/Ggginc
       DO Ngr=1,resParData%getNumSpinGroups()
@@ -288,7 +315,7 @@
                END IF
                Acn = Alj(Kn,Ngr)
                IF (Acn.NE.Zero) THEN
-                  Jxn = Jxnnn (Ncx,Nc,Ngr)
+                  Jxn = Jxnnn (Ncx,Nc,Ngr, Ntotc, Ngrps)
                   Acn = Aaa*Acn
                   Llnx = chanNcx%getL()
                   Elnx = dFLOAT(Llnx)
@@ -299,8 +326,7 @@
                      call resParData%getSpinGroupInfo(spinMgr, Mgr)
                      call resParData%getSpinGroup(spinGroupMgr, spinMgr)
                      isoMgr = spinMgr%getIsotopeIndex() 
-                     IF (Numiso.LE.0 .OR. IsoNgr.EQ.IsoMgr)   &
-                        THEN
+                     IF (IsoNgr.EQ.IsoMgr)  THEN
                         Spinm = dABS(spinGroupMgr%getJ())
                         Acn2 = Acn
 !v                        IF (Mgr.NE.Ngr) THEN
@@ -327,7 +353,7 @@
       call spinMgr%getChannelInfo(chanInfoMcx, Mcx)
       call resParData%getChannel(chanMcx,chanInfoMcx)                                      
       IF (dABS(dABS(ChanMcx%getSch())-Cjnx).LE.0.001d0) THEN
-         Jxm = Jxmmm (Mcx,Mc,Mgr,Jxn)
+         Jxm = Jxmmm (Mcx,Mc,Mgr,Jxn, Ntotc, Lllmax)
          Llmx = chanMcx%getL()
          Elmx = dFLOAT(Llmx)
          Lminx = IABS(Llnx-Llmx) + 1
@@ -424,12 +450,12 @@
             call resParData%getSpinGroupInfo(spinNgr, Ngr)           
             DO Nc=1,spinNgr%getNumChannels()              
                DO Ncx=1,spinNgr%getNumEntryChannels()                 
-                  Jxn = Jxnnn (Ncx,Nc,Ngr)
+                  Jxn = Jxnnn (Ncx,Nc,Ngr, Ntotc, Ngrps)
                   DO Mgr=1,resParData%getNumSpinGroups()
                      call resParData%getSpinGroupInfo(spinMgr, Mgr)       
                      DO Mc=1,spinMgr%getNumChannels()        
                         DO Mcx=1,spinMgr%getNumEntryChannels()      
-                           Jxm = Jxmmm (Mcx,Mc,Mgr,Jxn)
+                           Jxm = Jxmmm (Mcx,Mc,Mgr,Jxn, Ntotc, Lllmax)
                            DO L=1,Lllmax
                               Jx = Jxm + L
                               IF (Kxlmn(Kountr+1).EQ.Jx) THEN
@@ -455,8 +481,10 @@
 ! ------------------------------------------------------------------
 !
       SUBROUTINE Zzz (X, Lm, Lllmmm)
-      DOUBLE PRECISION X(*)
-      DATA Zero /0.0d0/
+      real(kind=8):: X(*)
+      real(kind=8),parameter::Zero=0.0d0
+      integer::Lm, L, Lllmmm
+
       LM = Lllmmm
       DO L=1,Lllmmm
          IF (X(Lm).NE.Zero) RETURN
@@ -474,8 +502,8 @@
 ! ***   in a fashion which requires as few as possible multiplications
 ! ***   and which is not likely to produce over- or under-flow problems
 !
-      DOUBLE PRECISION A, B, C, D, E, F, G, H, Q, X, Y, Z, One, Zero
-      DATA Zero /0.0d0/, One /1.0d0/
+      DOUBLE PRECISION A, B, C, D, E, F, G, H, Q, X, Y, Z
+      real(kind=8),parameter::One=1.0d0, Zero=0.0d0
 !
       Q = 0.05d0
 ! ***         q is used to be sure our real numbers are "almost integers"
@@ -539,8 +567,11 @@
 ! ***            where 2m = L1 + L2 + Ll
 ! ***            using as few as possible multiplications
 !
-      DOUBLE PRECISION H, One
-      DATA One /1.0d0/
+
+      integer::L1, L2, Ll
+      DOUBLE PRECISION H
+      integer::M, Mg, M1, M2, M3, J, K
+      real(kind=8),parameter:: One=1.0d0
 !
       M = (L1+L2+Ll)/2
       Mg = M
@@ -594,7 +625,7 @@
 ! --------------------------------------------------------------------------
 !
       DOUBLE PRECISION FUNCTION Ww (EL1, Aj1, EL2, Aj2, CJ, EL)
-      DIMENSION NN(7)
+      integer:: NN(7)
 !
 ! *** purpose -- generate the following expression:
 !
@@ -609,8 +640,9 @@
 !                                   (L2+J1+Ll+j-n)! ]
 !
       DOUBLE PRECISION EL1, Aj1, EL2, Aj2, CJ, EL
-      DOUBLE PRECISION G, H, Q, Sign, Zero, One
-      DATA Q /0.01d0/, Zero /0.0d0/, One /1.0d0/
+      DOUBLE PRECISION G, H, Sign
+      integer::I1,I2, M1, M2, K1, K2, K3, N, I, Maxz, Minz, Nm
+      real(kind=8),parameter:: Q=0.01d0, Zero=0.0d0, One=1.0d0
 !
 !
       H = Zero
@@ -672,7 +704,8 @@
 ! ----------------------------------------------------------------------
 !
       SUBROUTINE Sortt (Nn)
-      DIMENSION Nn(7)
+      integer:: Nn(7)
+      integer::I,IS, K, Mm
 !
 ! *** Purpose -- Order Nn(I) from high to low
 !
@@ -690,3 +723,4 @@
       END DO
       RETURN
       END
+end module mxct25_m
diff --git a/sammy/src/xct/mxct26.f90 b/sammy/src/xct/mxct26.f90
index cc810c29f..55ce15f84 100644
--- a/sammy/src/xct/mxct26.f90
+++ b/sammy/src/xct/mxct26.f90
@@ -1,29 +1,32 @@
+module mxct26_m
+contains
 !
 !
 ! --------------------------------------------------------------
 !
-      Subroutine Find_If_Coulomb (Zeta, IfCoul, Ifdif)
-      use  EndfData_common_m
+      Subroutine Find_If_Coulomb (calc, IfCoul, Ifdif)
+      use CrossSectionCalculator_M
       use SammySpinGroupInfo_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      IMPLICIT none
+      class(CrossSectionCalculator)::calc
       type(SammySpinGroupInfo)::spinInfo
-      DIMENSION Zeta(*)
+
+      integer::IfCoul, Ifdif
+      integer::I, Nent, Nn
+      logical::hasCoulomb
 ! *** On output, IfCoul = Maximum number of entrance channels which
 ! ***    require Coulomb
+       IfCoul = 0
       IF (Ifdif.EQ.1) THEN
          Nn = 0
-         DO I=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo(spinInfo, I)
+         hasCoulomb = .false.
+         DO I=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinInfo, I)
             nent = spinInfo%getNumEntryChannels()
             IF (Nent.GT.Nn) Nn = Nent
+            if (calc%Zeta(1,I).ne.0.0d0) hasCoulomb = .true.
          END DO
-         IF (Zeta(1).NE.0.0d0) THEN
-            IfCoul = Nn
-         ELSE
-            IfCoul = 0
-         END IF
-      ELSE
-         IfCoul = 0
+         if (hasCoulomb)  IfCoul = Nn
       END IF
       RETURN
       END
@@ -32,15 +35,16 @@
 ! ______________________________________________________________________
 !
       SUBROUTINE Start_Coul (Zke, Ccoulx)
-      use fixedi_m
-      use fixedr_m
+      use fixedi_m, only  : Ntotc
       use  EndfData_common_m
       use SammySpinGroupInfo_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      DIMENSION Zke(Ntotc,*),   &
-         Ccoulx(Ntotc,*)
+      IMPLICIT None
+      real(kind=8)::Zke(Ntotc,*), Ccoulx(Ntotc,*)
       type(SammySpinGroupInfo)::spinInfo
-      DATA Hth/0.01d0/
+      real(kind=8)::A
+      integer::Igroup, Ich, Nenti, Ntoti
+      real(kind=8),parameter::Hth=0.01d0
+
       DO Igroup=1,resParData%getNumSpinGroups()
          call resParData%getSpinGroupInfo(spinInfo, Igroup)
          A = Hth*   &
@@ -59,13 +63,16 @@
 ! ______________________________________________________________________
 !
       SUBROUTINE Store_Coul (Ccoul, Dcoul, Crssx, Derivx, Ccoulx, Jdat)
-      use fixedi_m
-      use ifwrit_m
+      use fixedi_m, only : Ntotc, Ngroup
+      use ifwrit_m, only : Nnpar
       use  EndfData_common_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      DIMENSION Ccoul(2,Ntotc,Ngroup,*), Crssx(2,Ntotc,Ntotc,*),   &
-         Dcoul(2,Ntotc,Nnpar,Ngroup,*),   &
-         Derivx(2,Ntotc,Ntotc,Nnpar,*), Ccoulx(Ntotc,*)
+      IMPLICIT NONE
+      integer::Jdat
+      real(kind=8)::Ccoul(2,Ntotc,Ngroup,*), Crssx(2,Ntotc,Ntotc,*),   &
+             Dcoul(2,Ntotc,Nnpar,Ngroup,*),   &
+             Derivx(2,Ntotc,Ntotc,Nnpar,*), Ccoulx(Ntotc,*)
+      integer::igroup, Nn, Ix, Iipar
+
       DO Igroup=1,resParData%getNumSpinGroups()
          DO Nn=1,Ntotc
             DO Ix=1,2
@@ -87,8 +94,11 @@
 ! --------------------------------------------------------------
 !
       SUBROUTINE Get_Coul_Phase (Cr, Ci, Lspin, Echan, Zeta, Su)
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      DATA Zero /0.0d0/, One /1.0d0/
+      IMPLICIT NONE
+      real(kind=8)::Aa, Cr, Ci, Cc, Ccx, Ccy, Eta,Echan, Zeta, Su
+      integer::L, Ll, Lspin
+      real(kind=8)::Ss, Ssx, Ssy
+      real(kind=8),parameter::Zero=0.0d0, One=1.0d0
 ! *** Purpose -- Calculate {Cr + i Ci} = exp (i omega_{l,alpha} )
 ! ***            where omega = Sum (n=1 to Lspin) tan^{-1} (eta/n)
       Eta = Zeta/Dsqrt(Su-Echan)
@@ -127,3 +137,4 @@
       END IF
       RETURN
       END
+end module mxct26_m
diff --git a/sammy/src/xct/mxct30.f90 b/sammy/src/xct/mxct30.f90
index 5534d19c1..c07c9646e 100644
--- a/sammy/src/xct/mxct30.f90
+++ b/sammy/src/xct/mxct30.f90
@@ -1,8 +1,19 @@
+module mxct30_m
+implicit none
+
+public Clbsch_Slow
+
+contains
+integer function Jxlmn(Mb,Ma,M, Ntotc)
+   integer::Mb,Ma,M, Ntotc
+   Jxlmn = ((M-1)*Ntotc+Ma-1)*Ntotc + Mb
+end function
+
 !
 !
 ! ------------------------------------------------------------------
 !
-      SUBROUTINE Clbsch_Slow (Xlmn, Alj, Kxlmn)
+      SUBROUTINE Clbsch_Slow (resparData, Xlmn, Alj, Kxlmn)
 !
 ! *** Purpose -- Generate portions of Clebsch-Gordan-ry for
 ! ***            combining two spin groups to give differential
@@ -10,27 +21,34 @@
 ! ***            and groups to use the regular version
 ! *** Output  -- Array Xlmn
 !
-      use fixedi_m
-      use fixedr_m
-      use lbro_common_m
-      use EndfData_common_m
+      use fixedi_m, only : Kkxlmn, Lllmax, Numiso, Ntotc
+      use fixedr_m, only : Spninc
+      use lbro_common_m, only : Debug
       use SammySpinGroupInfo_M
       use RMatResonanceParam_M    
       use SammyChannelInfo_M
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-!
-      EXTERNAL Del2, Ww, Gll    
-      DIMENSION Alj(Ntriag,*)
-      DIMENSION Xlmn(*), Kxlmn(Kkxlmn,*)
+      use SammyRMatrixParameters_M
+      use mxct25_m
+!  
+      real(kind=8):: Alj(:,:)
+      real(kind=8)::Xlmn(*)
+      integer::Kxlmn(Kkxlmn,*)
+      type(SammyRMatrixParameters)::resParData
       type(SammySpinGroupInfo)::spinMgr, spinNgr
       type(RMatSpinGroup)::spinGroupMgr, spinGroupNgr    
       type(RMatChannelParams)::chanNc, chanNcx
       type(RMatChannelParams)::chanMc, chanMcx
       type(SammyChannelInfo)::chanInfoNc, chanInfoNcx
       type(SammyChannelInfo)::chanInfoMc, chanInfoMcx
-      DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
-      Jxlmn (Mb,Ma,M) = ((M-1)*Ntotc+Ma-1)*Ntotc + Mb
-!
+      real(kind=8),parameter:: Zero=0.0d0, One=1.0d0, Two = 2.0d0
+      real(kind=8)::A, Aaa, Acn, Acn2, B, Bx, C, Cjn, Cjnx, Cx, D
+      real(kind=8)::Dx, EL, Elm, Elmx, Eln, Elnx, Ggginc, Acnm
+      real(kind=8)::Spinm, Spinn
+      integer::isoMgr, isoNgr, Jn, Jnx, Jxm, K, Kk, Km, Kn, Kountr
+      integer::L, Ll, Llm, Jxn, Llmx, Lln, Lmax, Lmaxx, Lmin, Lminx
+      integer::Mc, Mcx, Mgr, Nc, Ncx, Ngr, Nmlmax, Nmlmin, Llnx
+
+
       CALL Zero_Array (Xlmn, Kkxlmn)
       CALL Zero_Integer (Kxlmn, 3*Kkxlmn)
 !
@@ -63,7 +81,7 @@
                END IF
                Acn = Alj(Kn,Ngr)
                IF (Acn.NE.Zero) THEN
-                  Jxn = Jxlmn (Ncx,Nc,Ngr)
+                  Jxn = Jxlmn (Ncx,Nc,Ngr, Ntotc)
                   Acn = Aaa*Acn
                   Llnx = chanNcx%getL()
                   Elnx = dFLOAT(Llnx)
@@ -74,9 +92,8 @@
                      call resParData%getSpinGroupInfo(spinMgr, Mgr)
                      call resParData%getSpinGroup(spinGroupMgr, spinMgr)
                      isoMgr = spinMgr%getIsotopeIndex()
-                     IF (Numiso.LE.0 .OR. IsoNgr.EQ.IsoMgr)   &
+                     IF (IsoNgr.EQ.IsoMgr)   &
                         THEN
-                        Spinm = dABS(spinGroupMgr%getJ())
                         Acn2 = Acn
 !v                        IF (Mgr.NE.Ngr) THEN
 !v                           Acn2 = Two*Acn
@@ -102,7 +119,7 @@
       call spinMgr%getChannelInfo(chanInfoMcx, Mcx)
       call resParData%getChannel(chanMcx,chanInfoMcx)                                       
       IF (dABS(dABS(chanMcx%getSch())-Cjnx).LE.0.001d0) THEN
-         Jxm = Jxlmn (Mcx,Mc,Mgr)
+         Jxm = Jxlmn (Mcx,Mc,Mgr, Ntotc)
          Llmx = chanMcx%getL()
          Elmx = dFLOAT(Llmx)
          Lminx = IABS(Llnx-Llmx) + 1
@@ -200,12 +217,12 @@
             call resParData%getSpinGroupInfo(spinNgr, Ngr)   
             DO Nc=1,spinNgr%getNumChannels()
                DO Ncx=1,spinNgr%getNumEntryChannels()  
-                  Jxn = Jxlmn (Ncx,Nc,Ngr)
+                  Jxn = Jxlmn (Ncx,Nc,Ngr, Ntotc)
                   DO Mgr=1,resParData%getNumSpinGroups()
                      call resParData%getSpinGroupInfo(spinMgr, Mgr)    
                      DO Mc=1,spinMgr%getNumChannels()
                         DO Mcx=1,spinMgr%getNumEntryChannels() 
-                           Jxm = Jxlmn (Mcx,Mc,Mgr)
+                           Jxm = Jxlmn (Mcx,Mc,Mgr, Ntotc)
                            DO L=1,Lllmax
                               IF (Kxlmn(Kountr+1,1).EQ.L   .AND.   &
                                   Kxlmn(Kountr+1,2).EQ.Jxm .AND.   &
@@ -226,3 +243,4 @@
       END IF
       RETURN
       END
+end module mxct30_m
diff --git a/sammy/src/xct/mxct31.f90 b/sammy/src/xct/mxct31.f90
index 8c4be5a02..9409f2309 100644
--- a/sammy/src/xct/mxct31.f90
+++ b/sammy/src/xct/mxct31.f90
@@ -1,27 +1,47 @@
+module mxct31_m
+use XctCrossCalc_M
+implicit none
+public Setleg_Slow, Find_Kountr_Jx_Slow
+contains
+
+integer function Jxlmn(Mb,Ma,M, Ntotc)
+   integer::Mb,Ma,M, Ntotc
+   Jxlmn = ((M-1)*Ntotc+Ma-1)*Ntotc + Mb
+end function
+
 !
 !
 ! --------------------------------------------------------------
 !
-      SUBROUTINE Setleg_Slow (Sigxxx,   &
+      SUBROUTINE Setleg_Slow (calc, Sigxxx,   &
          Ccclll, Xlmn, Kxlmn, Crssx, Echan, Cmlab, Iso_Qv, Lllmmm, Eb)
 !
 ! *** Purpose -- set Ccclll(L,Iso) = coefficient of Legendre polynomial
 ! ***                         P-sub-(L-1) for Isotope Iso
 !
-      use fixedi_m
-      use ifwrit_m
-      use lbro_common_m
-      use EndfData_common_m
+      !use fixedi_m
+      !use ifwrit_m
+      !use lbro_common_m
+      !use EndfData_common_m
+      use fixedi_m, only : Nnnsig, Kkxlmn, Ntotc, Iq_Iso, Iq_Val, Lllmax, &
+                                 Nnniso, Numiso
       use SammySpinGroupInfo_M
       use mdat9_m
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
 !
-      DIMENSION    &
-         Sigxxx(Nnnsig,*), Ccclll(Lllmmm,*), Xlmn(*), Kxlmn(Kkxlmn,3),    &
-         Crssx(2,Ntotc,Ntotc,*), Echan(Ntotc,*), Cmlab(3,*), Iso_Qv(*)
+      class(XctCrossCalc)::calc
+      integer::Lllmmm
+      real(kind=8)::Eb
+      real(kind=8)::    &
+         Sigxxx(Nnnsig,*), Ccclll(Lllmmm,*), Xlmn(*),    &
+         Crssx(2,Ntotc,Ntotc,*), Echan(Ntotc,*), Cmlab(3,*)
+      integer::     Kxlmn(Kkxlmn,3),Iso_Qv(*)
       type(SammySpinGroupInfo)::spinMgr, spinNgr
-      DATA Zero /0.0d0/
-      Jxlmn (Mb,Ma,M) = ((M-1)*Ntotc+Ma-1)*Ntotc + Mb
+      real(kind=8),parameter:: Zero = 0.0d0
+      real(kind=8)::Ai, Ar, Br, C2
+      integer::Iq, Iso, isoMgr, isoNgr, Jxm, Jxn
+      integer::Klmn, Kountr, L, Mchan, Mchanx, Mgr, Ngr
+      integer::Nchan, Nchanx
 !
 !
       CALL Zero_Array (Ccclll, Iq_Iso*Lllmax)
@@ -35,8 +55,8 @@
             Iso = Iq
             C2 = Zero
          END IF
-         DO Ngr=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo(spinNgr, Ngr) 
+         DO Ngr=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinNgr, Ngr)
             IF (spinNgr%getIncludeInCalc()) THEN
                isoNgr = spinNgr%getIsotopeIndex()
                IF (IsoNgr.EQ.Iso .OR. Nnniso.NE.Numiso) THEN
@@ -53,9 +73,9 @@
 ! ***                         reaction under consideration
                            IF (Ar.NE.Zero .OR. Ai.NE.Zero) THEN
 ! --------------------------------------------------------------------
-         Jxn = Jxlmn (Nchanx,Nchan,Ngr)
-         DO Mgr=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo(spinMgr, Mgr) 
+         Jxn = Jxlmn (Nchanx,Nchan,Ngr, Ntotc)
+         DO Mgr=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinMgr, Mgr)
             IF (spinMgr%getIncludeInCalc()) THEN               
                isoMgr = spinMgr%getIsotopeIndex()
                IF (IsoMgr.EQ.Iso .OR. Nnniso.NE.Numiso) THEN
@@ -68,7 +88,7 @@
                            Br = Ar*Crssx(1,Mchanx,Mchan,Mgr) +   &
                                 Ai*Crssx(2,Mchanx,Mchan,Mgr)
                            IF (Br.NE.Zero) THEN
-                              Jxm = Jxlmn (Mchanx,Mchan,Mgr)
+                              Jxm = Jxlmn (Mchanx,Mchan,Mgr, Ntotc)
                               DO L=1,Lllmax
                                  CALL Find_Kountr_Jx_Slow (Kxlmn,   &
                                     Kkxlmn, L, Jxm, Jxn, Kountr, Klmn)
@@ -118,7 +138,9 @@
 ! ***               but the Jx do not necessarily occur in order.
 ! ***            Ergo, must search through all of Kxlmn each time.
 !
-      DIMENSION Kxlmn(Kkxlmn,3)
+      integer::L, Jxm, Jxn, Kountr, K, Kkxlmn
+      integer:: Kxlmn(Kkxlmn,3)
+      integer::Kk, Lw, Mw1, Mw2, Nw1, Nw2
       Kountr = 0
       Kk = K
       CALL Wherek (Kxlmn(1,3), Jxn, Kkxlmn, Kk, Nw1, Nw2)
@@ -138,7 +160,9 @@
       SUBROUTINE Wherek (Kx, Jx, Nn, Kk, Nw1, Nw2)
       use mdat9_m
 ! *** Purpose -- Locate the first and last positions for which Kx = Jx
-      DIMENSION Kx(*)
+      integer::Kx(*)
+      integer::Jx, Nn, Kk, Nw1, Nw2
+      integer::I, Iz, K, Km, M, Mz, N, Nz
       Nw1 = 0
       Nw2 = 0
       IF (Kx(1 ).GT.Jx) RETURN
@@ -234,7 +258,9 @@
       SUBROUTINE Wherel (Kx, Jx, Nn, Lw)
       use mdat9_m
 ! *** Purpose -- Locate position Lw for which Kx(Lw) = Jx
-      DIMENSION Kx(*)
+      integer:: Kx(*)
+      integer::Jx, Nn, Lw
+      integer::I,K, M, N
       Lw = 0
       IF (Kx(1 ).GT.Jx) RETURN
       IF (Kx(Nn).LT.Jx) RETURN
@@ -260,3 +286,4 @@
       END DO
       RETURN
       END
+end module mxct31_m
diff --git a/sammy/src/xct/mxct32.f90 b/sammy/src/xct/mxct32.f90
index 478eebe77..cfc6a5a5f 100644
--- a/sammy/src/xct/mxct32.f90
+++ b/sammy/src/xct/mxct32.f90
@@ -1,8 +1,20 @@
+module mxct32_m
+use XctCrossCalc_M
+IMPLICIT none
+
+public Derleg_Slow
+
+contains
+
+integer function Jxlmn(Mb,Ma,M, Ntotc)
+   integer::Mb,Ma,M, Ntotc
+   Jxlmn = ((M-1)*Ntotc+Ma-1)*Ntotc + Mb
+end function
 !
 !
 ! ______________________________________________________________________
 !
-      SUBROUTINE Derleg_Slow (   &
+      SUBROUTINE Derleg_Slow ( calc,    &
          Sigxxx, Dasigx, Ccclll, Dddlll, Xlmn, Kxlmn, Crssx,   &
          Derivx, Isopar, Echan, Cmlab, Iso_Qv, Lllmmm, Eb)
 !
@@ -11,25 +23,33 @@
 ! ***
 ! *** Note    -- Enter this routine only if Ndasig > 0 and Kslow=1
 !
-      use fixedi_m
-      use ifwrit_m
-      use lbro_common_m
-      use EndfData_common_m
+      use fixedi_m, only : Nnnsig, Iq_Iso, Iq_val, &
+                           Lllmax, Ndasig, Nfpiso, Nnniso, Numiso, Ntotc, &
+                           Kkxlmn
+      use ifwrit_m, only : Nnpar
       use SammySpinGroupInfo_M
       use SammyIsoInfo_M
       use mdat9_m
+      use  mxct31_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 !
-      DIMENSION    &
+      class(XctCrossCalc)::calc
+      real(kind=8)::Eb
+      integer::Lllmmm
+      real(kind=8)::    &
          Sigxxx(Nnnsig,*),   &
          Dasigx(Nnnsig,*), Ccclll(Lllmmm,*), Dddlll(Lllmmm,*),   &
-         Xlmn(*), Kxlmn(Kkxlmn,*), Crssx(2,Ntotc,Ntotc,*),   &
-         Derivx(2,Ntotc,Ntotc,Nnpar,*), Isopar(*),   &
-         Echan(Ntotc,*), Cmlab(3,*), Iso_Qv(*)
+         Xlmn(*), Crssx(2,Ntotc,Ntotc,*),   &
+         Derivx(2,Ntotc,Ntotc,Nnpar,*),     &
+         Echan(Ntotc,*), Cmlab(3,*)
+      integer:: Kxlmn(Kkxlmn,*), Isopar(*), Iso_Qv(*)
       type(SammySpinGroupInfo)::spinNgr, spinMgr
       type(SammyIsoInfo)::isoInfo
-      DATA Zero /0.0d0/
-      Jxlmn (Mb,Ma,M) = ((M-1)*Ntotc+Ma-1)*Ntotc + Mb
+      real(kind=8),parameter :: Zero = 0.0d0
+      real(kind=8)::val
+      real(kind=8)::Ai, Ar, Bi, Br, C2, Dai, Dar, Dbi, Dbr, Dr
+      integer::Ifl, Iipar, Iq, Iso, isoMgr, isoNgr, Jxm, Jxn,  Klmn, Kountr
+       integer::L,Mchan, Mchanx, Mgr, Nchan, Nchanx, Ngr
 !
       CALL Zero_Array (Dddlll, Lllmmm*Ndasig)
       CALL Findpr (Kkxlmn, Klmn)
@@ -42,8 +62,8 @@
             Iso = Iq
             C2 = Zero
          END IF
-         DO Ngr=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo(spinNgr, Ngr) 
+         DO Ngr=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinNgr, Ngr)
             IF (spinNgr%getIncludeInCalc()) THEN               
                isoNgr = spinNgr%getIsotopeIndex()
                IF (Numiso.LE.0 .OR. IsoNgr.EQ.Iso .OR.   &
@@ -58,9 +78,9 @@
                            Ai = Crssx(2,Nchanx,Nchan,Ngr)
                            IF (Ar.NE.Zero .OR. Ai.NE.Zero) THEN
 ! ----------------------------------------------------------
-         Jxn = Jxlmn (Nchanx,Nchan,Ngr)
-         DO Mgr=1,resParData%getNumSpinGroups()
-            call resParData%getSpinGroupInfo(spinMgr, Mgr) 
+         Jxn = Jxlmn (Nchanx,Nchan,Ngr, Ntotc)
+         DO Mgr=1,calc%resData%getNumSpinGroups()
+            call calc%resData%getSpinGroupInfo(spinMgr, Mgr)
             IF (spinMgr%getIncludeInCalc()) THEN
                isoMgr = spinMgr%getIsotopeIndex() 
                IF (Numiso.LE.0 .OR. IsoMgr.EQ.Iso .OR.   &
@@ -73,7 +93,7 @@
                         DO Mchanx=1,spinMgr%getNumEntryChannels()
                            Br = Crssx(1,Mchanx,Mchan,Mgr)
                            Bi = Crssx(2,Mchanx,Mchan,Mgr)
-                           Jxm = Jxlmn (Mchanx,Mchan,Mgr)
+                           Jxm = Jxlmn (Mchanx,Mchan,Mgr, Ntotc)
 ! ----------------------------------------------------------
          DO Iipar=1,Ndasig
             Dar = Derivx(1,Nchanx,Nchan,Iipar,Ngr)
@@ -124,14 +144,14 @@
 ! ##################### maybe NOT CORRECT YET FOR Iq_Val>0
       IF (Nfpiso.GT.0) THEN
          DO Iso=1,Numiso
-            call resParData%getIsoInfo(isoInfo, Iso)
+            call calc%resData%getIsoInfo(isoInfo, Iso)
             Ifl = isoInfo%getFitOption()
             IF (Ifl.GT.0) THEN
                Ifl = Ifl
                Isopar(Ifl) = Iso
                DO L=1,Lllmax
                   Dddlll(L,Ifl) = Ccclll(L,Iso)/   &
-                resParData%getAbundanceByIsotope(Iso)
+                calc%resData%getAbundanceByIsotope(Iso)
                END DO
             END IF
          END DO
@@ -146,3 +166,4 @@
       END IF
       RETURN
       END
+end module mxct32_m
-- 
GitLab