anchor.diff:                  u = calc_u( dir_values(i), spd_values(k) )
                   v = calc_v( dir_values(i), spd_values(k) )
   
  -                ! if (this_image()==1) print*, i,j,k, u,v, dir_values(i), spd_values(k)
                   debug_print = .False.
                   if ((spd_values(k)==15).and.(abs(u-15) < 1).and.(j==15)) debug_print=.True.
                   ! calculate the linear wind field for the current u and v values
  
positive.code:             return
          endif
  
          nx = size(z_top,1)
          ny = size(z_top,2)
  
          ! print*, this_image(), "minval(z_bottom), maxval(z_top)", minval(z_bottom), maxval(z_top)
  
          start_z = minval(z_bottom(2:nx-1,2:ny-1))
          end_z = maxval(z_top(2:nx-1,2:ny-1))
          ! print*, this_image(), start_z, end_z
          allocate(internal_z_top(size(lt_data%u_perturb,1),size(lt_data%u_perturb,2)), source=maxval(z_top))
          ! internal_z_top = maxval(z_top(2:nx-1,2:ny-1))
          ! internal_z_top(buffer+1:buffer+size(z_top,1)-2, buffer+1:buffer+size(z_top,2)-2) = z_top(2:nx-1,2:ny-1)
          internal_z_top(buffer:buffer+size(z_top,1)-1, buffer:buffer+size(z_top,2)-1) = z_top
  
          allocate(internal_z_bottom(size(lt_data%u_perturb,1),size(lt_data%u_perturb,2)), source=minval(z_bottom))
          ! internal_z_bottom(buffer+1:buffer+size(z_bottom,1)-2, buffer+1:buffer+size(z_bottom,2)-2) = z_bottom(2:nx-1,2:ny-1)
          internal_z_bottom(buffer:buffer+size(z_bottom,1)-1, buffer:buffer+size(z_bottom,2)-1) = z_bottom
  
          allocate(layer_count(size(lt_data%u_perturb,1), size(lt_data%u_perturb,2)))
  
negative.code:         allocate(layer_fraction, source=layer_count)
  
  
          step_size = min(minimum_step, minval(z_top - z_bottom))
          if (step_size < 5) then
              if (this_image()==1) write(*,*) "WARNING: Very small vertical step size in linear winds, check layer thicknesses in output"
          endif
  
          current_z = start_z + step_size/2 ! we want the value in the middle of each theoretical layer
  
          ! sum up u/v perturbations over all layers evaluated
          lt_data%u_accumulator = 0
          lt_data%v_accumulator = 0
  
          do while (current_z < end_z)
  
              call linear_perturbation_at_height(U,V,Nsq,current_z, fourier_terrain, lt_data)
  
              layer_fraction = max(0.0,                                                                                       &
                                  min(step_size/2, current_z - internal_z_bottom) + min(0.0, internal_z_top - current_z)      &
                                + min(step_size/2, internal_z_top - current_z)    + min(0.0, current_z - internal_z_bottom)   ) / step_size
  
              layer_count = layer_count + layer_fraction
              lt_data%u_accumulator = lt_data%u_accumulator + lt_data%u_perturb * layer_fraction
              lt_data%v_accumulator = lt_data%v_accumulator + lt_data%v_perturb * layer_fraction
  
              current_z = current_z + step_size
          enddo
  
          lt_data%u_perturb = lt_data%u_accumulator / layer_count
          lt_data%v_perturb = lt_data%v_accumulator / layer_count
  