anchor.diff:    do iiter=0,lastiter-1
       write(24,1155,advance="no") iiter
       do ifunct=1,num_ufunct
  -      write(24,1156,advance="no") err_ufunctional(ifunct,1,iiter)
  +      write(24,1156,advance="no") err_ufunctional(ifunct,iiter)
       enddo
       write(24,1154)
     enddo
     close(unit=24)    
   
     end subroutine iter_calcerr_print
  -
  -
  -
  -  subroutine spatial_calcerr_print
  -  !This subroutrine calculates the error (based off the most converged level)
  -  !for each functional and the log-log slopes (convergence rates),
  -  !and prints each of these values to a file.
  -  use MLMCvars, only: Q_ufunctional, spatial_Level, num_ufunct, def_ufunct, &
  -                      numcellsLevel0, nextLevelFactor
  -
  -  integer :: ilevel, ifunct, ibase, igap
  -  real(8) :: slope
  -  real(8), allocatable :: err_ufunctional(:,:,:)
  -
  -  if(.not.allocated(err_ufunctional)) allocate(err_ufunctional(num_ufunct,1,0:spatial_Level-1))
  -  err_ufunctional = 0.0d0
  -
  -  do ilevel=0,spatial_Level-1
  -    do ifunct=1,num_ufunct
  -      err_ufunctional(ifunct,1,ilevel) = &
  -                  abs( Q_ufunctional(ifunct,1,spatial_Level)-Q_ufunctional(ifunct,1,ilevel) ) / &
  -                       Q_ufunctional(ifunct,1,spatial_Level)
  -    enddo
  -  enddo
  -
  -  !Here I calculate the log-log slopes and print in each variation possible
  -  !for each functional chosen.  I calculate slope between the first and second
  -  !data points, first and third, etc, then first and third, first and fourth, etc.
  -  !Each of these slope calculations ought to be about the same, but this is a sanity
  -  !check that they are!
  -  call system("test -e plots/MLMCfuncts/spatialslopes.out && rm plots/MLMCfuncts/spatialslopes.out")
  -
  -  open(unit=24, file="plots/MLMCfuncts/spatialslopes.out")
  -
  -  1130 format(f15.7)
  -  1131 format("slope for ")
  -  do ifunct=1,num_ufunct
  -    write(24,1131,advance="no")
  -    if(def_ufunct(ifunct,3)==1) then
  -      write(24,1121) def_ufunct(ifunct,1),def_ufunct(ifunct,2)
  -    elseif(def_ufunct(ifunct,3)==2) then
  -      write(24,1122) def_ufunct(ifunct,1),def_ufunct(ifunct,2)
  -    elseif(def_ufunct(ifunct,3)==3) then
  -      write(24,1123) def_ufunct(ifunct,1)
  -    endif
  -    do igap=1,spatial_Level-1
  -      do ibase=0,spatial_Level-igap-1
  -        slope = spatiallogslope(err_ufunctional,ifunct,ibase,ibase+igap)
  -        write(24,1130,advance="no") slope
  -      enddo
  -      write(24,1124)
  -    enddo
  -  enddo
  -
  -  close(unit=24)      
  -
  -
  -  !Here we print the actual functional error data so that it can be examined
  -  !by hand and/or plotted.
  -  call system("test -e plots/MLMCfuncts/spatialconv.out && rm plots/MLMCfuncts/spatialconv.out")
  -
  -  open(unit=24, file="plots/MLMCfuncts/spatialconv.out")
  -  1120 format("#ilevel     num of cells   ")
  -  1121 format("L1 cell",i5," to",i5,"  ")
  -  1122 format("L2 cell",i5," to",i5,"  ")
  -  1123 format(" center cell",i5,"     ")
  -  1124 format(" ")
  -  write(24,1120,advance="no")
  -  do ifunct=1,num_ufunct
  -    if(def_ufunct(ifunct,3)==1) then
  -      write(24,1121,advance="no") def_ufunct(ifunct,1),def_ufunct(ifunct,2)
  -    elseif(def_ufunct(ifunct,3)==2) then
  -      write(24,1122,advance="no") def_ufunct(ifunct,1),def_ufunct(ifunct,2)
  -    elseif(def_ufunct(ifunct,3)==3) then
  -      write(24,1123,advance="no") def_ufunct(ifunct,1)
  -    endif
  -  enddo
  -  write(24,1124)
  -
  -  1125 format(i5,"     ",i13,"    ")
  -  1126 format(es14.7,"        ")
  -  do ilevel=0,spatial_Level-1
  -    write(24,1125,advance="no") ilevel,numcellsLevel0*nextLevelFactor**ilevel
  -    do ifunct=1,num_ufunct
  -      write(24,1126,advance="no") err_ufunctional(ifunct,1,ilevel)
  -    enddo
  -    write(24,1124)
  -  enddo
  -  close(unit=24)    
  -
  -  end subroutine spatial_calcerr_print
  -
  -
  -  function spatiallogslope(err_ufunctional2,ifunct,ilevel1,ilevel2)
  -  !This function finds the log-log slope of spatial convergence of functionals
  -  !err_ufunctional has a 2 because it is passed by reference, and the new array
  -  !is starts at 1, not 0.  The '2' denotes this difference for understandability.
  -  !Each 'ilevel#' has a plus one to compensate for this offset.
  -  use MLMCvars, only: numcellsLevel0, nextLevelFactor
  -  real(8) :: err_ufunctional2(:,:,:)
  -  real(8) :: spatiallogslope
  -  integer :: ifunct,ilevel1,ilevel2
  -
  -  spatiallogslope = ( log(err_ufunctional2(ifunct,1,ilevel1+1))                  - &
  -                      log(err_ufunctional2(ifunct,1,ilevel2+1))                ) / &
  -                    ( log(real(numcellsLevel0*nextLevelFactor**(ilevel1+1),8))   - &
  -                      log(real(numcellsLevel0*nextLevelFactor**(ilevel2+1),8)) )
  -  end function spatiallogslope
  -
  -
  -
  -  subroutine UQ_iterconv( icase )
  -  !This subroutine performs an iterative convergence study for all chosen functionals in an
  -  !effort to determine the iterative convergence parameter for each problem, solve, and QoI
  -  use genRealzvars, only: s
  -  use MLMCvars, only: Q_ufunctional, num_ufunct, ncellwidth, numcellsLevel0
  -  use FEDiffSn, only: setflvarspassedtrue, a, fliterstudy, max_iter, flnewiter, fliterstudy
  -  integer :: icase, iiter
  -  real(8), allocatable :: trarray3(:,:,:)
  -
  -  if(allocated(Q_ufunctional)) deallocate(Q_ufunctional)
  -  allocate(Q_ufunctional(num_ufunct,1,0:0))
  -  Q_ufunctional = 0.0d0
  -  if(allocated(ncellwidth)) deallocate(ncellwidth)
  -  allocate(ncellwidth(0:0))
  -  ncellwidth = 0.0d0
  -  ncellwidth = s/real(numcellsLevel0,8)
  -  call setflvarspassedtrue                             !tell FE mod to accept input from here
  -  a             = s
  -
  -  max_iter = 0
  -  flnewiter = .true.
  -  fliterstudy = .true.
  -  do while(flnewiter)                       !solve until FEDiffSn solver says it's converged
  -    call sampleconvInput( 0 )               !samples average values, ilevel=0 so to set # of cells
  -    call solveSamples( max_iter,1,icase )   !solves QoIs, ilevel=max_iter, isamp=1
  -
  -    if(flnewiter) then
  -      !add new max_iter level to Q_ufunctional and initialize it
  -      call move_alloc(Q_ufunctional,trarray3)
  -      allocate(Q_ufunctional(size(trarray3(:,1,1)),size(trarray3(1,:,1)),0:size(trarray3(1,1,:))))
  -      Q_ufunctional = 0.0d0
  -      Q_ufunctional(:,:,0:size(trarray3(1,1,:))-1) = trarray3
  -      deallocate(trarray3)
  -    endif
  -
  -    max_iter = max_iter + 1
  -  enddo
  -  call iter_calcerr_print                   !calc and print errs for functs
  -  end subroutine UQ_iterconv
  -
  -
  -
   
   
   
  
positive.code: 
    integer :: iiter, ifunct, ibase, igap, lastiter
    real(8) :: Rval, Rvalsum
    real(8), allocatable :: err_ufunctional(:,:,:)
  
    lastiter = size(Q_ufunctional(1,1,:))-1
    if(.not.allocated(err_ufunctional)) allocate(err_ufunctional(num_ufunct,1,0:lastiter))
    err_ufunctional = 0.0d0
  
    do iiter=0,lastiter
      do ifunct=1,num_ufunct
        err_ufunctional(ifunct,1,iiter) = &
                    abs( Q_ufunctional(ifunct,1,lastiter)-Q_ufunctional(ifunct,1,iiter) ) / &
                         Q_ufunctional(ifunct,1,lastiter)
      enddo
  
negative.code:     endif
      Rvalsum = 0.0d0
      do iiter=0,lastiter-3
        Rval = log(err_ufunctional(ifunct,iiter+1)/err_ufunctional(ifunct,iiter  )) / &
               log(err_ufunctional(ifunct,iiter+2)/err_ufunctional(ifunct,iiter+1))
        Rvalsum = Rvalsum + Rval
        write(24,1140,advance="no") Rval
      enddo
  