anchor.diff:  		if (linear_contribution/=1) then
   			write(*,*) "Using a fraction of the linear perturbation:",linear_contribution
   		endif
  +
  +		! set up linear_mask variable
  +		nx=size(domain%terrain,1)
  +		ny=size(domain%terrain,2)
  +		if ( (.not.reverse) .and. (.not.allocated(linear_mask)) ) then
  +			allocate(linear_mask(nx,ny))
  +			linear_mask=linear_contribution
  +	
  +			if (options%linear_mask) then
  +				print*, "Reading Linear Mask"
  +				print*, "  from file: "//trim(options%linear_mask_file)
  +				print*, "  varname: "//trim(options%linear_mask_var)
  +				call io_read2d(options%linear_mask_file,options%linear_mask_var,domain%linear_mask)
  +				linear_mask=(1-(1-domain%linear_mask)*0.8) * linear_contribution
  +			endif
  +		endif
   		
   		if (options%spatial_linear_fields) then
   			if ((.not.allocated(hi_u_LUT) .and. (.not.reverse)) .or. ((.not.allocated(rev_u_LUT)) .and. reverse)) then
   				
  -				nx=size(domain%terrain,1)+1
  -				ny=size(domain%terrain,2)+1
  -				allocate(linear_mask(nx,ny))
  -				linear_mask=1.0
  -				if (.not.reverse) then
  -					if (options%linear_mask) then
  -						print*, "Reading Linear Mask"
  -						print*, "  from file: "//trim(options%linear_mask_file)
  -						print*, "  varname: "//trim(options%linear_mask_var)
  -						call io_read2d(options%linear_mask_file,options%linear_mask_var,domain%linear_mask)
  -						linear_mask(:nx-1,:ny-1)=1-(1-domain%linear_mask)*0.8
  -						linear_mask(nx,:ny-1)=linear_mask(nx-1,:ny-1)
  -						linear_mask(:nx-1,ny)=linear_mask(:nx-1,ny-1)
  -						linear_mask(nx,ny)=linear_mask(nx-1,ny-1)
  -					endif
  -				endif
  -			
   				write(*,*) "Generating a spatially variable linear perturbation look up table"
   				call initialize_spatial_winds(domain,options,reverse,useDensity)
  -				deallocate(linear_mask)
   			else
   				write(*,*) "Skipping spatial wind field for presumed domain repeat"
   			endif
  
positive.code: 	            m = sqrt(msq)         ! # % vertical wave number, hydrostatic
  				where(sig<0) m=m*(-1)   ! equivilant to m=m*sign(sig)
  	    		where(real(msq)<0) m=mimag
  			
  	            ineta=j*domain%fzs*exp(j*m* &
  						sum(domain%z(:,z,:)-domain%z(:,1,:)) &
  						/(realnx*realny))
  				!  what=sig*ineta
  
  				!# with coriolis : [+/-]j*[l/k]*f
  
negative.code:         complex,parameter :: j= (0,1)
          real::gain,offset !used in setting up k and l arrays
          integer::nx,ny,nz,i,midpoint,z,realnx,realny,realnx_u,realny_v,x,y,bottom,top
  		logical :: staggered, zaxis_is_third
          real::U,V
  		real,dimension(:),allocatable::U_layers,V_layers,preU_layers,preV_layers
  		real,parameter::pi=3.1415927
  