anchor.diff:          class(linearizable_type),intent(inout)::domain
   		integer :: nx,ny
   		
  -		ny=size(domain%terrain,1)
  -		nx=size(domain%terrain,2)
  +		nx=size(domain%terrain,1)
  +		ny=size(domain%terrain,2)
   ! 		compute u at x mid-point
  -		domain%u(:,:,1:nx-1)=(domain%u(:,:,1:nx-1)+domain%u(:,:,2:nx))/2
  +		domain%u(1:nx-1,:,:)=(domain%u(1:nx-1,:,:)+domain%u(2:nx,:,:))/2
   ! 		compute v at y mid-point
  -		domain%v(1:ny-1,:,:)=(domain%v(1:ny-1,:,:)+domain%v(2:ny,:,:))/2
  +		domain%v(:,:,1:ny-1)=(domain%v(:,:,1:ny-1)+domain%v(:,:,2:ny))/2
   		
   	end subroutine stagger_winds
   	
  
positive.code: 		real, dimension(:,:),allocatable :: real_terrain
  		integer::nx,ny,i,pos
  		real::weight
  ! 		WARNING x,y dimensions are swapped this doesn't affect output but may confuse readability
  		ny=size(terrain,1)+buffer*2
  		nx=size(terrain,2)+buffer*2
  		allocate(buffer_topo(ny,nx))
  		buffer_topo=minval(terrain)
  		
  		buffer_topo(1+buffer:ny-buffer,1+buffer:nx-buffer)=terrain
  		do i=1,buffer
  			weight=i/(real(buffer)*2)
  			pos=buffer-i
  			buffer_topo(pos+1,1+buffer:nx-buffer)  =terrain(1,1:nx-buffer*2)*(1-weight)+terrain(ny-buffer*2,1:nx-buffer*2)*weight
  			buffer_topo(ny-pos-1,1+buffer:nx-buffer) =terrain(1,1:nx-buffer*2)*(weight)+terrain(ny-buffer*2,1:nx-buffer*2)*(1-weight)
  			buffer_topo(1+buffer:ny-buffer,pos+1)  =terrain(1:ny-buffer*2,1)*(1-weight)+terrain(1:ny-buffer*2,nx-buffer*2)*weight
  			buffer_topo(1+buffer:ny-buffer,nx-pos-1) =terrain(1:ny-buffer*2,1)*(weight)+terrain(1:ny-buffer*2,nx-buffer*2)*(1-weight)
  		enddo
  		print*, "There may be a problem in linear_winds-add_buffer_topo"
  ! 		clearly something isn't quiet right here...
  		buffer_topo(ny,:)=buffer_topo(ny-1,:)
  		buffer_topo(:,nx)=buffer_topo(:,nx-1)
  		allocate(real_terrain(ny,nx))
  		real_terrain=buffer_topo
  		call io_write2d("complex_terrain.nc","HGT",real_terrain)
  		deallocate(real_terrain)
  
negative.code: 	        allocate(v_hat(nx,ny))
  	        allocate(m(nx,ny))
  	        allocate(msq(nx,ny))
  	        allocate(mimag(nx,ny))
  	        allocate(ineta(nx,ny))
          
  	!         # % Compute 2D k and l wavenumber fields (could be stored in options or domain or something)
  