anchor.diff:          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
  +        ! 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
  
positive.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*2+1)>nx).and.(ydim==3)) then
              write(*,*) "WARNING can not operate if windowsize*2+1 is larger than nx"
              write(*,*) "NX         = ", nx
              write(*,*) "windowsize = ", windowsize
              stop
          endif
  
  
negative.code:                 if (ydim==3) then
                      ! if we are pinned to the top edge
                      if ((k-windowsize)<=1) then
                          starty  = 1
                          endy    = min(nz, k+windowsize)
                          rowsums = rowsums - inputwind(1:nx, j, starty) + inputwind(1:nx,j,endy)
  
                      ! if we are pinned to the bottom edge
                      else if ((k+windowsize)>nz) then
                          starty  = max(1, k-windowsize)
                          endy    = nz
                          rowsums = rowsums - inputwind(1:nx, j, starty-1) + inputwind(1:nx, j, endy)
  
                      ! if we are in the middle (this is the most common)
                      else
                          starty  = max(1, k - windowsize)
                          endy    = min(nz, k + windowsize)
                          rowsums = rowsums - inputwind(1:nx, j, starty-1) + inputwind(1:nx, j, endy)
                      endif
                      rowmeans = rowsums / nrows
                      cursum = sum(rowmeans(1:windowsize)) + rowmeans(1) * (windowsize+1)
                  endif
  
                  do i=1,nx
                      if (ydim==3) then
                          ! if we are pinned to the left edge
                          if ((i-windowsize)<=1) then
                              startx = 1
                              endx   = min(nx, i + windowsize)
                              cursum = cursum - rowmeans(startx) + rowmeans(endx)
                          ! if we are pinned to the right edge
                          else if ((i+windowsize)>nx) then
                              startx = max(2, i - windowsize)
                              endx   = nx
                              cursum = cursum - rowmeans(startx-1) + rowmeans(endx)
                          ! if we are in the middle (this is the most common)
                          else
                              startx = max(2, i - windowsize)
                              endx   = min(nx, i + windowsize)
                              cursum = cursum - rowmeans(startx-1) + rowmeans(endx)
                          endif
  
                          wind(i,j,k) = cursum / ncols
  
                      else ! ydim==2
                          ! ydim=2 for the input data which is a small grid, thus cheap so we still use the slow method
                          !first find the current window bounds
                          startx=max(1, i - windowsize)
                          endx  =min(nx,i + windowsize)
                          starty=max(1, j - windowsize)
                          endy  =min(ny,j + windowsize)
                          ! then compute the mean within that window (sum/n)
                          ! note, artifacts near the borders (mentioned in ydim==3) don't affect this
                          ! because the borders *should* be well away from the domain
                          ! if the forcing data are not much larger than the model domain this *could* create issues
                          wind(i,j,k) = sum(inputwind(startx:endx,starty:endy,k)) &
                                         / ((endx-startx+1)*(endy-starty+1))
                      endif
                  enddo
              enddo
  