anchor.diff:  
          
         ! External Hydrodynamic Side Loads
  -      DO i =1,N+1    ! loop through member elements    
  -         h_c = 0.0_ReKi !TODO: Determine actual h_c based on water level, etc. GJH 3/23/20 see table in Section 7.1.1
  - !TODO:        call GetExtHydroForceMomentMultipliers( i, N+1, Za, Zb, Zi, h_c, deltal )
  +      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)
  +            else
  +               ! This node is a fully submerged interior node
  +               deltal = mem%dl
  +               h_c    = 0.0_ReKi
  +            end if
  +            
  +         end if
  +
            if (i == 1) then
               dRdl_p  = abs(mem%dRdl_mg(i))
               dRdl_pp = mem%dRdl_mg(i)   
  
positive.code:                          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 DO ! i =1,N+1    ! loop through member elements       
        
        
        ! Any end plate loads that are modeled on a per-member basis
  
negative.code:             Fl      = -mem%FillDens * g * pi *mem%Rin(1)**2*mem%l_fill*cosPhi
              Moment  =  mem%FillDens * g * pi *0.25*mem%Rin(1)**4*sinPhi
              call AddEndLoad(Fl, Moment, sinPhi1, cosPhi1, sinBeta1, cosBeta1, m%F_BF_End(:, mem%NodeIndx(1)))
           else
              ! no load if member is not flooded at all
           end if
           
        elseif ( mem%i_floor < mem%NElements ) then ! upper node is still above the seabed
           if (mem%z_overfill > 0) then 
              Fl      =   mem%FillDens * g * pi *mem%Rin(N+1)**2* (mem%z_overfill + max(z1-z2, 0.0_ReKi))
              Moment  =  -mem%FillDens * g * pi *0.25*mem%Rin(N+1)**4*sinPhi            
              call AddEndLoad(Fl, Moment, sinPhi2, cosPhi2, sinBeta2, cosBeta2, m%F_BF_End(:, mem%NodeIndx(N+1)))
           end if
           
        else    
           ! no loads because both end nodes are below seabed
        end if
  
        ! --- no inertia loads from water ballast modeled on ends
  
        ! --- external buoyancy loads: ends ---
  
        if ( .not. mem%PropPot ) then
           if (mem%i_floor == 0) then  ! both ends above or at seabed
              if (z2<= 0.0_ReKi) then
                 ! Compute loads on both ends
                 Fl      = -p%WtrDens * g * pi *mem%RMG(1)**2*z1
  