anchor.diff:  
   
   
  +  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(allocated(err_ufunctional)) deallocate(err_ufunctional)
  +  allocate(err_ufunctional(num_ufunct,0:spatial_Level-1))
  +  err_ufunctional = 0.0d0
  +
  +  do ilevel=0,spatial_Level-1
  +    do ifunct=1,num_ufunct
  +      err_ufunctional(ifunct,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,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,ilevel1+1))                  - &
  +                      log(err_ufunctional2(ifunct,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
  +
  +
  +
     subroutine iter_calcerr_print
     !This subroutrine calculates the error (based off the most converged iteration)
     !for each functional and iterative convergence values and prints each to a file.
  
positive.code:     enddo
  
    enddo
  !  call benchmark_calcerr_print                   !calc and print errs for functs
    end subroutine UQ_benchmark
  
  
    subroutine UQ_spatialconv( icase )
  
negative.code: 
    integer :: iiter, ifunct, ibase, igap, lastiter
    real(8) :: Rval, Rvalsum
    real(8), allocatable :: err_ufunctional(:,:)
  
    lastiter = size(Q_ufunctional(1,1,:))-1
    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,iiter) = &
                    abs( Q_ufunctional(ifunct,1,lastiter)-Q_ufunctional(ifunct,1,iiter) ) / &
                         Q_ufunctional(ifunct,1,lastiter)
      enddo
  