anchor.diff:    C_w  = max( abs(Gave(Level  ))/(ncellwidth(Level  )**qq*(nextLevelFactor**qq-1)), &
                 abs(Gave(Level-1))/(ncellwidth(Level-1)**qq*(nextLevelFactor**qq-1))    )
     err1 = C_w * ncellwidth(Level)**qq
  -print *,"err1:",err1
  +print *,"err1, disc err:",err1
     Vsum = 0.0d0
     do ilevel=0,Level
       Vsum = Vsum + abs(Gvar(ilevel)/M_optsamps(1,ilevel))
  
positive.code:   subroutine solveSamples( ilevel,isamp )
    !This subroutine drives the solver for a new response function flux for a new sample
    !collects solution, and deallocated module variables from FEDiffSn.
    use MLMCvars, only: Q_ufunctional, G_ufunctional, nextLevelFactor, numcellsLevel0
    use FEDiffSn, only: FEmain,FEDiffSn_externaldeallocate, &
                        solve,phidiff,phiSnl,phiSnr,phiDSAl,phiDSAr
  
    integer :: ilevel, isamp, i, icell, first, last
    real(8), allocatable :: flux(:)
  
    call FEmain
  
  
negative.code:       flMLMC = .false.
      else
        Level = Level + 1
  !if(Level==4) flMLMC = .false.
      endif
    enddo !loops over realizations
  
  