anchor.diff:                          !print *, "turbI, turb: ", turbI, turb
                           ! calculate TI value at each turbine
   !                         print *, "turb, turbI: ", turb, turbI
  +                        
  +                        ! save ti_area_ratio and ti_dst to new memory locations to avoid 
  +                        ! aliasing during differentiation
  +                        TI_area_ratio_temp = TI_area_ratio
  +                        TI_dst_tmp = TIturbs(turbI)
  +                        
                           call added_ti_func(TI, Ct(turb), x, ky_local, rotorDiameter(turb), & 
                                              & rotorDiameter(turbI), deltay, turbineZ(turb), &
                                              & turbineZ(turbI), TIturbs(turb), &
  -                                           & TI_calculation_method, &
  -                                           & ti_area_ratio, TIturbs(turbI))
  +                                           & TI_calculation_method, TI_area_ratio_temp, &
  +                                           & TI_dst_tmp, TI_area_ratio, TIturbs(turbI))
                                              
                           !print *, "rotorDiameter after TI calcs", rotorDiameter
                       end if
  
positive.code:         do, p=1, nRotorPoints
          
              ! initialize the TI_area_ratio to 0.0 for each turbine
              ti_area_ratio = 0.0_dp
      
              ! initialize deficit summation term to zero
              deficit_sum = 0.0_dp
  
negative.code:                         ! vertical spread at discontinuity point
                          call sigmaz_func(kz_local(turb), deltax0_dp, rotorDiameter(turb), sigmaz_dp)
  
                          ! print *, "inputs in parent: ", deltay, deltaz, Ct(turb), yaw(turb), sigmay_dp, sigmaz_dp, &
  !                                          & rotorDiameter(turb), x, discontinuity_point, sigmay_dp, sigmaz_dp, &
  !                                          & wake_model_version, kz_local, x0, &
  !                                          & opt_exp_fac
  
                          ! velocity deficit in the nearwake (linear model)
                          call deltav_near_wake_lin_func(deltay, deltaz, &
  