diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index a88200f2ac..1192036db1 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -2335,7 +2335,9 @@ SUBROUTINE AllocateNodeLoadVariables(InitInp, p, m, NNodes, errStat, errMsg ) ! Initialize errStat errStat = ErrID_None errMsg = "" - + + call AllocAry( m%DispNodePosHdn, 3, NNodes , 'm%DispNodePosHdn', errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) + call AllocAry( m%DispNodePosHst, 3, NNodes , 'm%DispNodePosHst', errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) call AllocAry( m%nodeInWater , NNodes , 'm%nodeInWater' , errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) call AllocAry( m%vrel , 3, NNodes , 'm%vrel' , errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) call AllocAry( m%FV , 3, NNodes , 'm%FV' , errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) @@ -2366,6 +2368,8 @@ SUBROUTINE AllocateNodeLoadVariables(InitInp, p, m, NNodes, errStat, errMsg ) if (errStat >= AbortErrLev) return + m%DispNodePosHdn = 0.0_ReKi + m%DispNodePosHst = 0.0_ReKi m%nodeInWater = 0 m%vrel = 0.0_ReKi m%FV = 0.0_ReKi @@ -2571,42 +2575,21 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, errMsg = "" Imat = 0.0_ReKi g = p%Gravity - WtrDpth = p%WtrDpth + p%MSL2SWL ! Water depth measured from the free surface + WtrDpth = p%WtrDpth + p%MSL2SWL ! Water depth measured from the still water level + + !=============================================================================================== + ! Get displaced positions of the hydrodynamic nodes + CALL GetDisplacedNodePosition( .FALSE., m%DispNodePosHdn ) ! For hydrodynamic loads; depends on WaveDisp and WaveStMod + CALL GetDisplacedNodePosition( .TRUE. , m%DispNodePosHst ) ! For hydrostatic loads; always use actual displaced position - !InterpolationSlope = GetInterpolationSlope(Time, p, m, IntWrapIndx) - !=============================================================================================== ! Calculate the fluid kinematics at all mesh nodes and store for use in the equations below - + CALL WaveField_GetWaveKin( p%WaveField, Time, m%DispNodePosHdn, .FALSE., m%nodeInWater, m%WaveElev1, m%WaveElev2, m%WaveElev, m%FDynP, m%FV, m%FA, m%FAMCF, ErrStat2, ErrMsg2 ) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + ! Compute fluid velocity relative to the structure DO j = 1, p%NNodes - IF (p%WaveDisp == 0 ) THEN - ! use the initial X,Y location - pos1(1) = u%Mesh%Position(1,j) - pos1(2) = u%Mesh%Position(2,j) - ELSE - ! Use current X,Y location - pos1(1) = u%Mesh%TranslationDisp(1,j) + u%Mesh%Position(1,j) - pos1(2) = u%Mesh%TranslationDisp(2,j) + u%Mesh%Position(2,j) - END IF - - IF (p%WaveStMod > 0 .AND. p%WaveDisp /= 0) THEN ! Wave stretching enabled - pos1(3) = u%Mesh%Position(3,j) + u%Mesh%TranslationDisp(3,j) - p%MSL2SWL ! Use the current Z location. - ELSE ! Wave stretching disabled - pos1(3) = u%Mesh%Position(3,j) - p%MSL2SWL ! We are intentionally using the undisplaced Z position of the node. - END IF - - ! Get the wave elevation and wave kinematics at each node - CALL WaveField_GetWaveKin( p%WaveField, Time, pos1, m%nodeInWater(j), m%WaveElev1(j), m%WaveElev2(j), m%WaveElev(j), FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - m%FDynP(j) = REAL(FDynP,ReKi) - m%FV(:, j) = REAL(FV, ReKi) - m%FA(:, j) = REAL(FA, ReKi) - IF (ALLOCATED(m%FAMCF)) THEN - m%FAMCF(:,j) = REAL(FAMCF,ReKi) - END IF m%vrel(:,j) = ( m%FV(:,j) - u%Mesh%TranslationVel(:,j) ) * m%nodeInWater(j) - - END DO ! j = 1, p%NNodes + END DO ! ============================================================================================== ! Calculate instantaneous loads on each member except for the hydrodynamic loads on member ends. @@ -2668,12 +2651,9 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, DO i = max(mem%i_floor,1), N ! loop through member elements that are not completely buried in the seabed ! calculate instantaneous incline angle and heading, and related trig values - ! the first and last NodeIndx values point to the corresponding Joint nodes idices which are at the start of the Mesh - - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i)) + u%Mesh%Position(:, mem%NodeIndx(i)) - pos1(3) = pos1(3) - p%MSL2SWL - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i+1)) + u%Mesh%Position(:, mem%NodeIndx(i+1)) - pos2(3) = pos2(3) - p%MSL2SWL + ! the first and last NodeIndx values point to the corresponding Joint nodes indices which are at the start of the Mesh + pos1 = m%DispNodePosHst(:, mem%NodeIndx(i )) + pos2 = m%DispNodePosHst(:, mem%NodeIndx(i+1)) call GetOrientationAngles( pos1, pos2, phi, sinPhi, cosPhi, tanPhi, sinBeta, cosBeta, k_hat, errStat2, errMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -2802,12 +2782,9 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, DO i = max(mem%i_floor,1), N ! loop through member elements that are not completely buried in the seabed ! calculate instantaneous incline angle and heading, and related trig values - ! the first and last NodeIndx values point to the corresponding Joint nodes idices which are at the start of the Mesh - - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i)) + u%Mesh%Position(:, mem%NodeIndx(i)) - pos1(3) = pos1(3) - p%MSL2SWL - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i+1)) + u%Mesh%Position(:, mem%NodeIndx(i+1)) - pos2(3) = pos2(3) - p%MSL2SWL + ! the first and last NodeIndx values point to the corresponding Joint nodes indices which are at the start of the Mesh + pos1 = m%DispNodePosHst(:,mem%NodeIndx(i )) + pos2 = m%DispNodePosHst(:,mem%NodeIndx(i+1)) call GetOrientationAngles( pos1, pos2, phi, sinPhi, cosPhi, tanPhi, sinBeta, cosBeta, k_hat, errStat2, errMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -2929,16 +2906,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, DO i = mem%i_floor+1,N ! loop through member nodes starting from the first node above seabed, but skip the last node which should not be submerged anyways ! Get positions of node i and i+1 - IF (p%WaveDisp /= 0) THEN ! Use current position - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i)) + u%Mesh%Position(:, mem%NodeIndx(i)) - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i+1)) + u%Mesh%Position(:, mem%NodeIndx(i+1)) - ELSE ! Use initial position - pos1 = u%Mesh%Position(:, mem%NodeIndx(i)) - pos2 = u%Mesh%Position(:, mem%NodeIndx(i+1)) - END if - ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - pos1(3) = pos1(3) - p%MSL2SWL - pos2(3) = pos2(3) - p%MSL2SWL + pos1 = m%DispNodePosHdn(:,mem%NodeIndx(i )) + pos2 = m%DispNodePosHdn(:,mem%NodeIndx(i+1)) ! Free surface elevation above or below node i and i+1 Zeta1 = m%WaveElev(mem%NodeIndx(i)) @@ -3048,8 +3017,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !----------------------------------------------------------------------------------------------------! ! Compute the distributed loads at the point of intersection between the member and the free surface ! !----------------------------------------------------------------------------------------------------! - ! Get wave dynamics at the free surface intersection - CALL WaveField_GetWaveKin( p%WaveField, Time, FSInt, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) + ! Get wave kinematics at the free-surface intersection. Set forceNodeInWater=.TRUE. to guarantee the free-surface intersection is in water. + CALL WaveField_GetNodeWaveKin( p%WaveField, Time, FSInt, .TRUE., nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FDynPFSInt = REAL(FDynP,ReKi) FVFSInt = REAL(FV, ReKi) @@ -3193,15 +3162,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ELSE IF ( MemSubStat .NE. 3_IntKi) THEN ! Skip members with centerline completely out of water !----------------------------No load smoothing----------------------------! DO i = mem%i_floor+1,N+1 ! loop through member nodes starting from the first node above seabed - ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - ! Using the initial z-position to be consistent with the evaluation of wave kinematics - IF (p%WaveStMod > 0 .AND. p%WaveDisp /= 0) THEN - ! Use current z-position - z1 = u%Mesh%Position(3, mem%NodeIndx(i)) + u%Mesh%TranslationDisp(3, mem%NodeIndx(i)) - p%MSL2SWL - ELSE - ! Use initial z-position - z1 = u%Mesh%Position(3, mem%NodeIndx(i)) - p%MSL2SWL - END IF + z1 = m%DispNodePosHdn(3, mem%NodeIndx(i)) !---------------------------------------------Compute deltal and h_c------------------------------------------! ! Cannot make any assumption about WaveStMod and member orientation IF ( m%NodeInWater(mem%NodeIndx(i)) .EQ. 0_IntKi ) THEN ! Node is out of water @@ -3219,11 +3180,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF ( m%NodeInWater(mem%NodeIndx(i-1)) .EQ. 1_IntKi ) THEN ! Node to the left is submerged deltalLeft = 0.5_ReKi * mem%dl ELSE ! Element i-1 crosses the free surface - IF (p%WaveStMod > 0 .AND. p%WaveDisp /= 0) THEN ! Use current z-position - z2 = u%Mesh%Position(3, mem%NodeIndx(i-1)) + u%Mesh%TranslationDisp(3, mem%NodeIndx(i-1)) - p%MSL2SWL - ELSE ! Use initial z-position - z2 = u%Mesh%Position(3, mem%NodeIndx(i-1)) - p%MSL2SWL - END IF + z2 = m%DispNodePosHdn(3, mem%NodeIndx(i-1)) IF ( p%WaveStMod > 0_IntKi ) THEN ! Wave stretching enabled zeta1 = m%WaveElev(mem%NodeIndx(i )) zeta2 = m%WaveElev(mem%NodeIndx(i-1)) @@ -3242,11 +3199,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF ( m%NodeInWater(mem%NodeIndx(i+1)) .EQ. 1_IntKi ) THEN ! Node to the right is submerged deltalRight = 0.5_ReKi * mem%dl ELSE ! Element i crosses the free surface - IF (p%WaveStMod > 0 .AND. p%WaveDisp /= 0) THEN ! Use current z-position - z2 = u%Mesh%Position(3, mem%NodeIndx(i+1)) + u%Mesh%TranslationDisp(3, mem%NodeIndx(i+1)) - p%MSL2SWL - ELSE ! Use initial z-position - z2 = u%Mesh%Position(3, mem%NodeIndx(i+1)) - p%MSL2SWL - END IF + z2 = m%DispNodePosHdn(3, mem%NodeIndx(i+1)) IF ( p%WaveStMod > 0_IntKi ) THEN ! Wave stretching enabled zeta1 = m%WaveElev(mem%NodeIndx(i )) zeta2 = m%WaveElev(mem%NodeIndx(i+1)) @@ -3348,12 +3301,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !-----------------------------------------------------------------------------------------------------! ! reassign convenience variables to correspond to member ends ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(1)) + u%Mesh%Position(:, mem%NodeIndx(1)) - pos1(3) = pos1(3) - p%MSL2SWL - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(2)) + u%Mesh%Position(:, mem%NodeIndx(2)) - pos2(3) = pos2(3) - p%MSL2SWL - z1 = pos1(3) - + pos1 = m%DispNodePosHst(:,mem%NodeIndx(1)) + pos2 = m%DispNodePosHst(:,mem%NodeIndx(2)) call GetOrientationAngles( pos1, pos2, phi1, sinPhi1, cosPhi1, tanPhi, sinBeta1, cosBeta1, k_hat1, errStat2, errMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) if ( N == 1 ) then ! Only one element in member @@ -3364,18 +3313,14 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, k_hat2 = k_hat1 else ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(N)) + u%Mesh%Position(:, mem%NodeIndx(N)) - pos1(3) = pos1(3) - p%MSL2SWL - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(N+1)) + u%Mesh%Position(:, mem%NodeIndx(N+1)) - pos2(3) = pos2(3) - p%MSL2SWL + pos1 = m%DispNodePosHst(:, mem%NodeIndx(N )) + pos2 = m%DispNodePosHst(:, mem%NodeIndx(N+1)) call GetOrientationAngles( pos1, pos2, phi2, sinPhi2, cosPhi2, tanPhi, sinBeta2, cosBeta2, k_hat2, errStat2, errMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) end if - - ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(N+1)) + u%Mesh%Position(:, mem%NodeIndx(N+1)) - pos2(3) = pos2(3) - p%MSL2SWL - z2 = pos2(3) + ! z-coordinates of the two ends of the member + z1 = m%DispNodePosHst(3,mem%NodeIndx( 1)) + z2 = m%DispNodePosHst(3,mem%NodeIndx(N+1)) !----------------------------------- filled buoyancy loads: starts -----------------------------------! !TODO: Do the equations below still work if z1 > z2 ? @@ -3420,13 +3365,10 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !---------------------------------- external buoyancy loads: starts ----------------------------------! if ( (.not. mem%PropPot) .AND. (mem%MHstLMod /= 0) ) then - - ! Get positions of member end nodes - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(1 )) + u%Mesh%Position(:, mem%NodeIndx(1 )) - pos1(3) = pos1(3) - p%MSL2SWL - r1 = mem%RMGB(1 ) - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(N+1)) + u%Mesh%Position(:, mem%NodeIndx(N+1)) - pos2(3) = pos2(3) - p%MSL2SWL + ! Get positions and scaled radii of member end nodes + pos1 = m%DispNodePosHst(:,mem%NodeIndx( 1)) + pos2 = m%DispNodePosHst(:,mem%NodeIndx(N+1)) + r1 = mem%RMGB( 1) r2 = mem%RMGB(N+1) if (mem%i_floor == 0) then ! both ends above or at seabed ! Compute loads on the end plate of node 1 @@ -3576,10 +3518,28 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !---------------------------------------------------------------------------------------------------------------! ! Map calculated results into the y%WriteOutput Array CALL MrsnOut_MapOutputs(y, p, u, m) - - CONTAINS + + SUBROUTINE GetDisplacedNodePosition( forceDisplaced, pos ) + LOGICAL, INTENT( IN ) :: forceDisplaced ! Set to true to return the exact displaced position no matter WaveDisp or WaveStMod + REAL(ReKi), INTENT( OUT ) :: pos(:,:) ! Displaced node positions + + ! Undisplaced node position + pos = u%Mesh%Position + pos(3,:) = pos(3,:) - p%MSL2SWL ! Z position measured from the SWL + IF ( (p%WaveDisp /= 0) .OR. forceDisplaced ) THEN + ! Use displaced X and Y position + pos(1,:) = pos(1,:) + u%Mesh%TranslationDisp(1,:) + pos(2,:) = pos(2,:) + u%Mesh%TranslationDisp(2,:) + IF ( (p%WaveStMod > 0) .OR. forceDisplaced ) THEN + ! Use displaced Z position only when wave stretching is enabled + pos(3,:) = pos(3,:) + u%Mesh%TranslationDisp(3,:) + END IF + END IF + + END SUBROUTINE GetDisplacedNodePosition + SUBROUTINE GetTotalWaveElev( Time, pos, Zeta, ErrStat, ErrMsg ) REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. @@ -3592,7 +3552,7 @@ SUBROUTINE GetTotalWaveElev( Time, pos, Zeta, ErrStat, ErrMsg ) ErrStat = ErrID_None ErrMsg = "" - Zeta = WaveField_GetTotalWaveElev( p%WaveField, Time, pos, ErrStat2, ErrMsg2 ) + Zeta = WaveField_GetNodeTotalWaveElev( p%WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END SUBROUTINE GetTotalWaveElev @@ -3610,7 +3570,7 @@ SUBROUTINE GetFreeSurfaceNormal( Time, pos, r, n, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = "" - CALL WaveField_GetWaveNormal( p%WaveField, Time, pos, r, n, ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveNormal( p%WaveField, Time, pos, r, n, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END SUBROUTINE GetFreeSurfaceNormal @@ -4229,7 +4189,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat END IF ! Get fluid velocity at the joint - CALL WaveField_GetWaveVel( p%WaveField, Time, pos, nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveVel( p%WaveField, Time, pos, .FALSE., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FV = REAL(FVTmp, ReKi) vrel = ( FV - u%Mesh%TranslationVel(:,J) ) * nodeInWater diff --git a/modules/hydrodyn/src/Morison.txt b/modules/hydrodyn/src/Morison.txt index a12a4a0cde..41be828b76 100644 --- a/modules/hydrodyn/src/Morison.txt +++ b/modules/hydrodyn/src/Morison.txt @@ -306,7 +306,9 @@ typedef ^ OtherStateType IntKi # ..... Misc/Optimization variables................................................................................................. # Define any data that are used only for efficiency purposes (these variables are not associated with time): # e.g. indices for searching in an array, large arrays that are local variables in any routine called multiple times, etc. -typedef ^ MiscVarType ReKi FV {:}{:} - - "Fluid velocity at line element node at time t, which may not correspond to the WaveTime array of times" - +typedef ^ MiscVarType ReKi DispNodePosHdn {:}{:} - - "Instantaneous displaced position of the line element nodes at time t for hydrodynamic load calculation" (m) +typedef ^ ^ ReKi DispNodePosHst {:}{:} - - "Instantaneous displaced position of the line element nodes at time t for hydrostatic and other load calcuation" (m) +typedef ^ ^ ReKi FV {:}{:} - - "Fluid velocity at line element node at time t, which may not correspond to the WaveTime array of times" - typedef ^ ^ ReKi FA {:}{:} - - "Fluid acceleration at line element node at time t, which may not correspond to the WaveTime array of times" - typedef ^ ^ ReKi FAMCF {:}{:} - - "Fluid acceleration at line element node at time t, which may not correspond to the WaveTime array of times" - typedef ^ ^ ReKi FDynP {:} - - "Fluid dynamic pressure at line element node at time t, which may not correspond to the WaveTime array of times" - diff --git a/modules/hydrodyn/src/Morison_Types.f90 b/modules/hydrodyn/src/Morison_Types.f90 index ac212b7386..165080b7af 100644 --- a/modules/hydrodyn/src/Morison_Types.f90 +++ b/modules/hydrodyn/src/Morison_Types.f90 @@ -369,6 +369,8 @@ MODULE Morison_Types ! ======================= ! ========= Morison_MiscVarType ======= TYPE, PUBLIC :: Morison_MiscVarType + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: DispNodePosHdn !< Instantaneous displaced position of the line element nodes at time t for hydrodynamic load calculation [(m)] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: DispNodePosHst !< Instantaneous displaced position of the line element nodes at time t for hydrostatic and other load calcuation [(m)] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: FV !< Fluid velocity at line element node at time t, which may not correspond to the WaveTime array of times [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: FA !< Fluid acceleration at line element node at time t, which may not correspond to the WaveTime array of times [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: FAMCF !< Fluid acceleration at line element node at time t, which may not correspond to the WaveTime array of times [-] @@ -9017,6 +9019,34 @@ SUBROUTINE Morison_CopyMisc( SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg ! ErrStat = ErrID_None ErrMsg = "" +IF (ALLOCATED(SrcMiscData%DispNodePosHdn)) THEN + i1_l = LBOUND(SrcMiscData%DispNodePosHdn,1) + i1_u = UBOUND(SrcMiscData%DispNodePosHdn,1) + i2_l = LBOUND(SrcMiscData%DispNodePosHdn,2) + i2_u = UBOUND(SrcMiscData%DispNodePosHdn,2) + IF (.NOT. ALLOCATED(DstMiscData%DispNodePosHdn)) THEN + ALLOCATE(DstMiscData%DispNodePosHdn(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%DispNodePosHdn.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + END IF + DstMiscData%DispNodePosHdn = SrcMiscData%DispNodePosHdn +ENDIF +IF (ALLOCATED(SrcMiscData%DispNodePosHst)) THEN + i1_l = LBOUND(SrcMiscData%DispNodePosHst,1) + i1_u = UBOUND(SrcMiscData%DispNodePosHst,1) + i2_l = LBOUND(SrcMiscData%DispNodePosHst,2) + i2_u = UBOUND(SrcMiscData%DispNodePosHst,2) + IF (.NOT. ALLOCATED(DstMiscData%DispNodePosHst)) THEN + ALLOCATE(DstMiscData%DispNodePosHst(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%DispNodePosHst.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + END IF + DstMiscData%DispNodePosHst = SrcMiscData%DispNodePosHst +ENDIF IF (ALLOCATED(SrcMiscData%FV)) THEN i1_l = LBOUND(SrcMiscData%FV,1) i1_u = UBOUND(SrcMiscData%FV,1) @@ -9273,6 +9303,12 @@ SUBROUTINE Morison_DestroyMisc( MiscData, ErrStat, ErrMsg ) ErrStat = ErrID_None ErrMsg = "" +IF (ALLOCATED(MiscData%DispNodePosHdn)) THEN + DEALLOCATE(MiscData%DispNodePosHdn) +ENDIF +IF (ALLOCATED(MiscData%DispNodePosHst)) THEN + DEALLOCATE(MiscData%DispNodePosHst) +ENDIF IF (ALLOCATED(MiscData%FV)) THEN DEALLOCATE(MiscData%FV) ENDIF @@ -9368,6 +9404,16 @@ SUBROUTINE Morison_PackMisc( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg Re_BufSz = 0 Db_BufSz = 0 Int_BufSz = 0 + Int_BufSz = Int_BufSz + 1 ! DispNodePosHdn allocated yes/no + IF ( ALLOCATED(InData%DispNodePosHdn) ) THEN + Int_BufSz = Int_BufSz + 2*2 ! DispNodePosHdn upper/lower bounds for each dimension + Re_BufSz = Re_BufSz + SIZE(InData%DispNodePosHdn) ! DispNodePosHdn + END IF + Int_BufSz = Int_BufSz + 1 ! DispNodePosHst allocated yes/no + IF ( ALLOCATED(InData%DispNodePosHst) ) THEN + Int_BufSz = Int_BufSz + 2*2 ! DispNodePosHst upper/lower bounds for each dimension + Re_BufSz = Re_BufSz + SIZE(InData%DispNodePosHst) ! DispNodePosHst + END IF Int_BufSz = Int_BufSz + 1 ! FV allocated yes/no IF ( ALLOCATED(InData%FV) ) THEN Int_BufSz = Int_BufSz + 2*2 ! FV upper/lower bounds for each dimension @@ -9505,6 +9551,46 @@ SUBROUTINE Morison_PackMisc( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg Db_Xferred = 1 Int_Xferred = 1 + IF ( .NOT. ALLOCATED(InData%DispNodePosHdn) ) 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%DispNodePosHdn,1) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%DispNodePosHdn,1) + Int_Xferred = Int_Xferred + 2 + IntKiBuf( Int_Xferred ) = LBOUND(InData%DispNodePosHdn,2) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%DispNodePosHdn,2) + Int_Xferred = Int_Xferred + 2 + + DO i2 = LBOUND(InData%DispNodePosHdn,2), UBOUND(InData%DispNodePosHdn,2) + DO i1 = LBOUND(InData%DispNodePosHdn,1), UBOUND(InData%DispNodePosHdn,1) + ReKiBuf(Re_Xferred) = InData%DispNodePosHdn(i1,i2) + Re_Xferred = Re_Xferred + 1 + END DO + END DO + END IF + IF ( .NOT. ALLOCATED(InData%DispNodePosHst) ) 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%DispNodePosHst,1) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%DispNodePosHst,1) + Int_Xferred = Int_Xferred + 2 + IntKiBuf( Int_Xferred ) = LBOUND(InData%DispNodePosHst,2) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%DispNodePosHst,2) + Int_Xferred = Int_Xferred + 2 + + DO i2 = LBOUND(InData%DispNodePosHst,2), UBOUND(InData%DispNodePosHst,2) + DO i1 = LBOUND(InData%DispNodePosHst,1), UBOUND(InData%DispNodePosHst,1) + ReKiBuf(Re_Xferred) = InData%DispNodePosHst(i1,i2) + Re_Xferred = Re_Xferred + 1 + END DO + END DO + END IF IF ( .NOT. ALLOCATED(InData%FV) ) THEN IntKiBuf( Int_Xferred ) = 0 Int_Xferred = Int_Xferred + 1 @@ -9883,6 +9969,52 @@ SUBROUTINE Morison_UnPackMisc( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, Err Re_Xferred = 1 Db_Xferred = 1 Int_Xferred = 1 + IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! DispNodePosHdn 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 + i2_l = IntKiBuf( Int_Xferred ) + i2_u = IntKiBuf( Int_Xferred + 1) + Int_Xferred = Int_Xferred + 2 + IF (ALLOCATED(OutData%DispNodePosHdn)) DEALLOCATE(OutData%DispNodePosHdn) + ALLOCATE(OutData%DispNodePosHdn(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%DispNodePosHdn.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + DO i2 = LBOUND(OutData%DispNodePosHdn,2), UBOUND(OutData%DispNodePosHdn,2) + DO i1 = LBOUND(OutData%DispNodePosHdn,1), UBOUND(OutData%DispNodePosHdn,1) + OutData%DispNodePosHdn(i1,i2) = ReKiBuf(Re_Xferred) + Re_Xferred = Re_Xferred + 1 + END DO + END DO + END IF + IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! DispNodePosHst 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 + i2_l = IntKiBuf( Int_Xferred ) + i2_u = IntKiBuf( Int_Xferred + 1) + Int_Xferred = Int_Xferred + 2 + IF (ALLOCATED(OutData%DispNodePosHst)) DEALLOCATE(OutData%DispNodePosHst) + ALLOCATE(OutData%DispNodePosHst(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%DispNodePosHst.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + DO i2 = LBOUND(OutData%DispNodePosHst,2), UBOUND(OutData%DispNodePosHst,2) + DO i1 = LBOUND(OutData%DispNodePosHst,1), UBOUND(OutData%DispNodePosHst,1) + OutData%DispNodePosHst(i1,i2) = ReKiBuf(Re_Xferred) + Re_Xferred = Re_Xferred + 1 + END DO + END DO + END IF IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! FV not allocated Int_Xferred = Int_Xferred + 1 ELSE diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index 4d6efea5ca..c33c053aaf 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -8,27 +8,29 @@ MODULE SeaSt_WaveField PRIVATE ! Public functions and subroutines -PUBLIC WaveField_GetWaveElev1 -PUBLIC WaveField_GetWaveElev2 -PUBLIC WaveField_GetTotalWaveElev -PUBLIC WaveField_GetWaveNormal +PUBLIC WaveField_GetNodeWaveElev1 +PUBLIC WaveField_GetNodeWaveElev2 +PUBLIC WaveField_GetNodeTotalWaveElev +PUBLIC WaveField_GetNodeWaveNormal +PUBLIC WaveField_GetNodeWaveKin +PUBLIC WaveField_GetNodeWaveVel + PUBLIC WaveField_GetWaveKin -PUBLIC WaveField_GetWaveVel CONTAINS !-------------------- Subroutine for wave elevation ------------------! -FUNCTION WaveField_GetWaveElev1( WaveField, Time, pos, ErrStat, ErrMsg ) +FUNCTION WaveField_GetNodeWaveElev1( WaveField, Time, pos, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None - REAL(SiKi) :: WaveField_GetWaveElev1 + REAL(SiKi) :: WaveField_GetNodeWaveElev1 REAL(SiKi) :: Zeta LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveElev1' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveElev1' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 @@ -42,21 +44,21 @@ FUNCTION WaveField_GetWaveElev1( WaveField, Time, pos, ErrStat, ErrMsg ) Zeta = 0.0_SiKi END IF - WaveField_GetWaveElev1 = Zeta + WaveField_GetNodeWaveElev1 = Zeta -END FUNCTION WaveField_GetWaveElev1 +END FUNCTION WaveField_GetNodeWaveElev1 -FUNCTION WaveField_GetWaveElev2( WaveField, Time, pos, ErrStat, ErrMsg ) +FUNCTION WaveField_GetNodeWaveElev2( WaveField, Time, pos, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None - REAL(SiKi) :: WaveField_GetWaveElev2 + REAL(SiKi) :: WaveField_GetNodeWaveElev2 REAL(SiKi) :: Zeta LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveElev2' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveElev2' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 @@ -70,36 +72,36 @@ FUNCTION WaveField_GetWaveElev2( WaveField, Time, pos, ErrStat, ErrMsg ) Zeta = 0.0_SiKi END IF - WaveField_GetWaveElev2 = Zeta + WaveField_GetNodeWaveElev2 = Zeta -END FUNCTION WaveField_GetWaveElev2 +END FUNCTION WaveField_GetNodeWaveElev2 -FUNCTION WaveField_GetTotalWaveElev( WaveField, Time, pos, ErrStat, ErrMsg ) +FUNCTION WaveField_GetNodeTotalWaveElev( WaveField, Time, pos, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None - REAL(SiKi) :: WaveField_GetTotalWaveElev + REAL(SiKi) :: WaveField_GetNodeTotalWaveElev REAL(SiKi) :: Zeta1, Zeta2 LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetTotalWaveElev' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeTotalWaveElev' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None ErrMsg = "" - Zeta1 = WaveField_GetWaveElev1( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + Zeta1 = WaveField_GetNodeWaveElev1( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - Zeta2 = WaveField_GetWaveElev2( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + Zeta2 = WaveField_GetNodeWaveElev2( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - WaveField_GetTotalWaveElev = Zeta1 + Zeta2 + WaveField_GetNodeTotalWaveElev = Zeta1 + Zeta2 -END FUNCTION WaveField_GetTotalWaveElev +END FUNCTION WaveField_GetNodeTotalWaveElev -SUBROUTINE WaveField_GetWaveNormal( WaveField, Time, pos, r, n, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, Time, pos, r, n, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time @@ -110,7 +112,7 @@ SUBROUTINE WaveField_GetWaveNormal( WaveField, Time, pos, r, n, ErrStat, ErrMsg CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None REAL(SiKi) :: ZetaP,ZetaM REAL(ReKi) :: r1,dZetadx,dZetady - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetFreeSurfaceNormal' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveNormal' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None @@ -118,28 +120,29 @@ SUBROUTINE WaveField_GetWaveNormal( WaveField, Time, pos, r, n, ErrStat, ErrMsg r1 = MAX(r,1.0e-6) ! In case r is zero - ZetaP = WaveField_GetTotalWaveElev( WaveField, Time, (/pos(1)+r1,pos(2)/), ErrStat2, ErrMsg2 ) + ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, Time, (/pos(1)+r1,pos(2)/), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - ZetaM = WaveField_GetTotalWaveElev( WaveField, Time, (/pos(1)-r1,pos(2)/), ErrStat2, ErrMsg2 ) + ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, Time, (/pos(1)-r1,pos(2)/), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) dZetadx = REAL(ZetaP-ZetaM,ReKi)/(2.0_ReKi*r1) - ZetaP = WaveField_GetTotalWaveElev( WaveField, Time, (/pos(1),pos(2)+r1/), ErrStat2, ErrMsg2 ) + ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, Time, (/pos(1),pos(2)+r1/), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - ZetaM = WaveField_GetTotalWaveElev( WaveField, Time, (/pos(1),pos(2)-r1/), ErrStat2, ErrMsg2 ) + ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, Time, (/pos(1),pos(2)-r1/), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) dZetady = REAL(ZetaP-ZetaM,ReKi)/(2.0_ReKi*r1) n = (/-dZetadx,-dZetady,1.0_ReKi/) n = n / SQRT(Dot_Product(n,n)) -END SUBROUTINE WaveField_GetWaveNormal +END SUBROUTINE WaveField_GetNodeWaveNormal !-------------------- Subroutine for full wave field kinematics --------------------! -SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetNodeWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(3) + LOGICAL, INTENT( IN ) :: forceNodeInWater REAL(SiKi), INTENT( OUT ) :: WaveElev1 REAL(SiKi), INTENT( OUT ) :: WaveElev2 REAL(SiKi), INTENT( OUT ) :: WaveElev @@ -155,7 +158,7 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, W REAL(ReKi) :: posXY(2), posPrime(3), posXY0(3) TYPE(SeaSt_Interp_MiscVarType) :: SeaSt_Interp_m LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveKin' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveKin' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 @@ -167,9 +170,9 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, W FAMCF(:) = 0.0 ! Wave elevation - WaveElev1 = WaveField_GetWaveElev1( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + WaveElev1 = WaveField_GetNodeWaveElev1( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - WaveElev2 = WaveField_GetWaveElev2( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + WaveElev2 = WaveField_GetNodeWaveElev2( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) WaveElev = WaveElev1 + WaveElev2 @@ -200,7 +203,7 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, W ELSE ! Wave stretching enabled - IF ( pos(3) <= WaveElev ) THEN ! Node is submerged + IF ( (pos(3) <= WaveElev) .OR. forceNodeInWater ) THEN ! Node is submerged nodeInWater = 1_IntKi @@ -259,6 +262,7 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, W ! 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 SeaSt_Interp_Setup( Time, posPrime, WaveField%seast_interp_p, seast_interp_m, ErrStat2, ErrMsg2 ) @@ -287,13 +291,14 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, W END IF ! If wave stretching is on or off -END SUBROUTINE WaveField_GetWaveKin +END SUBROUTINE WaveField_GetNodeWaveKin !-------------------- Subroutine for wave field velocity only --------------------! -SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, nodeInWater, FV, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetNodeWaveVel( WaveField, Time, pos, forceNodeInWater, nodeInWater, FV, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(3) + LOGICAL, INTENT( IN ) :: forceNodeInWater INTEGER(IntKi), INTENT( OUT ) :: nodeInWater REAL(SiKi), INTENT( OUT ) :: FV(3) INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation @@ -303,7 +308,7 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, nodeInWater, FV, ErrStat, REAL(ReKi) :: posXY(2), posPrime(3), posXY0(3) TYPE(SeaSt_Interp_MiscVarType) :: SeaSt_Interp_m LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveVel' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveVel' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 @@ -314,7 +319,7 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, nodeInWater, FV, ErrStat, posXY0 = (/pos(1),pos(2),0.0_ReKi/) ! Wave elevation - WaveElev = WaveField_GetTotalWaveElev( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + WaveElev = WaveField_GetNodeTotalWaveElev( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) IF (WaveField%WaveStMod == 0) THEN ! No wave stretching @@ -333,7 +338,7 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, nodeInWater, FV, ErrStat, ELSE ! Wave stretching enabled - IF ( pos(3) <= WaveElev ) THEN ! Node is submerged + IF ( (pos(3) <= WaveElev) .OR. forceNodeInWater ) THEN ! Node is submerged nodeInWater = 1_IntKi @@ -368,6 +373,7 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, nodeInWater, FV, ErrStat, ! 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 SeaSt_Interp_Setup( Time, posPrime, WaveField%seast_interp_p, seast_interp_m, ErrStat2, ErrMsg2 ) @@ -386,6 +392,46 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, nodeInWater, FV, ErrStat, END IF ! If wave stretching is on or off -END SUBROUTINE WaveField_GetWaveVel +END SUBROUTINE WaveField_GetNodeWaveVel + +SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) + TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField + REAL(DbKi), INTENT( IN ) :: Time + REAL(ReKi), INTENT( IN ) :: pos(:,:) + LOGICAL, INTENT( IN ) :: forceNodeInWater + REAL(SiKi), INTENT( OUT ) :: WaveElev1(:) + REAL(SiKi), INTENT( OUT ) :: WaveElev2(:) + REAL(SiKi), INTENT( OUT ) :: WaveElev(:) + REAL(ReKi), INTENT( OUT ) :: FV(:,:) + REAL(ReKi), INTENT( OUT ) :: FA(:,:) + REAL(ReKi), INTENT( OUT ) :: FAMCF(:,:) + REAL(ReKi), 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 + + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveKin' + INTEGER(IntKi) :: errStat2 + CHARACTER(ErrMsgLen) :: errMsg2 + + INTEGER(IntKi) :: NumPoints, i + REAL(SiKi) :: FDynP_node, FV_node(3), FA_node(3), FAMCF_node(3) + + ErrStat = ErrID_None + ErrMsg = "" + + NumPoints = size(pos, dim=2) + DO i = 1, NumPoints + CALL WaveField_GetNodeWaveKin( WaveField, Time, pos(:,i), forceNodeInWater, nodeInWater(i), WaveElev1(i), WaveElev2(i), WaveElev(i), FDynP_node, FV_node, FA_node, FAMCF_node, ErrStat2, ErrMsg2 ) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + FDynP(i) = REAL(FDynP_node,ReKi) + FV(:, i) = REAL(FV_node, ReKi) + FA(:, i) = REAL(FA_node, ReKi) + IF (ALLOCATED(WaveField%WaveAccMCF)) THEN + FAMCF(:,i) = REAL(FAMCF_node,ReKi) + END IF + END DO + +END SUBROUTINE WaveField_GetWaveKin END MODULE SeaSt_WaveField diff --git a/modules/seastate/src/SeaState.f90 b/modules/seastate/src/SeaState.f90 index 91dd0e054c..0855e36a45 100644 --- a/modules/seastate/src/SeaState.f90 +++ b/modules/seastate/src/SeaState.f90 @@ -857,7 +857,7 @@ SUBROUTINE SeaSt_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, Er DO i = 1, p%NWaveKin positionXYZ = (/p%WaveKinxi(i),p%WaveKinyi(i),p%WaveKinzi(i)/) - CALL WaveField_GetWaveKin( p%WaveField, Time, positionXYZ, nodeInWater, zeta1, zeta2, zeta, WaveDynP(i), WaveVel(:,i), WaveAcc(:,i), WaveAccMCF(:,i), ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveKin( p%WaveField, Time, positionXYZ, .FALSE., nodeInWater, zeta1, zeta2, zeta, WaveDynP(i), WaveVel(:,i), WaveAcc(:,i), WaveAccMCF(:,i), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END DO @@ -865,9 +865,9 @@ SUBROUTINE SeaSt_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, Er DO i = 1, p%NWaveElev positionXY = (/p%WaveElevxi(i),p%WaveElevyi(i)/) - WaveElev1(i) = WaveField_GetWaveElev1( p%WaveField, Time, positionXY, ErrStat2, ErrMsg2 ) + WaveElev1(i) = WaveField_GetNodeWaveElev1( p%WaveField, Time, positionXY, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - WaveElev2(i) = WaveField_GetWaveElev2( p%WaveField, Time, positionXY, ErrStat2, ErrMsg2 ) + WaveElev2(i) = WaveField_GetNodeWaveElev2( p%WaveField, Time, positionXY, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) WaveElev(i) = WaveElev1(i) + WaveElev2(i) END DO