anchor.diff:                                              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) = domain%nsquared(i,j,k) * nsq_calibration(i,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
  
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:     type(C_PTR) :: uh_aligned_data, u_h_aligned_data, vh_aligned_data, v_h_aligned_data
      type(C_PTR), allocatable :: uplans(:), vplans(:)
      
      real, allocatable, dimension(:) :: dir_values, nsq_values, spd_values
      ! Look Up Tables for linear perturbation are nspd x n_dir_values x n_nsq_values x nx x nz x ny
      real, allocatable, target, dimension(:,:,:,:,:,:) :: hi_u_LUT, hi_v_LUT, rev_u_LUT, rev_v_LUT
  