anchor.diff:                  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
  +                
  +                where(nsq_calibration<1) nsq_calibration = 1+1/((1-1/nsq_calibration)/100)
  +                where(nsq_calibration>1) nsq_calibration = 1+(nsq_calibration-1)/100
               endif
           endif
           
  
positive.code:     ! 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
      ! 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, 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(1e-4)
  !     real, parameter :: nsqmin=log(3e-5)
      
  
negative.code:                                             domain%z(i,bottom,k),  domain%z(i,top,k),   &
                                              domain%qv(i,bottom,k), domain%qv(i,top,k),  &
                                              domain%cloud(i,j,k)+domain%ice(i,j,k)+domain%qrain(i,j,k)+domain%qsnow(i,j,k))
                      domain%nsquared(i,j,k) = max(min(domain%nsquared(i,j,k) * nsq_calibration(i,k), max_stability), min_stability)
                  end do
              end do
          end do
  