anchor.diff:  
           nxu = size(domain%u,1)
           nyv = size(domain%v,3)
  -        
  +
           fftnx = size(domain%fzs,1)
           fftny = size(domain%fzs,2)
  -        ! note: 
  +        ! note:
           ! buffer = (fftnx - nx)/2
   
           ! default assumes no errors in reading the LUT
           error = 0
  -        
  +
           ! store to make it easy to check dim sizes in read_LUT
           LUT_dims(:,1) = [nxu,nz,ny]
           LUT_dims(:,2) = [nx,nz,nyv]
  
positive.code:                 write(*,*) "  from file: " // trim(options%linear_mask_file)
                  write(*,*) "  with var: "  // trim(options%linear_mask_var)
                  call io_read2d(options%linear_mask_file, options%linear_mask_var, domain%linear_mask)
                  
                  linear_mask = domain%linear_mask * linear_contribution
              endif
  
              ! Stupidly simple adjustment to the linear wind field to account for using density. 
              ! If we are using density in the advection calculations, modify the linear perturbation
              ! to get the vertical velocities closer to what they would be without density (boussinesq)
              ! need to check if this makes the most sense when close to the surface
  
negative.code:             ! calculate the mean stability over the entire profile
              BV_freq=sum(vertical_N)/nz
              ! impose limits so the linear solution doesn't go crazy in really stable or unstable air
              BV_freq=min(max(BV_freq, min_stability), max_stability)
          else
              ! or just use the supplied BV frequency
              BV_freq=N_squared
  