anchor.diff:          ! initialize deficit summation term to zero
           deficit_sum = 0.0_dp
           
  -!         print *, turbineXw(turbI)
  -        
           do, u=1, nTurbines ! at turbineX-locations
               
               ! get index of upstream turbine
  
positive.code:                 end if
                  
              end if
          end do
          
          ! final velocity calculation for turbine turbI
          wtVelocity(turbI) = wind_speed - deficit_sum
  
      !! make sure turbine inflow velocity is non-negative
  !             if (wtVelocity(turbI) .lt. 0.0_dp) then 
  !                 wtVelocity(turbI) = 0.0_dp
  !             end if
      !print *, "fortran"
      
      end do
      
      
  
  end subroutine porteagel_analyze
  
  
negative.code:                 ! calculate wake spreading parameter at each turbine if desired
                  if (calc_k_star .eqv. .true.) then
                      call k_star_func(TIturbs(turb), k_star(turb))
                      ky_local = k_star(turb)
                      kz_local = k_star(turb)
                  end if
  
                  ! determine the initial wake angle at the onset of far wake
                  call theta_c_0_func(yaw(turb), Ct(turb), theta_c_0)
          
                  ! horizontal spread
                  call sigmay_func(ky_local, deltax0, rotorDiameter(turb), yaw(turb), sigmay)
                  !print *, "sigmay ", sigmay
                  ! vertical spread
                  call sigmaz_func(kz_local, deltax0, rotorDiameter(turb), sigmaz)
                  !print *, "sigmaz ", sigmaz
                  ! horizontal cross-wind wake displacement from hub
                  call wake_offset_func(rotorDiameter(turb), theta_c_0, x0, yaw(turb), &
                                       & ky_local, kz_local, Ct(turb), sigmay, sigmaz, wake_offset)
                  !print *, "wake_offset ", wake_offset                 
                  ! cross wind distance from downstream hub location to wake center
                  deltay = turbineYw(turbI) - (turbineYw(turb) + wake_offset)
              
                  ! cross wind distance from hub height to height of point of interest
                  deltaz = turbineZ(turbI) - turbineZ(turb)
                  
                  ! far wake region
                  if (x >= x0) then
      
                      ! velocity difference in the wake
                      call deltav_func(deltay, deltaz, Ct(turb), & 
                                       & yaw(turb), sigmay, sigmaz, rotorDiameter(turb), deltav)  
  
                  ! near wake region (linearized)
                  else
                  
                      ! horizontal spread at far wake onset
                      call sigmay_func(ky_local, 0.0_dp, rotorDiameter(turb), yaw(turb), sigmay0)
  