From ce8351ea2fbb356ae105d6e38ae00a630268c305 Mon Sep 17 00:00:00 2001
From: andrew-platt <andy.platt@nrel.gov>
Date: Mon, 25 Nov 2024 13:05:51 -0700
Subject: [PATCH] IfW: use single point for disk average if no radius given

---
 modules/inflowwind/src/InflowWind_Subs.f90 | 39 +++++++++++++---------
 1 file changed, 23 insertions(+), 16 deletions(-)

diff --git a/modules/inflowwind/src/InflowWind_Subs.f90 b/modules/inflowwind/src/InflowWind_Subs.f90
index 8b5a18b0c..c8349ffbd 100644
--- a/modules/inflowwind/src/InflowWind_Subs.f90
+++ b/modules/inflowwind/src/InflowWind_Subs.f90
@@ -967,6 +967,7 @@ SUBROUTINE InflowWind_SetParameters( InitInp, InputFileData, p, m, ErrStat, ErrM
       ! Temporary variables
    INTEGER(IntKi)                                     :: TmpErrStat           !< Temporary error status  for subroutine and function calls
    CHARACTER(ErrMsgLen)                               :: TmpErrMsg            !< Temporary error message for subroutine and function calls
+   integer(IntKi)                                     :: NumPtsAvg            !< Number of points to use for disk average vel (1 if no radius)
 
       ! Local variables
    INTEGER(IntKi)                                     :: I                    !< Generic counter
@@ -1029,28 +1030,34 @@ SUBROUTINE InflowWind_SetParameters( InitInp, InputFileData, p, m, ErrStat, ErrM
    CALL AllocAry( m%y_Hub%VelocityUVW, 3, 1, "Array of velocities for hub values", TmpErrStat, TmpErrMsg )
       CALL SetErrStat(TmpErrStat,TmpErrMsg,ErrStat,ErrMsg,RoutineName)
       
-   CALL AllocAry( m%u_Avg%PositionXYZ, 3, IfW_NumPtsAvg, "Array of positions for rotor-averaged values", TmpErrStat, TmpErrMsg )
+   ! Disk Velocity calculations
+   if (InitInp%RadAvg < 0.0_ReKi)  then
+      NumPtsAvg = 1_IntKi              ! Use only hub point
+   else
+      NumPtsAvg = IfW_NumPtsAvg        ! Use a field of points for disk average calculations
+   endif
+
+   CALL AllocAry( m%u_Avg%PositionXYZ, 3, NumPtsAvg, "Array of positions for rotor-averaged values", TmpErrStat, TmpErrMsg )
       CALL SetErrStat(TmpErrStat,TmpErrMsg,ErrStat,ErrMsg,RoutineName)
-   CALL AllocAry( m%y_Avg%VelocityUVW, 3, IfW_NumPtsAvg, "Array of velocities for rotor-averaged values", TmpErrStat, TmpErrMsg )
+   CALL AllocAry( m%y_Avg%VelocityUVW, 3, NumPtsAvg, "Array of velocities for rotor-averaged values", TmpErrStat, TmpErrMsg )
       CALL SetErrStat(TmpErrStat,TmpErrMsg,ErrStat,ErrMsg,RoutineName)
-   CALL AllocAry( p%PositionAvg, 3, IfW_NumPtsAvg, "Array of positions for computing average wind speed", TmpErrStat, TmpErrMsg )
+   CALL AllocAry( p%PositionAvg, 3, NumPtsAvg, "Array of positions for computing average wind speed", TmpErrStat, TmpErrMsg )
       CALL SetErrStat(TmpErrStat,TmpErrMsg,ErrStat,ErrMsg,RoutineName)
       IF ( ErrStat>= AbortErrLev ) RETURN
       
       
    if (InitInp%RadAvg < 0.0_ReKi)  then
-      R  =  max(1.0_ReKi, InputFileData%Uniform_RefLength)/2.0_ReKi    ! We'll use this as a guess for the rotor radius
+      p%PositionAvg = 0.0_ReKi         ! Use only hub
    else
-      R = InitInp%RadAvg
+      ! Calculate a ring of points at 70 rotor radius
+      R = InitInp%RadAvg * 0.7_ReKi !70% radius
+      do i=1,NumPtsAvg
+         theta = pi +(i-1)*TwoPi/NumPtsAvg
+         p%PositionAvg(1,i) = 0.0_ReKi          ! Hub X (perpindicular to rotor plane)
+         p%PositionAvg(2,i) = R*cos(theta)      ! Hub Y
+         p%PositionAvg(3,i) = R*sin(theta)      ! Hub Z (in vertical plane when azimuth=0)
+      end do
    end if
-   R = R * 0.7_ReKi !70% radius
-   
-   do i=1,IfW_NumPtsAvg
-      theta = pi +(i-1)*TwoPi/IfW_NumPtsAvg
-      p%PositionAvg(1,i) = 0.0_ReKi          ! Hub X (perpindicular to rotor plane)
-      p%PositionAvg(2,i) = R*cos(theta)      ! Hub Y
-      p%PositionAvg(3,i) = R*sin(theta)      ! Hub Z (in vertical plane when azimuth=0)
-   end do
    
    p%OutputAccel = InitInp%OutputAccel
 
@@ -1709,17 +1716,17 @@ SUBROUTINE InflowWind_GetSpatialAverage( Time, InputData, p, x, xd, z, OtherStat
    m%u_Avg%HubPosition    = InputData%HubPosition
    m%u_Avg%HubOrientation = InputData%HubOrientation
    
-   do i=1,IfW_NumPtsAvg
+   do i=1,size(m%u_Avg%PositionXYZ,DIM=2)
       m%u_Avg%PositionXYZ(:,i) = matmul(InputData%HubOrientation,p%PositionAvg(:,i)) + InputData%HubPosition
    end do
 
    CALL CalculateOutput( Time, m%u_Avg, p, x, xd, z, OtherStates, m%y_Avg, m, FillWrOut, ErrStat2, ErrMsg2 )
       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
 
-   do i=1,IfW_NumPtsAvg
+   do i=1,size(m%u_Avg%PositionXYZ,DIM=2)
       MeanVelocity = MeanVelocity + m%y_Avg%VelocityUVW(:,i)
    end do
-   MeanVelocity = MeanVelocity / REAL(IfW_NumPtsAvg,ReKi)
+   MeanVelocity = MeanVelocity / REAL(size(m%u_Avg%PositionXYZ,DIM=2),ReKi)
         
    
    RETURN