anchor.diff:  		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,allocatable,dimension(:) :: rowsums,rowmeans
   		real :: cursum
  -		integer :: cur_n,curcol,ncols
  +		integer :: cur_n,curcol,ncols,nrows
   		
   		ncols=windowsize*2+1
   		nx=size(wind,1)
  -		ny=size(wind,2) !note, this could be the Y or Z dimension depending on ydim
  +		ny=size(wind,2) !note, this is Z for the high-res domain (ydim=3)
   		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), &
  +		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
  +		
  +		!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
  +		!$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
   		!$omp do schedule(static)
  -		do k=1,nz
  -			do j=1,ny
  +		do j=1,ny
  +			
  +			! so we pre-compute the sum over rows for each column in the current window
  +			if (ydim==3) then
  +				rowsums=0
  +				do i=1,windowsize
  +					rowsums=rowsums+inputwind(:,j,i)
  +				enddo
  +				nrows=windowsize
  +			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
   				! 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
  +					! if we are pinned to the top edge
  +					if ((k-windowsize)<=1) then
  +						starty=1
  +						endy  =k+windowsize
  +						nrows =nrows+1
  +						rowsums=rowsums+inputwind(:,j,endy)
  +					! if we are pinned to the bottom edge
  +					else if ((k+windowsize)>nz) then
  +						starty=k-windowsize
  +						endy  =nz
  +						nrows=nrows-1
  +						rowsums=rowsums-inputwind(:,j,starty-1)
  +					! if we are in the middle (this is the most common)
  +					else
  +						starty=k-windowsize
  +						endy  =k+windowsize
  +						rowsums=rowsums-inputwind(:,j,starty-1)+inputwind(:,j,endy)
  +					endif
  +					if ((endy-starty+1)/=nrows) then
  +						print*, "Y windowing error"
  +						print*,endy,starty,nrows
  +						print*,i,j,k
  +					endif
  +					rowmeans=rowsums/nrows
  +					cursum=sum(rowmeans(1:windowsize))
  +					ncols=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)
  +						! if we are pinned to the left edge
  +						if ((i-windowsize)<=1) then
  +							startx=1
  +							endx  =i+windowsize
  +							ncols =ncols+1
  +							cursum=cursum+rowmeans(endx)
  +						! if we are pinned to the right edge
  +						else if ((i+windowsize)>nx) then
  +							startx=i-windowsize
  +							endx  =nx
  +							ncols =ncols-1
  +							cursum=cursum-rowmeans(startx-1)
  +						! if we are in the middle (this is the most common)
   						else
  -							rowsums(curcol)=0 ! this shouldn't be necessary as we shouldn't get back to it
  +							startx=i-windowsize
  +							endx  =i+windowsize
  +							cursum=cursum-rowmeans(startx-1)+rowmeans(endx)
   						endif
   						
  -						wind(i,j,k)=cursum/cur_n
  +						if ((endx-startx+1)/=ncols) then
  +							print*, "X windowing error"
  +							print*,endx,startx,ncols
  +							print*,i,j,k
  +						endif
  +						wind(i,j,k)=cursum/ncols
  +						
   						! old SLOW way
   						! wind(i,j,k)=sum(inputwind(startx:endx,j,starty:endy)) &
   						! 			/ ((endx-startx+1)*(endy-starty+1))
  
positive.code: ! 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
  	subroutine smooth_wind(wind,windowsize,ydim)
  		implicit none
  		real, intent(inout), dimension(:,:,:):: wind	! 3 dimensional wind field to be smoothed
  
negative.code: 			enddo
  		enddo
  		!$omp end do
  		deallocate(rowmeans,rowsums)
  		!$omp end parallel
  		
  		deallocate(inputwind)
  