anchor.diff:                          if (this_image()==1) write(*,"(f5.1,A)") loops_completed/real(nz*(stop_pos-start_pos+1))*100," %"
                           !$omp end critical (print_lock)
                       endif
  -                    ! if (reverse) then
  -                    !     layer_height        = domain%z(1,1,z) - domain%terrain%data_2d(1,1)
  -                    !     layer_height_bottom = layer_height - (options%parameters%dz_levels(z) / 2)
  -                    !     layer_height_top    = layer_height + (options%parameters%dz_levels(z) / 2)
  -                    ! else
  +
  +                    if (options%parameters%space_varying_dz) then
  +
  +                        call linear_perturbation(u, v, exp(nsq_values(j)),                                                                          &
  +                                                 domain%z_interface%data_3d(:,z,:) - domain%terrain%data_2d,                                        &
  +                                                 domain%z_interface%data_3d(:,z,:) - domain%terrain%data_2d + domain%dz_interface%data_3d(:,z,:),   &
  +                                                 minimum_layer_size, domain%terrain_frequency, lt_data_m)
  +                    else
                           layer_height = domain%z%data_3d(ims,z,jms) - domain%terrain%data_2d(ims,jms)
                           layer_height_bottom = layer_height - (options%parameters%dz_levels(z) / 2)
                           layer_height_top    = layer_height + (options%parameters%dz_levels(z) / 2)
  -                    ! endif
  -
  -                    call linear_perturbation(u, v, exp(nsq_values(j)),                                  &
  -                                             layer_height_bottom, layer_height_top, minimum_layer_size, &
  -                                             domain%terrain_frequency, lt_data_m)
  -
  +
  +                        call linear_perturbation(u, v, exp(nsq_values(j)),                                  &
  +                                                 layer_height_bottom, layer_height_top, minimum_layer_size, &
  +                                                 domain%terrain_frequency, lt_data_m)
  +                    endif
                       ! need to handle stagger (nxu /= nx) and the buffer around edges of the domain
                       if (nxu /= nx) then
                           temporary_u(:,:) = real( real(                                              &
  
positive.code: 
          lt_data%u_perturb = lt_data%u_accumulator / n_steps
          lt_data%v_perturb = lt_data%v_accumulator / n_steps
      end subroutine linear_perturbation
  
      !>----------------------------------------------------------
      !! Add a smoothed buffer around the edge of the terrain to prevent crazy wrap around effects
      !! in the FFT due to discontinuities between the left and right (or top and bottom) edges of the domain
  
negative.code:         ! returns u_perturb and v_perturb
      end subroutine linear_perturbation_at_height
  
      subroutine linear_perturbation_constz(U, V, Nsq, z_bottom, z_top, minimum_step, fourier_terrain, lt_data)
          use, intrinsic :: iso_c_binding
          implicit none
          real,                     intent(in)    :: U, V             ! U and V components of background wind
  