diff --git a/sammy/src/amr/mamr.f b/sammy/src/amr/mamr.f
index 2d9ac08b566162701e3a0bbee64fa4c6d63d2330..c9445bb7f2ba8de6761b3c18f63b166a3069dec0 100644
--- a/sammy/src/amr/mamr.f
+++ b/sammy/src/amr/mamr.f
@@ -42,7 +42,6 @@ C
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
       
       real(kind=8) ::  A(5000000)
-      real(kind=8),allocatable,dimension(:)::tmpUCov
 C                       Msize below must be the same as this number
 C
       CHARACTER*3 Action, Type
@@ -78,7 +77,6 @@ C
       Iu62 = 62
 C *** READ FIRST TWO LINES OF COVARIANCE FILE TO GET DIMENSIONS
       CALL If_Read_Covar (Iu62, Ifold, 1, Nvpall_Old)
-      allocate(tmpUCov(Nvpall_Old))
       CALL Startx
 C
 C *** Fix array sizes
@@ -91,7 +89,6 @@ C
       IF (K.EQ.0) K = 1
       Nt = (K*(K+1))/2
       Ivrpr  = Idimen (Nt, 1, 'Vrpr    Nt, 1')
-      Ivvv   = Idimen (K , 1, 'Vvv     K , 1')
       Kompct = 0
       Iuncer = Idimen (K, 1, 'K, 1')
 C
@@ -141,7 +138,7 @@ C ***         VALUES, ETC.
      2      A_Iprbgf , I_Iflbgf , A_Indbgf , A_Ibgfmi , A_Ibgfma ,
      3      A_Itexbg , A_Iteabg , A_Iprdtp , I_Ifldtp ,
      4      A_Iprusd , I_Iflusd , A_Iprbag , I_Iflbag ,
-     *      tmpUCov  , A(Ivrpr ), A(Iallva), A(Ivvv  ), A(Ipoint),
+     *      A(Ivrpr ), A(Iallva), A(Ipoint),
      *      A(Idum1 ), A(Idum2 ), A(Idiot1), A(Idiot2),
      *      A(Ivsqua), A(Ieival), A(Ieivec), A(Izeiv ), A(Ideiv ),
      *      A(Iiord ), A(Ikord ),
@@ -162,11 +159,10 @@ C ***         VALUES, ETC.
      *      A_Iprbgf , I_Iflbgf , A_Indbgf , A_Ibgfmi , A_Ibgfma ,
      *      A_Iprdtp , I_Ifldtp ,
      *      A_Iprusd , I_Iflusd , A_Iprbag , I_Iflbag ,
-     *      tmpUCov  , A(Ia    ), A(Ia    ), A(Ipoint),
+     *      A(Ia    ), A(Ia    ), A(Ipoint),
      *      A(Isingl), Krext, Nowrt, Nsingl, Nvpall_Old)
 C
       END IF
-      deallocate(tmpUCov)
 C
 C
 C *** get user-supplied parameter information
diff --git a/sammy/src/amx/mamx.f b/sammy/src/amx/mamx.f
index d91529c283d3ab556f574c9912de81be287e43fc..86c182066cc28e05ad71250d8a249561d4911062 100644
--- a/sammy/src/amx/mamx.f
+++ b/sammy/src/amx/mamx.f
@@ -88,7 +88,6 @@ C
       Nt = (K*(K+1))/2
       Ivrpr  = Idimen (Nt, 1, 'Nt, 1')
       Iallva = Idimen (Nt, 1, 'Nt, 1')
-      Ivvv   = Idimen (K , 1, 'K , 1')
       Iuncer = Idimen (K , 1, 'K , 1')
       Kompct = 0
 C
@@ -139,7 +138,7 @@ C &&&    from old/mold2.f
      *      A_Iprbgf , I_Iflbgf , A_Indbgf , A_Ibgfmi , A_Ibgfma ,
      *      A_Itexbg , A_Iteabg , A_Iprdtp , I_Ifldtp ,
      *      A_Iprusd , I_Iflusd , A_Iprbag , I_Iflbag ,
-     *      A(Iu    ), A(Ivrpr ), A(Iallva), A(Ivvv  ), A(Ipoint),
+     *      A(Ivrpr ), A(Iallva), A(Ipoint),
      *      A(Idum1 ), A(Idum2 ), A(Idiot1), A(Idiot2),
      *      A(Ivsqua), A(Ieival), A(Ieivec), A(Izeiv ), A(Ideiv ),
      *      A(Iiord ), A(Ikord ),
@@ -161,7 +160,7 @@ c &&&    from old/mold5.f
      *      A_Iprbgf , I_Iflbgf , A_Indbgf , A_Ibgfmi , A_Ibgfma ,
      *      A_Iprdtp , I_Ifldtp ,
      *      A_Iprusd , I_Iflusd , A_Iprbag , I_Iflbag ,
-     *      A(Iu)    , A(Ia    ), A(Ia    ), A(Ipoint),
+     *      A(Ia    ), A(Ia    ), A(Ipoint),
      *      A(Isingl), Krext, Nowrt, Nsingl, Nvpall_Old)
 C
       END IF
diff --git a/sammy/src/odf/modf0.f b/sammy/src/odf/modf0.f
index 0a0eb366609b2552eab76d03647fc0259d7a28bf..d53be5d121feaf99938468a9f4ee55dd0579813e 100644
--- a/sammy/src/odf/modf0.f
+++ b/sammy/src/odf/modf0.f
@@ -13,9 +13,10 @@ C
       use namfil_common_m
       use aaaodf_common_m
       use ffnnnn_common_m
+      use AllocateFunctions_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      EXTERNAL Idimen
       DIMENSION A(-Msize:Msize)
+      real(kind=8),allocatable,dimension(:)::A_Ie, A_id, A_Iu, A_Idum
 C
 C
       WRITE (6,99999)
@@ -40,15 +41,18 @@ C *** Reads file SAM18.DAT for real, also scans the data file
 C
       Nn = Ndat
       IF (Nangle.GT.0) Nn = Ndat*Nangle
-      I  = Idimen (1, 1, '1, 1')
-      Ie = Idimen (Ndat, 1, 'Ndat, 1')
-      Id = Idimen (Nn, 1, 'Nn, 1')
-      Iu = Idimen (Nn, 1, 'Nn, 1')
+      call allocate_real_data(A_Ie,  Ndat)
+      call allocate_real_data(A_Id,  Nn)
+      call allocate_real_data(A_Iu,  Nn)
       N = (Ndat+1)
-      Idum = Idimen (N, 1, 'N, 1')
-      CALL Fixodf (A(Ie), A(Id), A(Iu), A(Idum), Ktzero, Nmax)
+      call allocate_real_data(A_Idum,N)
+      CALL Fixodf (A_Ie, A_Id, A_Iu, A_Idum, Ktzero, Nmax)
 C *** Reads data files (Unit 13) and generates ODF files (unit 72 etc)
-      I = Idimen (Ie, -1, 'Ie, -1')
+
+      deallocate(A_Ie)
+      deallocate(A_Id)
+      deallocate(A_Iu)
+      deallocate(A_Idum)
 C
       CLOSE (UNIT=18)
       Ix = 0
diff --git a/sammy/src/old/ReadCovarianceParts.f90 b/sammy/src/old/ReadCovarianceParts.f90
index a42901a40b07827e13fee228ec993819f2a7f5b2..c61efe1e3fe5768badcfe990a01bd4f48b6f9cd3 100644
--- a/sammy/src/old/ReadCovarianceParts.f90
+++ b/sammy/src/old/ReadCovarianceParts.f90
@@ -15,19 +15,43 @@ contains
             
               
   !  Drop correlation if uncertainty is large             
-  subroutine Xdrop2Modern (cov, diag)
-     real(kind=8)::Diag(*), Droppp
-     type(ResonanceCovariance)::cov     
-          
-     integer::npar, irow, icol     
+  subroutine Xdrop2Modern (cov, Droppp)
+     use EndfData_common_m
+     real(kind=8)::Droppp
+     type(ResonanceCovariance)::cov
+
+     integer::npar, irow, icol
+     type(ResonanceCovariance)::uCov
+     real(kind=8)::Dr, Di, V
           
      npar = cov%getNumRows()
+     call covData%getUCovariance(uCov)
      if (npar.le.1) return
      do irow = 1, npar
-        do icol = 1, irow - 1
-          if( diag(irow).eq.0.0d0.or.diag(icol).eq.0.0d0) then
-              call cov%setCovariance(irow, icol, 0.0d0)
-           end if
+         Dr =  uCov%getCovariance(irow, irow)
+         Dr =  dsqrt(Dr)
+         V = dAbs(covData%getUParamValue(irow))
+         if (V.ne.0.0) then
+            Dr  = Dr/V
+         else
+             Dr = 0.0d0
+         end if
+         if( Dr.GT.Droppp) Dr = 0.0d0   ! Dr was previously set in Xdrop1
+
+         do icol = 1, irow - 1
+            Di =  uCov%getCovariance(icol, icol)
+            Di =  dsqrt(Di)
+            V = dAbs(covData%getUParamValue(icol))
+            if (V.ne.0.0) then
+                Di  = Di/V
+            else
+                Di = 0.0d0
+            end if
+            if( Di.GT.Droppp) Dr = 0.0d0 ! Di was previously set in Xdrop1
+
+            if( Dr.eq.0.0d0.or.Di.eq.0.0d0) then
+                call covData%setCovarianceData( ucov, irow, icol, 0.0d0)
+            end if
         end do
      end do
    end subroutine Xdrop2Modern
@@ -98,14 +122,12 @@ contains
   !! @param old if true data are in single precision
   !! @param Nvpall_Old parameters to consider
   !! @param Allvr covariance matrix for all physical parameters (p-parameters) 
-  !! @param U   values for the u parameters
-  !! @param Vvv   Array to write correlation matrix (of Allvr) if Kompct.eq.2
   !! @param Vrpr covariance matrix for U parameters 
   !! @param Nowrt   flag  to indicate whether data in the correlation matrix should be dropped 
   !! @param idropp if dropping correlation the lowest correlation to drop (to be converted to float by multiplying by 0,01
   !! @param Kompct  option to print and compact correlation matrix
   !!
-  subroutine readHeaderInfo(old, Nvpall_Old, Allvr, U, Vvv,   Vrpr, Nowrt, idropp, Kompct)
+  subroutine readHeaderInfo(old, Nvpall_Old, Allvr, Vrpr, Nowrt, idropp, Kompct)
       use fixedi_m
       use fixedr_m
       use samxxx_common_m
@@ -114,13 +136,12 @@ contains
       integer::Nsingl, Nvpall_Old, nowrt, idropp, Kompct
       logical::old
     
-      integer::Ii,Ij,  I, IPar, JPar, nzero
-      real(kind=8)::Allvr(*), U(*), Vvv(*), Vrpr(*)
+      integer::Ij,  I, IPar, JPar, nzero
+      real(kind=8)::Allvr(*), Vrpr(*)
+      real(kind=8),allocatable, dimension(:)::Dummy
       real(kind=8)::Droppp, Zero = 0.0d0
-      logical::changed
       type(ResonanceCovariance)::physCov, uCov   
 
-      Ii = 0
       if( old) then
          if (allocated(SingleArray).and.size(SingleArray).lt.Nvpall_Old) then
             deallocate(SingleArray)
@@ -129,82 +150,83 @@ contains
             allocate(SingleArray(Nvpall_Old))
          end if
       end if
-                 
+
+      allocate(Dummy(Nvpall_Old))
+
+      ! read physical parameter covariance
+      call covData%getCovariance(physCov)
       DO Ipar=1,Nvpall_Old        
-         call readArrayData(old, Allvr(Ii+1:Ii+Ipar), Ipar, iu62)
-         Ii = Ii + Ipar     
-      end do           
-      call covData%getCovariance(physCov)         
-      call setCovMatrix(Nvpall_Old, physCov, Allvr)
+         call readArrayData(old, Dummy, Ipar, iu62)
+         do Ij = 1, Ipar
+            call covData%setCovarianceData(physCov, Ipar, Ij, Dummy(Ij))
+         end do
+      end do
+
+      ! read u parameter covariance
+      Nzero = 0
+      call covData%getUCovariance(uCov)
+      DO Ipar=1,Nvpall_Old
+         call readArrayData(old, Dummy, Ipar, iu62)
+         do Ij = 1, Ipar
+            call covData%setCovarianceData(uCov, Ipar, Ij, Dummy(Ij))
+         end do
+      end do
+
+      ! read u parameters
+      call readArrayData(old, Dummy, Nvpall_old, iu62)
+      do i = 1, Nvpall_old
+         call covData%setUParamValue(i, Dummy(I))
+      end do
+
+     deallocate(Dummy)
+
+      Droppp = dFLOAT(Idropp)*0.01d0
+      IF (Kdropp.EQ.1) THEN
+         WRITE ( 6,10000) Droppp
+         WRITE (21,10000) Droppp
+10000    FORMAT (' Drop terms in correlation matrix smaller than', F7.4)
+      END IF
 
       if (nowrt.ne.0) then
-         changed = .false.     
          IF (Kdropp.EQ.1) THEN
-            Droppp = dFLOAT(Idropp)*0.01d0
-            call covData%dropCorrelation(physCov, Droppp)       
-            changed = .true.  
+            call covData%dropCorrelation(physCov, Droppp)
          ELSE IF (Kdropp.EQ.2) THEN
-            call Xdrop2Modern (physCov, U)     
-            changed = .true.    
+            call Xdrop2Modern(physCov, Droppp)
          END IF
          IF (Fcovin.NE.0.0d0) THEN
-            call covData%scaleCovariance(physCov, Fcovin)              
-            changed = .true.     
-         end if
-              
-        if( changed) then
-           call setCovArray(physCov, Allvr)           
+            call covData%scaleCovariance(physCov, Fcovin)
          end if
       end if          
           
       if (Nowrt.NE.0.or.old) then
          IF (Kompct.EQ.2) THEN
-             CALL Covnd3 (vvv, Allvr, Nvpall_old, Idropp)   
+             CALL Covnd3 (Nvpall_old, Idropp)
          end if
       end if
-          
-      Ii = 0
-      Nzero = 0
-      DO Ipar=1,Nvpall_Old       
-         call readArrayData(old, Vrpr(Ii+1:Ii+Ipar), Ipar, iu62)      
-         Ii = Ii + Ipar     
-      end do        
-      call covData%getUCovariance(uCov)       
-      call setCovMatrix(Nvpall_Old, uCov, Vrpr)    
-          
+
       if( .not.old) then
           if (nowrt.ne.0) then
-             changed = .false.     
              IF (Kdropp.EQ.1) THEN
-                 Droppp = dFLOAT(Idropp)*0.01d0
-                 call covData%dropCorrelation(uCov, Droppp)                   
-                 changed = .true.  
+                 call covData%dropCorrelation(uCov, Droppp)
              ELSE IF (Kdropp.EQ.2) THEN
-                call Xdrop2Modern (uCov, U)     
-                changed = .true.    
+                call Xdrop2Modern(uCov, Droppp)
              END IF
              IF (Fcovin.NE.0.0d0) THEN
-                call covData%scaleCovariance(uCov, Fcovin)                       
-                changed = .true.     
+                call covData%scaleCovariance(uCov, Fcovin)
              end if
-              
-            if( changed) then
-               call setCovArray(uCov, Vrpr)                        
-            end if
          end if
          nzero = countNonZeroOffDiag(uCov)
       end if
-         
+
+      call setCovArray(physCov, Allvr)
+      call setCovArray(uCov, Vrpr)
+
       if( .not.old) then
            WRITE ( 6,10100) Nzero
            WRITE (21,10100) Nzero
 10100      FORMAT (' Number of non-zero off-diagonal cov matrix elements is', I8)
        end if
 
-       call readArrayData(old, U, Nvpall_old, iu62)
-       do i = 1, Nvpall_old
-          call covData%setUParamValue(i, U(I))
-       end do
      end subroutine readHeaderInfo
 
      !!
diff --git a/sammy/src/old/mold0.f b/sammy/src/old/mold0.f
index b868ebcdfe4d08ef6ea0f28ea2fa8af5e322f002..a049b1f74e61e4f6c36f938bac7830337da9d9fa 100644
--- a/sammy/src/old/mold0.f
+++ b/sammy/src/old/mold0.f
@@ -40,7 +40,6 @@ C
       Nt = (K*(K+1))/2
       Ivrpr  = Idimen (Nt, 1, 'Vrpr   Nt, 1')
       Iallvr = Idimen (Nt, 1, 'Allvr  Nt, 1')     
-      Ivvv   = Idimen (K , 1, 'Vvv    K , 1')
 
             N = K**2
             M = (K+1)/2
@@ -89,7 +88,7 @@ C ***    Read the COVariance file, extract information
      *      A_Iprbgf , I_Iflbgf , A_Indbgf , A_Ibgfmi , A_Ibgfma ,
      *      A_Itexbg , A_Iteabg , A_Iprdtp , I_Ifldtp ,
      *      A_Iprusd , I_Iflusd , A_Iprbag , I_Iflbag ,
-     *      A(Iu    ), A(Ivrpr ), A(Iallvr), A(Ivvv  ), A(Ipoint),
+     *      A(Ivrpr ), A(Iallvr), A(Ipoint),
      *      A(Idum1 ), A(Idum2 ), A(Idiot1), A(Idiot2),
      *      A(Ivsqua), A(Ieival), A(Ieivec), A(Izeiv ), A(Ideiv ),
      *      A(Iiord ), A(Ikord ),
@@ -115,7 +114,7 @@ C ***    Read the very old COVariance file, extract information
      *      A_Iprbgf , I_Iflbgf , A_Indbgf , A_Ibgfmi , A_Ibgfma ,
      *      A_Iprdtp , I_Ifldtp ,
      *      A_Iprusd , I_Iflusd , A_Iprbag , I_Iflbag ,
-     *      A(Iu    ), A(Ivrpr ), A(Iallvr), A(Ipoint),
+     *      A(Ivrpr ), A(Iallvr), A(Ipoint),
      *      A(Isingl), Krext, Nowrt, Nsingl, Nvpall_Old)
       END IF
 C
diff --git a/sammy/src/old/mold2.f b/sammy/src/old/mold2.f
index c643d3e7d8d1c6d508988168834f8f05091e8634..998ab51ce56c97332fd4d36ad7526b4ef010511a 100644
--- a/sammy/src/old/mold2.f
+++ b/sammy/src/old/mold2.f
@@ -108,7 +108,7 @@ C
      *   Parrpi, Iflrpi, Parudr, Ifludr, Parnbk, Iflnbk,
      *   Parbgf, Iflbgf, kndbgf, Bgfmin, Bgfmax, Texbgf, Teabgf,
      *   Pardtp, Ifldtp,
-     *   Parusd, Iflusd, Parbag, Iflbag, U, Vrpr, Allvr, Vvv   ,
+     *   Parusd, Iflusd, Parbag, Iflbag, Vrpr, Allvr,
      *   Point , Dum1  , Dum2  , Idiot1, Idiot2,
      *   Vsqua , Eival , Eivec , Zeiv  , Deiv  , Iord  , Kord,
      *   Krext , Nowrt , Nvpall_Old)
@@ -158,7 +158,7 @@ C#      CHARACTER*6 What
      *   Parbgf(*), Iflbgf(*), Kndbgf(*), Bgfmin(*), Bgfmax(*),
      *   Texbgf(Ntepnt,*),     Teabgf(Ntepnt,*),
      *   Pardtp(*), Ifldtp(*),            Parusd(*), Iflusd(*),
-     *   Parbag(*), Iflbag(*), U(*), Vrpr(*), Allvr(*), Vvv(*),
+     *   Parbag(*), Iflbag(*), Vrpr(*), Allvr(*),
      *   Point( *), Dum1(*), Dum2(*), Idiot1(*), Idiot2(*),
      *   Vsqua(Nvpall,*), Eival(Nvpall,*), Eivec(Nvpall,*),
      *   Zeiv (Nvpall,*), Deiv(*), Iord(*), Kord(*)
@@ -187,20 +187,9 @@ C
       DATA Zero /0.0d0/
 C
 C
-      IF (Kdropp.EQ.1) THEN
-         Droppp = dFLOAT(Idropp)*0.01d0
-         WRITE ( 6,10000) Droppp
-         WRITE (21,10000) Droppp
-10000    FORMAT (' Drop terms in correlation matrix smaller than', F7.4)
-      ELSE IF (Kdropp.EQ.2) THEN
-C ***      Note that Allvr and VVV are dummys here.
-           CALL Xdrop1 (Allvr, U, Vvv, Nvpall, Idropp, Iu62)
-C ***      Output has U = 0 if drop off-diagonal covariances for that par
-      END IF
 C
-
-      call readHeaderInfo(.false., Nvpall_old, Allvr, U,
-     *                     Vvv,  Vrpr, Nowrt, idropp, Kompct)
+      call readHeaderInfo(.false., Nvpall_old, Allvr,
+     *                    Vrpr, Nowrt, idropp, Kompct)
 
 C
 C     *** Nuclide abundances etc.      
diff --git a/sammy/src/old/mold5.f b/sammy/src/old/mold5.f
index 4c9795e3732a355c6605333f279877a9478a5ed2..44f7e40fd79807a234b8aebf4f4afee85ecf8800 100644
--- a/sammy/src/old/mold5.f
+++ b/sammy/src/old/mold5.f
@@ -12,7 +12,7 @@ C
      *                     Sesese, Eseses, Sigdts,
      *     Parrpi, Iflrpi, Parnbk, Iflnbk,
      *     Parbgf, Iflbgf, Kndbgf, Bgfmin, Bgfmax, Pardtp, Ifldtp,
-     *     Parusd, Iflusd, Parbag, Iflbag, U, Vrpr, Allvr,
+     *     Parusd, Iflusd, Parbag, Iflbag, Vrpr, Allvr,
      *     Point, Single, Krext, Nowrt, Nsingl, Nvpall_Old)
 C
 C       Here everything must be read in Single-precision and converted to
@@ -54,12 +54,11 @@ C
      *   Parnbk(*), Iflnbk(*),
      *   Parbgf(*), Iflbgf(*), Kndbgf(*), Bgfmin(*), Bgfmax(*),
      *   Pardtp(*), Ifldtp(*), Parusd(*), Iflusd(*),
-     *   Parbag(*), Iflbag(*), U(*), Vrpr(*), Allvr(*), Point(*)
+     *   Parbag(*), Iflbag(*), Vrpr(*), Allvr(*), Point(*)
 
       type(SammySpinGroupInfo)::spinInfo
       type(SammyResonanceInfo)::resInfo, resNew
       type(RMatResonance)::resonance, resonanceNew
-      real(kind=8),allocatable::diag(:)
       logical(C_BOOL)::reduced, inc
 C     
 C
@@ -67,15 +66,8 @@ C
 C
       IF (Nsingl.LT.Nvpall) CALL Stopss (Nvpall)
 C
-      if( Kompct.EQ.2) then
-         allocate(diag(Nvpall))
-      end if
-      call readHeaderInfo(.true., Nvpall, Allvr, U,
-     *                     Diag,  Vrpr, 0, idropp, Kompct)
-      if( Kompct.EQ.2) then
-         deallocate(diag)
-      end if
-
+      call readHeaderInfo(.true., Nvpall, Allvr,
+     *                    Vrpr, 0, idropp, Kompct)
 C
 C *** RESONANCE PARAMETERS
       IF (Nsingl.LT.Nres*ntriag) CALL Stopss (Nres*ntriag)
diff --git a/sammy/src/old/mold6.f b/sammy/src/old/mold6.f
index fdb2cc7d82d7ed735531ab189a67ccc279e17b64..429bb98754f094216f36b162bcb08d7e2adee91c 100755
--- a/sammy/src/old/mold6.f
+++ b/sammy/src/old/mold6.f
@@ -2,7 +2,7 @@ C
 C
 C --------------------------------------------------------------
 C
-      SUBROUTINE Covnd3 (Diag, Allvr, Nnpar, Idropp)
+      SUBROUTINE Covnd3 (Nnpar, Idropp)
 C
 C *** Purpose -- Write the initial parameter correlation matrix in compact 
 C ***               format (with Ndigit = 2 only, apparently)
@@ -10,11 +10,13 @@ C
       use EndfData_common_m
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
       CHARACTER*1 Ch(3,18), Bl, Minus
-      DIMENSION Diag(*), Allvr(*)
+      type(ResonanceCovariance)::physCov
+      real(kind=8),allocatable,dimension(:)::diag
       DATA Line /18/
       DATA Bl /' '/, Minus /'-'/
       DATA Zero/0.0d0/, One /1.0d0/, Hundred /100.0d0/
 C
+      call covData%getCovariance(physCov)
       IF (Idropp.NE.0) THEN
          Small = dfloat(Idropp)
       ELSE
@@ -34,13 +36,14 @@ C
 C
       Nvpall = Nnpar
       call covData%setNumParameters(Nvpall)
-      Ii = 0
+      allocate(diag(Nvpall))
       DO Ipar=1,Nvpall
-         Ii = Ii + Ipar
-         Diag(Ipar) = dSQRT(Allvr(Ii))
+         cov = physCov%getCovariance(Ipar, Ipar)
+         Diag(Ipar) = dSQRT(cov)
       END DO
       WRITE (77,10300) (Diag(Ipar),Ipar=1,Nvpall)
 10300 FORMAT ((6(1PG11.5)))
+      deallocate(diag)
 C
 C
       Ii = 1
@@ -48,13 +51,16 @@ C
       DO Ipar=2,Nvpall
          Ij = Ii
          Ii = Ii + Ipar
-         Di = Diag(Ipar)
+         Di = physCov%getCovariance(Ipar, Ipar)
+         Di = dSQRT(Di)
          Im = 0
          DO Jpar=1,Ipar-1
             Ij = Ij + 1
             IF (Ij.GT.Im) THEN
-               Dj = Diag(Jpar)
-               Vij = Allvr(Ij)/(Di*Dj) * Hundred
+               Dj = physCov%getCovariance(Jpar, Jpar)
+               Dj = dSQRT(Dj)
+               Vij = physCov%getCovariance(Ipar, Jpar)
+               Vij = Vij/(Di*Dj) * Hundred
                IF (dABS(Vij).GT.Small) THEN
 C
                   Im = Ij + Line - 1
@@ -71,16 +77,18 @@ C
                   DO Kpar=Jpar,Mk
                      Kk = Kk + 1
                      Ik = Ik + 1
-                     Dk = Diag(Kpar)
-                     Vik = dABS(Allvr(Ik))/(Di*Dk) * Hundred
+                     Dk = physCov%getCovariance(Kpar, Kpar)
+                     Dk = dSQRT(Dk)
+                     VikO = physCov%getCovariance(Ipar, Kpar)
+                     Vik = dABS(VikO)/(Di*Dk) * Hundred
                      IF (Vik.GT.Small) THEN
                         Mm = Vik
                         CALL Fixmm (Mm, Ch(2,Kk), 2)
                         IF (Ch(3,Kk).NE.Bl) THEN
                            IF (Ch(2,Kk).NE.Bl) THEN
-                              IF (Allvr(Ik).LT.Zero) Ch(1,Kk) = Minus
+                              IF (VikO.LT.Zero) Ch(1,Kk) = Minus
                            ELSE
-                              IF (Allvr(Ik).LT.Zero) Ch(2,Kk) = Minus
+                              IF (VikO.LT.Zero) Ch(2,Kk) = Minus
                            END IF
                         END IF
                      END IF
@@ -95,70 +103,4 @@ C
       END DO
       CLOSE (UNIT=77)
       RETURN
-      END
-C
-C
-C --------------------------------------------------------------
-C
-      SUBROUTINE Xdrop1 (X, U, Vvv, Nvpall_old, Idropp, Iu62)
-C
-C *** Purpose -- Prepare to drop correlations when uncertainty is large
-C
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      CHARACTER*3 Alfnm3
-      DIMENSION X(*), U(*), Vvv(*)
-      DATA Zero /0.0d0/, Hundth /0.01d0/
-C
-      Droppp = dFLOAT(Idropp)*Hundth
-      WRITE (6,10000) Idropp
-      WRITE (21,10000) Idropp
-10000 FORMAT (/, ' *** Dropping off-diagonal covariances when uncertaint
-     *y is greater',
-     */, 5X, 'than ', I2, ' percent of parameter value.')
-C
-C *** READ Allvr, ignore
-      DO Ipar=1,Nvpall_Old
-         READ (Iu62) (X(Jpar),Jpar=1,Ipar)
-      END DO
-C
-C *** Read Vrpr, store diagonals
-      DO Ipar=1,Nvpall_Old
-         READ (Iu62) (X(Jpar),Jpar=1,Ipar)
-         Vvv(Ipar) = dSQRT(X(Ipar))
-      END DO
-C
-C *** Read U, decide whether to keep off-diagonal covariances, store flag
-C ***    in U array
-      READ (Iu62) (U(Ipar),Ipar=1,Nvpall_Old)
-      DO Ipar=1,Nvpall_Old
-         IF (Vvv(Ipar)/dABS(U(Ipar)).GT.Droppp) U(Ipar) = Zero
-      END DO
-C
-C *** get back to original site in Iu62
-      REWIND Iu62
-      READ (Iu62) Alfnm3
-      READ (Iu62) Llf
-      RETURN
-      END
-C
-C
-C --------------------------------------------------------------
-C
-      SUBROUTINE Xdrop2 (X, U, Nvpall, Ipar)
-C
-C *** Purpose -- Drop correlation if uncertainty is large
-C
-      IMPLICIT DOUBLE PRECISION (a-h,o-z)
-      DIMENSION X(*), U(*)
-      DATA Zero /0.0d0/
-C
-      IF (Ipar.EQ.1) RETURN
-      Di = U(Ipar)
-      DO Jpar=1,Ipar-1
-         Dj = U(Jpar)
-         IF (Di.EQ.Zero .OR. Dj.EQ.Zero) X(Jpar) = Zero
-      END DO
-C *** From Xdrop1, U(Ipar)=Zero means this parameter has very large
-C ***    uncertainty, so we're dropping off-diagonal correlations
-      RETURN
-      END
+      END
\ No newline at end of file