anchor.diff:         
         ! External Hydrodynamic Side Loads
         DO i =1,N+1    ! loop through member nodes
  -         ! TODO: Note that for computational efficiency, we could precompute h_c and deltal for each element when we are NOT using wave stretching
  -         ! We would still need to test at time marching for nodes just below the free surface because that uses the current locations not the reference locations
  -         ! see table in Section 7.1.1
  -         if ( i == 1 ) then
  -            deltal = mem%dl/2.0_ReKi
  -            h_c    = mem%dl/4.0_ReKi
  -         elseif (i == N+1) then
  -            deltal =  mem%dl/2.0_ReKi
  -            h_c    = -mem%dl/4.0_ReKi
  -         elseif ( mem%i_floor == i ) then ! This node is the upper node of an element which crosses the seabed
  -            deltal = mem%dl/2.0_ReKi - mem%h_floor  ! TODO: h_floor is negative valued, should we be subrtracting it from dl/2? GJH
  -            h_c    = 0.5_ReKi*(mem%dl/2.0_ReKi + mem%h_floor)
  -         else
  -            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))
  -            if (pos1(3) <= 0.0 .and. 0.0 < pos2(3) ) then ! This node is just below the free surface !TODO: Needs to be augmented for wave stretching
  -               !TODO: Fix this one
  -               pos1 =  u%Mesh%Position(:, mem%NodeIndx(i)) ! use reference position for following equation
  -               h = (  pos1(3) ) / mem%cosPhi_ref !TODO: Needs to be augmented for wave stretching
  -               deltal = mem%dl/2.0 + h
  -               h_c    = 0.5*(h-mem%dl/2.0)
  +         z1 = u%Mesh%TranslationDisp(3, mem%NodeIndx(i))   + u%Mesh%Position(3, mem%NodeIndx(i))
  +         if ( i > mem%i_floor .and. z1 <= 0.0 ) then  ! node is above (or at? TODO: check) seabed and below or at free-surface)
  +            ! TODO: Note that for computational efficiency, we could precompute h_c and deltal for each element when we are NOT using wave stretching
  +            ! We would still need to test at time marching for nodes just below the free surface because that uses the current locations not the reference locations
  +            ! see table in Section 7.1.1
  +            if ( i == 1 ) then
  +               deltal = mem%dl/2.0_ReKi
  +               h_c    = mem%dl/4.0_ReKi
  +            elseif (i == N+1) then
  +               deltal =  mem%dl/2.0_ReKi
  +               h_c    = -mem%dl/4.0_ReKi
  +            elseif ( mem%i_floor == i+1 ) then ! This node is the upper node of an element which crosses the seabed
  +               deltal = mem%dl/2.0_ReKi - mem%h_floor  ! TODO: h_floor is negative valued, should we be subrtracting it from dl/2? GJH
  +               h_c    = 0.5_ReKi*(mem%dl/2.0_ReKi + mem%h_floor)
               else
  -               ! This node is a fully submerged interior node
  -               deltal = mem%dl
  -               h_c    = 0.0_ReKi
  +               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))
  +               if (pos1(3) <= 0.0 .and. 0.0 < pos2(3) ) then ! This node is just below the free surface !TODO: Needs to be augmented for wave stretching
  +                  !TODO: Fix this one
  +                  pos1 =  u%Mesh%Position(:, mem%NodeIndx(i)) ! use reference position for following equation
  +                  h = (  pos1(3) ) / mem%cosPhi_ref !TODO: Needs to be augmented for wave stretching
  +                  deltal = mem%dl/2.0 + h
  +                  h_c    = 0.5*(h-mem%dl/2.0)
  +               else
  +                  ! This node is a fully submerged interior node
  +                  deltal = mem%dl
  +                  h_c    = 0.0_ReKi
  +               end if
  +            
               end if
  -            
  -         end if
  -
  -         if (i == 1) then
  -            dRdl_p  = abs(mem%dRdl_mg(i))
  -            dRdl_pp = mem%dRdl_mg(i)   
  -         elseif ( i > 1 .and. i < (N+1)) then
  -            dRdl_p  = 0.5*( abs(mem%dRdl_mg(i-1)) + abs(mem%dRdl_mg(i)) )
  -            dRdl_pp = 0.5*( mem%dRdl_mg(i-1) + mem%dRdl_mg(i) )
  -         else 
  -            dRdl_p  = abs(mem%dRdl_mg(N))
  -            dRdl_pp = mem%dRdl_mg(N)
  -         end if
  -         
  -         ! ------------------- hydrodynamic drag loads: sides: Section 7.1.2 ------------------------ 
  -         vec = matmul( mem%Ak,m%vrel(:,mem%NodeIndx(i)) )
  -         f_hydro = mem%Cd(i)*p%WtrDens*mem%RMG(i)*TwoNorm(vec)*vec  +  &
  -                   0.5*mem%AxCd(i)*p%WtrDens*pi*mem%RMG(i)*dRdl_p * matmul( dot_product( mem%k, m%vrel(:,mem%NodeIndx(i)) )*mem%kkt, m%vrel(:,mem%NodeIndx(i)) )
  -         call LumpDistrHydroLoads( f_hydro, mem%k, mem%dl, h_c, m%F_D(:, mem%NodeIndx(i)) )
  -         
  -         if ( .not. mem%PropPot ) then
  -            ! ------------------- hydrodynamic added mass loads: sides: Section 7.1.3 ------------------------
  -            Am = mem%Ca(i)*p%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*mem%Ak + 2.0*mem%AxCa(i)*p%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p*mem%kkt
  -            f_hydro = -matmul( Am, u%Mesh%TranslationAcc(:,mem%NodeIndx(i)) )
  -            call LumpDistrHydroLoads( f_hydro, mem%k, mem%dl, h_c, m%F_A(:, mem%NodeIndx(i)) )
  -         
  -            ! ------------------- hydrodynamic inertia loads: sides: Section 7.1.4 ------------------------
  -            f_hydro=mem%Ca(i)*mem%Cp(i)*p%WtrDens*pi*mem%RMG(i)*mem%RMG(i)        * matmul( mem%Ak,  m%FA(:,mem%NodeIndx(i)) ) + &
  -                         2.0*mem%AxCa(i)*p%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p * matmul( mem%kkt, m%FA(:,mem%NodeIndx(i)) ) + &
  -                         2.0*m%FDynP(mem%NodeIndx(i))*mem%AxCp(i)*pi*mem%RMG(i)*dRdl_pp*mem%k 
  -            call LumpDistrHydroLoads( f_hydro, mem%k, mem%dl, h_c, m%F_I(:, mem%NodeIndx(i)) )
  -         end if
  +
  +            if (i == 1) then
  +               dRdl_p  = abs(mem%dRdl_mg(i))
  +               dRdl_pp = mem%dRdl_mg(i)   
  +            elseif ( i > 1 .and. i < (N+1)) then
  +               dRdl_p  = 0.5*( abs(mem%dRdl_mg(i-1)) + abs(mem%dRdl_mg(i)) )
  +               dRdl_pp = 0.5*( mem%dRdl_mg(i-1) + mem%dRdl_mg(i) )
  +            else 
  +               dRdl_p  = abs(mem%dRdl_mg(N))
  +               dRdl_pp = mem%dRdl_mg(N)
  +            end if
  +         
  +            ! ------------------- hydrodynamic drag loads: sides: Section 7.1.2 ------------------------ 
  +            vec = matmul( mem%Ak,m%vrel(:,mem%NodeIndx(i)) )
  +            f_hydro = mem%Cd(i)*p%WtrDens*mem%RMG(i)*TwoNorm(vec)*vec  +  &
  +                      0.5*mem%AxCd(i)*p%WtrDens*pi*mem%RMG(i)*dRdl_p * matmul( dot_product( mem%k, m%vrel(:,mem%NodeIndx(i)) )*mem%kkt, m%vrel(:,mem%NodeIndx(i)) )
  +            call LumpDistrHydroLoads( f_hydro, mem%k, mem%dl, h_c, m%F_D(:, mem%NodeIndx(i)) )
  +         
  +            if ( .not. mem%PropPot ) then
  +               ! ------------------- hydrodynamic added mass loads: sides: Section 7.1.3 ------------------------
  +               Am = mem%Ca(i)*p%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*mem%Ak + 2.0*mem%AxCa(i)*p%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p*mem%kkt
  +               f_hydro = -matmul( Am, u%Mesh%TranslationAcc(:,mem%NodeIndx(i)) )
  +               call LumpDistrHydroLoads( f_hydro, mem%k, mem%dl, h_c, m%F_A(:, mem%NodeIndx(i)) )
  +         
  +               ! ------------------- hydrodynamic inertia loads: sides: Section 7.1.4 ------------------------
  +               f_hydro=mem%Ca(i)*mem%Cp(i)*p%WtrDens*pi*mem%RMG(i)*mem%RMG(i)        * matmul( mem%Ak,  m%FA(:,mem%NodeIndx(i)) ) + &
  +                            2.0*mem%AxCa(i)*p%WtrDens*pi*mem%RMG(i)*mem%RMG(i)*dRdl_p * matmul( mem%kkt, m%FA(:,mem%NodeIndx(i)) ) + &
  +                            2.0*m%FDynP(mem%NodeIndx(i))*mem%AxCp(i)*pi*mem%RMG(i)*dRdl_pp*mem%k 
  +               call LumpDistrHydroLoads( f_hydro, mem%k, mem%dl, h_c, m%F_I(:, mem%NodeIndx(i)) )
  +            end if
  +         end if ! ( i > mem%i_floor .and. Zi <= 0.0 )
  +         
         END DO ! i =1,N+1    ! loop through member nodes       
         
         
  
positive.code:             exit
           end if
        end do
     end if
  
  
  
negative.code: 
     ! calculate h_floor if seabed-piercing
     member%h_floor = 0.0_ReKi
     member%i_floor = member%NElements+1  ! Default to entire member is below the seabed
     if (Za < -WtrDepth) then
        do i= 2, member%NElements+1
           Za = InitInp%Nodes(member%NodeIndx(i))%Position(3)
  