anchor.diff:          !$omp parallel firstprivate(windowsize,nx,ny,nz,ydim), &
           !$omp private(i,j,k,startx,endx,starty,endy, rowsums,rowmeans,nrows,ncols,cursum), &
           !$omp shared(wind,inputwind)
  +
           allocate(rowsums(nx)) !this is only used when ydim=3, so nz is really ny
           allocate(rowmeans(nx)) !this is only used when ydim=3, so nz is really ny
  -        nrows=windowsize*2+1
  -        ncols=windowsize*2+1
  +        nrows = windowsize * 2 + 1
  +        ncols = windowsize * 2 + 1
  +
           !$omp do schedule(static)
           do j=1,ny
   
               ! so we pre-compute the sum over rows for each column in the current window
               if (ydim==3) then
  -                rowsums=inputwind(1:nx,j,1)*(windowsize+1)
  -                do i=1,windowsize
  -                    rowsums=rowsums+inputwind(1:nx,j,i)
  +                ! start by adding the outer row n+1 times
  +                rowsums = inputwind(1:nx,j,1) * (windowsize+1)
  +                do i=1, windowsize
  +                    rowsums = rowsums + inputwind(1:nx, j, i)
                   enddo
               endif
  +
               ! don't parallelize over this loop because it is much more efficient to be able to assume
               ! that you ran the previous row in serial for the slow (ydim=3) case
               do k=1,nz ! note this is y for ydim=3
  
positive.code:         allocate(inputtemp(nx,ny))
  
          inputtemp = input !make a copy so we always use the unsmoothed data when computing the smoothed data
          if ((windowsize*2+1)>nx) then
              write(*,*) "WARNING can not operate if windowsize*2+1 is larger than nx"
              write(*,*) "NX         = ", nx
              write(*,*) "windowsize = ", windowsize
              stop
          endif
  
          !parallelize over a slower dimension (not the slowest because it is MUCH easier this way)
          ! as long as the inner loops (the array operations) are over the fastest dimension we are mostly OK
  
negative.code: 
          allocate(inputwind(nx,ny,nz)) ! Can't be module level because nx,ny,nz could change between calls,
  
          inputwind = wind !make a copy so we always use the unsmoothed data when computing the smoothed data
          if (((windowsize+1)>nx).and.(ydim==3)) then
              write(*,*) "WARNING can not operate if windowsize+1 is larger than nx"
              write(*,*) "NX         = ", nx
              write(*,*) "windowsize = ", windowsize
              print*, "Image = ",this_image()
              stop
          endif
  
  