anchor.diff:              y = geolut%y(1,i,j)
               output = input_data(:, x, y)
           end select
  -
  -        
       end subroutine read_point
       
       subroutine transform_data(transform_type, pred_data,  &
  
positive.code: 
          endif
          
          do i = 1, n
              if (options%pure_analog) then
                  call downscale_pure_analog(predictor(i,:), atm, obs, output_coeff(:,i),     &
                                             output(i), errors(i), logistic(i),               &
                                             options, timers)
                  
              elseif (options%analog_regression) then
                  call System_Clock(timeone)
                  call find_analogs(analogs, predictor(i,:), atm, n_analogs, analog_threshold)
                  real_analogs = size(analogs)
  
                  if (real_analogs /= n_analogs) then
                      if (allocated(obs_analogs)) then
                          if (size(obs_analogs) /= real_analogs) then
                              deallocate(regression_data, obs_analogs)
                              allocate(regression_data(real_analogs, nvars))
                              allocate(obs_analogs(real_analogs))
                              if (logistic_threshold/=kFILL_VALUE) then
                                  allocate(threshold_packing(real_analogs))
                              endif
                          endif
                      else
                          allocate(regression_data(real_analogs, nvars))
                          allocate(obs_analogs(real_analogs))
                          if (logistic_threshold/=kFILL_VALUE) then
                              allocate(threshold_packing(real_analogs))
                          endif
                      endif
                  endif
                  do a=1,real_analogs
                      obs_analogs(a) = obs(analogs(a))
                      regression_data(a,:) = atm(analogs(a),:)
                  enddo
                  call System_Clock(timetwo)
                  timers(1) = timers(1) + (timetwo-timeone)
  
                  call System_Clock(timeone)
                  if (logistic_threshold==kFILL_VALUE) then
                      output(i) = compute_regression(predictor(i,:), regression_data, obs_analogs, coefficients, errors(i))
                  else
                      threshold_packing = obs_analogs > logistic_threshold
                      n_packed = count(threshold_packing)
                      
                      if (n_packed > nvars) then
                          if (allocated(threshold_atm)) deallocate(threshold_atm)
                          allocate(threshold_atm(n_packed, nvars))
                          do j=1,nvars
                              threshold_atm(:,j) = pack(regression_data(:,j), threshold_packing)
                          enddo
                          
                          if (allocated(threshold_obs)) deallocate(threshold_obs)
                          allocate(threshold_obs(n_packed))
                          threshold_obs = pack(obs_analogs, threshold_packing)
                          
                          output(i) = compute_regression(predictor(i,:), threshold_atm, threshold_obs, coefficients, errors(i))
                          
                      elseif (n_packed>0) then
                          output(i) = sum(pack(obs_analogs, threshold_packing))/n_packed
                          
                      else
                          output(i) = compute_regression(predictor(i,:), regression_data, obs_analogs, coefficients, errors(i))
                      endif
                      
                  endif
                  if (options%debug) then
                      output_coeff(1:nvars,i) = coefficients
                  endif
                  call System_Clock(timetwo)
                  timers(2) = timers(2) + (timetwo-timeone)
  
                  call System_Clock(timeone)
                  if (logistic_threshold/=kFILL_VALUE) then
                      if (options%logistic_from_analog_exceedance) then
                          
                          ! if the user specified using fewer analogs specifically for the logistic component, find them here
                          if ((options%n_log_analogs /= real_analogs).and.(options%n_log_analogs > 0)) then
                              
                              ! if the number of logistic analogs is greater than the real number of analogs, then reallocate analogs
                              if (real_analogs < options%n_log_analogs) then
                                  deallocate(analogs)
                                  allocate(analogs(options%n_log_analogs))
                              endif
                              
                              ! find the logistic analogs
                              call find_analogs(analogs, predictor(i,:), atm, options%n_log_analogs, -1.0)
                              call System_Clock(timetwo)
                              timers(9) = timers(9) + (timetwo-timeone)
                              
                              ! compute the logistic as the probability of threshold exceedance in the analog population
                              call System_Clock(timeone)
                              logistic(i) = compute_analog_exceedance(obs, analogs(1:options%n_log_analogs), logistic_threshold)
                          else
                              ! compute the logistic as the probability of threshold exceedance in the analog population
                              logistic(i) = compute_analog_exceedance(obs, analogs, logistic_threshold)
                          endif
                          
                      else
                          logistic(i) = compute_logistic_regression(predictor(i,:), regression_data, obs_analogs, coefficients, logistic_threshold)
                          if (options%debug) then
                              do j = 1,nvars
                                  output_coeff(j+nvars,i) = coefficients(j)
                              enddo
                          endif
                      endif
                  endif
                  call System_Clock(timetwo)
                  timers(3) = timers(3) + (timetwo-timeone)
                  
              elseif (options%pure_regression) then
                  ! to test matmul(predictor, output_coeff(:,1)) should provide all time values efficiently?
                  call apply_pure_regression(output(i), predictor(i,:), output_coeff(:,1), logistic(i), logistic_threshold, timers)
  
negative.code: 
      end function downscale_point
  
      subroutine downscale_analog_regression(x, atm, obs, output_coeff, output, error, logistic, options, timers)
          implicit none
          real,           intent(in),     dimension(:,:)  :: atm
          real,           intent(in),     dimension(:)    :: x, obs
          real,           intent(inout),  dimension(:)    :: output_coeff
          real,           intent(inout)                   :: output, logistic, error
          type(config),   intent(in)                      :: options
          integer*8,      intent(inout),  dimension(:)    :: timers
          
          real,    dimension(:),   allocatable :: obs_analogs
          real,    dimension(:,:), allocatable :: regression_data
          integer, dimension(:),   allocatable :: analogs
          real(8), dimension(:),   allocatable :: coefficients
          integer*8   :: timeone, timetwo
          real        :: analog_threshold, logistic_threshold
          integer     :: real_analogs, selected_analog, nvars, n_analogs
          integer     :: j, a
          
          ! variables for handling packing e.g. precip to only use positive values in computing the amount
          real,    dimension(:,:),allocatable :: threshold_atm 
          real,    dimension(:),  allocatable :: threshold_obs
          logical, dimension(:),  allocatable :: threshold_packing
          integer :: n_packed
          
          n_analogs           = options%n_analogs
          analog_threshold    = options%analog_threshold
          logistic_threshold  = options%logistic_threshold
          nvars               = size(x)
                  
          call System_Clock(timeone)
          call find_analogs(analogs, x, atm, n_analogs, analog_threshold)
          real_analogs = size(analogs)
  
          allocate(coefficients(nvars))
          allocate(regression_data(real_analogs, nvars))
          allocate(obs_analogs(real_analogs))
          if (logistic_threshold/=kFILL_VALUE) then
              allocate(threshold_packing(real_analogs))
          endif
          
          do a=1,real_analogs
              obs_analogs(a) = obs(analogs(a))
              regression_data(a,:) = atm(analogs(a),:)
          enddo
          call System_Clock(timetwo)
          timers(1) = timers(1) + (timetwo-timeone)
  
          call System_Clock(timeone)
          if (logistic_threshold==kFILL_VALUE) then
              output = compute_regression(x, regression_data, obs_analogs, coefficients, error)
          else
              
              threshold_packing = obs_analogs > logistic_threshold
              n_packed = count(threshold_packing)
              
              if (n_packed > nvars) then
                  allocate(threshold_atm(n_packed, nvars))
                  allocate(threshold_obs(n_packed))
                  do j=1,nvars
                      threshold_atm(:,j) = pack(regression_data(:,j), threshold_packing)
                  enddo
                  threshold_obs = pack(obs_analogs, threshold_packing)
                  
                  output = compute_regression(x, threshold_atm, threshold_obs, coefficients, error)
                  
              elseif (n_packed > 0) then
                  output = sum(pack(obs_analogs, threshold_packing)) / n_packed
              else
                  output = compute_regression(x, regression_data, obs_analogs, coefficients, error)
              endif
              
          endif
          if (options%debug) then
              output_coeff(1:nvars) = coefficients
          endif
          call System_Clock(timetwo)
          timers(2) = timers(2) + (timetwo-timeone)
  
          call System_Clock(timeone)
          if (logistic_threshold/=kFILL_VALUE) then
              if (options%logistic_from_analog_exceedance) then
                  
                  ! if the user specified using fewer analogs specifically for the logistic component, find them here
                  if ((options%n_log_analogs /= real_analogs).and.(options%n_log_analogs > 0)) then
                      
                      ! find the logistic analogs
                      call find_analogs(analogs, x, atm, options%n_log_analogs, -1.0)
                      call System_Clock(timetwo)
                      timers(9) = timers(9) + (timetwo-timeone)
                      
                      ! compute the logistic as the probability of threshold exceedance in the analog population
                      call System_Clock(timeone)
                      logistic = compute_analog_exceedance(obs, analogs(1:options%n_log_analogs), logistic_threshold)
                  else
                      ! compute the logistic as the probability of threshold exceedance in the analog population
                      logistic = compute_analog_exceedance(obs, analogs, logistic_threshold)
                  endif
                  
              else
                  logistic = compute_logistic_regression(x, regression_data, obs_analogs, coefficients, logistic_threshold)
                  if (options%debug) then
                      do j = 1,nvars
                          output_coeff(j+nvars) = coefficients(j)
                      enddo
                  endif
              endif
          endif
          
          call System_Clock(timetwo)
          timers(3) = timers(3) + (timetwo-timeone)
      end subroutine downscale_analog_regression
  
      subroutine downscale_pure_analog(x, atm, obs, output_coeff, output, error, logistic, options, timers)
          implicit none
  