anchor.diff:  
     integer :: iiter, ifunct, ibase, igap, lastiter
     real(8) :: Rval, Rvalsum
  -  real(8), allocatable :: err_ufunctional(:,:,:)
  +  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))
  +  if(allocated(err_ufunctional)) deallocate(err_ufunctional)
  +  allocate(err_ufunctional(num_ufunct,0:lastiter))
     err_ufunctional = 0.0d0
   
     do iiter=0,lastiter
       do ifunct=1,num_ufunct
  -      err_ufunctional(ifunct,1,iiter) = &
  +      err_ufunctional(ifunct,iiter) = &
                     abs( Q_ufunctional(ifunct,1,lastiter)-Q_ufunctional(ifunct,1,iiter) ) / &
                          Q_ufunctional(ifunct,1,lastiter)
       enddo
  
positive.code:   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)
      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
  
  
  
  
  
  
  
negative.code:     enddo
  
    enddo
    call benchmark_calcerr_print                    !calc and print errs for functs
    end subroutine UQ_benchmark
  
  
  
    subroutine benchmark_calcerr_print
    !This subroutrine calculates the error of the ensemble averaged functional
    !values (based off the most converged level) for each functional along with
    !confidence bars produced by propagating SEM at chosen CI through error calculation,
    !and prints each of these values to a file, along with another file which is only
    !of the most converged level which can be used to compare functional profiles to MLMC.
    use MLMCvars, only: spatial_Level, num_ufunct, def_ufunct, &
                        numcellsLevel0, nextLevelFactor, Gave, Gvar
  
    integer :: ilevel, ifunct, ibase, igap
    real(8) :: slope
    real(8), allocatable :: err_ufunctional(:,:), err_ufunctSEM(:,:)
  
    if(allocated(err_ufunctional)) deallocate(err_ufunctional)
    allocate(err_ufunctional(num_ufunct,0:spatial_Level-1))
    err_ufunctional = 0.0d0
    if(allocated(err_ufunctSEM)) deallocate(err_ufunctSEM)
    allocate(err_ufunctSEM(num_ufunct,0:spatial_Level-1))
    err_ufunctSEM = 0.0d0
  
    do ilevel=0,spatial_Level-1
      do ifunct=1,num_ufunct
        err_ufunctional(ifunct,ilevel) = &
            abs( Gave(ifunct,spatial_Level)-Gave(ifunct,ilevel) ) / &
                 Gave(ifunct,spatial_Level)
        err_ufunctSEM(ifunct,ilevel) = &
             sqrt(                  (1.0d0/Gave(ifunct,spatial_Level)**2 * Gvar(ifunct,ilevel)**2 + &
         (2.0d0*Gave(ifunct,ilevel)/(Gave(ifunct,spatial_Level)**2)))**2 * Gvar(ifunct,spatial_Level)**2 )
      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 benchmark_calcerr_print
  
  
  
  
  
  
    subroutine UQ_spatialconv( icase )
  