Add SeaSt_C_GetFluidDynP routine to SeaState_C_Binding interface

modules/seastate/src/SeaSt_WaveField.f90

@@ -17,6 +17,7 @@ MODULE SeaSt_WaveField
 PUBLIC WaveField_GetNodeWaveKin
 PUBLIC WaveField_GetNodeWaveVelAcc
 PUBLIC WaveField_GetWaveKin
+PUBLIC WaveField_GetDynP
 
 CONTAINS
 
@@ -310,6 +311,103 @@ logical function Failed()
 END SUBROUTINE WaveField_GetNodeWaveKin
 
 
+!-------------------- Subroutine for dynamic pressure --------------------!
+!NOTE: this is a subset of WaveField_GetNodeWaveKin
+SUBROUTINE WaveField_GetDynP( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, FDynP, ErrStat, ErrMsg )
+   type(SeaSt_WaveFieldType),          intent(in   ) :: WaveField
+   type(GridInterp_MiscVarType),       intent(inout) :: WaveField_m
+   real(DbKi),                         intent(in   ) :: Time
+   real(ReKi),                         intent(in   ) :: pos(3)
+   logical,                            intent(in   ) :: forceNodeInWater
+   real(SiKi),                         intent(  out) :: FDynP
+   integer(IntKi),                     intent(  out) :: nodeInWater
+   integer(IntKi),                     intent(  out) :: ErrStat ! Error status of the operation
+   character(*),                       intent(  out) :: ErrMsg  ! Error message if errStat /= ErrID_None
+
+   real(ReKi)                                        :: posXY(2), posPrime(3), posXY0(3)
+   character(*),                       parameter     :: RoutineName = 'WaveField_GetDynP'
+   integer(IntKi)                                    :: errStat2
+   character(ErrMsgLen)                              :: errMsg2
+
+   ! Temporary vars not kept
+   real(SiKi)                                        :: WaveElev1
+   real(SiKi)                                        :: WaveElev2
+   real(SiKi)                                        :: WaveElev
+
+   ErrStat   = ErrID_None
+   ErrMsg    = ""
+
+   posXY    = pos(1:2)
+   posXY0   = (/pos(1),pos(2),0.0_ReKi/)
+
+   ! Wave elevation
+   WaveElev1 = WaveField_GetNodeWaveElev1( WaveField, WaveField_m, Time, pos, ErrStat2, ErrMsg2 ); if (Failed()) return;
+   WaveElev2 = WaveField_GetNodeWaveElev2( WaveField, WaveField_m, Time, pos, ErrStat2, ErrMsg2 ); if (Failed()) return;
+   WaveElev  = WaveElev1 + WaveElev2
+
+   IF (WaveField%WaveStMod == 0) THEN ! No wave stretching
+
+      IF ( pos(3) <= 0.0_ReKi) THEN ! Node is at or below the SWL
+         nodeInWater = 1_IntKi
+         ! Use location to obtain interpolated values of kinematics
+         CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, pos, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
+         FDynP = GridInterp4D   ( WaveField%WaveDynP, WaveField_m )
+      ELSE ! Node is above the SWL
+         nodeInWater = 0_IntKi
+         FDynP       = 0.0
+      END IF
+
+   ELSE ! Wave stretching enabled
+      IF ( (pos(3) <= WaveElev) .OR. forceNodeInWater ) THEN ! Node is submerged
+         nodeInWater = 1_IntKi
+         IF ( WaveField%WaveStMod < 3 ) THEN ! Vertical or extrapolated wave stretching
+
+            IF ( pos(3) <= 0.0_SiKi) THEN ! Node is below the SWL - evaluate wave dynamics as usual
+               ! Use location to obtain interpolated values of kinematics
+               CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, pos, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
+               FDynP = GridInterp4D   ( WaveField%WaveDynP, WaveField_m )
+            ELSE ! Node is above SWL - need wave stretching
+
+               ! Vertical wave stretching
+               CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
+               FDynP = GridInterp4D   ( WaveField%WaveDynP, WaveField_m )
+
+               ! Extrapoled wave stretching
+               IF (WaveField%WaveStMod == 2) THEN
+                  CALL WaveField_Interp_Setup3D( Time+WaveField%WaveTimeShift, posXY, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
+                  FDynP = FDynP + GridInterp3D   ( WaveField%PWaveDynP0, WaveField_m ) * pos(3)
+               END IF
+
+            END IF ! Node is submerged
+
+         ELSE ! Wheeler stretching - no need to check whether the node is above or below SWL
+
+            ! Map the node z-position linearly from [-EffWtrDpth,m%WaveElev(j)] to [-EffWtrDpth,0]
+            posPrime    = pos
+            posPrime(3) = WaveField%EffWtrDpth*(WaveField%EffWtrDpth+pos(3))/(WaveField%EffWtrDpth+WaveElev)-WaveField%EffWtrDpth
+            posPrime(3) = MIN( posPrime(3), 0.0_ReKi) ! Clamp z-position to zero. Needed when forceNodeInWater=.TRUE.
+
+            ! Obtain the wave-field variables by interpolation with the mapped position.
+            CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, posPrime, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
+            FDynP = GridInterp4D   ( WaveField%WaveDynP, WaveField_m )
+         END IF
+
+      ELSE ! Node is out of water - zero-out all wave dynamics
+
+         nodeInWater = 0_IntKi
+         FDynP       = 0.0
+
+      END IF ! If node is in or out of water
+   END IF ! If wave stretching is on or off
+
+contains
+   logical function Failed()
+      call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
+      Failed = ErrStat >= AbortErrLev
+   end function
+END SUBROUTINE WaveField_GetDynP
+
+
 !-------------------- Subroutine for wave field velocity only --------------------!
 SUBROUTINE WaveField_GetNodeWaveVelAcc( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, FV, FA, ErrStat, ErrMsg )
    type(SeaSt_WaveFieldType),          intent(in   ) :: WaveField

modules/seastate/src/SeaState_C_Binding.f90

@@ -44,6 +44,7 @@ MODULE SeaState_C_Binding
    PUBLIC :: SeaSt_C_GetDens
    PUBLIC :: SeaSt_C_GetDpth
    PUBLIC :: SeaSt_C_GetMSL2SWL
+   PUBLIC :: SeaSt_C_GetDynPressure
 
    !------------------------------------------------------------------------------------
    !  Debugging: DebugLevel -- passed at PreInit
@@ -898,6 +899,73 @@ subroutine CheckWaveFieldPtr(callingRoutine,valid,ErrStat,ErrMsg)
 end subroutine
 
 
+!> return the surface elevation and normal vector at a point.
+subroutine SeaSt_C_GetDynPressure(Time_C, Pos_C, DynP_C, ErrStat_C,ErrMsg_C) BIND (C, NAME='SeaSt_C_GetDynPressure')
+#ifndef IMPLICIT_DLLEXPORT
+!DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetDynPressure
+!GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetDynPressure
+#endif
+   real(c_double),            intent(in   ) :: Time_C
+   real(c_float),             intent(in   ) :: Pos_c(3)
+   real(c_float),             intent(  out) :: DynP_C
+   integer(c_int),            intent(  out) :: ErrStat_C
+   character(kind=c_char),    intent(  out) :: ErrMsg_C(ErrMsgLen_C)
+
+   real(DbKi)                 :: Time
+   real(ReKi)                 :: Pos(3)
+   real(SiKi)                 :: FDynP
+   integer                    :: ErrStat
+   character(ErrMsgLen)       :: ErrMsg
+   character(*), parameter    :: RoutineName = 'SeaSt_C_GetDynPressure'
+   logical                    :: valid
+   ! Temporary vars for calling GetNodeWaveKin -- remove if new interface developed
+   integer(IntKi)             :: nodeInWater
+
+   if (DebugLevel > 1) call ShowPassedData()
+
+   ! convert position and time to fortran types
+   Time = real(Time_C, DbKi)
+   Pos = real(Pos_C(1:3), ReKi)
+
+   ! verify there is actually wavefield data
+   call CheckWaveFieldPtr(RoutineName, valid, ErrStat, ErrMsg)
+   if (.not. valid) then
+      call Cleanup()
+      return
+   endif
+
+   ! get wave elevation (total combined first and second order)
+   ! Notes:
+   !     - if position is outside the wave field boundary, it will simply return boundary edge value
+   !     - time must be positive or a fatal error occurs
+   !     - forceNodeInWater = .false. -- perhaps update this later to allow throug interface?
+   call WaveField_GetDynP( p%WaveField, m%WaveField_m, Time, pos, .false., nodeInWater, FDynP, ErrStat, ErrMsg )
+
+
+   ! Store resulting elevation as C type
+   DynP_C = real(FDynP,c_float)
+
+   if (DebugLevel > 1) call ShowReturnData()
+   call Cleanup()
+   return
+contains
+   subroutine Cleanup()    ! NOTE: we are ignoring any error reporting from here
+      CALL SetErrStat_F2C(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C)
+   end subroutine Cleanup
+   subroutine ShowPassedData()
+      call WrScr("-----------------------------------------------------------")
+      call WrScr("Interface debugging:  SeaSt_C_GetDynPressure")
+      call WrScr("   --------------------------------------------------------")
+      call WrScr("   Time_C                 -> "//trim(Num2LStr(Time_C)))
+      call WrScr("   Pos_C                  -> ("//trim(Num2LStr(Pos_C(1)))//","//trim(Num2LStr(Pos_C(2)))//","//trim(Num2LStr(Pos_C(3)))//")")
+   end subroutine ShowPassedData
+   subroutine ShowReturnData()
+      call WrScr("   DynP_C                 <- "//trim(Num2LStr(DynP_C)))
+      call WrScr("-----------------------------------------------------------")
+   end subroutine ShowReturnData
+end subroutine SeaSt_C_GetDynPressure
+
+
 !FIXME: the following visualization writer should be merged into the library vtk.f90
 !        this is a modified duplicate of the routine from FAST_Subs by the same name.
 !FIXME: this routine currently only grabs the closest timestep and does not do any interpolation