anchor.diff:      ! store the linear perturbation so we can update it slightly each time step
       ! this permits the linear field to build up over time. 
       real, allocatable, dimension(:,:,:) :: u_perturbation, v_perturbation
  +    
  +    integer :: buffer, original_buffer=50 ! number of grid cells to buffer around the domain MUST be >=1
  +    integer,parameter :: stability_window_size=2
  +    
  +    real, parameter :: max_stability= 6e-4 ! limits on the calculated Brunt Vaisala Frequency
  +    real, parameter :: min_stability= 1e-7 ! these may need to be a little narrower. 
  +    real :: linear_contribution = 1.0 ! multiplier on uhat,vhat before adding to u,v
       ! controls the rate at which the linearfield updates (should be calculated as f(in_dt))
       ! new/current perturbation is multiplited by linear_update and added to (1-linear_update) * the previous combined perturbation
  -    real :: linear_update_fraction = 0.5
  +    real :: linear_update_fraction = 1.0
       
       real, parameter :: dirmax=2*pi
       real, parameter :: dirmin=0
       real, parameter :: spdmax=30
       real, parameter :: spdmin=0
  -    real, parameter :: nsqmax=log(6e-4)
  -    real, parameter :: nsqmin=log(1e-7)
  +    real, parameter :: nsqmax=log(max_stability)
  +    real, parameter :: nsqmin=log(min_stability)
   !     real, parameter :: nsqmax=log(1e-4)
   !     real, parameter :: nsqmin=log(3e-5)
       
  
positive.code:                 print*, "  varname: "//trim(options%nsq_calibration_var)
                  call io_read2d(options%nsq_calibration_file,options%nsq_calibration_var,domain%nsq_calibration)
                  nsq_calibration=domain%nsq_calibration
              endif
          endif
          
  
negative.code:                     top   = min(j+winsz,nz)
                      south = max(k-winsz,1)
                      north = min(k+winsz,ny)
                      n = (((west-east)+1) * ((top-bottom)+1) * ((north-south)+1))
                      
                      curnsq = sum(log(domain%nsquared(east:west, bottom:top,south:north)))
                      curnsq = curnsq / n
  
                      npos=1
  