anchor.diff:  ! Smooth an array (written for wind but will work for anything)
   ! only smooths over the first (x) and second or third (y) dimension
   ! ydim can be specified to allow working with (x,y,z) data or (x,z,y) data
  +! WARNING: this is a moderately complex setup to be efficient for the ydim=3 (SLOW) case
  +! be careful when editting
   	subroutine smooth_wind(wind,windowsize,ydim)
   		implicit none
   		real, intent(inout), dimension(:,:,:):: wind	! 3 dimensional wind field to be smoothed
  
positive.code: 		real,allocatable,dimension(:,:,:)::inputwind    ! temporary array to store the input data in
  		integer::i,j,k,nx,ny,nz,startx,endx,starty,endy ! various array indices/bounds
  		! intermediate sums to speed up the computation
  		real,dimension(windowsize*2+1) :: rowsums
  		real :: cursum
  		integer :: cur_n,curcol,ncols
  		
  		ncols=windowsize*2+1
  		nx=size(wind,1)
  		ny=size(wind,2) !note, this could be the Y or Z dimension depending on ydim
  		nz=size(wind,3) !note, this could be the Y or Z dimension depending on ydim
  		! this is ~20MB that has to be allocated deallocated every call...
  		allocate(inputwind(nx,ny,nz)) ! Can't be module level because nx,ny,nz could change between calls, 
  									  ! could be part of a "smoothable" object to avoid allocate-deallocating constantly
  		
  		inputwind=wind !make a copy so we always use the unsmoothed data when computing the smoothed data
  		
  		!parallelize over the slowest dimension
  		!$omp parallel firstprivate(windowsize,nx,ny,nz,ydim,ncols), &
  		!$omp private(i,j,k,startx,endx,starty,endy, rowsums,cursum,cur_n,curcol), &
  		!$omp shared(wind,inputwind)
  		!$omp do schedule(static)
  		do k=1,nz
  			do j=1,ny
  				! ydim=3 for the main model grid which is large and takes a long time
  				! so we pre-compute the sum over rows for each column in the current window
  				if (ydim==3) then
  					rowsums=0
  					starty=max(1, k-windowsize)
  					endy  =min(nz,k+windowsize)
  ! 					do i=1,windowsize+1
  ! 						rowsums(i)=sum(inputwind(i,j,starty:endy))/(endy-starty+1)
  ! 					enddo
  					do i=starty,endy
  						rowsums(1:windowsize+1)=rowsums(1:windowsize+1)+inputwind(1:windowsize+1,j,i)
  					enddo
  					rowsums=rowsums/(endy-starty+1)
  					cursum=sum(rowsums(1:windowsize))
  					cur_n=windowsize
  					curcol=windowsize
  				endif
  				
  				do i=1,nx
  					if (ydim==3) then
  						! find the current window bounds
  						startx=max(1, i-windowsize)
  						endx  =min(nx,i+windowsize)
  						curcol=mod(curcol+1,ncols)
  						if (curcol==0) curcol=ncols
  						
  						cursum=cursum-rowsums(curcol)
  						
  						! general window size
  						cur_n=endx-startx+1
  						! if the window is not pinned to the right edge, add the new column
  						if ((i+windowsize)<=nx) then
  							! this is an internal window set the "old" column to its new value
  							rowsums(curcol)=sum(inputwind(endx,j,starty:endy))/(endy-starty+1)
  							! then compute the mean within that window (sum/n)
  							cursum=cursum+rowsums(curcol)
  						else
  							rowsums(curcol)=0 ! this shouldn't be necessary as we shouldn't get back to it
  						endif
  						
  						wind(i,j,k)=cursum/cur_n
  						! old SLOW way
  						! wind(i,j,k)=sum(inputwind(startx:endx,j,starty:endy)) &
  						! 			/ ((endx-startx+1)*(endy-starty+1))
  
negative.code: 			enddo
  		enddo
  		!$omp end do
  		deallocate(rowmeans,rowsums)
  		!$omp end parallel
  		
  		deallocate(inputwind)
  