diff --git a/ci/linux/core_linux.sh b/ci/linux/core_linux.sh
index ca0d32b6f245fa90ced213c23e24acffc5b76918..8cbc5c05155e8ab85620a41327c17c2dcf468653 100755
--- a/ci/linux/core_linux.sh
+++ b/ci/linux/core_linux.sh
@@ -94,6 +94,7 @@ echo "...CMAKE_BUILD_TYPE is $CMAKE_BUILD_TYPE"
 export SCALE_PACKAGES="-DSCALE_ENABLE_ALL_OPTIONAL_PACKAGES=OFF  \
                        -DSCALE_ENABLE_ScaleUtilsEndfLib=ON       \
                        -DSCALE_ENABLE_ScaleUtilsAmpxLib=ON       \
+                       -DSCALE_ENABLE_ScaleData=ON               \
                        -DSCALE_ENABLE_ScaleUtilsResonance=ON     \
                        -DSCALE_ENABLE_ScaleExternalGoogletest=ON"
 
diff --git a/sammy/cmake/TPLs/FindTPLSCALE.cmake b/sammy/cmake/TPLs/FindTPLSCALE.cmake
index 43568494a0ac6ebe844da7afd2bda8b0a8bf2c38..e4c4238450015ab2ad2902d191efdeec7d8b14d7 100755
--- a/sammy/cmake/TPLs/FindTPLSCALE.cmake
+++ b/sammy/cmake/TPLs/FindTPLSCALE.cmake
@@ -70,12 +70,15 @@ SET( REQ_SCALE_HEADERS
         ScaleUtils/Resonance/ResonanceCalc.h
         ScaleUtils/Resonance/ReichMoore/CalculatedParams.h
         ScaleUtils/Math/Lapack.h
+        ScaleData/Core/io/CoverxFile.h
    )
 SET( REQ_SCALE_LIBS
         ScaleUtilsEndfLib
         ScaleUtilsResonance
         ScaleUtilsMath
         Interface
+        ScaleData
+        AmpxLibraryUtilsLib
    )
 
 #TODO: If we set the upper level directory correctly, we can use the relative path
diff --git a/sammy/src/acs/macs0.f90 b/sammy/src/acs/macs0.f90
index c852d19dc5beeb08cad4fcfc5234ea640183ae59..16817cbe489a46ce1a7b92e27d863e8d0c202e81 100644
--- a/sammy/src/acs/macs0.f90
+++ b/sammy/src/acs/macs0.f90
@@ -113,7 +113,8 @@ module acs_m
            A_Kktheo , A_Kkthel , A_Kkdthe , A_Kkksiz )
 !
       ELSE IF (Kpntws.NE.0) THEN
-         CALL Make_File_3_Urr (I_Kkktyp , I_Kkkknt , A_Keeset , A_Kktheo )
+         CALL Make_File_3_Urr (I_Kkktyp , I_Kkkknt , A_Keeset , A_Kktheo, &
+                               A_Kkdthe,A_Kkkkex)
          WRITE (6,10100)
 10100    FORMAT ('Normal end of SAMMY-URR')
          STOP
diff --git a/sammy/src/acs/macs2.f90 b/sammy/src/acs/macs2.f90
index 789cd3d408ea534afdeb1d07c87f2bfb6c126180..6a91f972ad104cb9edacb142f8dd4ca3cef27e6c 100644
--- a/sammy/src/acs/macs2.f90
+++ b/sammy/src/acs/macs2.f90
@@ -20,7 +20,7 @@ module acs2_m
       use xxx2
       use xxx5
       use xxx9
-      use EndfData_common_m, only : radFitFlags
+      use EndfData_common_m, only : radFitFlags,expData
 
       real(8), intent(in)  :: Edirin(*),Sdirin(*),Edircp(*),Sdircp(*),Streng(Numelv,*),Distnt(Numelv,*),Engurr(*),Eeeset(*)
       real(8), intent(out) :: Theorl(Numelv,*),Dtheor(Nvpall,*),Dth(*),Cl(*)
@@ -55,7 +55,7 @@ module acs2_m
          Energy = Eeeset(Ienerg)
          CALL Find_Urr (Engurr, Energy, Kumurr, Numurr, Iset)
          Rho = Rk1*dSQRT(Energy)
-         IF (Ksolve.NE.2) CALL Zero_Array (Dth, Nvpall)
+         IF (Ksolve.NE.2 .or. Kpntws.ne.0 ) CALL Zero_Array (Dth, Nvpall)
 !
          Ipar = 0
 ! ***    Sum OVER PARTIAL WAVES
@@ -94,19 +94,19 @@ module acs2_m
             Bbr = Cosphi*Aar + Sinphi*Aai
             Bbi = Cosphi*Aai - Sinphi*Aar
 !
-            IF (Ksolve.NE.2 .AND. Iflstr(L,Kumurr).NE.0) THEN
+            IF ( (Ksolve.NE.2 .or. Kpntws.ne.0) .AND. Iflstr(L,Kumurr).NE.0) THEN
                Ipar = Iflstr(L,Kumurr)
                Dth(Ipar) = Ce*PL*Bbr*PI*(Streng(L,Kumurr)/Rk1)
 ! ***          Remember that it is ln(Streng) which is the
 ! ***             "u-parameter"; hence use Dth=Streng(L,Kumurr)*whatever
             END IF
 !
-            IF (Ksolve.NE.2 .AND. Ifldst(L,Kumurr).NE.0) THEN
+            IF ( (Ksolve.NE.2 .or. Kpntws.ne.0) .AND. Ifldst(L,Kumurr).NE.0) THEN
                Ipar = Ifldst(L,Kumurr)
                Dth(Ipar) = Two * Ce*PL*Bbi
             END IF
 !
-            IF (Ksolve.NE.2 .AND. L.LE.2 .AND. Iflggg(L,Kumurr).NE.0) THEN
+            IF ( (Ksolve.NE.2 .or. Kpntws.ne.0) .AND. L.LE.2 .AND. Iflggg(L,Kumurr).NE.0) THEN
                Ipar = Iflggg(L,Kumurr)
                Dth(Ipar) = Zero
             END IF
@@ -133,7 +133,7 @@ module acs2_m
    10       CONTINUE
             Theory(Ienerg) = Theory(Ienerg) + A
          END IF
-         IF (Ksolve.NE.2) THEN
+         IF (Ksolve.NE.2 .or. Kpntws.ne.0 ) THEN
             DO Ipar=1,Nvpall
                Dtheor(Ipar,Ienerg) = Dth(Ipar)
             END DO
@@ -157,7 +157,7 @@ module acs2_m
    20       CONTINUE
             Theory(Ienerg) = Theory(Ienerg) + A
          END IF
-         IF (Ksolve.NE.2) THEN
+         IF (Ksolve.NE.2 .or. Kpntws.ne.0 ) THEN
             DO Ipar=1,Nvpall
                Dtheor(Ipar,Ienerg) = Dth(Ipar)
             END DO
diff --git a/sammy/src/acs/macs3.f90 b/sammy/src/acs/macs3.f90
index 1b2ecdad952416615746e8cd43a6ac4a5d29275d..2cea94d89e63a0ae920796fbfa006a0d934435f6 100644
--- a/sammy/src/acs/macs3.f90
+++ b/sammy/src/acs/macs3.f90
@@ -77,7 +77,7 @@ module acs3_m
                END DO
             END IF
 !
-            IF (Ksolve.NE.2) CALL Zero_Array (Dth, Nvpall)
+            IF (Ksolve.NE.2 .or. Kpntws.ne.0 ) CALL Zero_Array (Dth, Nvpall)
 !
             Ce = Ca/Energy
             Rho = Rk1*dSQRT(Energy)
@@ -182,7 +182,7 @@ module acs3_m
                            END IF
 !
 ! ***                      Derivatives Calculation
-                           IF (Ksolve.NE.2) THEN
+                           IF ( Ksolve.NE.2 .or. Kpntws.ne.0 ) THEN
                               CALL Deriva (Gf(1,Kumjjv,Kumurr),       &
                                  Iflstr(1,Kumurr), Ifldst(1,Kumurr),  &
                                  Iflggg(1,Kumurr),                    &
@@ -254,7 +254,7 @@ module acs3_m
 !         # END IF (Type.EQ.Capt .AND. Ndircp.GT.1)
             END IF
 !
-            IF (Ksolve.NE.2) THEN
+            IF ( Ksolve.NE.2 .or. Kpntws.ne.0 ) THEN
                DO Ipar=1,Nvpall
                   Dtheor(Ipar,Ienerg) = Dth(Ipar) + Dtheor(Ipar,Ienerg)
                END DO
@@ -382,7 +382,7 @@ module acs3_m
                Ainela(7,Iz) = Ainela(7,Iz) + Term
             END IF
 !
-            IF (Ksolve.NE.2) THEN
+            IF ( Ksolve.NE.2 .or. Kpntws.ne.0 ) THEN
 ! ***          elastic -- DERIVATIVES WITH RESPECT TO Bn(M) = Aelast(2,M)
                IF (Jelast.GT.0) THEN
                   DO M=1,Jelast
@@ -447,7 +447,7 @@ module acs3_m
       END DO
       Dof_Fis = AnfO
 !
-      IF (Ksolve.NE.2) THEN
+      IF ( Ksolve.NE.2 .or. Kpntws.ne.0 ) THEN
          Der_Fis = Der_Fis*Bbb_Fis
          DO N=1,Jelast
             Aelast(6,N) = Aelast(6,N)*Aelast(2,N)
diff --git a/sammy/src/acs/macs5.f90 b/sammy/src/acs/macs5.f90
index b4faac9f3b46264a81c26e98857cd10781b7fd14..3d8f5ea9466c564104877ec4913041306377ae84 100644
--- a/sammy/src/acs/macs5.f90
+++ b/sammy/src/acs/macs5.f90
@@ -214,21 +214,32 @@ module acs5_m
 !
 ! --------------------------------------------------------------
 !
-      subroutine Make_File_3_Urr (Kktype, Kkkntx, Eeeset, Theory)
+      SUBROUTINE Make_File_3_Urr (Kktype, Kkkntx, Eeeset, Theory, Dtheor, Ex)
 !
       use fixedi_m
       use ifwrit_m
       use samxxx_common_m
       use namfil_common_m
+      use EndfData_common_m
+      use CoverxIO_m
+      use Hdf5IO_m
 
+      ! --- input variables -----------
       integer, intent(in) :: Kktype(*),Kkkntx(*)
-      real(8), intent(in) :: Eeeset(Kntmax,*),Theory(Kntmax,Kdtset,*)
-
+      real(8), intent(in) :: Eeeset(Kntmax,*),Theory(Kntmax,Kdtset,*), &
+                             Dtheor(Nvpall,Kntmax,*),Ex(*)
+      ! --- local variables -----------
       real(8) :: Awr,Defunc,Eee,Za
       integer :: I,Int(1),K,Ktype(5),Kz,Lrf,Mat,Mf,Mt,Nbt(1),Ndum1,Ndum2,Ndum3, &
-                 Ndum4,Nro,Ns
+                 Ndum4,Nro,Ns,energInd,theoInd,iE,iPar
       real(8), parameter :: Zero = 0.0d0
       integer, parameter :: L0 = 0, L1 = 1, Iu = 59
+      character(len=80) :: h5filename,coverxfilename
+      logical :: mt51 = .true.
+      ! --- class types ---------------
+      type(GridData) :: grid
+      type(CoverxIO) :: coverx
+      type(Hdf5IO)   :: h5
 !
       Krdndf = 1
       CALL Filopn (10, Fpntws, 0)
@@ -236,6 +247,22 @@ module acs5_m
                      Ndum3, Ndum4)
 ! *** Note that Eee, Defunc, Lrf, Kz, Ndum. are ignored here
       CLOSE (UNIT=10)
+
+      write(6,*) "*******************************************************"
+      write(6,*) "If Coverx for File 33 is wanted, type name of H5 parameter"
+      write(6,*) "covariance file: (else press enter)"
+      read(5,'(A80)') h5filename 
+      write(6,*) "Reading file: ", trim(h5filename)
+      write(6,*) "Please enter name of ouput Coverx file"
+      read(5,'(A80)') coverxfilename
+      write(6,*) "Writing file: ", trim(coverxfilename)
+      write(6,*) "*******************************************************"
+      ! if( filename has values ) else( dont make file 33 )
+      ! read covariance from h5 file into covData object
+      call h5%initialize()
+      call coverx%initialize()
+      call h5%readRPCM(covData,trim(h5filename))
+
 !
       CALL Newopn (59, Sam57x, 0)
 !
@@ -243,14 +270,37 @@ module acs5_m
       DO K=1,5
          Ktype(K) = 0
       END DO
+      energInd = 0
       DO I=1,Kdtset
          IF (Kktype(I).EQ.1) Ktype(1) = I
          IF (Kktype(I).EQ.2) Ktype(3) = I
          IF (Kktype(I).EQ.3) Ktype(4) = I
          IF (Kktype(I).EQ.4) Ktype(5) = I
          IF (Kktype(I).EQ.5) Ktype(2) = I
+         ! --------------------
+         ! Set up C++ API memory
+         ! TODO: should use GridData for all grid memory in acs/fff
+         ! --------------------
+         call grid%initialize()
+         if(I.eq.1) then
+            energInd = 1  ! 0 since angle-integ., no time-shift, no data
+         else
+            energInd = 0  ! save space, only store energies once
+         end if
+         call grid%setDataIndex(energInd+1)
+         call expData%addGrid(grid)
+         do iE=1,Kntmax
+            if(I.eq.1) call grid%addData(iE,energInd,Eeeset(iE,I))
+            call grid%addData(iE,energInd+1,Theory(iE,I,1))
+            do Ipar=1,Nvpall
+               call grid%addData(iE,energInd+1+Ipar,Dtheor(Ipar,iE,I))
+            end do
+         end do
       END DO
 !
+!     if( 2nd inelastic channel is energetically available )
+      if( Eeeset(Kntmax,1)>Ex(3) ) mt51 = .false.
+!
 ! *** Do-loop over cross section types
       DO K=1,5
          I = Ktype(K)
@@ -261,7 +311,11 @@ module acs5_m
                WRITE (59,10100)
 10100          FORMAT ('####' /, '#### Total cross section ####')
             ELSE IF (Kktype(I).EQ.2) THEN
-               Mt = 4
+               if( mt51 ) then
+                  Mt = 51
+               else 
+                  Mt = 4
+               end if
                WRITE (59,10200)
 10200          FORMAT ('####' /, '#### Inelastic cross section ####')
             ELSE IF (Kktype(I).EQ.3) THEN
@@ -281,6 +335,13 @@ module acs5_m
             Nro= 0
             Mf = 3
             Ns = 1
+
+            ! ------------------------------------
+            ! Set MAT and MT ID numbers for the coverx 
+            ! file calculations
+            ! ------------------------------------
+            call coverx%setMatMt(Mat,Mt)
+
 ! ***       first entry in ENDF/B-VI file -- HEAD record
             WRITE (59,10600)
 10600       FORMAT ('####' /, '####  Z_A   Weight')
@@ -293,6 +354,7 @@ module acs5_m
             CALL Tab1x (Zero, Zero, L0, L0, L1, Kkkntx(I), Nbt, Int, &
                Eeeset(1,I), Theory(1,I,1), Mat, Mf, Mt, Ns, Iu)
 !
+!                section end
             CALL Send_Ndf (Mat, Mf, Ns, Iu)
          END IF
       END DO
@@ -300,6 +362,16 @@ module acs5_m
 !
       CALL Rewrit_Ndf (Iu, Ns)
       CLOSE (UNIT=55)
+
+      !-----------------------------------------------------
+      ! Option to dump Coverx file with pointwise xs covariance
+      !-----------------------------------------------------
+      !
+      
+      call coverx%setCoverxFile(expData,covData)
+      call coverx%writeCoverxFile(trim(coverxfilename))
+      call coverx%writeFile33(Za,Awr)
+
       RETURN
 !
   200 WRITE (6,99996) Summry
diff --git a/sammy/src/blk/ExpPars_common.f90 b/sammy/src/blk/ExpPars_common.f90
index 80d9b9dff81e081b5511814126e9ae3ab3a9aa5a..4030a7a0a92cec48cf8210a3702a07423ea712c2 100644
--- a/sammy/src/blk/ExpPars_common.f90
+++ b/sammy/src/blk/ExpPars_common.f90
@@ -1,5 +1,5 @@
 ! /**********************************************************
-! * This common block style module is a temporary solution
+! * TODO: This common block style module is a temporary solution
 ! * used to pass the variables contained in C++ classes until
 ! * the code has been fully modernized. 
 ! **********************************************************/
diff --git a/sammy/src/endf/CovarianceData.h b/sammy/src/endf/CovarianceData.h
index dd88b1e742d24311730ae4ac28a27af56844d3be..0fed6c12aa4212bcf5f4431c3fb4ecdf8fb61bb6 100644
--- a/sammy/src/endf/CovarianceData.h
+++ b/sammy/src/endf/CovarianceData.h
@@ -176,7 +176,11 @@ namespace sammy{
         * Get the number of total parameters: The pup'd and non-pup'd
         */
        int getNumTotalParam() const{
-           return (int)paramIndex.size();
+           int np = (int)paramIndex.size();
+           if( np != 0 ) return np;
+           else{
+               return numParameters;
+           }
        }
 
        void resetTotalParam(int size){
diff --git a/sammy/src/endf/EndfData.h b/sammy/src/endf/EndfData.h
index 6469a33a15a2649a82e6f8f3bd4f055a1e57ed32..e43e31b8e877c02f8af6a5f204e1c4db19f9fbc3 100644
--- a/sammy/src/endf/EndfData.h
+++ b/sammy/src/endf/EndfData.h
@@ -1,5 +1,5 @@
 #ifndef SAMMY_ENDF_ENDF_DATA_H
-#define	SAMMY_ENDF_ENDF_DATA_H
+#define SAMMY_ENDF_ENDF_DATA_H
 
 #include "ScaleUtils/EndfLib/Tab1Container.h"
 #include "ScaleUtils/EndfLib/ResonanceParameters.h"
@@ -10,93 +10,92 @@
 
 
 namespace endf{
-  class EndfData{
+    class EndfData{
     public:
-       EndfData():mat(-1){}
-       EndfData(const EndfData & orig);
-       virtual ~EndfData(){}
-
-       std::unique_ptr<endf::ResonanceInfo> & getResonanceDataPtr();
-       endf::ResonanceInfo * getResonanceData(){
-	 return getResonanceDataPtr().get();
-       }
-       void setResonanceData(endf::ResonanceInfo * info){
-       	 resonanceData = std::unique_ptr<endf::ResonanceInfo>(info);
-       }
-
-       std::unique_ptr<endf::CovarianceContainer> & getCovarianceDataPtr();
-       endf::CovarianceContainer *  getCovarianceData(){
-	 return getCovarianceDataPtr().get();
-       }
-       void setCovarianceData(endf::CovarianceContainer * info){
-       	 covarianceData = std::unique_ptr<endf::CovarianceContainer>(info);
-       }
+        EndfData():mat(-1){}
+        EndfData(const EndfData & orig);
+        virtual ~EndfData(){}
+
+        std::unique_ptr<endf::ResonanceInfo> & getResonanceDataPtr();
+        endf::ResonanceInfo * getResonanceData(){
+            return getResonanceDataPtr().get();
+        }
+        void setResonanceData(endf::ResonanceInfo * info){
+            resonanceData = std::unique_ptr<endf::ResonanceInfo>(info);
+        }
+
+        std::unique_ptr<endf::CovarianceContainer> & getCovarianceDataPtr();
+        endf::CovarianceContainer *  getCovarianceData(){
+            return getCovarianceDataPtr().get();
+        }
+        void setCovarianceData(endf::CovarianceContainer * info){
+            covarianceData = std::unique_ptr<endf::CovarianceContainer>(info);
+        }
        
 
-       std::unique_ptr<endf::Tab1Container> & getPointwiseDataPtr();
-       endf::Tab1Container *  getPointwiseData(){
-	 return getPointwiseDataPtr().get();
-       }
-       void setPointwiseData(endf::Tab1Container * info){
-       	 pointwiseData = std::unique_ptr<endf::Tab1Container>(info);
-       }
+        std::unique_ptr<endf::Tab1Container> & getPointwiseDataPtr();
+        endf::Tab1Container *  getPointwiseData(){
+            return getPointwiseDataPtr().get();
+        }
+        void setPointwiseData(endf::Tab1Container * info){
+            pointwiseData = std::unique_ptr<endf::Tab1Container>(info);
+        }
        
 
-       std::unique_ptr<endf::EvaluationInfo> & getEvalInfoPtr();
-       endf::EvaluationInfo *  getEvalInfo(){
-	 return getEvalInfoPtr().get();
-       }
-       void setEvalInfo(endf::EvaluationInfo * info){
-       	   evalInfo = std::unique_ptr<endf::EvaluationInfo>(info);
-       }
+        std::unique_ptr<endf::EvaluationInfo> & getEvalInfoPtr();
+        endf::EvaluationInfo *  getEvalInfo(){
+            return getEvalInfoPtr().get();
+        }
+        void setEvalInfo(endf::EvaluationInfo * info){
+            evalInfo = std::unique_ptr<endf::EvaluationInfo>(info);
+        }
 
       
-       void setFileName(const std::string & file){
-	 fileName.clear();
-	 fileName.append(file);
-       }
+        void setFileName(const std::string & file){
+            fileName.clear();
+            fileName.append(file);
+        }
 
-       std::string getFileName() const{
-	 return fileName;
-       }
+        std::string getFileName() const{
+            return fileName;
+        }
 
-       void setMat(int m){
-	 mat = m;	
-       }
+        void setMat(int m){
+            mat = m;   
+        }
 
-       double getMat() const{
-	 return mat;
-       }
+        double getMat() const{
+            return mat;
+        }
 
-       void getMatFromEndf();
+        void getMatFromEndf();
 
-       void convertToRMat(endf::ResonanceInfo & info);
+        void convertToRMat(endf::ResonanceInfo & info);
        
-       const endf::Tab1 * getTab1ByMt(int mt);
+        const endf::Tab1 * getTab1ByMt(int mt);
 
 
-       const endf::Tab1 & crossAddedToStellarAverage();
+        const endf::Tab1 & crossAddedToStellarAverage();
     private:
-       std::string fileName;
+        std::string fileName;
 
-       int mat;
+        int mat;
        
-      /** The resonance data but not covariance data */
-      std::unique_ptr<endf::ResonanceInfo> resonanceData;
+        /** The resonance data but not covariance data */
+        std::unique_ptr<endf::ResonanceInfo> resonanceData;
 
-       /** All covariance data  */
-      std::unique_ptr<endf::CovarianceContainer> covarianceData;
+        /** All covariance data  */
+        std::unique_ptr<endf::CovarianceContainer> covarianceData;
 
-      /** All file 3 data */
-      std::unique_ptr<endf::Tab1Container> pointwiseData;
+        /** All file 3 data */
+        std::unique_ptr<endf::Tab1Container> pointwiseData;
 
-      /** File 1 information */
-      std::unique_ptr<endf::EvaluationInfo> evalInfo;
+        /** File 1 information */
+        std::unique_ptr<endf::EvaluationInfo> evalInfo;
 
-
-      /** Tab1 data to add in the case if 0K data if stellar averages are calculated */
-      endf::Tab1 tab1Addcross;
-  };
+        /** Tab1 data to add in the case if 0K data if stellar averages are calculated */
+        endf::Tab1 tab1Addcross;
+    };
 }
 
 #endif /* SAMMY_ENDF_ENDF_DATA_H */
diff --git a/sammy/src/fff/mfff0.f90 b/sammy/src/fff/mfff0.f90
index 7bcb4348f1cd434dd876a9bce8e5cb71e57ffa4c..223117fbf251728ab409ca9fa584bf84608690af 100644
--- a/sammy/src/fff/mfff0.f90
+++ b/sammy/src/fff/mfff0.f90
@@ -27,8 +27,11 @@ module fff_m
       use fff2_m
       use fff3_m
       use fff9_m
+      use lru_m
       use EndfData_common_m, only : radFitFlags
+
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
+      
       real(kind=8),allocatable,dimension(:)::A_Idum
       real(kind=8),allocatable,dimension(:)::A_Nkdnrm
       integer,allocatable,dimension(:)::I_Nkknrm
diff --git a/sammy/src/fff/mfff6.f90 b/sammy/src/fff/mfff6.f90
index f4296f951ea4da66152c863f84d5771793037104..62c23f9ff723b0150c7d821e6f16335ac3314b6e 100644
--- a/sammy/src/fff/mfff6.f90
+++ b/sammy/src/fff/mfff6.f90
@@ -43,10 +43,11 @@ module fff6_m
       IF (Numelv.NE.Numelv1) CALL Xwrong ('Numelv', Numelv, Numelv1)
       IF (Numjjv.NE.Numjjv1) CALL Xwrong ('Numjjv', Numjjv, Numjjv1)
 !      
-      READ (32) (Vrpr(N),N=1,Nnn1)     
+      READ (32) (Vrpr(N),N=1,Nnn1)
       
-      call covData%getUCovariance(uCov)    
-      call setCovMatrix(Nvpall1, uCov, Vrpr)    
+      call covData%getUCovariance(uCov)
+      call covData%resetTotalParam(Nvpall1) ! no PUPs in URR
+      call setCovMatrix(Nvpall1, uCov, Vrpr)
 !
       READ (32) (A,N=1,Numexc), (Dum(N),N=1,Numexc), (A,N=1,Numexc)
       DO N=1,Numexc
diff --git a/sammy/src/fff/mfff9.f90 b/sammy/src/fff/mfff9.f90
index 1fb99ea4ded824cbc39692ce62751d06a86c03f6..15296f07dcd91216d48478731b7be57c4adb5afb 100644
--- a/sammy/src/fff/mfff9.f90
+++ b/sammy/src/fff/mfff9.f90
@@ -39,15 +39,19 @@ module fff9_m
       Iu17 = 15
       IF (Kexptd.EQ.1) Iu17 = 17
       numDataSets = Kdtset
+!     Kpntws > 0 if we want an ENDF/B File 3 output
       IF (Kpntws.NE.0) numDataSets = Kpntws
       DO I=1,numDataSets
          ! ----------------------------------
          ! add grids to the gridDataList
          ! ----------------------------------
-         call expGrid%initialize()
-         dataInd = 2  ! 2 since angle-integrated and no time-shift (URR)
-         call expGrid%setDataIndex(dataInd)
-         call expData%addGrid(expGrid)
+!        if( we don't want pointwise xs for file 3 )
+         if(Kpntws.eq.0) then
+            call expGrid%initialize()
+            dataInd = 2  ! 2 since angle-integrated and no time-shift (URR)
+            call expGrid%setDataIndex(dataInd)
+            call expData%addGrid(expGrid)
+         end if
          ! ----------------------------------
          IF (Kexptd.EQ.1) THEN
             IF (Fdatax(I).EQ.Fblank) GO TO 30
@@ -97,6 +101,7 @@ module fff9_m
       IF (Emax.GE.Emaxxx) Emax = Emaxa
 !
       IF (Kpntws.EQ.0) RETURN
+!     ELSE we want ENDF/B File 3 output
       Igrid = 0
       DO I=1,Kpntws
          IF (Kktype(I).EQ.6) THEN
diff --git a/sammy/src/fit/mfit0.f b/sammy/src/fit/mfit0.f
index 94bd63a4b733f12fcca03ab8f4ec549d8534ecf5..0294519fbd32be65ba390635db199af70e482c08 100644
--- a/sammy/src/fit/mfit0.f
+++ b/sammy/src/fit/mfit0.f
@@ -17,6 +17,7 @@ C
       use AllocateFunctions_m
       use SammyLptPrinting_m
       use fit1
+      use fit3_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
       real(kind=8),allocatable,dimension(:)::A_Icorrx
       integer,allocatable,dimension(:)::I_Itabcx
diff --git a/sammy/src/fit/mfit3.f b/sammy/src/fit/mfit3.f90
similarity index 64%
rename from sammy/src/fit/mfit3.f
rename to sammy/src/fit/mfit3.f90
index 77c3b23445bb2297e1105ae21456976333976677..e4872063319a21ba9406df2374a8a2b1d0326c07 100644
--- a/sammy/src/fit/mfit3.f
+++ b/sammy/src/fit/mfit3.f90
@@ -1,14 +1,16 @@
-C
-C
-C ______________________________________________________________________
-C
-      SUBROUTINE Writep ( Anrm, Dnrm, Inrm, Kktype, Edirin, Sdirin,
-     *   Edircp, Sdircp, Streng, Dstren, Distnt, Ddistn, Gg, Dlgg,
-     *   Bethed, Gf, Fnu, Ethr, Wthr, Dlgf, Ex, Spin, Pty, Jnl, Jxl,
-     *   J2n, Bindee, Pairee, Engurr, Aold, Iflold, Ktold)
-C
-C *** Purpose -- Write parameters in the same format as parameter file
-C
+!
+module fit3_m
+   contains
+!
+! ______________________________________________________________________
+!
+      SUBROUTINE Writep ( Anrm, Dnrm, Inrm, Kktype, Edirin, Sdirin, &
+         Edircp, Sdircp, Streng, Dstren, Distnt, Ddistn, Gg, Dlgg, &
+         Bethed, Gf, Fnu, Ethr, Wthr, Dlgf, Ex, Spin, Pty, Jnl, Jxl, &
+         J2n, Bindee, Pairee, Engurr, Aold, Iflold, Ktold)
+!
+! *** Purpose -- Write parameters in the same format as parameter file
+!
       use fixedi_m
       use ifwrit_m
       use samxxx_common_m
@@ -19,40 +21,40 @@ C
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
       CHARACTER*8 Type(5)
       CHARACTER*40 Title(5)
-      DIMENSION Anrm(3,*), Dnrm(3,*), Inrm(*), Kktype(*),
-     * Edirin(*), Sdirin(*), Edircp(*), Sdircp(*), Streng(Numelv,*),
-     * Dstren(Numelv,*), Distnt(Numelv,*), Ddistn(Numelv,*),
-     * Gg(Numelv,*), Dlgg(Numelv,*), Bethed(Numelv,*),
-     * Gf(Numelv,Numjjv,*), Fnu(Numelv,*), Ethr(Numelv,Numjjv,*),
-     * Wthr(Numelv,Numjjv,*), Dlgf(Numelv,Numjjv,*), Ex(*), Spin(*),
-     * Pty(*), Jnl(*), Jxl(*), J2n(*), Bindee(*), Pairee(*), Engurr(*),
-     * Aold(3,*), Iflold(3,*), Ktold(*)
+      DIMENSION Anrm(3,*), Dnrm(3,*), Inrm(*), Kktype(*), &
+       Edirin(*), Sdirin(*), Edircp(*), Sdircp(*), Streng(Numelv,*), &
+       Dstren(Numelv,*), Distnt(Numelv,*), Ddistn(Numelv,*), &
+       Gg(Numelv,*), Dlgg(Numelv,*), Bethed(Numelv,*), &
+       Gf(Numelv,Numjjv,*), Fnu(Numelv,*), Ethr(Numelv,Numjjv,*), &
+       Wthr(Numelv,Numjjv,*), Dlgf(Numelv,Numjjv,*), Ex(*), Spin(*), &
+       Pty(*), Jnl(*), Jxl(*), J2n(*), Bindee(*), Pairee(*), Engurr(*), &
+       Aold(3,*), Iflold(3,*), Ktold(*)
       type(ResonanceCovariance)::uCov
       DIMENSION Ddddum(3)
-C
+!
       DATA Title(1) /'Covariance matrix from old run is used  '/
       DATA Title(2) /'Output file from SAMMY/FITACS run       '/
       DATA Title(3) /'                                        '/
       DATA Title(4) /'                                        '/
       DATA Title(5) /'----------------------------------------'/
-      DATA Type /'total   ', 'inelastc', 'fission ', 'capture ',
-     *   'elastic '/
-C
+      DATA Type /'total   ', 'inelastc', 'fission ', 'capture ', &
+         'elastic '/
+!
       DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/, Amillion /1.0d6/
-C
-C
+!
+!
       CALL Newopn (32, Samppp, 0)
       WRITE (32,10000) Title
 10000 FORMAT (A40)
-C
+!
       IF (Kannot.EQ.0 .AND. Kannox.EQ.0) THEN
          Zitmax = dFLOAT(Itmax+1)
          WRITE (32,10100) Zitmax, Tolera
 10100    FORMAT (F10.1, F10.6)
          Bex = Bindee(1)/Amillion
          Pex = Pairee(1)/Amillion
-         WRITE (32,10200) Aaawww,
-     *         radFitFlags%getMatchingRadius(), Bex, Pex
+         WRITE (32,10200) Aaawww, &
+               radFitFlags%getMatchingRadius(), Bex, Pex
 10200    FORMAT (F10.4, 3F10.6)
       ELSE
          Itmaxx = Itmax + 1
@@ -70,18 +72,18 @@ C
          WRITE (32,10240)
 10240    FORMAT ('Elastic and Inelastic States (with energy in eV)')
       END IF
-C
+!
       DO N=1,Numexc
          Exa = Ex(N)/(One+A_Mass_Small/Aaawww)
          WRITE (32,10300) Exa, Spin(N), Pty(N)
 10300    FORMAT (3F10.1)
       END DO
       WRITE (32,99999)
-C
-C
-C *** For each energy region
+!
+!
+! *** For each energy region
       DO Kumurr=1,Numurr
-C
+!
          IF (Kannot.EQ.1 .OR. Kannox.EQ.1) THEN
             Bex = Bindee(Kumurr)/Amillion
             Pex = Pairee(Kumurr)/Amillion
@@ -93,31 +95,31 @@ C
             WRITE (32,99999)
             WRITE (32,10220)
             WRITE (32,10350)
-10350       FORMAT ('Streng    Del_S', 8X,
-     *         'Distnt   Del_D   Gamma width Del_g      Bethed')
+10350       FORMAT ('Streng    Del_S', 8X, &
+               'Distnt   Del_D   Gamma width Del_g      Bethed')
          END IF
-C
-         WRITE (32,10400) Streng(1,Kumurr), Dstren(1,Kumurr),
-     *      Distnt(1,Kumurr), Ddistn(1,Kumurr),
-     *      GG(1,Kumurr), Dlgg(1,Kumurr), Bethed(1,Kumurr)
+!
+         WRITE (32,10400) Streng(1,Kumurr), Dstren(1,Kumurr), &
+            Distnt(1,Kumurr), Ddistn(1,Kumurr), &
+            GG(1,Kumurr), Dlgg(1,Kumurr), Bethed(1,Kumurr)
 10400    FORMAT (2F10.8, 2F10.7, 2F10.7, F10.4)
          IF (Numelv.GE.2) THEN
             DO L=2,Numelv
-               WRITE (32,10400) Streng(L,Kumurr), Dstren(L,Kumurr),
-     *            Distnt(L,Kumurr), Ddistn(L,Kumurr), GG(L,Kumurr),
-     *            Dlgg(L,Kumurr)
+               WRITE (32,10400) Streng(L,Kumurr), Dstren(L,Kumurr), &
+                  Distnt(L,Kumurr), Ddistn(L,Kumurr), GG(L,Kumurr), &
+                  Dlgg(L,Kumurr)
             END DO
          END IF
          WRITE (32,99999)
-C
+!
          IF (Kannot.EQ.1 .OR. Kannox.EQ.1) THEN
             WRITE (32,10220)
             WRITE (32,10450)
-10450       FORMAT ('Fission Width   Fnu    Ethr        Wthr', 4x,
-     *         'del_F_W         J    L')
+10450       FORMAT ('Fission Width   Fnu    Ethr        Wthr', 4x, &
+               'del_F_W         J    L')
          END IF
-C
-C ***    WRITE FISSION PARAMETERS
+!
+! ***    WRITE FISSION PARAMETERS
          DO L=1,Numelv
             Jlo = Jnl(L)
             Jhi = Jxl(L)
@@ -125,16 +127,16 @@ C ***    WRITE FISSION PARAMETERS
             DO J=Jlo,Jhi
                IF (Gf(L,J,Kumurr).EQ.Zero) GO TO 100
                Lll = L - 1
-               WRITE (32,10500) Gf(L,J,Kumurr), Fnu(L,J),
-     *            Ethr(L,J,Kumurr), Wthr(L,J,Kumurr), Dlgf(L,J,Kumurr),
-     *            Ajj, L-1
+               WRITE (32,10500) Gf(L,J,Kumurr), Fnu(L,J), &
+                  Ethr(L,J,Kumurr), Wthr(L,J,Kumurr), Dlgf(L,J,Kumurr), &
+                  Ajj, L-1
 10500          FORMAT (F10.7, F10.2, 2F10.0, F10.7, F10.1, I5)
                Ajj = Ajj + One
             END DO
          END DO
   100    CONTINUE
          WRITE (32,99999)
-C
+!
          IF (Kannot.EQ.1 .OR. Kannox.EQ.1) THEN
             WRITE (32,10220)
             IF (Kumurr.EQ.1) THEN
@@ -147,17 +149,17 @@ C
             WRITE (32,10310) Engurr(Kumurr)
 10310       FORMAT ('Energy Maximum in eV =', 1PG20.10)
          END IF
-C
+!
       END DO
-C *** End of energy-regions
-C
+! *** End of energy-regions
+!
       IF (Kannot.EQ.1 .OR. Kannox.EQ.1) THEN
          WRITE (32,10000) Title(4)
          WRITE (32,10550)
 10550    FORMAT ('End of resonance parameter description')
          WRITE (32,10000) Title(4)
       END IF
-C
+!
       IF (Kdtset.GT.0) THEN
          WRITE (32,10600)
 10600    FORMAT ('Earlier Normalizations')
@@ -181,25 +183,25 @@ C
                      Ddddum(J) = dSQRT(val)
                   END IF
                END DO
-               WRITE (32,10750) Iset, Type(Ktold(Iset)), 
-     *            (Aold(J,Iset),Ddddum(J),J=1,3)
+               WRITE (32,10750) Iset, Type(Ktold(Iset)),  &
+                  (Aold(J,Iset),Ddddum(J),J=1,3)
             END DO
          END IF
          IF (Inrm(1).EQ.0) THEN
             DO Iset=1,Kdtset
                Is = Iset + Kdtold
-               WRITE (32,10750) Is, Type(Kktype(Iset)),
-     *            (Anrm(J,Iset),Dnrm(J,Iset),J=1,3)
+               WRITE (32,10750) Is, Type(Kktype(Iset)), &
+                  (Anrm(J,Iset),Dnrm(J,Iset),J=1,3)
             END DO
          END IF
 10750    FORMAT ('# for data set number', I3, /, '#', A8, 2X, 6F10.6)
          WRITE (32,10760)
 10760    FORMAT ('#')
       END IF
-C
+!
       IF (Ndirin.GT.0) THEN
-10800 FORMAT ('Direct Inelastic contribution (energy in eV, cross section
-     *n in barns)')
+10800 FORMAT ('Direct Inelastic contribution (energy in eV, cross section &
+      n in barns)')
          WRITE (32,10800)
          DO I=1,Ndirin
             WRITE (32,10900) Edirin(I), Sdirin(I)
@@ -207,33 +209,33 @@ C
 10900    FORMAT ('E=', F20.8, '   Sigma=', F15.10)
          WRITE (32,99999)
       END IF
-C
+!
       IF (Ndircp.GT.0) THEN
-11000 FORMAT ('Direct Capture contribution (energy in eV, cross section i
-     *in barns)')
+11000 FORMAT ('Direct Capture contribution (energy in eV, cross section i &
+      in barns)')
          WRITE (32,11000)
          DO I=1,Ndircp
             WRITE (32,10900) Edircp(I), Sdircp(I)
          END DO
          WRITE (32,99999)
       END IF
-C
+!
       CLOSE (UNIT=32)
-C
+!
       RETURN
 99999 FORMAT (' ')
       END
-C
-C
-C ______________________________________________________________________
-C
-      SUBROUTINE Write_C ( Anrm, Inrm, Kkktyp, Streng, Distnt, Gg,
-     *   Bethed, Gf, Fnu, Ethr, Wthr, Ex, Spin, Pty, Jnl, Jxl,
-     *   Iflstr, Ifldst, Iflggg, Iflgff, Iflnrm, Bindee, Pairee,
-     *   Engurr, Aold, Iflold, Ktold)
-C
-C *** Purpose -- Create URR covariance file
-C
+!
+!
+! ______________________________________________________________________
+!
+      SUBROUTINE Write_C ( Anrm, Inrm, Kkktyp, Streng, Distnt, Gg, &
+         Bethed, Gf, Fnu, Ethr, Wthr, Ex, Spin, Pty, Jnl, Jxl, &
+         Iflstr, Ifldst, Iflggg, Iflgff, Iflnrm, Bindee, Pairee, &
+         Engurr, Aold, Iflold, Ktold)
+!
+! *** Purpose -- Create URR covariance file
+!
       use fixedi_m
       use ifwrit_m
       use samxxx_common_m
@@ -241,27 +243,29 @@ C
       use constn_common_m
       use EndfData_common_m
       use ResonanceCovariance_M
+      use Hdf5IO_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      DIMENSION Anrm(3,*), Inrm(*), Kkktyp(*), Streng(Numelv,*),
-     *   Distnt(Numelv,*), Gg(Numelv,*), Bethed(Numelv,*),
-     *   Gf(Numelv,Numjjv,*), Fnu(Numelv,*), Ethr(Numelv,Numjjv,*),
-     *   Wthr(Numelv,Numjjv,*), Ex(*), Spin(*), Pty(*), Jnl(*), Jxl(*),
-     *   Iflstr(Numelv,*), Ifldst(Numelv,*), Iflggg(Numelv,*),
-     *   Iflgff(Numelv,Numjjv,*), Iflnrm(3,*), Bindee(*), Pairee(*),
-     *   Engurr(*), Aold(3,*), Iflold(3,*), Ktold(*)
+      DIMENSION Anrm(3,*), Inrm(*), Kkktyp(*), Streng(Numelv,*), &
+         Distnt(Numelv,*), Gg(Numelv,*), Bethed(Numelv,*), &
+         Gf(Numelv,Numjjv,*), Fnu(Numelv,*), Ethr(Numelv,Numjjv,*), &
+         Wthr(Numelv,Numjjv,*), Ex(*), Spin(*), Pty(*), Jnl(*), Jxl(*), &
+         Iflstr(Numelv,*), Ifldst(Numelv,*), Iflggg(Numelv,*), &
+         Iflgff(Numelv,Numjjv,*), Iflnrm(3,*), Bindee(*), Pairee(*), &
+         Engurr(*), Aold(3,*), Iflold(3,*), Ktold(*)
       type(ResonanceCovariance)::uCov, invCov
-C
-C
+      type(Hdf5IO)::h5
+!
+!
       CALL Newopn (62, Samcov, 1)
       Nnn = (Nvpall*(Nvpall+1))/2
       K = Kdtold
       IF (Inrm(1).EQ.0) K = K + Kdtset
       WRITE (62) Nvpall, Nnn, Nvpthe, Numexc, Numurr, Numelv, Numjjv, K
       call covData%getUCovariance(uCov)
-      WRITE (62) ((uCov%getCovariance(ipar, i), i = 1, ipar),
-     *            ipar=1,Nvpall)
-      WRITE (62) (Ex(N),N=1,Numexc), (Spin(N),N=1,Numexc),
-     *          (Pty(N),N=1,Numexc)
+      WRITE (62) ((uCov%getCovariance(ipar, i), i = 1, ipar), &
+                  ipar=1,Nvpall)
+      WRITE (62) (Ex(N),N=1,Numexc), (Spin(N),N=1,Numexc), &
+                (Pty(N),N=1,Numexc)
       WRITE (62) (Bindee(N),N=1,Numurr)
       WRITE (62) (Pairee(N),N=1,Numurr)
       WRITE (62) (Engurr(N),N=1,Numurr)
@@ -279,31 +283,40 @@ C
       WRITE (62) (((Wthr  (L,J,N),L=1,Numelv),J=1,Numjjv),N=1,Numurr)
       WRITE (62) (((Iflgff(L,J,N),L=1,Numelv),J=1,Numjjv),N=1,Numurr)
       IF (Kdtold.EQ.0) THEN
-         WRITE (62) ((Anrm(J,Iset),J=1,3),Iset=1,Kdtset),
-     *            ((Iflnrm(J,Iset),J=1,3),Iset=1,Kdtset),
-     *             (Kkktyp(  Iset),       Iset=1,Kdtset)
+         WRITE (62) ((Anrm(J,Iset),J=1,3),Iset=1,Kdtset), &
+                  ((Iflnrm(J,Iset),J=1,3),Iset=1,Kdtset), &
+                   (Kkktyp(  Iset),       Iset=1,Kdtset)
       ELSE IF (Inrm(1).EQ.0) THEN
-         WRITE (62) ((Aold(J,Iset),J=1,3),Iset=1,Kdtold),
-     *            ((  Anrm(J,Iset),J=1,3),Iset=1,Kdtset),
-     *            ((Iflold(J,Iset),J=1,3),Iset=1,Kdtold),
-     *            ((Iflnrm(J,Iset),J=1,3),Iset=1,Kdtset),
-     *             ( Ktold(  Iset),       Iset=1,Kdtold),
-     *             (Kkktyp(  Iset),       Iset=1,Kdtset)
+         WRITE (62) ((Aold(J,Iset),J=1,3),Iset=1,Kdtold), &
+                  ((  Anrm(J,Iset),J=1,3),Iset=1,Kdtset), &
+                  ((Iflold(J,Iset),J=1,3),Iset=1,Kdtold), &
+                  ((Iflnrm(J,Iset),J=1,3),Iset=1,Kdtset), &
+                   ( Ktold(  Iset),       Iset=1,Kdtold), &
+                   (Kkktyp(  Iset),       Iset=1,Kdtset)
       ELSE
-         WRITE (62)(( Aold(J,Iset),J=1,3),Iset=1,Kdtold),
-     *            ((Iflold(J,Iset),J=1,3),Iset=1,Kdtold),
-     *             ( Ktold(  Iset),       Iset=1,Kdtold)
+         WRITE (62)(( Aold(J,Iset),J=1,3),Iset=1,Kdtold), &
+                  ((Iflold(J,Iset),J=1,3),Iset=1,Kdtold), &
+                   ( Ktold(  Iset),       Iset=1,Kdtold)
       END IF
-C
-C *** Read and write the parameter values U    
-      WRITE (62) (covData%getUParamValue(I),
-     *              i = 1, covData%getNumParam())
+!
+! *** Read and write the parameter values U    
+      WRITE (62) (covData%getUParamValue(I), &
+                    i = 1, covData%getNumParam())
   
       call covData%getInverseUCovariance(invCov)     
-      WRITE (62) ( (invCov%getCovariance(ipar, i), i = 1,ipar), 
-     *              ipar=1,Nvpall)  
+      WRITE (62) ( (invCov%getCovariance(ipar, i), i = 1,ipar),  &
+                    ipar=1,Nvpall)  
       
-C     
+!     
       CLOSE (UNIT=62)
+
+      ! -----------------------------------------
+      ! Write simple res param covariance matrix
+      ! to HDF5 file (u-param & u-covar)
+      ! -----------------------------------------
+      call h5%writeRPCM(covData,trim(Samcov)//".h5")
+
       RETURN
       END
+
+end module fit3_m
diff --git a/sammy/src/io/CMakeLists.txt b/sammy/src/io/CMakeLists.txt
index d5af49a65b26d7d5553f8d94f530ca1a67561397..2fc423405c8f8f11a0c40ff0d5f16d9a41b95cb7 100644
--- a/sammy/src/io/CMakeLists.txt
+++ b/sammy/src/io/CMakeLists.txt
@@ -10,20 +10,30 @@ INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR}
                     ${HDF5_INCLUDE_DIRS}) # from find_package()
 
 # Set headers
-SET(HEADERS ExperimentalParametersIO.h
-            OdfIO.h
+SET(HEADERS CoverxIO.h
+            interface/cpp/CoverxIOInterface.h
+            ExperimentalParametersIO.h
+            interface/cpp/ExperimentalParametersIOInterface.h
             Hdf5IO.h
-            interface/cpp/ExperimentalParametersIOInterface.h 
+            interface/cpp/Hdf5IOInterface.h
+            OdfIO.h
 )
 
 APPEND_SET(IO_SOURCES
 
+        CoverxIO.cpp
+        interface/cpp/CoverxIOInterface.cpp
+        interface/fortran/CoverxIO_M.f90
+        interface/fortran/CoverxIO_I.f90
         ExperimentalParametersIO.cpp
-        OdfIO.cpp
-        Hdf5IO.cpp
         interface/cpp/ExperimentalParametersIOInterface.cpp
-        interface/fortran/ExperimentalParametersIO_M.f90   
-        interface/fortran/ExperimentalParametersIO_I.f90      
+        interface/fortran/ExperimentalParametersIO_M.f90
+        interface/fortran/ExperimentalParametersIO_I.f90
+        Hdf5IO.cpp
+        interface/cpp/Hdf5IOInterface.cpp
+        interface/fortran/Hdf5IO_M.f90
+        interface/fortran/Hdf5IO_I.f90
+        OdfIO.cpp
 )
 
 
diff --git a/sammy/src/io/CoverxIO.cpp b/sammy/src/io/CoverxIO.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eb3bf9ad74c67e51b1b362e3be7c8f4cf7ae665a
--- /dev/null
+++ b/sammy/src/io/CoverxIO.cpp
@@ -0,0 +1,143 @@
+#include "CoverxIO.h"
+
+
+namespace sammy{
+
+    CoverxIO::CoverxIO() {}
+
+    void CoverxIO::setCoverxFile(GridDataList & gdList, CovarianceData & cov)
+    {
+        // variable setup
+        int numGrids = gdList.getLength();
+        int numPar = cov.getNumTotalParam();
+        std::cerr << "Num par = " << numPar << std::endl;
+        endf::ResonanceCovariance * parCov = cov.getUCovariance();
+
+
+
+        covFile.setDescription("This Coverx file generated in SAMMY.");
+        covFile.setHname("coverx");
+
+        if( matmt.size() < numGrids ){
+            std::cerr << "---------------------------------" << std::endl;
+            std::cerr << "Warning: the number of MAT/MT IDs" << std::endl
+                      << "is less than the number of grids " << std::endl
+                      << "being written to CoverxFile.     " << std::endl
+                      << "MAT/MT IDs should match ordering " << std::endl
+                      << "of the GridDataList.             " << std::endl;
+            std::cerr << "---------------------------------" << std::endl;
+            throw std::runtime_error("Please set the MAT/MT numbers.");
+        }
+
+        // cycle through the grid list
+        for( int i=0; i<numGrids; ++i ){
+            // set up input info
+            bool mt1NonZero = false;
+            int numEnergy = gdList.getGrid(i)->getLength();
+            int dataColI = gdList.getGrid(i)->getDataColumn();
+            int derCol1I = dataColI+1;
+
+            // set up numGrids many cross sections
+            coverx::CrossSectionData * cs = new coverx::CrossSectionData(numEnergy);
+            cs->changeRelative(false); // absolute error on cs will be given
+            cs->setMat(matmt[i].first); 
+            cs->setMt(matmt[i].second);
+
+            if( i==0 ){
+                covFile.setNeutGroups(numEnergy); // add 1 for upper/lower bounds
+                // check if data column == 1, if so we should have energy info
+                if( dataColI < 1 ){
+                    std::cerr << "Uh oh, data column is equal to: " << dataColI
+                              << std::endl;
+                    throw std::runtime_error("No energy data found in grid # 1!");
+                }
+            }
+            for( int j=0; j<numGrids; ++j ){
+                bool mt2NonZero = false;
+                if( numEnergy != gdList.getGrid(j)->getLength() ){
+                    throw std::runtime_error("Energy grids don't match!");
+                }
+                int dataColJ = gdList.getGrid(j)->getDataColumn();
+                int derCol1J = dataColJ+1;
+
+                // set up output info
+                coverx::CovMatrix * csCov = new coverx::CovMatrix(2) ; // type = 2 = absolute covar
+                csCov->setSize(numEnergy);
+                csCov->setMatValues(matmt[i].first,matmt[i].second,
+                                    matmt[j].first,matmt[j].second);
+
+                // matrix multiplication loop -------------------------------------
+                for( int iE=0; iE<numEnergy; ++iE ){
+                    if( i==j ){
+                        double xs = gdList.getGrid(i)->getData(iE,dataColI);
+                        if( xs > 0.0 ) mt1NonZero = true;
+                        cs->setCrossAt(iE,xs);
+                        if( i==0 ){
+                            float E = (float)gdList.getGrid(i)->getData(iE,dataColI-1);
+                            covFile.setNeutGroupBounds(iE,E);
+                        }
+                    }
+                    for( int jE=0; jE<numEnergy; ++jE ){
+                        double val = 0.0;
+                        double xs = gdList.getGrid(j)->getData(jE,dataColI);
+                        if( xs > 0.0 ) mt2NonZero = true;
+                        for( int iPar=0; iPar<numPar; ++iPar ){
+                            for( int jPar=0; jPar<numPar; ++jPar ){
+                                val += gdList.getGrid(i)->getData(iE,iPar+derCol1I) * 
+                                       gdList.getGrid(j)->getData(jE,jPar+derCol1J) *
+                                       parCov->getCovariance(iPar,jPar);
+                            }
+                        }
+                        csCov->setValue(iE,jE,val);
+                        if( i==j && iE==jE && val>0.0 ){
+                            // default is zero
+                            cs->setErrorAt( iE,std::sqrt(val) );
+                        }
+                    }
+                }
+                if( i==j && i==0 ){
+                    double maxE = covFile.getNeutGroupBounds(numEnergy-1);
+                    double plus1Perc = maxE + maxE*0.01;
+                    covFile.setNeutGroupBounds(numEnergy,plus1Perc);
+                }
+                if( mt1NonZero && mt2NonZero ) covFile.addMatrix( csCov );
+                // ----------------------------------------------------------------
+            }
+            if( mt1NonZero ) covFile.addCross( cs );
+        } // end grid loops
+    }
+
+    void CoverxIO::writeCoverxFile(const std::string & filename)
+    {
+        coverx::CoverxWriter cw;
+        int type = 1; // types should be relative
+        bool swap = false; // should we swap little/big endian?
+        cw.writeCoverxFile( filename.c_str(), covFile, swap, type);
+    }
+
+    void CoverxIO::readCoverxFile(const std::string & filename)
+    {
+        coverx::CoverxReader cr;
+        bool swap = false; // swap little/big endian?
+        cr.readCoverxFile(filename.c_str(),covFile,swap);
+    }
+
+    void CoverxIO::writeFile33(double za, double awr)
+    {
+        if( covFile.getNumCross() < 1 ){
+            std::string error;
+            error.append("No cross sections found in CoverxFile.\n");
+            error.append("I'm guessing you never set them.\n");
+            throw std::runtime_error(error);
+        }
+
+        CovToolBox tools;
+        endfendf::File33 f33;
+        std::ofstream outFile("sammy.f33");
+        tools.fillCrossCovFromCoverx(covFile,ccList);
+        f33.writeData(outFile,ccList,za,awr);
+    }
+}
+
+
+
diff --git a/sammy/src/io/CoverxIO.h b/sammy/src/io/CoverxIO.h
new file mode 100644
index 0000000000000000000000000000000000000000..bfc2ff9ef9dc2eb6fe7c2be665016c0f55c0ddd0
--- /dev/null
+++ b/sammy/src/io/CoverxIO.h
@@ -0,0 +1,101 @@
+#ifndef SAMMY_COVERXIO_H
+#define SAMMY_COVERXIO_H
+
+#include <vector>
+#include <fstream> // ifstream
+#include <iostream>
+#include <map>
+#include <cmath>
+#include "ScaleData/Core/io/CoverxFile.h"
+#include "ScaleData/Core/io/CoverxWriter.h"
+#include "ScaleData/Core/io/CoverxReader.h"
+#include "ScaleUtils/EndfLib/CovarianceContainer.h"
+#include "../salmon/CovToolBox.h"
+#include "ScaleUtils/EndfLib/endf/File33.h"
+#include "../salmon/GridData.h"
+#include "../endf/CovarianceData.h"
+
+
+namespace sammy{
+
+    /*****************************************************************
+     * The CoverxIO class provides methods to deal with input and 
+     * output for the Coverx file format, which contains cross section
+     * covariance information.
+     ****************************************************************/
+    class CoverxIO {
+    
+    public:
+        CoverxIO();
+        virtual ~CoverxIO() {}
+
+        /***********************************************************************
+         * Set all the information up for the Coverx file-writing. Only non-zero
+         * cross sections will be written, and have associated matrices written.
+         * 
+         * @param gdList a list of GridData's containing XS and derivatives
+         * @param parCov the resonance u-covariance and u-parameters
+         **********************************************************************/
+        void setCoverxFile(GridDataList & gdList, CovarianceData & parCov);
+
+        /***********************************************************************
+         * Set the MAT/MT material/reaction ID numbers, one at a time. 
+         * These should be in the same order as the `gdList` input 
+         * GridDataList parameter given in `setCoverxFile()`
+         * 
+         * @param mat the material ID number
+         * @param mt the reaction ID number
+         **********************************************************************/
+        void setMatMt(int mat,int mt){
+            std::pair<int,int> temp = {mat,mt};
+            matmt.push_back(temp);
+        }
+
+        /***********************************************************************
+         * Write the CoverxFile info to an actual file on disk
+         * !! WARNING !!: this function changes the values in the `CoverxFile`
+         *                object to "relative" values (correlations)
+         * 
+         * @param filename the name of the file
+         **********************************************************************/
+        void writeCoverxFile(const std::string & filename);
+
+        /***********************************************************************
+         * Read the Coverx formatted file into the CoverxFile object
+         * 
+         * @param filename the name of the file
+         **********************************************************************/
+        void readCoverxFile(const std::string & filename);
+
+        /***********************************************************************
+         * Write the pointwise cross section covariance in ENDF/B File 33 format
+         * 
+         * @param za the ENDF/B ZA parameter: (protons*1000)+(neutrons+protons)
+         * @param awr the mass of the nucleus divide by the neutron mass
+         **********************************************************************/
+        void writeFile33(double za, double awr);
+
+        const coverx::CovMatrix & getMatrix(int mat1, int mt1, int mat2, int mt2){
+            const coverx::CovMatrix * cm = covFile.getMatrix(mat1,mt1,mat2,mt2);
+            if( cm == nullptr ) throw std::runtime_error("matrix pointer null");
+            return *cm;
+        }
+        const coverx::CrossSectionData & getCross(int mat, int mt){
+            const coverx::CrossSectionData * cs = covFile.getCross(mat,mt);
+            if( cs == nullptr ) throw std::runtime_error("cross section pointer null");
+            return *cs;
+        }
+
+    private:
+        // The object containing info necessary to write file
+        coverx::CoverxFile covFile;
+        // list of MAT/MT material/reaction ID numbers corresponding to
+        // grids in `gdList`
+        std::vector<std::pair<int,int>> matmt;
+
+        endf::CrossCovList ccList;
+        
+    };
+}
+
+#endif /* SAMMY_COVERXIO_H */
\ No newline at end of file
diff --git a/sammy/src/io/Hdf5IO.cpp b/sammy/src/io/Hdf5IO.cpp
index b81cc466aee0927a331f377c9230ffe978e4bb1e..baac8ff0699cb55b6a8d990f0f71c196215b8821 100644
--- a/sammy/src/io/Hdf5IO.cpp
+++ b/sammy/src/io/Hdf5IO.cpp
@@ -1,4 +1,5 @@
 #include "Hdf5IO.h"
+#include <cmath>
 
 namespace sammy {
     
@@ -101,4 +102,95 @@ namespace sammy {
         }
     }
 
+    void Hdf5IO::writeRPCM(CovarianceData & cov,std::string hdf5Filename){
+        std::cout << "--------------------------------------------" << std::endl;
+        std::cout << "writing resonance covariance to HDF5 file..." << std::endl;
+        std::cout << "--------------------------------------------" << std::endl;
+
+        // get parameter covariance out
+        endf::ResonanceCovariance * parCov = cov.getUCovariance();
+        int numParam = cov.getNumTotalParam();
+        int numRows = parCov->getNumRows();
+        if( numParam != numRows ) {
+            std::cerr << "numParam != numRows of covariance!" << std::endl;
+        }
+
+        // get parameters out
+        std::vector<double> param;
+        for( int i=0; i<numParam; ++i ) param.push_back(cov.getUParamValue(i));
+
+        double data[numRows][numRows];
+        for( int i=0; i<numRows; ++i ){
+            for( int j=0; j<numRows; ++j ){
+                data[i][j] = parCov->getCovariance(i,j);
+            }
+        }
+
+        // ----------------------------------------------------------
+        // write 2 datasets: u-parameters and u-covariance
+        // ----------------------------------------------------------
+
+        // --- covar space -----
+        hsize_t cdim[2];
+        cdim[0] = numRows;
+        cdim[1] = cdim[0];
+        int crank = sizeof(cdim) / sizeof(hsize_t);
+        std::string covarname("/uCovar/");
+
+        // --- param space -----
+        hsize_t pdim[1];
+        pdim[0] = numParam;
+        int prank = sizeof(pdim) / sizeof(hsize_t);
+        std::string paramname("/uParam/");
+
+        H5::DataSpace cspace(crank,cdim);
+        H5::DataSpace pspace(prank,pdim);
+        H5::H5File  file = H5::H5File(hdf5Filename,H5F_ACC_TRUNC);
+        H5::DataSet covdataset = H5::DataSet(file.createDataSet(covarname,H5::PredType::NATIVE_DOUBLE,cspace));
+        H5::DataSet pardataset = H5::DataSet(file.createDataSet(paramname,H5::PredType::NATIVE_DOUBLE,pspace));
+        // write param covariance to HDF5 file
+        covdataset.write( data,         H5::PredType::NATIVE_DOUBLE );
+        // write parameters to HDF5 file
+        pardataset.write( param.data(), H5::PredType::NATIVE_DOUBLE );
+
+    }
+
+    void Hdf5IO::readRPCM(CovarianceData & cov,std::string hdf5Filename){
+        // Keep in mind that this depends on reading the file we wrote
+        // with the function Hdf5IO::writeRPCM(). Otherwise we'd have 
+        // to pass more information into the function.
+
+        try{
+            H5::H5File file(hdf5Filename, H5F_ACC_RDONLY);
+            H5::DataSet h5covSet = file.openDataSet( "uCovar" );
+            H5::DataSet h5parSet = file.openDataSet( "uParam" );
+            H5::DataSpace h5covSpace = h5covSet.getSpace();
+            H5::DataSpace h5parSpace = h5parSet.getSpace();
+            // --- get dimensions of data to size local arrays
+            hsize_t cov_dims[2];
+            hsize_t par_dims[1];
+            int ndims = h5covSpace.getSimpleExtentDims( cov_dims, NULL);
+            ndims     = h5parSpace.getSimpleExtentDims( par_dims, NULL);
+            double h5cov[cov_dims[0]][cov_dims[1]];
+            double h5par[par_dims[0]];
+            // --- read data into local arrays
+            h5covSet.read(h5cov,H5::PredType::NATIVE_DOUBLE);
+            h5parSet.read(h5par,H5::PredType::NATIVE_DOUBLE);
+            // --- load into CovarianceData object
+            cov.makeNewUCovariance();
+            endf::ResonanceCovariance * parCov = cov.getUCovariance();
+            cov.clearUParameterValues();
+            for( int i=0; i<cov_dims[0]; ++i ){
+                cov.setUParamValue( i,h5par[i] );
+                for( int j=0; j<cov_dims[1]; ++j ){
+                    parCov->setCovariance( i,j,h5cov[i][j] );
+                }
+            }
+        }
+        catch(...)
+        {
+           throw std::runtime_error("HDF5 file reading failed in io/Hdf5IO.cpp");
+        }
+    }
+
 }
\ No newline at end of file
diff --git a/sammy/src/io/Hdf5IO.h b/sammy/src/io/Hdf5IO.h
index dd073a0f4a1b8785c94721d30ff260a961390bd0..396439e51f05215beccc63a9f55c1d7b0b5ed76d 100644
--- a/sammy/src/io/Hdf5IO.h
+++ b/sammy/src/io/Hdf5IO.h
@@ -8,6 +8,7 @@
 #include <stdexcept>
 #include "OdfIO.h"
 #include "H5Cpp.h"
+#include "../endf/CovarianceData.h"
 
 namespace sammy{
 
@@ -52,6 +53,24 @@ namespace sammy{
          ************************************************************/
         void fillHDF5header(std::map<std::string,int> &odfhead);
 
+        /*************************************************************
+         * Write the resonance parameter covariance matrix to an HDF5
+         * file
+         * 
+         * @param cov class that contains all param. covariance info
+         * @param hdf5Filename the string specifying the H5 file name
+         ************************************************************/
+        void writeRPCM(CovarianceData & cov,std::string hdf5Filename);
+
+        /*************************************************************
+         * Read the resonance parameter covariance matrix from an HDF5
+         * file which was written by Hdf5IO::writeRPCM()
+         * 
+         * @param cov class that we fill with param. covariance info
+         * @param hdf5filename the string specifying the H5 file name
+         ************************************************************/
+        void readRPCM(CovarianceData & cov,std::string hdf5Filename);
+
     private:
         std::vector<headKeypair> h5head;
     };
diff --git a/sammy/src/io/interface/cix/CoverxIO.cpp2f.xml b/sammy/src/io/interface/cix/CoverxIO.cpp2f.xml
new file mode 100644
index 0000000000000000000000000000000000000000..f0d626d991334fb6a35fbe7421368b1f84a5b8cb
--- /dev/null
+++ b/sammy/src/io/interface/cix/CoverxIO.cpp2f.xml
@@ -0,0 +1,23 @@
+<generate name="CoverxIO">
+    <include_relative name="../../CoverxIO.h"/>
+    <include_relative name="../../../salmon/GridData.h"/>
+    <using_namespace name="sammy"/>
+
+    <class name="CoverxIO">
+      <method name="setCoverxFile">
+         <param name="gdList" type="GridDataList"/>
+         <param name="parCov" type="CovarianceData"/>
+      </method>
+      <method name="setMatMt">
+        <param name="mat" type="int"/>
+        <param name="mt" type="int"/>
+      </method>
+      <method name="writeCoverxFile">
+        <param name="filename" type="char*" ctype="std::string"/>
+      </method>
+      <method name="writeFile33">
+        <param name="za" type="double"/>
+        <param name="awr" type="double"/>
+      </method>
+    </class>
+</generate>
\ No newline at end of file
diff --git a/sammy/src/io/interface/cix/Hdf5IO.cpp2f.xml b/sammy/src/io/interface/cix/Hdf5IO.cpp2f.xml
new file mode 100644
index 0000000000000000000000000000000000000000..87b2cf43ebe3681180867899b8f421e11904a5a5
--- /dev/null
+++ b/sammy/src/io/interface/cix/Hdf5IO.cpp2f.xml
@@ -0,0 +1,15 @@
+<generate name="Hdf5IO">
+    <include_relative name="../../Hdf5IO.h"/>
+    <using_namespace name="sammy"/>
+
+    <class name="Hdf5IO">
+      <method name="writeRPCM">
+         <param name="cov" type="CovarianceData"/>
+         <param name="hdf5Filename" type="char*" ctype="std::string"/>
+      </method>
+      <method name="readRPCM">
+         <param name="cov" type="CovarianceData"/>
+         <param name="hdf5Filename" type="char*" ctype="std::string"/>
+      </method>
+    </class>
+</generate>
\ No newline at end of file
diff --git a/sammy/src/io/interface/cpp/CoverxIOInterface.cpp b/sammy/src/io/interface/cpp/CoverxIOInterface.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d272a66e683aa39e7cf377da254709442db6480d
--- /dev/null
+++ b/sammy/src/io/interface/cpp/CoverxIOInterface.cpp
@@ -0,0 +1,41 @@
+/*!
+* This file has been dynamically generated by Class Interface Xml (CIX) 
+* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+* If changes need to occur, modify the appropriate CIX xml file
+* Date Generated: Fri Feb 25 11:18:44 EST 2022
+* If any issues are experiences with this generated file that cannot be fixed
+* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+*/
+#include <string.h>
+#include "CoverxIOInterface.h"
+using namespace sammy;
+void CoverxIO_setCoverxFile(void * CoverxIO_ptr,GridDataList * gdList,CovarianceData * parCov)
+{
+    ((CoverxIO*)CoverxIO_ptr)->setCoverxFile(*gdList,*parCov);
+}
+
+void CoverxIO_setMatMt(void * CoverxIO_ptr,int * mat,int * mt)
+{
+    ((CoverxIO*)CoverxIO_ptr)->setMatMt(*mat,*mt);
+}
+
+void CoverxIO_writeCoverxFile(void * CoverxIO_ptr,char* filename)
+{
+    ((CoverxIO*)CoverxIO_ptr)->writeCoverxFile( std::string(filename));
+}
+
+void CoverxIO_writeFile33(void * CoverxIO_ptr,double * za,double * awr)
+{
+    ((CoverxIO*)CoverxIO_ptr)->writeFile33(*za,*awr);
+}
+
+void* CoverxIO_initialize()
+{
+    return new CoverxIO();
+}
+
+void CoverxIO_destroy(void * CoverxIO_ptr)
+{
+    delete (CoverxIO*)CoverxIO_ptr;
+}
+
diff --git a/sammy/src/io/interface/cpp/CoverxIOInterface.h b/sammy/src/io/interface/cpp/CoverxIOInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..9954e7fc77e0a52cc4d2544d1622e6cb9e595c00
--- /dev/null
+++ b/sammy/src/io/interface/cpp/CoverxIOInterface.h
@@ -0,0 +1,26 @@
+/*!
+* This file has been dynamically generated by Class Interface Xml (CIX) 
+* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+* If changes need to occur, modify the appropriate CIX xml file
+* Date Generated: Fri Feb 25 11:18:44 EST 2022
+* If any issues are experiences with this generated file that cannot be fixed
+* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+*/
+#ifndef COVERXIOINTERFACE_H
+#define COVERXIOINTERFACE_H
+#include "../../CoverxIO.h"
+#include "../../../salmon/GridData.h"
+using namespace sammy;
+#ifdef __cplusplus
+extern "C" {
+#endif
+void CoverxIO_setCoverxFile(void * CoverxIO_ptr,GridDataList * gdList,CovarianceData * parCov);
+void CoverxIO_setMatMt(void * CoverxIO_ptr,int * mat,int * mt);
+void CoverxIO_writeCoverxFile(void * CoverxIO_ptr,char* filename);
+void CoverxIO_writeFile33(void * CoverxIO_ptr,double * za,double * awr);
+void* CoverxIO_initialize();
+void CoverxIO_destroy(void * CoverxIO_ptr);
+#ifdef __cplusplus
+}
+#endif
+#endif /* COVERXIOINTERFACE_H */
diff --git a/sammy/src/io/interface/cpp/Hdf5IOInterface.cpp b/sammy/src/io/interface/cpp/Hdf5IOInterface.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..419d9310b3cf55cdc2f98bf4b52bdac87843277d
--- /dev/null
+++ b/sammy/src/io/interface/cpp/Hdf5IOInterface.cpp
@@ -0,0 +1,31 @@
+/*!
+* This file has been dynamically generated by Class Interface Xml (CIX) 
+* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+* If changes need to occur, modify the appropriate CIX xml file
+* Date Generated: Wed Jan 26 13:57:17 EST 2022
+* If any issues are experiences with this generated file that cannot be fixed
+* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+*/
+#include <string.h>
+#include "Hdf5IOInterface.h"
+using namespace sammy;
+void Hdf5IO_writeRPCM(void * Hdf5IO_ptr,CovarianceData * cov,char* hdf5Filename)
+{
+    ((Hdf5IO*)Hdf5IO_ptr)->writeRPCM(*cov,std::string(hdf5Filename));
+}
+
+void Hdf5IO_readRPCM(void * Hdf5IO_ptr,CovarianceData * cov,char* hdf5Filename)
+{
+    ((Hdf5IO*)Hdf5IO_ptr)->readRPCM(*cov,std::string(hdf5Filename));
+}
+
+void* Hdf5IO_initialize()
+{
+    return new Hdf5IO();
+}
+
+void Hdf5IO_destroy(void * Hdf5IO_ptr)
+{
+    delete (Hdf5IO*)Hdf5IO_ptr;
+}
+
diff --git a/sammy/src/io/interface/cpp/Hdf5IOInterface.h b/sammy/src/io/interface/cpp/Hdf5IOInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2ea363683a1a15b17666673a77233eb067573dc
--- /dev/null
+++ b/sammy/src/io/interface/cpp/Hdf5IOInterface.h
@@ -0,0 +1,23 @@
+/*!
+* This file has been dynamically generated by Class Interface Xml (CIX) 
+* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+* If changes need to occur, modify the appropriate CIX xml file
+* Date Generated: Wed Jan 26 13:57:17 EST 2022
+* If any issues are experiences with this generated file that cannot be fixed
+* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+*/
+#ifndef HDF5IOINTERFACE_H
+#define HDF5IOINTERFACE_H
+#include "../../Hdf5IO.h"
+using namespace sammy;
+#ifdef __cplusplus
+extern "C" {
+#endif
+void Hdf5IO_writeRPCM(void * Hdf5IO_ptr,CovarianceData * cov,char* hdf5Filename);
+void Hdf5IO_readRPCM(void * Hdf5IO_ptr,CovarianceData * cov,char* hdf5Filename);
+void* Hdf5IO_initialize();
+void Hdf5IO_destroy(void * Hdf5IO_ptr);
+#ifdef __cplusplus
+}
+#endif
+#endif /* HDF5IOINTERFACE_H */
diff --git a/sammy/src/io/interface/fortran/CoverxIO_I.f90 b/sammy/src/io/interface/fortran/CoverxIO_I.f90
new file mode 100644
index 0000000000000000000000000000000000000000..431294919319debb1eb494acb14b05a0c88b8021
--- /dev/null
+++ b/sammy/src/io/interface/fortran/CoverxIO_I.f90
@@ -0,0 +1,50 @@
+!>
+!! This file has been dynamically generated by Class Interface Xml (CIX) 
+!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+!! If changes need to occur, modify the appropriate CIX xml file
+!! Date Generated: Fri Feb 25 11:18:44 EST 2022
+!! If any issues are experiences with this generated file that cannot be fixed
+!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+!!/
+module CoverxIO_I
+use, intrinsic :: ISO_C_BINDING
+interface
+subroutine f_CoverxIO_setCoverxFile(CoverxIO_ptr, gdList,parCov ) BIND(C,name="CoverxIO_setCoverxFile")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR), value :: CoverxIO_ptr;
+    type(C_PTR), value :: gdList;
+    type(C_PTR), value :: parCov;
+end subroutine
+subroutine f_CoverxIO_setMatMt(CoverxIO_ptr, mat,mt ) BIND(C,name="CoverxIO_setMatMt")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR), value :: CoverxIO_ptr;
+    integer(C_INT) :: mat;
+    integer(C_INT) :: mt;
+end subroutine
+subroutine f_CoverxIO_writeCoverxFile(CoverxIO_ptr, filename ) BIND(C,name="CoverxIO_writeCoverxFile")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR), value :: CoverxIO_ptr;
+    character(C_CHAR) :: filename;
+end subroutine
+subroutine f_CoverxIO_writeFile33(CoverxIO_ptr, za,awr ) BIND(C,name="CoverxIO_writeFile33")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR), value :: CoverxIO_ptr;
+    real(C_DOUBLE) :: za;
+    real(C_DOUBLE) :: awr;
+end subroutine
+type(C_PTR) function f_CoverxIO_initialize( )BIND(C,name="CoverxIO_initialize")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR) :: CoverxIO_ptr;
+end function
+subroutine f_CoverxIO_destroy(this) BIND(C,name="CoverxIO_destroy")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR), value :: this;
+end subroutine
+end interface
+end module CoverxIO_I
diff --git a/sammy/src/io/interface/fortran/CoverxIO_M.f90 b/sammy/src/io/interface/fortran/CoverxIO_M.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b168740a3e6512c690a8214cb44fd59f103075be
--- /dev/null
+++ b/sammy/src/io/interface/fortran/CoverxIO_M.f90
@@ -0,0 +1,63 @@
+!>
+!! This file has been dynamically generated by Class Interface Xml (CIX) 
+!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+!! If changes need to occur, modify the appropriate CIX xml file
+!! Date Generated: Fri Feb 25 11:18:44 EST 2022
+!! If any issues are experiences with this generated file that cannot be fixed
+!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+!!/
+module CoverxIO_M
+use, intrinsic :: ISO_C_BINDING
+use CoverxIO_I
+use CovarianceData_M
+use GridData_M
+type CoverxIO
+    type(C_PTR) :: instance_ptr=C_NULL_PTR
+    contains
+    procedure, pass(this) :: setCoverxFile => CoverxIO_setCoverxFile
+    procedure, pass(this) :: setMatMt => CoverxIO_setMatMt
+    procedure, pass(this) :: writeCoverxFile => CoverxIO_writeCoverxFile
+    procedure, pass(this) :: writeFile33 => CoverxIO_writeFile33
+    procedure, pass(this) :: initialize => CoverxIO_initialize
+    procedure, pass(this) :: destroy => CoverxIO_destroy
+end type CoverxIO
+contains
+subroutine CoverxIO_setCoverxFile(this, gdList, parCov)
+    implicit none
+    class(CoverxIO)::this
+    class(GridDataList)::gdList
+    class(CovarianceData)::parCov
+    call f_CoverxIO_setCoverxFile(this%instance_ptr, gdList%instance_ptr,parCov%instance_ptr)
+end subroutine
+subroutine CoverxIO_setMatMt(this, mat, mt)
+    implicit none
+    class(CoverxIO)::this
+    integer(C_INT)::mat
+    integer(C_INT)::mt
+    call f_CoverxIO_setMatMt(this%instance_ptr, mat,mt)
+end subroutine
+subroutine CoverxIO_writeCoverxFile(this, filename)
+    implicit none
+    class(CoverxIO)::this
+    character(len=*)::filename
+    call f_CoverxIO_writeCoverxFile(this%instance_ptr, filename//C_NULL_CHAR)
+end subroutine
+subroutine CoverxIO_writeFile33(this, za, awr)
+    implicit none
+    class(CoverxIO)::this
+    real(C_DOUBLE)::za
+    real(C_DOUBLE)::awr
+    call f_CoverxIO_writeFile33(this%instance_ptr, za,awr)
+end subroutine
+subroutine CoverxIO_initialize(this)
+    implicit none
+    class(CoverxIO) :: this
+    this%instance_ptr = f_CoverxIO_initialize()
+end subroutine
+subroutine CoverxIO_destroy(this)
+    implicit none
+    class(CoverxIO) :: this
+    call f_CoverxIO_destroy(this%instance_ptr)
+    this%instance_ptr = C_NULL_PTR
+end subroutine
+end module CoverxIO_M
diff --git a/sammy/src/io/interface/fortran/Hdf5IO_I.f90 b/sammy/src/io/interface/fortran/Hdf5IO_I.f90
new file mode 100644
index 0000000000000000000000000000000000000000..f9f92f7f4784f6beeca75ea27dbf08ff10a7c0cd
--- /dev/null
+++ b/sammy/src/io/interface/fortran/Hdf5IO_I.f90
@@ -0,0 +1,37 @@
+!>
+!! This file has been dynamically generated by Class Interface Xml (CIX) 
+!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+!! If changes need to occur, modify the appropriate CIX xml file
+!! Date Generated: Wed Jan 26 13:57:17 EST 2022
+!! If any issues are experiences with this generated file that cannot be fixed
+!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+!!/
+module Hdf5IO_I
+use, intrinsic :: ISO_C_BINDING
+interface
+subroutine f_Hdf5IO_writeRPCM(Hdf5IO_ptr, cov,hdf5Filename ) BIND(C,name="Hdf5IO_writeRPCM")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR), value :: Hdf5IO_ptr;
+    type(C_PTR), value :: cov;
+    character(C_CHAR) :: hdf5Filename;
+end subroutine
+subroutine f_Hdf5IO_readRPCM(Hdf5IO_ptr, cov,hdf5Filename ) BIND(C,name="Hdf5IO_readRPCM")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR), value :: Hdf5IO_ptr;
+    type(C_PTR), value :: cov;
+    character(C_CHAR) :: hdf5Filename;
+end subroutine
+type(C_PTR) function f_Hdf5IO_initialize( )BIND(C,name="Hdf5IO_initialize")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR) :: Hdf5IO_ptr;
+end function
+subroutine f_Hdf5IO_destroy(this) BIND(C,name="Hdf5IO_destroy")
+    use,intrinsic :: ISO_C_BINDING
+    implicit none
+    type(C_PTR), value :: this;
+end subroutine
+end interface
+end module Hdf5IO_I
diff --git a/sammy/src/io/interface/fortran/Hdf5IO_M.f90 b/sammy/src/io/interface/fortran/Hdf5IO_M.f90
new file mode 100644
index 0000000000000000000000000000000000000000..c6c492573d986c0996aca19b7600a8d90a8d4a44
--- /dev/null
+++ b/sammy/src/io/interface/fortran/Hdf5IO_M.f90
@@ -0,0 +1,47 @@
+!>
+!! This file has been dynamically generated by Class Interface Xml (CIX) 
+!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+!! If changes need to occur, modify the appropriate CIX xml file
+!! Date Generated: Wed Jan 26 13:57:17 EST 2022
+!! If any issues are experiences with this generated file that cannot be fixed
+!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+!!/
+module Hdf5IO_M
+use, intrinsic :: ISO_C_BINDING
+use Hdf5IO_I
+use CovarianceData_M
+type Hdf5IO
+    type(C_PTR) :: instance_ptr=C_NULL_PTR
+    contains
+    procedure, pass(this) :: writeRPCM => Hdf5IO_writeRPCM
+    procedure, pass(this) :: readRPCM => Hdf5IO_readRPCM
+    procedure, pass(this) :: initialize => Hdf5IO_initialize
+    procedure, pass(this) :: destroy => Hdf5IO_destroy
+end type Hdf5IO
+contains
+subroutine Hdf5IO_writeRPCM(this, cov, hdf5Filename)
+    implicit none
+    class(Hdf5IO)::this
+    class(CovarianceData)::cov
+    character(len=*)::hdf5Filename
+    call f_Hdf5IO_writeRPCM(this%instance_ptr, cov%instance_ptr,hdf5Filename//C_NULL_CHAR)
+end subroutine
+subroutine Hdf5IO_readRPCM(this, cov, hdf5Filename)
+    implicit none
+    class(Hdf5IO)::this
+    class(CovarianceData)::cov
+    character(len=*)::hdf5Filename
+    call f_Hdf5IO_readRPCM(this%instance_ptr, cov%instance_ptr,hdf5Filename//C_NULL_CHAR)
+end subroutine
+subroutine Hdf5IO_initialize(this)
+    implicit none
+    class(Hdf5IO) :: this
+    this%instance_ptr = f_Hdf5IO_initialize()
+end subroutine
+subroutine Hdf5IO_destroy(this)
+    implicit none
+    class(Hdf5IO) :: this
+    call f_Hdf5IO_destroy(this%instance_ptr)
+    this%instance_ptr = C_NULL_PTR
+end subroutine
+end module Hdf5IO_M
diff --git a/sammy/src/io/tests/CMakeLists.txt b/sammy/src/io/tests/CMakeLists.txt
index 2318d6e1225a364c123ef105208160163a043f15..2b21afbcde595a32401840fbfc93d98a445bece6 100644
--- a/sammy/src/io/tests/CMakeLists.txt
+++ b/sammy/src/io/tests/CMakeLists.txt
@@ -2,3 +2,5 @@ include(NemesisTest)
 
 ADD_NEMESIS_TEST(ExperimentalParametersIOTest.cpp TIMEOUT 60000 NP 1 )
 ADD_NEMESIS_TEST(OdfIOTest.cpp                    TIMEOUT 60000 NP 1 )
+ADD_NEMESIS_TEST(CoverxIOTest.cpp                 TIMEOUT 60000 NP 1 )
+ADD_NEMESIS_TEST(Hdf5IOTest.cpp                   TIMEOUT 60000 NP 1 )
diff --git a/sammy/src/io/tests/CoverxIOTest.cpp b/sammy/src/io/tests/CoverxIOTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1edc441bf1e75584a47f8b2ebdacdb353dd648b4
--- /dev/null
+++ b/sammy/src/io/tests/CoverxIOTest.cpp
@@ -0,0 +1,107 @@
+#include "Nemesis/gtest/Gtest_Functions.hh"
+#include "Nemesis/gtest/nemesis_gtest.hh"
+#include "../CoverxIO.h"
+
+/*******************************************************************************
+ * Test whether we can set the CoverxFile
+ ******************************************************************************/
+TEST(CoverxIO,setCoverxFile)
+{
+    sammy::CoverxIO coverx;
+    sammy::GridDataList gdList;
+    sammy::CovarianceData covData;
+    int mat = 7328; // ta-181
+    int mt = 1;
+
+    // ----- set up the resonance covariance ---------------
+    covData.makeNewUCovariance();
+    endf::ResonanceCovariance * cov = covData.getUCovariance();
+    covData.setNumParameters(2);
+    cov->setCovariance(0,0,2.0); cov->setCovariance(0,1,1.0);
+    cov->setCovariance(1,0,1.0); cov->setCovariance(1,1,4.0);
+    // ----- set up the derivatives ------------------------
+    sammy::GridDataList gridList;
+    std::unique_ptr<sammy::GridData> grid = std::make_unique<sammy::GridData>();
+    // add energy data in first column
+    grid->addData(0,0,1.0); grid->addData(1,0,2.0); grid->addData(2,0,3.0);
+    // add xs data in second column
+    grid->addData(0,1,0.5); grid->addData(1,1,0.5); grid->addData(2,1,0.5);
+    // add derivative 1 data in third column
+    grid->addData(0,2,2.0); grid->addData(1,2,3.0); grid->addData(2,2,4.0);
+    // add derivative 2 data in fourth column
+    grid->addData(0,3,5.0); grid->addData(1,3,6.0); grid->addData(2,3,7.0);
+    grid->setDataIndex(1);
+    gridList.addGrid(grid);
+    
+    // ----- call setCoverxFile() --------------------------
+    coverx.setMatMt(mat,mt);
+    coverx.setCoverxFile(gridList,covData);
+
+    //  /2 5/ /2 1/ /2 3 4/  =  /9  22/ /2 3 4/  =  /128 159 190/ 
+    //  /3 6/ /1 4/ /5 6 7/  =  /12 27/ /5 6 7/  =  /159 198 237/
+    //  /4 7/                   /15 32/             /190 237 284/
+    const coverx::CovMatrix & cm = coverx.getMatrix(mat,mt,mat,mt);
+    double gold00 = cm.getValue(0,0);
+    double gold01 = cm.getValue(0,1);
+    double gold02 = cm.getValue(0,2);
+    double gold11 = cm.getValue(1,1);
+    double gold12 = cm.getValue(1,2);
+    double gold22 = cm.getValue(2,2);
+    
+    ASSERT_EQ( 128.0, gold00 );
+    ASSERT_EQ( 159.0, gold01 );
+    ASSERT_EQ( 190.0, gold02 );
+    ASSERT_EQ( 198.0, gold11 );
+    ASSERT_EQ( 237.0, gold12 );
+    ASSERT_EQ( 284.0, gold22 );
+}
+
+/*******************************************************************************
+ * Test whether we read in the same values that we wrote out
+ * - NOTE: keep in mind that the write operation switches to correlations
+ ******************************************************************************/
+TEST(CoverxIO,readWriteCoverx)
+{
+    sammy::CoverxIO coverx;
+    sammy::CoverxIO coverxRead;
+    sammy::GridDataList gdList;
+    sammy::CovarianceData covData;
+    int mat = 7328; // ta-181
+    int mt = 1;
+
+    // ----- set up the resonance covariance ---------------
+    covData.makeNewUCovariance();
+    endf::ResonanceCovariance * cov = covData.getUCovariance();
+    covData.setNumParameters(2);
+    cov->setCovariance(0,0,2.0); cov->setCovariance(0,1,1.0);
+    cov->setCovariance(1,0,1.0); cov->setCovariance(1,1,4.0);
+    // ----- set up the derivatives ------------------------
+    sammy::GridDataList gridList;
+    std::unique_ptr<sammy::GridData> grid = std::make_unique<sammy::GridData>();
+    // add energy data in first column
+    grid->addData(0,0,1.0); grid->addData(1,0,2.0); grid->addData(2,0,3.0);
+    // add xs data in second column
+    grid->addData(0,1,0.5); grid->addData(1,1,0.5); grid->addData(2,1,0.5);
+    // add derivative 1 data in third column
+    grid->addData(0,2,2.0); grid->addData(1,2,3.0); grid->addData(2,2,4.0);
+    // add derivative 2 data in fourth column
+    grid->addData(0,3,5.0); grid->addData(1,3,6.0); grid->addData(2,3,7.0);
+    grid->setDataIndex(1);
+    gridList.addGrid(grid);
+    
+    // ----- call setCoverxFile() --------------------------
+    coverx.setMatMt(mat,mt);
+    coverx.setCoverxFile(gridList,covData);
+    coverx.writeCoverxFile("test.coverx");
+
+    coverxRead.readCoverxFile("test.coverx");
+
+    const coverx::CovMatrix & write = coverx.getMatrix(mat,mt,mat,mt);
+    const coverx::CovMatrix & read  = coverxRead.getMatrix(mat,mt,mat,mt);
+
+    for( int i=0;i<3;++i ){
+        for( int j=0;j<=i;++j){
+            ASSERT_EQ( write.getValue(i,j), read.getValue(i,j) );
+        }
+    }
+}
\ No newline at end of file
diff --git a/sammy/src/io/tests/Hdf5IOTest.cpp b/sammy/src/io/tests/Hdf5IOTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..53fb088ac48f40c42c7262776062a999c61968c1
--- /dev/null
+++ b/sammy/src/io/tests/Hdf5IOTest.cpp
@@ -0,0 +1,35 @@
+#include "Nemesis/gtest/Gtest_Functions.hh"
+#include "Nemesis/gtest/nemesis_gtest.hh"
+#include "../Hdf5IO.h"
+
+/*******************************************************************************
+ * Test whether we can read and write the res param cov matrix
+ ******************************************************************************/
+TEST(Hdf5IO,readWriteRPCM)
+{
+    sammy::CovarianceData covData;
+    sammy::CovarianceData covDataR;
+    sammy::Hdf5IO h5;
+    sammy::Hdf5IO h5R;
+
+    // ----- set up the resonance covariance ---------------
+    covData.makeNewUCovariance();
+    endf::ResonanceCovariance * cov = covData.getUCovariance();
+    covData.setNumParameters(2);
+    cov->setCovariance(0,0,2.0); cov->setCovariance(0,1,1.0);
+    cov->setCovariance(1,0,1.0); cov->setCovariance(1,1,4.0);
+
+    h5.writeRPCM(covData, "testreadwrite.h5");
+    h5R.readRPCM(covDataR,"testreadwrite.h5");
+
+    endf::ResonanceCovariance * covR = covDataR.getUCovariance();
+
+    double gold00 = 2.0;
+    double gold01 = 1.0;
+    double gold10 = 1.0;
+    double gold11 = 4.0;
+    ASSERT_EQ( covR->getCovariance(0,0), gold00 );
+    ASSERT_EQ( covR->getCovariance(0,1), gold01 );
+    ASSERT_EQ( covR->getCovariance(1,0), gold10 );
+    ASSERT_EQ( covR->getCovariance(1,1), gold11 );
+}
diff --git a/sammy/src/lru/mlru0.f b/sammy/src/lru/mlru0.f90
similarity index 71%
rename from sammy/src/lru/mlru0.f
rename to sammy/src/lru/mlru0.f90
index d3bf4861fda31377917dc543066a6ecea5685e34..f6add8ad5c5e1c9b9c5b8db9ccd9fcf341d706ed 100644
--- a/sammy/src/lru/mlru0.f
+++ b/sammy/src/lru/mlru0.f90
@@ -1,10 +1,11 @@
-C
-C
-C
+module lru_m
+    contains
+!
+!
       SUBROUTINE Samlru_0
-C
-C *** Purpose -- write ENDF file 2 and 32 for unresolved resonance region
-C
+!
+! *** Purpose -- write ENDF file 2 and 32 for unresolved resonance region
+!
       use oops_common_m
       use fixedi_m
       use ifwrit_m
@@ -14,14 +15,16 @@ C
       use exploc_urr_common_m
       use AllocateFunctions_m
       use SammyLptPrinting_m
+      use lru1_m
+      use lru3_m
       IMPLICIT DOUBLE PRECISION (A-h,o-z)
       real(kind=8),allocatable,dimension(:)::A_Idiag
       integer,allocatable,dimension(:)::I_Jderiv
       real(kind=8),allocatable,dimension(:)::A_Ixderi
       real(kind=8),allocatable,dimension(:)::A_Icovar
       character(len=80)::line
-C
-C
+!
+!
 
           WRITE (line,99999)
 99999     FORMAT (' *** SAMMY-LRU   13 Jul 06 ***')
@@ -30,43 +33,42 @@ C
       Segmen(2) = 'R'
       Segmen(3) = 'U'
       Nowwww = 0
-C
+!
       CALL Initil
       I = Idimen (1, 1, '1, 1')
       I = Idimen (I, -1, 'I, -1')
       I = Idimen (0, 0, '0, 0')
-C
+!
       Num_Par_Max = 5000
       If (Ndfcov.NE.0) THEN
          N = (Nvpall*(Nvpall+1))/2        
          call allocate_integer_data(I_Jderiv, 2* Num_Par_Max) ! now allocated as integer
          call allocate_real_data(A_Ixderi, Num_Par_Max*2)
       END IF
-C
-      CALL Urr (  A_Istren , A_Idistn , A_Kkkkgg , A_Kkkkgf , A_Kkkfnu ,
-     * A_Kkethr , A_Kkwthr , A_Ibethj , A_Kkkkex , A_Kkspin , A_Kkkpty ,
-     * I_Kkkjnl , I_Kkkjxl , I_Kkkj2n , A_Kbinde , A_Kpaire , A_Iengur ,
-     * A_Ialevl , I_Iflstr , I_Ifldst , I_Iflggg , I_Iflgff , I_Jderiv ,
-     * A_Ixderi , Num_Par, Lcomp, Mat, Ns)
+!
+      CALL Urr (  A_Istren , A_Idistn , A_Kkkkgg , A_Kkkkgf , A_Kkkfnu , &
+       A_Kkethr , A_Kkwthr , A_Ibethj , A_Kkkkex , A_Kkspin , A_Kkkpty , &
+       I_Kkkjnl , I_Kkkjxl , I_Kkkj2n , A_Kbinde , A_Kpaire , A_Iengur , &
+       A_Ialevl , I_Iflstr , I_Ifldst , I_Iflggg , I_Iflgff , I_Jderiv , &
+       A_Ixderi , Num_Par, Lcomp, Mat, Ns)
       IF (Num_Par.GT.Num_Par_Max) THEN
          WRITE (6,10000) Num_Par, Num_Par_Max
-10000    FORMAT ('Num_Par,Num_Par_Max =', 2I10, /,
-     *      'Ask NML to change Num_Par_Max to be > Num_Par')
+10000    FORMAT ('Num_Par,Num_Par_Max =', 2I10, /, &
+            'Ask NML to change Num_Par_Max to be > Num_Par')
          STOP '[STOP in Samlru_0 in lru/mlru0.f]'
       END IF
-C
+!
       IF (Ndfcov.NE.0) THEN
          call allocate_real_data(A_Icovar, (Num_Par*(Num_Par+1))/2)
-         CALL Fix_Covar (I_Jderiv, A_Ixderi,  A_Icovar,
-     *      Num_Par)
+         CALL Fix_Covar (I_Jderiv, A_Ixderi,  A_Icovar, &
+            Num_Par)
          IF (Lcomp.EQ.2) THEN
             allocate(A_Idiag(Num_par))
-            CALL Write_Cmpc_Cov_2 (A_Icovar, A_Idiag, Num_Par,
-     *         Mat, Ns)
+            CALL Write_Cmpc_Cov_2 (A_Icovar, A_Idiag, Num_Par, Mat, Ns)
             deallocate(A_Idiag)
          ELSE IF (Lcomp.EQ.1) THEN
-            CALL Write_Urr_Lcomp_1 (A_Icovar,
-     *         Nvpall, Num_Par, If_Fis, Mat, Ns)
+            CALL Write_Urr_Lcomp_1 (A_Icovar, &
+               Nvpall, Num_Par, If_Fis, Mat, Ns)
          ELSE
             WRITE (21,10500) Lcomp
             WRITE (6,10500) Lcomp
@@ -76,7 +78,7 @@ C
       if (allocated(I_Jderiv)) deallocate(I_Jderiv)
       if (allocated(A_Ixderi)) deallocate(A_Ixderi)
       if (allocated(A_Icovar)) deallocate(A_Icovar)
-C
+!
       CLOSE (UNIT=55)
       I = Idimen (0, 0, '0, 0')
       CALL Timer (0)
@@ -85,3 +87,6 @@ C
 10100 FORMAT (' Normal finish to samlru')
       STOP
       END
+
+end module lru_m
+
diff --git a/sammy/src/lru/mlru1.f b/sammy/src/lru/mlru1.f90
similarity index 59%
rename from sammy/src/lru/mlru1.f
rename to sammy/src/lru/mlru1.f90
index f9a58782b7664928487e0960bbb2f7746d1a58ec..8276872d9e413c6d9eaae9d27737193148e251bc 100644
--- a/sammy/src/lru/mlru1.f
+++ b/sammy/src/lru/mlru1.f90
@@ -1,12 +1,14 @@
-C
-C
-C ______________________________________________________________
-C
-      Subroutine Urr (Streng, Distnt, Gg, Gf, Fnu, Ethr, Wthr, Bethdj,
-     *   Ex, Spin, Pty, Jnl, Jxl, J2n, Bindee, Pairee, Engurr, Alevel,
-     *   Iflstr, Ifldst, Iflggg, Iflgff, Jderiv, Xderiv,  Num_Par,
-     *   Lcomp, Mat, Ns)
-C
+module lru1_m
+    contains
+!
+!
+! ______________________________________________________________
+!
+      Subroutine Urr (Streng, Distnt, Gg, Gf, Fnu, Ethr, Wthr, Bethdj, &
+         Ex, Spin, Pty, Jnl, Jxl, J2n, Bindee, Pairee, Engurr, Alevel, &
+         Iflstr, Ifldst, Iflggg, Iflgff, Jderiv, Xderiv,  Num_Par, &
+         Lcomp, Mat, Ns)
+!
       use fixedi_m
       use ifwrit_m
       use samxxx_common_m
@@ -15,31 +17,32 @@ C
       use constn_common_m
       use EndfData_common_m, only : radFitFlags
       use acs6_m
+      use lru2_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
-      DIMENSION Streng(Numelv,*), Distnt(Numelv,*), Gg(Numelv,*),
-     *   Gf(Numelv,Numjjv,*),Fnu(Numelv,*),
-     *   Ethr(Numelv,Numjjv,*), Wthr(Numelv,Numjjv,*), Bethdj(Numjjv,*),
-     *   Ex(*), Spin(*), Pty(*), Jnl(*), Jxl(*), J2n(*),
-     *   Bindee(*), Pairee(*), Engurr(*), Alevel(*),
-     *   Iflstr(Numelv,*), Ifldst(Numelv,*), Iflggg(Numelv,*),
-     *   Iflgff(Numelv,Numjjv,*), Jderiv(2,*), Xderiv(2,*)
+!
+      DIMENSION Streng(Numelv,*), Distnt(Numelv,*), Gg(Numelv,*), &
+         Gf(Numelv,Numjjv,*),Fnu(Numelv,*), &
+         Ethr(Numelv,Numjjv,*), Wthr(Numelv,Numjjv,*), Bethdj(Numjjv,*), &
+         Ex(*), Spin(*), Pty(*), Jnl(*), Jxl(*), J2n(*), &
+         Bindee(*), Pairee(*), Engurr(*), Alevel(*), &
+         Iflstr(Numelv,*), Ifldst(Numelv,*), Iflggg(Numelv,*), &
+         Iflgff(Numelv,Numjjv,*), Jderiv(2,*), Xderiv(2,*)
       DIMENSION Eeex(100)
       DATA One /1.0d0/, Two /2.0d0/, Ten /10.0d0/
       DATA Lrf /2/, Lx /0/, Mf /2/
-C     Lru = 2
-C ***       2 for Unresolved Resonance Region (URR)
-C     Lrf = 2
-C ***       2 for "everything is energy-dependent" more-or-less
-C ***         for the File 2 case.  For File 32, have to use only
-C ***         energy-independent stuff... bad approximation.
-C
-C
+!     Lru = 2
+! ***       2 for Unresolved Resonance Region (URR)
+!     Lrf = 2
+! ***       2 for "everything is energy-dependent" more-or-less
+! ***         for the File 2 case.  For File 32, have to use only
+! ***         energy-independent stuff... bad approximation.
+!
+!
       Mfx =32
       Mt = 151
       Ns = 1
       Nisotp = 1
-C
+!
       CALL Oldopn (16, Sam16x, 0)
       READ (16,89999,END=210,ERR=210) Summry
       READ (16,89999,END=210,ERR=210) Summry
@@ -48,17 +51,17 @@ C
 89999 FORMAT (A70)
       IF (Summry.EQ.Fblank) GO TO 200
       CLOSE (UNIT=16)
-C
+!
       Kz = 0
       CALL Filopn (10, Summry, 0)
       Krdndf = 1
       Kdebug = 1
-C
-      CALL Read_Ndf_Urr (Engurr, Eeex, Engur0, Za, Awr, Eee, Eminz,
-     *   Emaxz, Iwhich, Lrf, Mat, Kz, Num_Energy, Lcomp, Numurr,
-     *   Numex, Nume)
-C
-C *** done inputting data ... now can start writing:
+!
+      CALL Read_Ndf_Urr (Engurr, Eeex, Engur0, Za, Awr, Eee, Eminz, &
+         Emaxz, Iwhich, Lrf, Mat, Kz, Num_Energy, Lcomp, Numurr, &
+         Numex, Nume)
+!
+! *** done inputting data ... now can start writing:
       CALL Newopn (57, Sam57x, 0)
       Iu = 57
       IF (Ndfcov.NE.0) THEN
@@ -67,41 +70,41 @@ C *** done inputting data ... now can start writing:
       ELSE
          Iux = 0
       END IF
-C
-C *** first entry in ENDF/B-VI file -- for sample material -- HEAD record
+!
+! *** first entry in ENDF/B-VI file -- for sample material -- HEAD record
       CALL Cont_Ndf (Za, Awr, Lx, Lx, Nisotp, Lx, Mat, Mf, Mt, Ns, Iu)
-      IF (Iux.NE.0) CALL Cont_Ndf (Za, Awr, Lx, Lx, Nisotp, Lx, Mat,
-     *   Mfx, Mt, Ns, Iux)
-C
-C
-      Ap = radFitFlags%getMatchingRadius()/Ten *
-     *       (One - Distnt(1,1))
-C *** Note that Distnt can depend on energy-region, but ENDF does not
-C ***    permit that; hence we assume Distnt(1,1) not (1,Kumurr)
-C
-C *** ENTRANCE CHANNEL SPINS (TIMES 2)
+      IF (Iux.NE.0) CALL Cont_Ndf (Za, Awr, Lx, Lx, Nisotp, Lx, Mat, &
+         Mfx, Mt, Ns, Iux)
+!
+!
+      Ap = radFitFlags%getMatchingRadius()/Ten * &
+             (One - Distnt(1,1))
+! *** Note that Distnt can depend on energy-region, but ENDF does not
+! ***    permit that; hence we assume Distnt(1,1) not (1,Kumurr)
+!
+! *** ENTRANCE CHANNEL SPINS (TIMES 2)
       I2   = Spin(1)*Two + 0.001d0
       Is2n = iABS(I2-1)
       Is2x =      I2+1
-C
-C *** determine Lfw =1 if fission widths are present, =0 if not
-C *** Note:  ENDF manual says this is for UNresolved region only.  So
-C ***        ingore this in /ndf/ but get it right here.
+!
+! *** determine Lfw =1 if fission widths are present, =0 if not
+! *** Note:  ENDF manual says this is for UNresolved region only.  So
+! ***        ingore this in /ndf/ but get it right here.
       Lfw = If_Fis
       IF (Lfw.NE.1) Lfw = 0
-C
-C *** At the moment we allow only one isotope in sam-URR.  Keep this loop,
-C ***    to make it easier to generalize eventually.
-C
-C *** Begin do loop over isotopes
+!
+! *** At the moment we allow only one isotope in sam-URR.  Keep this loop,
+! ***    to make it easier to generalize eventually.
+!
+! *** Begin do loop over isotopes
       DO Misotp=1,Nisotp
          Mmsotp = Misotp
-C
-C ***    several entries in ENDF/B-VI file -- for each isotope in material
-         CALL Write_Urr_1 (Spin(1), Ap, Za, Eminz, Emaxz, Lcomp, Lfw,
-     *      Numelv, Mat, Iu, Iux)
-C
-C ---------------------------------
+!
+! ***    several entries in ENDF/B-VI file -- for each isotope in material
+         CALL Write_Urr_1 (Spin(1), Ap, Za, Eminz, Emaxz, Lcomp, Lfw, &
+            Numelv, Mat, Iu, Iux)
+!
+! ---------------------------------
          Ipar = 0
          Parity = - Pty(1)
          DO LL=1,Numelv
@@ -110,44 +113,44 @@ C ---------------------------------
             Jlo = Jnl(L)
             Jhi = Jxl(L)
             J2  = J2n(L)
-C
-C           Find values of Pj = J^pi for first J_value for this L
+!
+!           Find values of Pj = J^pi for first J_value for this L
             Qj = dFLOAT(J2)/2.0d0
             Pj = Qj*Parity
-C
+!
             Nj = Jhi - Jlo + 1
             CALL Write_Urr_2 (Awr, Lcomp, L, Nj, Mat, Iu, Iux)
-C
+!
             DO Jxxx=Jlo,Jhi
                Jxx = Jxxx
-C
+!
                CALL Write_Urr_3 (Qj, Nume, Numex, Lcomp, Mat, Iu, Iux)
-C
+!
                Amuf = Fnu(L,Jxx)
                Kk = Kwhich (J2, L, Amun, Is2n, Is2x)
                IF (Numexc.GT.1) Amux = One
                CALL Write_Urr_4 (Amux, Amun, Amuf, Lcomp, Mat, Iu, Iux)
-C
-               CALL Parout_Endf (Streng, Distnt, Gg, Gf, Ethr, Wthr,
-     *            Bethdj, Ex, Spin, Pty, Bindee, Pairee, Engurr, Alevel,
-     *            Iflstr, Ifldst, Iflggg, Iflgff, Jderiv, Xderiv, 
-     *            Eeex, Numex, Pj, L, Jxx, Num_Energy, Mat, Iu, Iux,
-     *            Num_Par, Ipar, Iwhich, Qj, Lcomp, J2, Parity)
-C
-C ---------    Now get ready for next J^pi value
+!
+               CALL Parout_Endf (Streng, Distnt, Gg, Gf, Ethr, Wthr, &
+                  Bethdj, Ex, Spin, Pty, Bindee, Pairee, Engurr, Alevel, &
+                  Iflstr, Ifldst, Iflggg, Iflgff, Jderiv, Xderiv,  &
+                  Eeex, Numex, Pj, L, Jxx, Num_Energy, Mat, Iu, Iux, &
+                  Num_Par, Ipar, Iwhich, Qj, Lcomp, J2, Parity)
+!
+! ---------    Now get ready for next J^pi value
                Qj = Qj + One
                Pj = Qj*Parity
                J2 = J2 + 2
-C
+!
             END DO
-C         # END DO Jxx=Jlo,Jhi
-C
+!         # END DO Jxx=Jlo,Jhi
+!
          END DO
-C      # END DO L=1,Numelv
-C
+!      # END DO L=1,Numelv
+!
       END DO
-C   # END DO Misotp=1,Nisotp
-C
+!   # END DO Misotp=1,Nisotp
+!
       CLOSE (UNIT=57)
       CALL Rewrit_Ndf (Iu, Ns)
       IF (Iux.NE.0) THEN
@@ -155,20 +158,20 @@ C
          CALL Rewrit_Ndf (Iux, Ns)
       END IF
       RETURN
-C
+!
   200 WRITE (6,11000) Summry
 11000 FORMAT (' ENDF file name is wrong: ##', A70, '##')
       STOP '[STOP in Urr in lru/mlru1.f    # 1]'
-C
+!
   210 WRITE (6,11100)
 11100 FORMAT (' Could not read SAM16.DAT to find name of file')
       STOP '[STOP in Urr in lru/mlru1.f    # 2]'
-C
+!
       END
-C
-C
-C ______________________________________________________________
-C
+!
+!
+! ______________________________________________________________
+!
       SUBROUTINE Read_Eee (Eeex, Eee, Numex)
       use fixedi_m
       use ifwrit_m
@@ -180,7 +183,7 @@ C
       DIMENSION Eeex(*)
       CHARACTER*1 Alpha(80)
       DATA Zero /0.0d0/
-C
+!
       IF (Eee.LT.Zero) THEN
          In = 0
       ELSE
@@ -216,34 +219,34 @@ C
             END IF
          END IF
          GO TO 10
-C
+!
    20 CONTINUE
       Numex = In
       RETURN
       END
-C
-C
-C ______________________________________________________________
-C
-      SUBROUTINE Read_Ndf_Urr (Engurr, Eeex, Engur0, Za, Awr, Eee,
-     *   Eminz, Emaxz, Iwhich, Lrf, Mat, Kz, Num_Energy, Lcomp,
-     *   Numurr, Numex, Nume)
+!
+!
+! ______________________________________________________________
+!
+      SUBROUTINE Read_Ndf_Urr (Engurr, Eeex, Engur0, Za, Awr, Eee, &
+         Eminz, Emaxz, Iwhich, Lrf, Mat, Kz, Num_Energy, Lcomp, &
+         Numurr, Numex, Nume)
       use ifwrit_m
       IMPLICIT DOUBLE PRECISION (A-h,o-z)
       DIMENSION Engurr(*), Eeex(*)
       DATA Zero /0.0d0/
-C
-C *** Initialize, then read the *.NDF file
+!
+! *** Initialize, then read the *.NDF file
       Eee = Zero
       Num_Energy = 0
       Lcomp = 2
-      CALL Read_Ndf (Za, Awr, Eee, Defunc, Lrf, Mat, Kz, Num_Energy,
-     *   Lcomp, Nd2, Ifflag)
+      CALL Read_Ndf (Za, Awr, Eee, Defunc, Lrf, Mat, Kz, Num_Energy, &
+         Lcomp, Nd2, Ifflag)
       IF (Num_Energy.EQ.0) Num_Energy = 1
-C
-C *** Fill in gaps not given in *.NDF file
+!
+! *** Fill in gaps not given in *.NDF file
       IF (Eee.EQ.Zero) THEN
-C        Eee = Zero means there's no info in Read_Ndf file re energies
+!        Eee = Zero means there's no info in Read_Ndf file re energies
          Iwhich = 1
          Eminz = Engur0
          Emaxz = Engurr(Numurr)
@@ -257,13 +260,13 @@ C        Eee = Zero means there's no info in Read_Ndf file re energies
             END IF
          ELSE
             Numex = Num_Energy*Numurr
-C                   Num_Energy points inside each of Numurr energy regions
+!                   Num_Energy points inside each of Numurr energy regions
             IF (Kkends.EQ.1) Numex = Numex + 2
-C                   2 for end-points Emin & Emax
+!                   2 for end-points Emin & Emax
             Nume = (Numex+1) * 6
          END IF
       ELSE
-C        Eee > Zero or Eee < Zero means read more energies
+!        Eee > Zero or Eee < Zero means read more energies
          CALL Read_Eee (Eeex, Eee, Numex)
          Nume = (Numex+1) * 6
          IF (Numex.GT.100) STOP '[STOP in Urr in lru/mlur1.f Numex>100]'
@@ -274,19 +277,19 @@ C        Eee > Zero or Eee < Zero means read more energies
       END IF
       RETURN
       END
-C
-C
-C ______________________________________________________________
-C
-      SUBROUTINE Write_Urr_1 (Spin, Ap, Za, Eminz, Emaxz, Lcomp, Lfw,
-     *   Numelv, Mat, Iu, Iux)
+!
+!
+! ______________________________________________________________
+!
+      SUBROUTINE Write_Urr_1 (Spin, Ap, Za, Eminz, Emaxz, Lcomp, Lfw, &
+         Numelv, Mat, Iu, Iux)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DATA Lru /2/, Lrf /2/, L1 /1/, Lx /0/, Mf /2/, Mfx /32/, Mt /151/
       DATA One /1.0d0/
       Ns = 1
-C
-C *** second entry in ENDF/B-VI file -- for each isotope in material
-C *** CONT record starting with (Z,A) for isotope
+!
+! *** second entry in ENDF/B-VI file -- for each isotope in material
+! *** CONT record starting with (Z,A) for isotope
       WRITE (Iu,10100)
 10100 FORMAT ('####' /, '#### Z_A')
       CALL Cont_Ndf (Za, One, Lx, Lfw, L1, Lx, Mat, Mf, Mt, Ns, Iu)
@@ -294,115 +297,117 @@ C *** CONT record starting with (Z,A) for isotope
          WRITE (Iux,10100)
          CALL Cont_Ndf (Za, One, Lx, Lfw, L1, Lx, Mat,Mfx,Mt,Ns,Iux)
       END IF
-C
-C *** third entry in ENDF/B-VI file -- for each isotope in material
-C *** CONT record starting with energy range for isotope
+!
+! *** third entry in ENDF/B-VI file -- for each isotope in material
+! *** CONT record starting with energy range for isotope
       WRITE (Iu, 10200) Lrf
-10200 FORMAT ('####' /, '#### Emin    Emax           Lru=2 => URR',
-     *      I4, '=LRF')
+10200 FORMAT ('####' /, '#### Emin    Emax           Lru=2 => URR', &
+            I4, '=LRF')
       Nro = 0
       Naps= 0
-      CALL Cont_Ndf (Eminz, Emaxz, Lru, Lrf, Nro, Naps,
-     *   Mat, Mf, Mt, Ns, Iu)
+      CALL Cont_Ndf (Eminz, Emaxz, Lru, Lrf, Nro, Naps, &
+         Mat, Mf, Mt, Ns, Iu)
       IF (Iux.NE.0) THEN
          Lrfx = Lrf
          IF (Lcomp.EQ.1) Lrfx = 1
          WRITE (Iux,10200) Lrfx
-         CALL Cont_Ndf (Eminz, Emaxz, Lru, Lrfx, Nro, Naps,
-     *      Mat, Mfx, Mt, Ns, Iux)
+         CALL Cont_Ndf (Eminz, Emaxz, Lru, Lrfx, Nro, Naps, &
+            Mat, Mfx, Mt, Ns, Iux)
       END IF
-C
-C *** Next entry in ENDF/B-VI file -- CONT record giving "number of
-C ***    values of orbital angular momentum" = Numelv
-10300 FORMAT ('####' /, '#### Spin    Radius', 13x,
-     *   ' Number of L values =', I2)
+!
+! *** Next entry in ENDF/B-VI file -- CONT record giving "number of
+! ***    values of orbital angular momentum" = Numelv
+10300 FORMAT ('####' /, '#### Spin    Radius', 13x, &
+         ' Number of L values =', I2)
       WRITE (Iu,10300) Numelv
-      CALL Cont_Ndf (Spin, Ap, Lx, Lx, Numelv, Lx, Mat, Mf,
-     *   Mt, Ns, Iu)
+      CALL Cont_Ndf (Spin, Ap, Lx, Lx, Numelv, Lx, Mat, Mf, &
+         Mt, Ns, Iu)
       IF (Iux.NE.0) THEN
          WRITE (Iux,10300) Numelv
-         CALL Cont_Ndf (Spin, Ap, Lx, Lx, Numelv,
-     *      Lx, Mat, Mfx, Mt, Ns, Iux)
+         CALL Cont_Ndf (Spin, Ap, Lx, Lx, Numelv, &
+            Lx, Mat, Mfx, Mt, Ns, Iux)
          IF (Lcomp.EQ.1) WRITE (Iux,10350)
-10350    FORMAT ('####' /, '#### E_Spacng  J       Gamma_n',
-     *      4X, 'Gamma_gam  Gamma_f    Gamma_x')
+10350    FORMAT ('####' /, '#### E_Spacng  J       Gamma_n', &
+            4X, 'Gamma_gam  Gamma_f    Gamma_x')
       END IF
-C
+!
       RETURN
       END
-C
-C
-C ______________________________________________________________
-C
+!
+!
+! ______________________________________________________________
+!
       SUBROUTINE Write_Urr_2 (Awr, Lcomp, L, Nj, Mat, Iu, Iux)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DATA Lx /0/, Mf /2/, Mfx /32/, Mt /151/
       DATA Zero /0.0d0/
       Ns = 1
       WRITE (Iu,10400) L-1, Nj
-10400 FORMAT ('####' /, '#### Awr', 21X, 'L=', I2,
-     *        ' Number of J values=', I2)
+10400 FORMAT ('####' /, '#### Awr', 21X, 'L=', I2, &
+              ' Number of J values=', I2)
       CALL Cont_Ndf (Awr, Zero, L-1, Lx, Nj, Lx, Mat, Mf, Mt, Ns, Iu)
       IF (Iux.NE.0) THEN
-C   # IF (want covariance matrix in File 32) then
+!   # IF (want covariance matrix in File 32) then
          IF (Lcomp.EQ.2) THEN
             WRITE (Iux,10400) L-1, Nj
-            CALL Cont_Ndf (Awr, Zero, L-1, Lx, Nj, Lx, Mat, Mfx, 
-     *         Mt, Ns, Iux)
+            CALL Cont_Ndf (Awr, Zero, L-1, Lx, Nj, Lx, Mat, Mfx,  &
+               Mt, Ns, Iux)
          ELSE IF (Lcomp.EQ.1) THEN
             WRITE (Iux,10500) L-1, Nj
-10500       FORMAT ('####' /, '#### Awr', 21X, 'L=', I2, 11X,
-     *              ' Number of J values=', I2)
+10500       FORMAT ('####' /, '#### Awr', 21X, 'L=', I2, 11X, &
+                    ' Number of J values=', I2)
             Nj6 = Nj*6
-            CALL Cont_Ndf (Awr, Zero, L-1, Lx, Nj6, Nj, Mat, Mfx,
-     *         Mt, Ns, Iux)
+            CALL Cont_Ndf (Awr, Zero, L-1, Lx, Nj6, Nj, Mat, Mfx, &
+               Mt, Ns, Iux)
          END IF
       END IF
-C   # END IF (Iux.NE.0) THEN
+!   # END IF (Iux.NE.0) THEN
       RETURN
       END
-C
-C
-C ______________________________________________________________
-C
+!
+!
+! ______________________________________________________________
+!
       SUBROUTINE Write_Urr_3 (Qj, Nume, Numex, Lcomp, Mat, Iu, Iux)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DATA Lx /0/, Mf /2/, Mfx /32/, Mt /151/, Int /2/
       DATA Zero /0.0d0/
       WRITE (Iu,10600) Numex
-10600 FORMAT ('####' /, '#### J^pi', 19X,
-     *        'INT=2            Number of E values=', I2)
+10600 FORMAT ('####' /, '#### J^pi', 19X, &
+              'INT=2            Number of E values=', I2)
       Ns = 1
       CALL Cont_Ndf (Qj, Zero, Int, Lx, Nume, Numex, Mat, Mf, Mt, Ns,Iu)
       IF (Iux.NE.0 .AND. Lcomp.EQ.2) THEN
          WRITE (Iux,10600) Numex
-         CALL Cont_Ndf (Qj, Zero, Int, Lx, Nume, Numex, Mat, Mfx, Mt,
-     *      Ns, Iux)
+         CALL Cont_Ndf (Qj, Zero, Int, Lx, Nume, Numex, Mat, Mfx, Mt, &
+            Ns, Iux)
       END IF
       RETURN
       END
-C
-C
-C ______________________________________________________________
-C
+!
+!
+! ______________________________________________________________
+!
       SUBROUTINE Write_Urr_4 (Amux, Amun, Amuf, Lcomp, Mat, Iu, Iux)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DATA Mf /2/, No /0/, Mfx /32/, Mt /151/
       DATA Zero /0.0d0/
       WRITE (Iu,10700)
-10700 FORMAT ('####' /, '#### Deg of Frdm for   inelastic  n',
-     *   10X, 'capture    fission')
-      CALL List_7x (Zero, Zero, Amux, Amun, Zero, Amuf,
-     *   Mat, Mf, Mt, Ns, No, No, No, No, No, No, Iu)
+10700 FORMAT ('####' /, '#### Deg of Frdm for   inelastic  n', &
+         10X, 'capture    fission')
+      CALL List_7x (Zero, Zero, Amux, Amun, Zero, Amuf, &
+         Mat, Mf, Mt, Ns, No, No, No, No, No, No, Iu)
       IF (Iux.GT.0 .AND. Lcomp.EQ.2) THEN
          WRITE (Iux,10700)
-         CALL List_7x (Zero, Zero, Amux, Amun, Zero, Amuf,
-     *      Mat, Mfx, Mt, Ns, No, No, No, No, No, No, Iux)
+         CALL List_7x (Zero, Zero, Amux, Amun, Zero, Amuf, &
+            Mat, Mfx, Mt, Ns, No, No, No, No, No, No, Iux)
       END IF
-C
+!
       WRITE (Iu,10800)
-10800 FORMAT ('####' /, '#### Energy Spacing    Gamma_x',
-     *   4X, 'Gamma_n    Gamma_gam  Gamma_f')
+10800 FORMAT ('####' /, '#### Energy Spacing    Gamma_x', &
+         4X, 'Gamma_n    Gamma_gam  Gamma_f')
       IF (Iux.GT.0 .AND. Lcomp.EQ.2) WRITE (Iux,10800)
       RETURN
       END
+
+end module lru1_m
\ No newline at end of file
diff --git a/sammy/src/lru/mlru2.f b/sammy/src/lru/mlru2.f90
similarity index 65%
rename from sammy/src/lru/mlru2.f
rename to sammy/src/lru/mlru2.f90
index 9b5af552ff36d0a8333b2c190e68ee5705f8f855..4946be48a4263553fe63b29e02ba9256e0a5aed3 100644
--- a/sammy/src/lru/mlru2.f
+++ b/sammy/src/lru/mlru2.f90
@@ -1,15 +1,16 @@
-C
-C
-C ______________________________________________________________________
-C
-      SUBROUTINE Parout_Endf (Streng, Distnt, Gg, Gf, Ethr, Wthr,
-     *   Bethdj, Ex, Spin, Pty, Bindee, Pairee, Engurr, Alevel, Iflstr,
-     *   Ifldst, Iflggg, Iflgff, Jderiv, Xderiv, Eeex, Numex, Pj,
-     *   L, Jxx, Num_Energy, Mat, Iu, Iux, Num_Par, Ipar, Iwhich, Qj,
-     *   Lcomp, J2, Parity)
-C
-C *** Purpose-- Find and print parameters into ENDF file
-C
+module lru2_m
+   contains
+!
+! ______________________________________________________________________
+!
+      SUBROUTINE Parout_Endf (Streng, Distnt, Gg, Gf, Ethr, Wthr, &
+         Bethdj, Ex, Spin, Pty, Bindee, Pairee, Engurr, Alevel, Iflstr, &
+         Ifldst, Iflggg, Iflgff, Jderiv, Xderiv, Eeex, Numex, Pj, &
+         L, Jxx, Num_Energy, Mat, Iu, Iux, Num_Par, Ipar, Iwhich, Qj, &
+         Lcomp, J2, Parity)
+!
+! *** Purpose-- Find and print parameters into ENDF file
+!
       use fixedi_m
       use ifwrit_m
       use fixedr_m
@@ -19,29 +20,29 @@ C
       use acs4_m
       use acs6_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
-      DIMENSION Streng(Numelv,*), Distnt(Numelv,*), Gg(Numelv,*),
-     *   Gf(Numelv,Numjjv,*), Ethr(Numelv,Numjjv,*),
-     *   Wthr(Numelv,Numjjv,*), Bethdj(Numjjv,*), Ex(*), Spin(*),
-     *   Pty(*), Bindee(*), Pairee(*), Engurr(*), Alevel(*),
-     *   Iflstr(Numelv,*), Ifldst(Numelv,*), Iflggg(Numelv,*),
-     *   Iflgff(Numelv,Numjjv,*), Jderiv(2,*), Xderiv(2,*), 
-     *   Eeex(*)
+!
+      DIMENSION Streng(Numelv,*), Distnt(Numelv,*), Gg(Numelv,*), &
+         Gf(Numelv,Numjjv,*), Ethr(Numelv,Numjjv,*), &
+         Wthr(Numelv,Numjjv,*), Bethdj(Numjjv,*), Ex(*), Spin(*), &
+         Pty(*), Bindee(*), Pairee(*), Engurr(*), Alevel(*), &
+         Iflstr(Numelv,*), Ifldst(Numelv,*), Iflggg(Numelv,*), &
+         Iflgff(Numelv,Numjjv,*), Jderiv(2,*), Xderiv(2,*),  &
+         Eeex(*)
       DATA Zero /0.0d0/, Two /2.0d0/
       DATA Unknown /0.00999d0/
-C
+!
       Rk1 = Twomhb*radFitFlags%getMatchingRadius()/Rrrmas
-C *** Rk1 is the conversion factor from E to Rho=ka (a=radius=Crfn),
-C ***      so Rk1 is what we call Zkte elsewhere in the code
-C
+! *** Rk1 is the conversion factor from E to Rho=ka (a=radius=Crfn),
+! ***      so Rk1 is what we call Zkte elsewhere in the code
+!
       D_Energy = Unknown
-C *** By fiat, energies have no uncertainties
-C
+! *** By fiat, energies have no uncertainties
+!
       Kumur0 = 0
       Gggf = Zero
       Gggx = Zero
       Itimes = 0
-C
+!
       Num_E = 0
       Del_E = Zero
       Emaxa = Engur0
@@ -52,9 +53,9 @@ C
          Num_E = 1
       END IF
       DO Num=1,Number
-C
-C *** -----------------------------------------------------------
-C ***    Locate the first energy and the parameter-set Kumurr
+!
+! *** -----------------------------------------------------------
+! ***    Locate the first energy and the parameter-set Kumurr
          IF (Iwhich.EQ.1) THEN
             Kumurr = Num
             Emina = Emaxa
@@ -71,7 +72,7 @@ C ***    Locate the first energy and the parameter-set Kumurr
                Energy = Emina + Del_E/Two
             END IF
          ELSE
-C        ELSE IF (Iwhich.EQ.2) THEN
+!        ELSE IF (Iwhich.EQ.2) THEN
             Energy = Eeex(Num)
             DO I=1,Numurr
                IF (Energy.LT.Engurr(I)) THEN
@@ -81,83 +82,83 @@ C        ELSE IF (Iwhich.EQ.2) THEN
             END DO
    10       CONTINUE
          END IF
-C *** -----------------------------------------------------------
-C
-C ***    Initialize for this parameter region
-         Density_1 = Get_Density (Bindee(Kumurr),
-     *                            Pairee(Kumurr), Alevel(Kumurr))
+! *** -----------------------------------------------------------
+!
+! ***    Initialize for this parameter region
+         Density_1 = Get_Density (Bindee(Kumurr), &
+                                  Pairee(Kumurr), Alevel(Kumurr))
          D1b = Density_1*Bethdj(Jxx,Kumurr)
-C
+!
          Kountx = 0
-C ***    Loop over the energies in this region
+! ***    Loop over the energies in this region
          DO I_Energy=1,Num_E
-C
-C ***       Level spacing Dd
-            CALL Urr_Dd (Bindee, Pairee, Alevel, Jderiv, Xderiv, Energy,
-     *         Dd, D_Dd, D1b, Jxx, Iux, Ipar, Lcomp, Itimes, Kountx,
-     *         Kumurr)
-C
-C ***       Neutron widths
-            CALL Urr_Gn (Streng, Distnt, Iflstr, Ifldst, Jderiv, Xderiv,
-     *         Energy, Gggn, D_Gggn, Rk1, Dd, L, Ipar, Lcomp,
-     *         Itimes, Kumurr, Iux)
-C
-C ***       Capture widths
-            CALL Endep (Bindee(Kumurr), Pairee(Kumurr),
-     *         Alevel(Kumurr), Energy, Edgg, Kumurr)
+!
+! ***       Level spacing Dd
+            CALL Urr_Dd (Bindee, Pairee, Alevel, Jderiv, Xderiv, Energy, &
+               Dd, D_Dd, D1b, Jxx, Iux, Ipar, Lcomp, Itimes, Kountx, &
+               Kumurr)
+!
+! ***       Neutron widths
+            CALL Urr_Gn (Streng, Distnt, Iflstr, Ifldst, Jderiv, Xderiv, &
+               Energy, Gggn, D_Gggn, Rk1, Dd, L, Ipar, Lcomp, &
+               Itimes, Kumurr, Iux)
+!
+! ***       Capture widths
+            CALL Endep (Bindee(Kumurr), Pairee(Kumurr), &
+               Alevel(Kumurr), Energy, Edgg, Kumurr)
             Gggg = Gg(L,Kumurr)*Edgg
             IF (Iux.NE.0 .AND. (Lcomp.EQ.2 .OR. Itimes.EQ.0)) THEN
-               CALL Urr_Gg_Cov (Iflggg, Jderiv, Xderiv, Gggg,
-     *            D_Gggg, L, Ipar, Lcomp, Itimes, Kumurr)
+               CALL Urr_Gg_Cov (Iflggg, Jderiv, Xderiv, Gggg, &
+                  D_Gggg, L, Ipar, Lcomp, Itimes, Kumurr)
             END IF
-C
-C ***       Fission widths
+!
+! ***       Fission widths
             IF (If_Fis.EQ.1) THEN
-               CALL Urr_Gf (Gf, Ethr, Wthr, Iflgff, Jderiv, Xderiv,
-     *            Energy, Gggf, D_Gggf, L, Jxx, Ipar, Lcomp,
-     *            Itimes, Kumurr, Iux)
+               CALL Urr_Gf (Gf, Ethr, Wthr, Iflgff, Jderiv, Xderiv, &
+                  Energy, Gggf, D_Gggf, L, Jxx, Ipar, Lcomp, &
+                  Itimes, Kumurr, Iux)
             END IF
-C
-C ***       Inelastic widths = "competitive" width
+!
+! ***       Inelastic widths = "competitive" width
             IF (Numexc.GT.1) THEN
-               CALL Urr_Gx (Streng, Distnt, Ex, Spin, Pty, Iflstr,
-     *            Ifldst, Jderiv, Xderiv, Energy, Gggx, D_Gggx,
-     *            Rk1, Dd, Ipar, Lcomp, Itimes, Kumurr, Iux, J2,
-     *            Parity)
+               CALL Urr_Gx (Streng, Distnt, Ex, Spin, Pty, Iflstr, &
+                  Ifldst, Jderiv, Xderiv, Energy, Gggx, D_Gggx, &
+                  Rk1, Dd, Ipar, Lcomp, Itimes, Kumurr, Iux, J2, &
+                  Parity)
             END IF
-C
-C ***       Insert results into ENDF File 2 (and maybe File 32)
-            CALL Write_Urr_5 (Energy, Dd, Gggx, Gggn, Gggg, Gggf, Qj,
-     *          D_Energy, D_Dd, D_Gggx, D_Gggn, D_Gggg, D_Gggf,
-     *          Itimes, Mat, Lcomp, Iu, Iux)
-C
-C *** -----------------------------------------------------------
-C ***       Figure the next energy
+!
+! ***       Insert results into ENDF File 2 (and maybe File 32)
+            CALL Write_Urr_5 (Energy, Dd, Gggx, Gggn, Gggg, Gggf, Qj, &
+                D_Energy, D_Dd, D_Gggx, D_Gggn, D_Gggg, D_Gggf, &
+                Itimes, Mat, Lcomp, Iu, Iux)
+!
+! *** -----------------------------------------------------------
+! ***       Figure the next energy
             IF (Iwhich.EQ.1) THEN
-               IF (Kumurr.EQ.1 .AND. Numurr.GT.1 .AND. I_Energy.EQ.1
-     *            .AND. Kkends.EQ.1) THEN
+               IF (Kumurr.EQ.1 .AND. Numurr.GT.1 .AND. I_Energy.EQ.1 &
+                  .AND. Kkends.EQ.1) THEN
                   Energy = Energy + Del_E/Two
-               ELSE IF (Kumurr.EQ.Numurr .AND. Numurr.GT.1 .AND.
-     *            I_Energy.EQ.Num_Energy .AND. Kkends.EQ.1) THEN
+               ELSE IF (Kumurr.EQ.Numurr .AND. Numurr.GT.1 .AND. &
+                  I_Energy.EQ.Num_Energy .AND. Kkends.EQ.1) THEN
                   Energy = Emaxa
-C                 Energy = Energy + Del_E/Two
+!                 Energy = Energy + Del_E/Two
                ELSE
                   Energy = Energy + Del_E
                END IF
             ELSE
             END IF
-C *** -----------------------------------------------------------
+! *** -----------------------------------------------------------
          END DO
-C      # END DO I_Energy=1,Num_E
-C
+!      # END DO I_Energy=1,Num_E
+!
       END DO
       Num_Par = Ipar
       RETURN
       END
-C
-C
-C ______________________________________________________________________
-C
+!
+!
+! ______________________________________________________________________
+!
       SUBROUTINE Zurr (Dd, Jderiv, Xderiv)
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
       DIMENSION Jderiv(2), Xderiv(2)
@@ -169,50 +170,50 @@ C
       Xderiv(2) = Zero
       RETURN
       END
-C
-C
-C ______________________________________________________________________
-C
-      SUBROUTINE Urr_Dd (Bindee, Pairee, Alevel, Jderiv, Xderiv, Energy,
-     *   Dd, D_Dd, D1b, Jxx, Iux, Ipar, Lcomp, Itimes, Kountx, Kumurr)
-C *** Purpose -- generate level density and uncertainty
-C
+!
+!
+! ______________________________________________________________________
+!
+      SUBROUTINE Urr_Dd (Bindee, Pairee, Alevel, Jderiv, Xderiv, Energy, &
+         Dd, D_Dd, D1b, Jxx, Iux, Ipar, Lcomp, Itimes, Kountx, Kumurr)
+! *** Purpose -- generate level density and uncertainty
+!
       use fixedi_m
       use ifwrit_m
       use fixedr_m
       use constn_common_m
       use acs6_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
+!
       DIMENSION Bindee(*), Pairee(*), Alevel(*), Jderiv(2,*),Xderiv(2,*)
       DATA Unknown /0.00999d0/
-C
-      Density_2 = Get_Density (Bindee(Kumurr)+Energy,
-     *                         Pairee(Kumurr), Alevel(Kumurr))
+!
+      Density_2 = Get_Density (Bindee(Kumurr)+Energy, &
+                               Pairee(Kumurr), Alevel(Kumurr))
       F = D1b/Density_2
-C       = ratio of level densities (with other term Bethdj)
+!       = ratio of level densities (with other term Bethdj)
       Fj = Fjx (Bindee(Kumurr), Pairee(Kumurr), Energy, Jxx, Kountx)
       Dd = Fj * F
-C *** Dd is avg level spacing for (this J^pi) and for this Energy
-C *** It appears that Dd is independent of L
+! *** Dd is avg level spacing for (this J^pi) and for this Energy
+! *** It appears that Dd is independent of L
       IF (Iux.NE.0 .AND. (Lcomp.EQ.2 .OR. Itimes.EQ.0)) THEN
          Ipar = Ipar + 1
          CALL Zurr (D_Dd, Jderiv(1,Ipar), Xderiv(1,Ipar))
          D_Dd = Unknown
-C ***    Level spacings have no uncertainties ?
+! ***    Level spacings have no uncertainties ?
       END IF
       RETURN
       END
-C
-C
-C ______________________________________________________________________
-C
-      SUBROUTINE Urr_Gn (Streng, Distnt, Iflstr, Ifldst, Jderiv, Xderiv,
-     *   Energy, Gggn, D_Gggn, Rk1, Dd, L, Ipar, Lcomp, Itimes,
-     *   Kumurr, Iux)
-C
-C *** Purpose -- Generate values and uncertainties for elastic width
-C
+!
+!
+! ______________________________________________________________________
+!
+      SUBROUTINE Urr_Gn (Streng, Distnt, Iflstr, Ifldst, Jderiv, Xderiv, &
+         Energy, Gggn, D_Gggn, Rk1, Dd, L, Ipar, Lcomp, Itimes, &
+         Kumurr, Iux)
+!
+! *** Purpose -- Generate values and uncertainties for elastic width
+!
       use fixedi_m
       use ifwrit_m
       use fixedr_m
@@ -221,16 +222,16 @@ C
       use ResonanceCovariance_M
       use xxx5
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
-      DIMENSION Streng(Numelv,*), Distnt(Numelv,*), Iflstr(Numelv,*),
-     *   Ifldst(Numelv,*), Jderiv(2,*), Xderiv(2,*)
+!
+      DIMENSION Streng(Numelv,*), Distnt(Numelv,*), Iflstr(Numelv,*), &
+         Ifldst(Numelv,*), Jderiv(2,*), Xderiv(2,*)
       DIMENSION Dsde(5), Drde(5)
       type(ResonanceCovariance)::uCov
       DATA Dsde /-2.2d-11, 2.0d-11, -2.8d-11, 1.8d-11, 0.0d0/
       DATA Drde / 1.8E-07, 0.0d0  ,  2.0d-07, 0.0d0  , 0.0d0/
       DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
       DATA Unknown /0.00999d0/
-C
+!
       Rho = Rk1*dSQRT(Energy)
       Bnd = -dFLOAT(L-1)
       PL  = Pf (Rho, L-1)
@@ -241,19 +242,19 @@ C
          Rbarr = Rbarr + Energy*Drde(L)
          Rbari = Rbari + Energy*Halfpi*Dsde(L)/Rk1
       END IF
-C
+!
       RL0r = Rbarr*SL - Rbari*PL
       RL0i = Rbarr*PL + Rbari*SL
       Denom = (One-RL0r)**2 + (RL0i)**2
-C ***       = denominator of equation (VIIID.5) Pg 624 SAMMY/R7 manual
-C
-C *** factor of PL omitted in numerator because it is the REDUCED neutron
-C ***    width which is wanted (Gamma_n = Gamma_n^0 * sqrt(E)/rho * PL)
-C ***    where rho/sqrt(E) = Rk1
+! ***       = denominator of equation (VIIID.5) Pg 624 SAMMY/R7 manual
+!
+! *** factor of PL omitted in numerator because it is the REDUCED neutron
+! ***    width which is wanted (Gamma_n = Gamma_n^0 * sqrt(E)/rho * PL)
+! ***    where rho/sqrt(E) = Rk1
       Gggn = Streng(L,Kumurr)*Dd/Denom
-c January 25, 2005 attempt
+! January 25, 2005 attempt
 	Gggn = Dd * Streng(L,Kumurr)
-C
+!
       call covData%getUCovariance(uCov)
       IF (Iux.NE.0 .AND. (Lcomp.EQ.2 .OR. Itimes.EQ.0)) THEN
          Ipar = Ipar + 1
@@ -262,10 +263,10 @@ C
          D = Zero
          Jpar = 0
          IF (Iflstr(L,Kumurr).NE.0) THEN
-            D = Dd/Denom + Gggn/Denom * Pi/Rk1 * (
-     *                                         (One-RL0r)*PL - RL0i*SL )
+            D = Dd/Denom + Gggn/Denom * Pi/Rk1 * ( &
+                                               (One-RL0r)*PL - RL0i*SL )
             D = D * Streng(L,Kumurr)
-C                 * Streng because Streng=exp(Strenl)
+!                 * Streng because Streng=exp(Strenl)
             Jpar = Iflstr(L,Kumurr)
             Jderiv(Ik,Ipar) = Jpar
             Xderiv(Ik,Ipar) = D
@@ -302,14 +303,14 @@ C                 * Streng because Streng=exp(Strenl)
       END IF
       RETURN
       END
-C
-C
-C ______________________________________________________________________
-C
-      SUBROUTINE Urr_Gg_Cov (Iflggg, Jderiv, Xderiv,  Gggg,
-     *            D_Gggg, L, Ipar, Lcomp, Itimes, Kumurr)
-C *** Purpose -- generate uncertainties for gamma (capture) width
-C
+!
+!
+! ______________________________________________________________________
+!
+      SUBROUTINE Urr_Gg_Cov (Iflggg, Jderiv, Xderiv,  Gggg, &
+                  D_Gggg, L, Ipar, Lcomp, Itimes, Kumurr)
+! *** Purpose -- generate uncertainties for gamma (capture) width
+!
       use fixedi_m
       use ifwrit_m
       use fixedr_m
@@ -317,12 +318,12 @@ C
       use EndfData_common_m
       use ResonanceCovariance_M
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
+!
       DIMENSION Iflggg(Numelv,*), Jderiv(2,*), Xderiv(2,*)
       type(ResonanceCovariance)::uCov
       DATA Zero /0.0d0/
       DATA Unknown /0.00999d0/
-C
+!
       call covData%getUCovariance(uCov)
       Ipar = Ipar + 1
       IF (Iflggg(L,Kumurr).NE.0) THEN
@@ -330,7 +331,7 @@ C
          Jderiv(1,Ipar) = Jpar
          Jderiv(2,Ipar) = 0
          Xderiv(1,Ipar) = Gggg
-C                         Gggg not Edgg because Gg=exp(Gglg)
+!                         Gggg not Edgg because Gg=exp(Gglg)
          val = uCov%getCovariance(Jpar, Jpar)
          D_Gggg = Gggg * dSQRT(val)
          Xderiv(2,Ipar) = Zero
@@ -343,15 +344,15 @@ C                         Gggg not Edgg because Gg=exp(Gglg)
       END IF
       RETURN
       END
-C
-C
-C ______________________________________________________________________
-C
-      SUBROUTINE Urr_Gf (Gf, Ethr, Wthr, Iflgff, Jderiv, Xderiv,
-     *   Energy, Gggf, D_Gggf, L, Jxx, Ipar, Lcomp, Itimes,
-     *   Kumurr, Iux)
-C *** Purpose -- generate uncertainties for fission width
-C
+!
+!
+! ______________________________________________________________________
+!
+      SUBROUTINE Urr_Gf (Gf, Ethr, Wthr, Iflgff, Jderiv, Xderiv, &
+         Energy, Gggf, D_Gggf, L, Jxx, Ipar, Lcomp, Itimes, &
+         Kumurr, Iux)
+! *** Purpose -- generate uncertainties for fission width
+!
       use fixedi_m
       use ifwrit_m
       use fixedr_m
@@ -359,14 +360,14 @@ C
       use EndfData_common_m
       use ResonanceCovariance_M
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
-      DIMENSION Gf(Numelv,Numjjv,*), Ethr(Numelv,Numjjv,*),
-     *   Wthr(Numelv,Numjjv,*), Iflgff(Numelv,Numjjv,*),
-     *   Jderiv(2,*), Xderiv(2,*)
+!
+      DIMENSION Gf(Numelv,Numjjv,*), Ethr(Numelv,Numjjv,*), &
+         Wthr(Numelv,Numjjv,*), Iflgff(Numelv,Numjjv,*), &
+         Jderiv(2,*), Xderiv(2,*)
       type(ResonanceCovariance)::uCov
       DATA Zero /0.0d0/, One /1.0d0/
       DATA Unknown /0.00999d0/
-C
+!
       IF (Gf(L,Jxx,Kumurr).EQ.Zero) THEN
          Gggf = Zero
          IF (Iux.NE.0 .AND. (Lcomp.EQ.2 .OR. Itimes.EQ.0)) THEN
@@ -376,10 +377,10 @@ C
          END IF
       ELSE
          Gggf = Gf(L,Jxx,Kumurr)
-         A_Gf = (One+dEXP(         Ethr(L,Jxx,Kumurr)  /
-     *                             Wthr(L,Jxx,Kumurr)))
-     *         /(One+dEXP(-(Energy-Ethr(L,Jxx,Kumurr)) /
-     *                             Wthr(L,Jxx,Kumurr)))
+         A_Gf = (One+dEXP(         Ethr(L,Jxx,Kumurr)  / &
+                                   Wthr(L,Jxx,Kumurr))) &
+               /(One+dEXP(-(Energy-Ethr(L,Jxx,Kumurr)) / &
+                                   Wthr(L,Jxx,Kumurr)))
          Gggf = Gggf * A_Gf
          IF (Iux.NE.0 .AND. (Lcomp.EQ.2 .OR. Itimes.EQ.0)) THEN
             Ipar = Ipar + 1
@@ -403,16 +404,16 @@ C
       END IF
       RETURN
       END
-C
-C
-C ______________________________________________________________________
-C
-      SUBROUTINE Urr_Gx (Streng, Distnt, Ex, Spin, Pty, Iflstr, Ifldst,
-     *   Jderiv, Xderiv, Energy, Gggx, D_Gggx, Rk1, Dd, Ipar,
-     *   Lcomp, Itimes, Kumurr, Iux, J2, Parity)
-C
-C *** Purpose -- generate uncertainties for inelastic width
-C
+!
+!
+! ______________________________________________________________________
+!
+      SUBROUTINE Urr_Gx (Streng, Distnt, Ex, Spin, Pty, Iflstr, Ifldst, &
+         Jderiv, Xderiv, Energy, Gggx, D_Gggx, Rk1, Dd, Ipar, &
+         Lcomp, Itimes, Kumurr, Iux, J2, Parity)
+!
+! *** Purpose -- generate uncertainties for inelastic width
+!
       use fixedi_m
       use ifwrit_m
       use fixedr_m
@@ -421,14 +422,14 @@ C
       use ResonanceCovariance_M
       use acs6_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
-      DIMENSION Streng(Numelv,*), Distnt(Numelv,*), Ex(*), Spin(*),
-     *   Pty(*), Iflstr(Numelv,*), Ifldst(Numelv,*), Jderiv(2,*),
-     *   Xderiv(2,*)
+!
+      DIMENSION Streng(Numelv,*), Distnt(Numelv,*), Ex(*), Spin(*), &
+         Pty(*), Iflstr(Numelv,*), Ifldst(Numelv,*), Jderiv(2,*), &
+         Xderiv(2,*)
       type(ResonanceCovariance)::uCov
       DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
       DATA Unknown /0.00999d0/
-C
+!
       call covData%getUCovariance(uCov)
       Gggx = Zero
       D_Gggx = Zero
@@ -444,23 +445,23 @@ C
             Is2fx =      I2f+1
             Llo = 1
             IF (Pty(Kumexc)*Parity.LT.Zero) Llo = 2
-C
+!
             Emex = Energy - Ex(Kumexc)
             Rhox = Rk1*dSQRT(Emex)
-C ***       LOOP OF OUTGOING (INELASTIC) PARTIAL WAVES:
-C ***          ALLOW FOR POSSIBLE PARITY CHANGE
+! ***       LOOP OF OUTGOING (INELASTIC) PARTIAL WAVES:
+! ***          ALLOW FOR POSSIBLE PARITY CHANGE
             IF (Llo.LE.Numelv) THEN
                DO Lx=Llo,Numelv,2
                   Dof = Zero
                   Kk = Kwhich (J2, Lx, Dof, Is2fn, Is2fx)
-C                 Kwhich determines value of Dof = degree-of-freedom
+!                 Kwhich determines value of Dof = degree-of-freedom
                   IF (Lx.GT.0 .AND. Dof.GT.Zero) THEN
                      LLL = Lx - 1
                      PL  = Pf (Rhox, LLL)
                      Gggxx = Dd * PL * Streng(Lx,Kumurr) / Rk1
                      Gggx  = Gggx + Gggxx*Dof
-C
-C ***             Derivative calculation is certainly not correct yet
+!
+! ***             Derivative calculation is certainly not correct yet
 	go to 99
                   IF (Iux.NE.0 .AND. (Lcomp.EQ.2 .OR. Itimes.EQ.0)) THEN
                      Ipar = Ipar + 1
@@ -468,10 +469,10 @@ C ***             Derivative calculation is certainly not correct yet
                      D_Gggxx = Zero
                      CALL Zurr (D_Gggxx, Jderiv(1,Ipar), Xderiv(1,Ipar))
                      IF (Iflstr(Lx,Kumurr).NE.0) THEN
-                        D = Dd/Denom + Gggxx/Denom * Pi/Rk1 * (
-     *                                         (One-RL0r)*PL - RL0i*SL )
+                        D = Dd/Denom + Gggxx/Denom * Pi/Rk1 * ( &
+                                               (One-RL0r)*PL - RL0i*SL )
                         D = D * Streng(Lx,Kumurr)
-C                             * Streng because Streng=exp(Strenl)
+!                             * Streng because Streng=exp(Strenl)
                         Jpar = Iflstr(Lx,Kumurr)
                         Jderiv(Ik,Ipar) = Jpar
                        Xd = D
@@ -497,7 +498,7 @@ C                             * Streng because Streng=exp(Strenl)
                            D_Gggxx = val*C*D + D_Gggxx
                         END IF
                      END IF
-C                  # END IF (Ifldst(Lx,Kumurr).NE.0) THEN
+!                  # END IF (Ifldst(Lx,Kumurr).NE.0) THEN
                      IF (D_Gggxx.GT.Zero) THEN
                         D_Gggxx = dSQRT(D_Gggxx) * PL
                      ELSE
@@ -509,32 +510,35 @@ C                  # END IF (Ifldst(Lx,Kumurr).NE.0) THEN
                   END IF
                END DO
             END IF
-C         # END IF (Llo.LE.Numelv) THEN
+!         # END IF (Llo.LE.Numelv) THEN
          END DO
       END IF
       RETURN
       END
-C
-C
-C ______________________________________________________________________
-C
-      SUBROUTINE Write_Urr_5 (Energy, Dd, Gggx, Gggn, Gggg, Gggf, Qj,
-     *   D_Energy, D_Dd, D_Gggx, D_Gggn, D_Gggg, D_Gggf, Itimes, Mat,
-     *   Lcomp, Iu, Iux)
+!
+!
+! ______________________________________________________________________
+!
+      SUBROUTINE Write_Urr_5 (Energy, Dd, Gggx, Gggn, Gggg, Gggf, Qj, &
+         D_Energy, D_Dd, D_Gggx, D_Gggn, D_Gggg, D_Gggf, Itimes, Mat, &
+         Lcomp, Iu, Iux)
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
       DATA Mf /2/, No /0/, Mfx /32/, Mt /151/
       Ns = 1
-      CALL List_7x (Energy, Dd, Gggx, Gggn, Gggg, Gggf,
-     *   Mat, Mf, Mt, Ns, No, No, No, No, No, No, Iu)
+      CALL List_7x (Energy, Dd, Gggx, Gggn, Gggg, Gggf, &
+         Mat, Mf, Mt, Ns, No, No, No, No, No, No, Iu)
       IF (Iux.NE.0) THEN
          IF (Lcomp.EQ.2) THEN
-            CALL List_7x (D_Energy, D_Dd, D_Gggx, D_Gggn, D_Gggg,
-     *         D_Gggf, Mat, Mfx, Mt, Ns, No, No, No, No, No, No, Iux)
+            CALL List_7x (D_Energy, D_Dd, D_Gggx, D_Gggn, D_Gggg, &
+               D_Gggf, Mat, Mfx, Mt, Ns, No, No, No, No, No, No, Iux)
          ELSE IF (Lcomp.EQ.1 .AND. Itimes.EQ.0) THEN
             Itimes = 1
-            CALL List_7x (Dd, Qj, Gggn, Gggg, Gggf, D_Gggx,
-     *         Mat, Mfx, Mt, Ns, No, No, No, No, No, No, Iux)
+            CALL List_7x (Dd, Qj, Gggn, Gggg, Gggf, D_Gggx, &
+               Mat, Mfx, Mt, Ns, No, No, No, No, No, No, Iux)
          END IF
       END IF
       RETURN
       END
+
+end module lru2_m
+
diff --git a/sammy/src/lru/mlru3.f b/sammy/src/lru/mlru3.f90
similarity index 75%
rename from sammy/src/lru/mlru3.f
rename to sammy/src/lru/mlru3.f90
index 32ffeee65d2f3cd059bf091b909f9339d524ad17..e11e295b55e64af2ea7738cec9ef537b98fafbc6 100755
--- a/sammy/src/lru/mlru3.f
+++ b/sammy/src/lru/mlru3.f90
@@ -1,18 +1,19 @@
-C
-C
-C ______________________________________________________________
-C
+module lru3_m
+   contains
+!
+! ______________________________________________________________
+!
       SUBROUTINE Fix_Covar (Jderiv, Xderiv, Covar, Num_Par)
-C
-C *** PURPOSE -- Determine Covar = covariance matrix for ENDF parameters
-C
+!
+! *** PURPOSE -- Determine Covar = covariance matrix for ENDF parameters
+!
       use EndfData_common_m
       use ResonanceCovariance_M
       IMPLICIT DOUBLE PRECISION (A-h,o-z)
       DIMENSION Jderiv(2,*), Xderiv(2,*), Covar(*)
       type(ResonanceCovariance)::uCov
       DATA Zero /0.0d0/
-C
+!
       call covData%getUCovariance(uCov)
       N = (Num_Par*(Num_Par+1))/2
       CALL Zero_Array (Covar, N)
@@ -54,29 +55,29 @@ C
       END DO
       RETURN
       END
-C
-C
-C --------------------------------------------------------------
-C
+!
+!
+! --------------------------------------------------------------
+!
       SUBROUTINE Write_Cmpc_Cov_2 (Covar, Diag, Num_Par, Mat, Ns)
-C
-C *** PURPOSE -- Write covariance matrix in compact format
-C *** INPUT   -- Diag   = dummy array of size (Num_Par)
-C ***         -- Num_Par= number of elements of full array
-C ***         -- Covar  = triangular array containing covariance matrix
-C ***                   1
-C ***                   2  3
-C ***                   4  5  6
-C ***                   7  8  9 10
-C ***                  11 12 13 14 15
-C ***                  16 17 18 19 20 21
-C ***                  ... up to (Num_Par*(Num_Par+1))/2
-C ***
-C *** OUTPUT  -- Diag contains sqrt of diagonal elements of Covar
-C ***         -- Output File contains abbreviated correlation matrix
-C ***
-C *** AUTHOR  -- N.M.Larson, July 2003
-C
+!
+! *** PURPOSE -- Write covariance matrix in compact format
+! *** INPUT   -- Diag   = dummy array of size (Num_Par)
+! ***         -- Num_Par= number of elements of full array
+! ***         -- Covar  = triangular array containing covariance matrix
+! ***                   1
+! ***                   2  3
+! ***                   4  5  6
+! ***                   7  8  9 10
+! ***                  11 12 13 14 15
+! ***                  16 17 18 19 20 21
+! ***                  ... up to (Num_Par*(Num_Par+1))/2
+! ***
+! *** OUTPUT  -- Diag contains sqrt of diagonal elements of Covar
+! ***         -- Output File contains abbreviated correlation matrix
+! ***
+! *** AUTHOR  -- N.M.Larson, July 2003
+!
       use samxxx_common_m
       use mold7_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
@@ -85,10 +86,10 @@ C
       DATA Bl /' '/, Minus /'-'/
       DATA Zero/0.0d0/, One /1.0d0/, Hundred /100.0d0/
       DATA Line /18/, Mfx /32/, Mt /151/
-C
+!
       Nn = Ns
       CALL Newopn (57, Sam57x, 0)
-C
+!
       Ii = 0
       DO Ipar=1,Num_Par
          Ii = Ii + Ipar
@@ -98,8 +99,8 @@ C
             Diag(Ipar) = Zero
          END IF
       END DO
-C
-C *** Convert from covariance to correlation*100
+!
+! *** Convert from covariance to correlation*100
       Iimax = 0
       Ij = 0
       DO Ipar=1,Num_Par
@@ -112,16 +113,15 @@ C *** Convert from covariance to correlation*100
                IF (Diag(Jpar).GT.Zero) THEN
                   Covar(Ij) = Covar(Ij)/(Di*Diag(Jpar)) * Hundred
                ELSE IF (Covar(Ij).NE.Zero) THEN
-                  WRITE (6,30000) Ipar, Jpar, Ij, Covar(Ij), Diag(Jpar),
-     *               Di
+                  WRITE (6,30000) Ipar, Jpar, Ij, Covar(Ij), Diag(Jpar), Di
 30000             FORMAT ('Ipar,Jpar,Ij,Covar=', 3I5, 1P5G14.6)
-cx                  STOP
+!x                  STOP
                END IF
             END IF
          END DO
       END DO
-C
-C
+!
+!
       Ii = 0
       Ij = 0
       DO Ipar=1,Num_Par
@@ -134,7 +134,7 @@ C
                IF (Jpar.LT.Ipar) THEN
                   Vij = dABS(Covar(Ij))
                   IF (Vij.GT.One) THEN
-C
+!
                      Im = Ij + Line - 1
                      IF (Im.GE.Ii) Im = Ii - 1
                      Mk = Jpar + Line - 1
@@ -168,8 +168,8 @@ C
          Nsw = Ns/100000
          Ns = Ns - 100000*Nsw
       END IF
-                     WRITE (57,10300) Ipar, Jpar,
-     *                  ((Ch(J,K),J=1,3),K=1,Line), Mat, Mfx, Mt, Ns
+                     WRITE (57,10300) Ipar, Jpar, &
+                        ((Ch(J,K),J=1,3),K=1,Line), Mat, Mfx, Mt, Ns
 10300                FORMAT (I5, I5, 1X, 54A1, 1X, I4, I2, I3, I5)
                      IF (Im.GE.Ii-1) GO TO 70
                   END IF
@@ -193,36 +193,36 @@ C
       CLOSE (UNIT=57)
       RETURN
       END
-C
-C
-C --------------------------------------------------------------
-C
-      SUBROUTINE Write_Urr_Lcomp_1 (Covar, Nnpar, Num_Par, If_Fis,
-     *   Mat, Ns)
-C
-C *** PURPOSE -- Write covariance matrix in Lcomp="undefined" format
-C *** INPUT   -- Nnpar = number of elements of arrays
-C ***                  = number of sammy-varied parameters
-C ***         -- Covar = triangular array containing RELATIVE covariance matrix
-C ***                   1
-C ***                   2  3
-C ***                   4  5  6
-C ***                   7  8  9 10
-C ***                  11 12 13 14 15
-C ***                  16 17 18 19 20 21
-C ***                  ... up to (Nnpar*(Nnpar+1))/2
-C ***            Note that ordering is opposite to that wanted in ENDF
-C ***
-C *** OUTPUT  -- Output File contains ENDF Lcomp=undefined covariance matrix
-C
-C *** AUTHOR  -- N.M.Larson
-C ***            This version prints correlation matrix, assuming energy-
-C ***               independence of all parameters.  Ergo, SAMMY just 
-C ***               takes the lowest-energy values
-C ***            Diagonal element = 0.000001 for unvaried parameters 
-C ***               (convention only, could be changed)
-C
-C
+!
+!
+! --------------------------------------------------------------
+!
+      SUBROUTINE Write_Urr_Lcomp_1 (Covar, Nnpar, Num_Par, If_Fis, &
+         Mat, Ns)
+!
+! *** PURPOSE -- Write covariance matrix in Lcomp="undefined" format
+! *** INPUT   -- Nnpar = number of elements of arrays
+! ***                  = number of sammy-varied parameters
+! ***         -- Covar = triangular array containing RELATIVE covariance matrix
+! ***                   1
+! ***                   2  3
+! ***                   4  5  6
+! ***                   7  8  9 10
+! ***                  11 12 13 14 15
+! ***                  16 17 18 19 20 21
+! ***                  ... up to (Nnpar*(Nnpar+1))/2
+! ***            Note that ordering is opposite to that wanted in ENDF
+! ***
+! *** OUTPUT  -- Output File contains ENDF Lcomp=undefined covariance matrix
+!
+! *** AUTHOR  -- N.M.Larson
+! ***            This version prints correlation matrix, assuming energy-
+! ***               independence of all parameters.  Ergo, SAMMY just 
+! ***               takes the lowest-energy values
+! ***            Diagonal element = 0.000001 for unvaried parameters 
+! ***               (convention only, could be changed)
+!
+!
       use samxxx_common_m
       use EndfData_common_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
@@ -231,17 +231,17 @@ C
       DIMENSION Dum(6)
       DATA Zero/0.0d0/
       DATA Mfx /32/, Mt /151/
-C
+!
       Nn = Ns
       CALL Newopn (57, Sam57x, 0)
-C
+!
       Nvpall = Nnpar
       call covData%setNumParameters(Nnpar)
-C
+!
       DO Idum=1,6
          Dum(Idum) = Zero
       END DO
-C
+!
       Ijmin = 0
       Idum = 0
       DO Ipar=1,Num_Par
@@ -274,16 +274,16 @@ C
       Nm = Ns - Nn
       X = Zero
       N = 0
-C
-C *** Write next line of File 32
+!
+! *** Write next line of File 32
       Npar = Num_Par
       Nx   = (Npar*(Npar+1))/2
       Mpar = 3
       IF (If_Fis.EQ.1) Mpar = 4
       WRITE (55,10400) X, X, Mpar, N, Nx, Npar, Mat, Mfx, Mt, Nn
 10400 FORMAT (2F11.1, 4I11, I4, I2, I3, I5)
-C
-C *** Read, massage, and rewrite the rest of the dummy file into File 32
+!
+! *** Read, massage, and rewrite the rest of the dummy file into File 32
       DO I=1,Nm
          READ (57,10500) Aaaa
 10500    FORMAT (78A1)
@@ -300,3 +300,5 @@ C *** Read, massage, and rewrite the rest of the dummy file into File 32
       CLOSE (UNIT=55)
       RETURN
       END
+
+end module lru3_m
diff --git a/sammy/src/salmon/CMakeLists.txt b/sammy/src/salmon/CMakeLists.txt
index 642c5f3730079d2210f162c06688a16b8910af9d..da0b358ac2bbb75eb9671ff5b790e4764c024b69 100644
--- a/sammy/src/salmon/CMakeLists.txt
+++ b/sammy/src/salmon/CMakeLists.txt
@@ -32,6 +32,8 @@ SET(HEADERS ExperimentalParameters.h
 
             DerivativeHandler.h
             interface/cpp/DerivativeHandlerInterface.h
+            
+            CovToolBox.h
 )
 
 APPEND_SET(SALMON_SOURCES
@@ -66,10 +68,14 @@ APPEND_SET(SALMON_SOURCES
         interface/fortran/DerivativeListHolder_M.f90
         interface/fortran/DerivativeListHolder_I.f90
 
+<<<<<<< HEAD
         DerivativeHandler.cpp
         interface/cpp/DerivativeHandlerInterface.cpp
         interface/fortran/DerivativeHandler_I.f90
         interface/fortran/DerivativeHandler_M.f90
+=======
+        CovToolBox.cpp
+>>>>>>> d5ff02b1... add cov memory management for file33 writing
 )
 
 
diff --git a/sammy/src/salmon/CovToolBox.cpp b/sammy/src/salmon/CovToolBox.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9fb03349714e6966c1003f405a2addb6c1a56776
--- /dev/null
+++ b/sammy/src/salmon/CovToolBox.cpp
@@ -0,0 +1,51 @@
+#include "CovToolBox.h"
+
+namespace sammy{
+
+    CovToolBox::CovToolBox() {}
+
+    /**
+     * Return whether integer `i` is less than integer `j`
+     * 
+     * @param i 1st test integer
+     * @param j 2nd test integer
+     * */
+    bool ascending(int i,int j) { return (i<j); }
+
+
+    void CovToolBox::fillCrossCovFromCoverx(coverx::CoverxFile & cf, 
+                                            endf::CrossCovList & cc){
+
+        int numCross = cf.getNumCross();
+        std::vector<int> mtList;
+        for( int i=0;i<numCross;++i ){ 
+            mtList.push_back(cf.getCrossAt(i)->getMt());
+        }
+        std::sort(mtList.begin(),mtList.end(),ascending);
+
+        // single MAT number, will have to cycle thru to expand to 
+        // multiple isotopes
+        int mat = cf.getCrossAt(0)->getMat();
+        int numE = cf.getCrossAt(0)->getSize();
+        std::vector<double> eGroups;
+        for( int i=0;i<=numE;++i ) eGroups.push_back(cf.getNeutGroupBounds(i));
+        // for every reaction
+        for( int i=0;i<numCross;++i ){
+            int mt1 = mtList[i];
+            // for every 2nd reaction (j<=i)
+            for( int j=0;j<numCross;++j){
+                if(i>j) continue;
+                int mt2 = mtList[j];
+                endf::CrossCovDirect * cov = new  endf::CrossCovDirect(mat,mt1,mat,mt2);
+                cov->setEnergy(numE,eGroups.data());
+                for( int iE=0;iE<numE;++iE ){
+                    for( int jE=0;jE<numE;++jE ){
+                        cov->setMatrixElement(iE,jE,
+                                cf.getMatrix(mat,mt1,mat,mt2)->getValue(iE,jE));
+                    }
+                }
+                cc.addCrossCovDirect(cov);
+            }
+        }
+    }
+}
\ No newline at end of file
diff --git a/sammy/src/salmon/CovToolBox.h b/sammy/src/salmon/CovToolBox.h
new file mode 100644
index 0000000000000000000000000000000000000000..7a17de76ff36be646853aba4abfb02ea467034f3
--- /dev/null
+++ b/sammy/src/salmon/CovToolBox.h
@@ -0,0 +1,38 @@
+
+
+#ifndef SAMMY_COVTOOLBOX_H
+#define SAMMY_COVTOOLBOX_H
+
+#include "ScaleData/Core/io/CoverxFile.h"
+#include "ScaleUtils/EndfLib/CrossSectionCovariance.h"
+#include <iostream>
+#include <algorithm>
+
+namespace sammy {
+
+    class CovToolBox{
+    public:
+        CovToolBox();
+        virtual ~CovToolBox() = default;
+
+        /***********************************************************************
+         * Pull data from CoverxFile object and fill a CrossCovDirect
+         * object. Helps separate EndfLib package from Coverx dependency.
+         * For now we only support one material (MAT), so we have MT^2
+         * many matrices in the CoverxFile. CoverxFile *must* have MAT and 
+         * reaction (MT) numbers associated with matrices, as we need to 
+         * retrieve matrices in specific MAT/MT order.
+         * 
+         * @param cf the object containing cross section covariance info
+         * @param cc the vector of objects we want to fill with cov info
+         **********************************************************************/
+        void fillCrossCovFromCoverx(coverx::CoverxFile & cf, 
+                                    endf::CrossCovList & cc);
+
+    private:
+
+    };
+}
+
+#endif  /* SAMMY_COVTOOLBOX_H */
+
diff --git a/sammy/src/salmon/GridData.h b/sammy/src/salmon/GridData.h
index f25064159532f5a387759165fa6e79a592591b16..5a26e7149b439fec645603508032cb17dcbc3fd0 100644
--- a/sammy/src/salmon/GridData.h
+++ b/sammy/src/salmon/GridData.h
@@ -81,6 +81,8 @@ namespace sammy{
 
       /**
        * If true, addData accumulates data instead of set data
+       * TODO: discuss whether we should break `addData()` into
+       * two functions: `addData()` and `addToCurrentData()`
        *
        * @param add set to true to accumulate data with addData
        */
diff --git a/sammy/src/sammy/CMakeLists.txt b/sammy/src/sammy/CMakeLists.txt
index bfb4d77a732b69e18e698c24434f6df2d5c6a54f..305fa27eeb711369ad28eed13ad73239edc58a30 100644
--- a/sammy/src/sammy/CMakeLists.txt
+++ b/sammy/src/sammy/CMakeLists.txt
@@ -154,7 +154,7 @@ APPEND_SET(SAMMY_SOURCES
             ../fit/mfit0.f
             ../fit/mfit1.f90
             ../fit/mfit2.f
-            ../fit/mfit3.f
+            ../fit/mfit3.f90
 
             ../fnc/aaaerf.f90
             ../fnc/bbberf.f90
@@ -224,10 +224,10 @@ APPEND_SET(SAMMY_SOURCES
             ../ipq/sammy_ipq_m.f90
             ../ipq/FillFitData_m.f90
 
-            ../lru/mlru0.f
-            ../lru/mlru1.f
-            ../lru/mlru2.f
-            ../lru/mlru3.f
+            ../lru/mlru0.f90
+            ../lru/mlru1.f90
+            ../lru/mlru2.f90
+            ../lru/mlru3.f90
 
             ../mas/mmas0.f90
             ../mas/mmas1.f90