Fix AeroDyn WriteOutput linearization (and cleanup some code)

modules/aerodyn/src/AeroAcoustics_TNO.f90

@@ -169,7 +169,6 @@ END FUNCTION f_int2
 FUNCTION Pressure(k1_in)
     ! Variables
    REAL(TNOKi)  :: a,b,answer
-   REAL(TNOKi)  :: omega
    REAL(TNOKi)  :: abserr,resabs,resasc
    REAL(TNOKi)  :: k1_in
    real(TNOKi)  :: Pressure
@@ -189,7 +188,7 @@ FUNCTION Pressure(k1_in)
 
     CALL slatec_qk61(f_int1,a,b,answer,abserr,resabs,resasc)
 
-    Pressure = 4.*rho**2*k1**2./(k1**2.+k3**2.)*answer
+    Pressure = 4.0_TNOKi*rho**2 * k1**2 / (k1**2 + k3**2)*answer
 
    RETURN
 END FUNCTION Pressure

modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90

@@ -112,11 +112,10 @@ MODULE AeroDyn_AllBldNdOuts_IO
 !----------------------------------------------------------------------------------------------------------------------------------
 !> This subroutine populates the headers with the blade node outputs.  The iteration cycle is blade:node:channel (channel iterated
 !! fastest).  If this iteration order is changed, it should be changed in the Calc_WriteAllBldNdOutput routine as well.
-SUBROUTINE AllBldNdOuts_InitOut( InitOut, p, p_AD, InputFileData, ErrStat, ErrMsg )
+SUBROUTINE AllBldNdOuts_InitOut( InitOut, p, InputFileData, ErrStat, ErrMsg )
 
    TYPE(RotInitOutputType),      INTENT(INOUT)  :: InitOut                          ! output data
    TYPE(RotParameterType),       INTENT(IN   )  :: p                                ! The rotor parameters
-   TYPE(AD_ParameterType),       INTENT(IN   )  :: p_AD                             ! The module parameters
    TYPE(RotInputFile),           INTENT(IN   )  :: InputFileData                    ! All the data in the AeroDyn input file (want Blade Span for channel name)
    INTEGER(IntKi),               INTENT(  OUT)  :: ErrStat                          ! The error status code
    CHARACTER(*),                 INTENT(  OUT)  :: ErrMsg                           ! The error message, if an error occurred
@@ -1253,7 +1252,7 @@ SUBROUTINE AllBldNdOuts_SetParameters( InputFileData, p, p_AD, ErrStat, ErrMsg )
    IF ( (InputFileData%BldNd_BladesOut < 0_IntKi) ) then
       p%BldNd_BladesOut = 0_IntKi
    ELSE IF ((InputFileData%BldNd_BladesOut > p%NumBlades) ) THEN
-      CALL SetErrStat( ErrID_Warn, " Number of blades to output data at all blade nodes (BldNd_BladesOut) must be less than "//TRIM(Num2LStr(p%NumBlades))//".", ErrStat, ErrMsg, RoutineName)
+      CALL SetErrStat( ErrID_Warn, " Number of blades to output data at all blade nodes (BldNd_BladesOut) must be no more than the total number of blades, "//TRIM(Num2LStr(p%NumBlades))//".", ErrStat, ErrMsg, RoutineName)
       p%BldNd_BladesOut = p%NumBlades ! NOTE: we are forgiving and plateau to numBlades
    ELSE
       p%BldNd_BladesOut = InputFileData%BldNd_BladesOut

modules/aerodyn/src/AeroDyn_Driver_Subs.f90

@@ -1396,8 +1396,8 @@ subroutine Dvr_ReadInputFile(fileName, dvr, errStat, errMsg )
          call ParseAry(FileInfo_In, CurLine, 'nacOrigin_t'//sWT        , wt%nac%origin_t, 3       , errStat2, errMsg2, unEc); if(Failed()) return
          call ParseAry(FileInfo_In, CurLine, 'hubOrigin_n'//sWT        , wt%hub%origin_n, 3       , errStat2, errMsg2, unEc); if(Failed()) return
          call ParseAry(FileInfo_In, CurLine, 'hubOrientation_n'//sWT   , wt%hub%orientation_n, 3  , errStat2, errMsg2, unEc); if(Failed()) return
-         wt%hub%orientation_n   = wt%hub%orientation_n*Pi/180_ReKi
-         wt%orientationInit     = wt%orientationInit*Pi/180_ReKi
+         wt%hub%orientation_n   = wt%hub%orientation_n*D2R
+         wt%orientationInit     = wt%orientationInit*D2R
          ! Blades
          call ParseCom(FileInfo_In, CurLine, Line, errStat2, errMsg2, unEc); if(Failed()) return
          call ParseVar(FileInfo_In, CurLine, 'numBlades'//sWT , wt%numBlades, errStat2, errMsg2, unEc); if(Failed()) return
@@ -1848,7 +1848,7 @@ subroutine Dvr_InitializeDriverOutputs(dvr, errStat, errMsg)
 
    dvr%out%WriteOutputHdr(j) = 'ShearExp'
    if (dvr%CompInflow==1) then
-      dvr%out%WriteOutputUnt(j) = '(NVALID)'; j=j+1
+      dvr%out%WriteOutputUnt(j) = '(INVALID)'; j=j+1
    else
       dvr%out%WriteOutputUnt(j) = '(-)'; j=j+1
    endif

modules/aerodyn/src/AeroDyn_IO.f90

@@ -1654,11 +1654,12 @@ END FUNCTION Calc_Chi0
 
 
 !----------------------------------------------------------------------------------------------------------------------------------
-SUBROUTINE Calc_WriteOutput( p, p_AD, u, m, m_AD, y, OtherState, xd, indx, iRot, ErrStat, ErrMsg )
+SUBROUTINE Calc_WriteOutput( p, p_AD, u, x, m, m_AD, y, OtherState, xd, indx, iRot, ErrStat, ErrMsg )
 
    TYPE(RotParameterType),       INTENT(IN   )  :: p                                 ! The rotor parameters
    TYPE(AD_ParameterType),       INTENT(IN   )  :: p_AD                              ! The module parameters
    TYPE(RotInputType),           INTENT(IN   )  :: u                                 ! inputs
+   TYPE(RotContinuousStateType), INTENT(IN   )  :: x                                 !< Continuous states at t
    TYPE(RotMiscVarType),         INTENT(INOUT)  :: m                                 ! misc variables
    TYPE(AD_MiscVarType),         INTENT(INOUT)  :: m_AD                              ! misc variables
    TYPE(RotOutputType),          INTENT(IN   )  :: y                                 ! outputs
@@ -2079,7 +2080,7 @@ SUBROUTINE ReadInputFiles( InputFileName, InputFileData, Default_DT, OutFileRoot
          iBld = iBld+1 ! Increment blade counter
       END DO
 
-   ENDDO ! Loop on rotors
+   end do ! loop on rotors
 
 
    CALL Cleanup ( )
@@ -2617,8 +2618,12 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg )
    CHARACTER(*), PARAMETER      :: FmtDatT   = '(A,T35,1(:,F13.8))'                ! Format for outputting time steps.
 
    CHARACTER(30)                :: OutPFmt                                         ! Format to print list of selected output channels to summary file
+   CHARACTER(30)                :: OutPFmtS                                        ! Format to print list of selected output channels to summary file
    CHARACTER(100)               :: Msg                                             ! temporary string for writing appropriate text to summary file
 
+   CHARACTER(ChanLen),PARAMETER :: TitleStr(2) = (/ 'Parameter', 'Units    ' /)
+   CHARACTER(ChanLen),PARAMETER :: TitleStrLines(2) = (/ '---------------', '---------------' /)
+
    ! Open the summary file and give it a heading.
 
    CALL GetNewUnit( UnSu, ErrStat, ErrMsg )
@@ -2803,7 +2808,7 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg )
             Msg = 'Stieg Oye dynamic stall model'
          case (UA_BV)
             Msg = 'Boeing-Vertol dynamic stall model (e.g. used in CACTUS)'
-         case default      
+         case default
             Msg = 'unknown'      
       end select
       WRITE (UnSu,Ec_IntFrmt) InputFileData%UAMod, 'UAMod', 'Unsteady Aero Model: '//TRIM(Msg)
@@ -2845,23 +2850,25 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg )
    end if
 
 
-   OutPFmt =  '( 15x, I4, 2X, A '//TRIM(Num2LStr(ChanLen))//',1 X, A'//TRIM(Num2LStr(ChanLen))//' )'
+   OutPFmt  = '( 15x, I4, 3X,A '//TRIM(Num2LStr(ChanLen))//',1 X, A'//TRIM(Num2LStr(ChanLen))//' )'
+   OutPFmtS = '( 15x, A4, 3X,A '//TRIM(Num2LStr(ChanLen))//',1 X, A'//TRIM(Num2LStr(ChanLen))//' )'
+   WRITE (UnSu,'(15x,A)')
    WRITE (UnSu,'(15x,A)')  'Requested Output Channels:'
-   WRITE (UnSu,'(15x,A)')  'Col   Parameter       Units'
-   WRITE (UnSu,'(15x,A)')  '----  --------------  -----'
-
+   WRITE (UnSu,OutPFmtS )  "Col", TitleStr
+   WRITE (UnSu,OutPFmtS )  "---", TitleStrLines
    DO I = 0,p%NumOuts
       WRITE (UnSu,OutPFmt)  I, p%OutParam(I)%Name, p%OutParam(I)%Units
    END DO             
 
    WRITE (UnSu,'(15x,A)')
    WRITE (UnSu,'(15x,A)')
    WRITE (UnSu,'(15x,A)')  'Requested Output Channels at each blade station:'
-   WRITE (UnSu,'(15x,A)')  'Col   Parameter       Units'
-   WRITE (UnSu,'(15x,A)')  '----  --------------  -----'
+   WRITE (UnSu,OutPFmtS )  "Col", TitleStr
+   WRITE (UnSu,OutPFmtS )  "---", TitleStrLines
    DO I = 1,p%BldNd_NumOuts
       WRITE (UnSu,OutPFmt)  I, p%BldNd_OutParam(I)%Name, p%BldNd_OutParam(I)%Units
-   END DO             
+   END DO
+
 
    CLOSE(UnSu)
 
@@ -3403,7 +3410,7 @@ SUBROUTINE SetOutParam(OutList, p, p_AD, ErrStat, ErrMsg )
 
          ! BNClrnc is set only when we're computing the tower influence
       do i = 1,size(BNClrnc,2)  ! all blades (need to do this in a loop because we need the index of InvalidOutput to be an array of rank one)
-         InvalidOutput( BNClrnc(:,i) ) = .true.  
+         InvalidOutput( BNClrnc(:,i) ) = .true.
       end do
 
    end if

modules/aerodyn/src/BEMT_Registry.txt

@@ -47,6 +47,7 @@ typedef   ^                            ^                             ReKi
 typedef   ^                            ^                             ReKi                     zLocal                          {:}{:}     -         -        "Distance to blade node, measured along the blade" m
 typedef   ^                            ^                             ReKi                     zTip                            {:}        -         -        "Distance to blade tip, measured along the blade" m
 typedef   ^                            ^                             ReKi                     rLocal                          {:}{:}     -         -        "Radial distance to blade node from the center of rotation, measured in the rotor plane, needed for DBEMT" m
+typedef   ^                            ^                             ReKi                     rTipFix                         {:}        -         -        "Nominally the coned rotor diameter (without prebend)" m
 typedef   ^                            ^                             INTEGER                  UAMod                           -          -         -        "Model for the dynamic stall equations [1 = Leishman/Beddoes, 2 = Gonzalez, 3 = Minnema]"    -
 typedef   ^                            ^                             LOGICAL                  UA_Flag                         -          -         -        "logical flag indicating whether to use UnsteadyAero"        -
 typedef   ^                            ^                             LOGICAL                  Flookup                         -          -         -        "Use table lookup for f' and f'' "    -

modules/aerodyn/src/BEMT_Types.f90

@@ -57,6 +57,7 @@ MODULE BEMT_Types
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: zLocal      !< Distance to blade node, measured along the blade [m]
     REAL(ReKi) , DIMENSION(:), ALLOCATABLE  :: zTip      !< Distance to blade tip, measured along the blade [m]
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: rLocal      !< Radial distance to blade node from the center of rotation, measured in the rotor plane, needed for DBEMT [m]
+    REAL(ReKi) , DIMENSION(:), ALLOCATABLE  :: rTipFix      !< Nominally the coned rotor diameter (without prebend) [m]
     INTEGER(IntKi)  :: UAMod      !< Model for the dynamic stall equations [1 = Leishman/Beddoes, 2 = Gonzalez, 3 = Minnema] [-]
     LOGICAL  :: UA_Flag      !< logical flag indicating whether to use UnsteadyAero [-]
     LOGICAL  :: Flookup      !< Use table lookup for f' and f''  [-]
@@ -296,6 +297,18 @@ SUBROUTINE BEMT_CopyInitInput( SrcInitInputData, DstInitInputData, CtrlCode, Err
     END IF
   END IF
     DstInitInputData%rLocal = SrcInitInputData%rLocal
+ENDIF
+IF (ALLOCATED(SrcInitInputData%rTipFix)) THEN
+  i1_l = LBOUND(SrcInitInputData%rTipFix,1)
+  i1_u = UBOUND(SrcInitInputData%rTipFix,1)
+  IF (.NOT. ALLOCATED(DstInitInputData%rTipFix)) THEN 
+    ALLOCATE(DstInitInputData%rTipFix(i1_l:i1_u),STAT=ErrStat2)
+    IF (ErrStat2 /= 0) THEN 
+      CALL SetErrStat(ErrID_Fatal, 'Error allocating DstInitInputData%rTipFix.', ErrStat, ErrMsg,RoutineName)
+      RETURN
+    END IF
+  END IF
+    DstInitInputData%rTipFix = SrcInitInputData%rTipFix
 ENDIF
     DstInitInputData%UAMod = SrcInitInputData%UAMod
     DstInitInputData%UA_Flag = SrcInitInputData%UA_Flag
@@ -359,6 +372,9 @@ SUBROUTINE BEMT_DestroyInitInput( InitInputData, ErrStat, ErrMsg )
 IF (ALLOCATED(InitInputData%rLocal)) THEN
   DEALLOCATE(InitInputData%rLocal)
 ENDIF
+IF (ALLOCATED(InitInputData%rTipFix)) THEN
+  DEALLOCATE(InitInputData%rTipFix)
+ENDIF
 IF (ALLOCATED(InitInputData%UAOff_innerNode)) THEN
   DEALLOCATE(InitInputData%UAOff_innerNode)
 ENDIF
@@ -445,6 +461,11 @@ SUBROUTINE BEMT_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrM
   IF ( ALLOCATED(InData%rLocal) ) THEN
     Int_BufSz   = Int_BufSz   + 2*2  ! rLocal upper/lower bounds for each dimension
       Re_BufSz   = Re_BufSz   + SIZE(InData%rLocal)  ! rLocal
+  END IF
+  Int_BufSz   = Int_BufSz   + 1     ! rTipFix allocated yes/no
+  IF ( ALLOCATED(InData%rTipFix) ) THEN
+    Int_BufSz   = Int_BufSz   + 2*1  ! rTipFix upper/lower bounds for each dimension
+      Re_BufSz   = Re_BufSz   + SIZE(InData%rTipFix)  ! rTipFix
   END IF
       Int_BufSz  = Int_BufSz  + 1  ! UAMod
       Int_BufSz  = Int_BufSz  + 1  ! UA_Flag
@@ -629,6 +650,21 @@ SUBROUTINE BEMT_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrM
           Re_Xferred = Re_Xferred + 1
         END DO
       END DO
+  END IF
+  IF ( .NOT. ALLOCATED(InData%rTipFix) ) THEN
+    IntKiBuf( Int_Xferred ) = 0
+    Int_Xferred = Int_Xferred + 1
+  ELSE
+    IntKiBuf( Int_Xferred ) = 1
+    Int_Xferred = Int_Xferred + 1
+    IntKiBuf( Int_Xferred    ) = LBOUND(InData%rTipFix,1)
+    IntKiBuf( Int_Xferred + 1) = UBOUND(InData%rTipFix,1)
+    Int_Xferred = Int_Xferred + 2
+
+      DO i1 = LBOUND(InData%rTipFix,1), UBOUND(InData%rTipFix,1)
+        ReKiBuf(Re_Xferred) = InData%rTipFix(i1)
+        Re_Xferred = Re_Xferred + 1
+      END DO
   END IF
     IntKiBuf(Int_Xferred) = InData%UAMod
     Int_Xferred = Int_Xferred + 1
@@ -866,6 +902,24 @@ SUBROUTINE BEMT_UnPackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, E
           Re_Xferred = Re_Xferred + 1
         END DO
       END DO
+  END IF
+  IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN  ! rTipFix not allocated
+    Int_Xferred = Int_Xferred + 1
+  ELSE
+    Int_Xferred = Int_Xferred + 1
+    i1_l = IntKiBuf( Int_Xferred    )
+    i1_u = IntKiBuf( Int_Xferred + 1)
+    Int_Xferred = Int_Xferred + 2
+    IF (ALLOCATED(OutData%rTipFix)) DEALLOCATE(OutData%rTipFix)
+    ALLOCATE(OutData%rTipFix(i1_l:i1_u),STAT=ErrStat2)
+    IF (ErrStat2 /= 0) THEN 
+       CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%rTipFix.', ErrStat, ErrMsg,RoutineName)
+       RETURN
+    END IF
+      DO i1 = LBOUND(OutData%rTipFix,1), UBOUND(OutData%rTipFix,1)
+        OutData%rTipFix(i1) = ReKiBuf(Re_Xferred)
+        Re_Xferred = Re_Xferred + 1
+      END DO
   END IF
     OutData%UAMod = IntKiBuf(Int_Xferred)
     Int_Xferred = Int_Xferred + 1