labelize_point Subroutine

   subroutine labelize_point(height, label, x, y)
   !================================================================================================
   !! Function to label a point as: peak, valley, saddle or nothing particular
   implicit none
   real   (kind=R8), intent(in ), dimension(1:3, 1:3) :: height      !! *nodal height values as a 2D array*
   character(len=1), intent(out)                      :: label       !! *kind of point*
   real   (kind=R8), intent(out), optional            :: x, y        !! *coordinates of the extremum found*

      real(kind=R8), dimension(1:9)      :: tab
      real(kind=R8), dimension(1:9, 1:6) :: derivatives

      real(kind=R8) :: h, dhdx, dhdy, d2hdx2, d2hdy2, d2hdxdy, delta, dx, dy
      real(kind=R8) :: xs, ys

      real(kind=R8), parameter :: eps = 1.0e-8_R8

      integer(kind=I4) :: k, status

      ! following QU9 node notation:
      tab(1:9) = [ height(1, 1), &  !  (3,1)----(3,2)----(3,3)
                   height(1, 3), &  !    |        |        |
                   height(3, 3), &  !    |        |        |
                   height(3, 1), &  !  (2,1)----(2,2)----(2,3)
                   height(1, 2), &  !    |        |        |
                   height(2, 3), &  !    |        |        |
                   height(3, 2), &  !  (1,1)----(1,2)----(1,3)
                   height(2, 1), &  !
                   height(2, 2) ]   !

      ! Newton Raphson to locate the null first derivatives, which means an extremum
      ! Intiate with the center of the square [-1,1]X[-1,1] and initalize the counting
      !
      xs = 0 ; ys = 0 ; k = 0
      do

         ! for the current coordinates, what are the surface height and its derivatives:
         call deriv_N(x = xs, y = ys, mat_d = derivatives(1:9, 1:6))

         dhdx    = sum( derivatives(1:9, 2) * tab(1:9) )
         dhdy    = sum( derivatives(1:9, 3) * tab(1:9) )
         d2hdx2  = sum( derivatives(1:9, 4) * tab(1:9) )
         d2hdxdy = sum( derivatives(1:9, 5) * tab(1:9) )
         d2hdy2  = sum( derivatives(1:9, 6) * tab(1:9) )

         ! jacobian denominator
         delta = d2hdx2 * d2hdy2 - d2hdxdy**2

         if ( abs(dhdx) < eps .and. abs(dhdy) < eps ) then             ! converge ok
            status = 0                                                 ! extremum found, whatever its kind
            exit                                                       ! nothing more to do
         endif

         if ( abs(xs) >= 5._R8 .or. abs(ys) >= 5._R8 ) then            ! during the convergence process, the point is far from the square, so exit
            status = 1                                                 ! with the appropriate status
            exit                                                       !
         endif                                                         !

         k = k + 1

         if (k > 1000) then   ! limit the number of iterations
            status = 1
            exit
         endif

         dx = (-1. / delta) * ( + d2hdy2  * dhdx - d2hdxdy * dhdy)
         dy = (-1. / delta) * ( - d2hdxdy * dhdx + d2hdx2  * dhdy)

         xs = xs + 0.9 * dx
         ys = ys + 0.9 * dy


      enddo

      ! outside the square [-0.5,0.5]X[-0.5,0.5] the extremum belongs to another node
      if ( abs(xs) > 0.5_R8 .or. abs(ys) > 0.5_R8 ) status = 1

      ! if derivatives are null, what kind of point is it ?
      !     Nodes order:
      !
      !     4---7---3
      !     |   |   |
      !     8---9-- 6
      !     |   |   |
      !     1---5---2

      if (status == 0) then

         call deriv_N(x = xs, y = ys, mat_d = derivatives(1:9, 1:6))

         ! height at the extremum point xs, ys
         h = sum( derivatives(1:9, 1) * tab(1:9) )

         if (     all(tab(1:8) >= h) ) then ; label = 'V' ! valley or pit
         elseif ( all(tab(1:8) <= h) ) then ; label = 'P' ! hill or peak
         else ;                               label = 'S' ! saddle
         endif

      else

         label = 'N' ! nothing particular

      endif

      if ( present(x) ) then
         x = xs
         y = ys
      endif

   return
   endsubroutine labelize_point