calcul_asfc_hermite Subroutine

public subroutine calcul_asfc_hermite(tab_in, scal, asfc_res, omp)

Return the asfc of a surface. The different grids are obtained by Hermite interpolation

Arguments

Type IntentOptional Attributes Name
real(kind=R8), intent(in), dimension(1:scal%xres, 1:scal%yres) :: tab_in

input surface

type(SCALE_SURF), intent(in) :: scal

surface characteristics

real(kind=R8), intent(out), dimension(1:2) :: asfc_res

result: asfc, adjustment factor

logical(kind=I4), intent(in) :: omp

with openmp ?


Calls

proc~~calcul_asfc_hermite~~CallsGraph proc~calcul_asfc_hermite calcul_asfc_hermite get_unit get_unit proc~calcul_asfc_hermite->get_unit lmder1 lmder1 proc~calcul_asfc_hermite->lmder1 omp_get_num_procs omp_get_num_procs proc~calcul_asfc_hermite->omp_get_num_procs proc~df_boltz df_boltz proc~calcul_asfc_hermite->proc~df_boltz proc~dnq_et_i dnq_et_i proc~calcul_asfc_hermite->proc~dnq_et_i proc~dnq_xi_i dnq_xi_i proc~calcul_asfc_hermite->proc~dnq_xi_i proc~f_boltz f_boltz proc~calcul_asfc_hermite->proc~f_boltz proc~init_beta_boltz init_beta_boltz proc~calcul_asfc_hermite->proc~init_beta_boltz proc~locate locate proc~calcul_asfc_hermite->proc~locate proc~locate2 locate2 proc~calcul_asfc_hermite->proc~locate2 proc~nq_i nq_i proc~calcul_asfc_hermite->proc~nq_i proc~dn_i dn_i proc~dnq_et_i->proc~dn_i proc~n_i n_i proc~dnq_et_i->proc~n_i proc~dnq_xi_i->proc~dn_i proc~dnq_xi_i->proc~n_i proc~nq_i->proc~n_i

Called by

proc~~calcul_asfc_hermite~~CalledByGraph proc~calcul_asfc_hermite calcul_asfc_hermite proc~calcul_asfc calcul_asfc proc~calcul_asfc->proc~calcul_asfc_hermite program~test_asfc test_asfc program~test_asfc->proc~calcul_asfc_hermite

Source Code

   subroutine calcul_asfc_hermite(tab_in, scal, asfc_res, omp)
   !================================================================================================
   !! Return the *asfc* of a surface. The different grids are obtained by Hermite interpolation
   implicit none
   type(SCALE_SURF), intent(in )                                      :: scal     !! *surface characteristics*
   real   (kind=R8), intent(in ), dimension(1:scal%xres, 1:scal%yres) :: tab_in   !! *input surface*
   real   (kind=R8), intent(out), dimension(1:2)                      :: asfc_res !! *result: asfc, adjustment factor*
   logical(kind=I4), intent(in )                                      :: omp      !! *with openmp ?*

      real(kind=R8), allocatable, dimension(:) :: x            ! x points in original grid
      real(kind=R8), allocatable, dimension(:) :: y            ! y points in original grid

      integer(kind=I4) :: long_new                             ! number of points in x dimension for new grid
      integer(kind=I4) :: larg_new                             ! number of points in y dimension for new grid

      real(kind=R8), allocatable, dimension(:)   :: x_new      ! new grid x points
      real(kind=R8), allocatable, dimension(:)   :: y_new      ! new grid y points
      real(kind=R8), allocatable, dimension(:,:) :: tab_ou     ! new grid function evaluations


      real(kind=R8), allocatable, dimension(:,:) :: tab_in_dx  ! new grid function evaluations
      real(kind=R8), allocatable, dimension(:,:) :: tab_in_dy  ! new grid function evaluations
      real(kind=R8), allocatable, dimension(:,:) :: tab_in_xy  ! new grid function evaluations

      real(kind=R8), allocatable, dimension(:) :: gx
      real(kind=R8), allocatable, dimension(:) :: gy
      real(kind=R8), allocatable, dimension(:) :: gw
      real(kind=R8), allocatable, dimension(:) :: tab_dnq

      real(kind=R8)    :: rr
      integer(kind=I4) :: i, ii, j, jj, k, long_tmp, larg_tmp

      real(kind=R8), dimension(1:nb_beta) :: beta
      real(kind=R8), dimension(1:npp)     :: vec_s

      real(kind=R8), dimension(1:npp) :: vec_x  ! points coordinates
      real(kind=R8), dimension(1:npp) :: vec_y  ! points coordinates

      real(kind=R8) :: asfc1, asfc2, aire, hx, hy, hhx, hhy, width, height


      real   (kind=R8) :: xi, yi, eps_x
      integer(kind=I4) :: it, ng, nbpt

      integer(kind=I4) :: long, larg

      integer(kind=I4) :: nb_th

      long = scal%xres
      larg = scal%yres

      if (out_her) call get_unit(unit_out_her)
      if (out_her) open(unit=unit_out_her, file="out/asfc_her_her_all.txt")

      width  = scal%lx
      height = scal%ly

      hx =  width/(long -1)
      hy = height/(larg -1)

      eps_x = min(hx/10**3, hy/10**3)

      ! définition d'abscisses pour l'interpolation par splines
      allocate( x(1:long), y(1:larg) )
      do i = 1, long
         x(i) = hx*real(i-1, kind=R8)
      enddo
      do j = 1, larg
         y(j) = hy*real(j-1, kind=R8)
      enddo


      allocate( x_new(1:long),               &  !
                y_new(1:larg),               &  !
               tab_ou(1:long, 1:larg),       &  !
               tab_in_dx(1:long, 1:larg),    &  !
               tab_in_dy(1:long, 1:larg),    &  !
               tab_in_xy(1:long, 1:larg) )      !

      rr =          (real(long0, kind=R8)/long)**(UN/npp)    ! facteur de réduction pour aller du maillage initial au maillage minimal avec npp points
      rr = max( rr, (real(larg0, kind=R8)/larg)**(UN/npp) )  ! facteur de réduction pour aller du maillage initial au maillage minimal avec npp points

      x_new(1:long)          = x(1:long)
      y_new(1:larg)          = y(1:larg)
      tab_ou(1:long, 1:larg) = tab_in(1:long, 1:larg)


      call init_aire_hermite(gx = gx, gy = gy, gw = gw, tab_dnq = tab_dnq, ng = ng)

      ! pour chaque réduction de maillage, calcul du maillage résultant et de l'aire relative associée
      !.............................................................
      long_new = long
      larg_new = larg
      hhx      = hx
      hhy      = hy


      call calcul_aire_hermite(   tab_in = tab_ou(1:long_new, 1:larg_new), &  !
                                    long = long_new,                       &  !
                                    larg = larg_new,                       &  !
                                      gw = gw(1:ng),                       &  !
                                 tab_dnq = tab_dnq(1:ng),                  &  !
                                      ng = ng,                             &  !
                                      hx = hhx,                            &  !
                                      hy = hhy,                            &  !
                                   width = width,                          &  !
                                  height = height,                         &  !
                                    aire = aire)                              !


      vec_s(1) = log( aire )
      vec_x(1) = log( (hhx*1e6)*(hhy*1e6)/2 )

      call calcul_tabd_hermite( tab_in = tab_in   (1:long, 1:larg),  &  !
                                tab_dx = tab_in_dx(1:long, 1:larg),  &  !
                                tab_dy = tab_in_dy(1:long, 1:larg),  &  !
                                tab_xy = tab_in_xy(1:long, 1:larg),  &  !
                                  long = long,                       &  !
                                  larg = larg,                       &  !
                                    hx = hx,                         &  !
                                    hy = hy )                           !

      nb_th = 1
      if (omp) then
         nb_th = omp_get_num_procs()
      endif


      !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nb_th) IF (omp)
      !$OMP DO ORDERED SCHEDULE (STATIC,1) PRIVATE(it, long_tmp, larg_tmp, long_new, larg_new, hhx, hhy, i, x_new, j, y_new, yi, jj, xi, ii, tab_ou, aire)


      do it = 1, npp - 1

         long_tmp = nint(long*(rr**it))   ! nb points en suite géométrique
         larg_tmp = nint(larg*(rr**it))

         if ( long_new == long_tmp .or. larg_new == larg_tmp) then
            vec_s(it + 1) = 0
            cycle ! à découper trop fin, on peut tomber sur les mêmes entiers
         endif

         long_new = long_tmp
         larg_new = larg_tmp

         hhx =  width/(long_new -1)
         hhy = height/(larg_new -1)

         vec_x(it + 1) = log( (hhx*1e6)*(hhy*1e6)/2 )


         deallocate( x_new, y_new, tab_ou )

         allocate( x_new(1:long_new),                 &  !
                   y_new(1:larg_new),                 &  ! nouvelles abscisses
                  tab_ou(1:long_new, 1:larg_new) )       !

         do i = 1, long_new
            x_new(i) = hhx*real(i-1, kind=R8)
         enddo

         do j = 1, larg_new
            y_new(j) = hhy*real(j-1, kind=R8)
         enddo

         do j = 1, larg_new

            yi = y_new(j)
            jj = locate2(n = larg, xx = y(1:larg), x = yi, eps = eps_x)
            yi = (yi -(y(jj)+y(jj+1))/2)/(hy/2)

            do i = 1, long_new

               xi = x_new(i)
               ii = locate2(n = long, xx = x(1:long), x = xi, eps = eps_x)
               xi = (xi -(x(ii)+x(ii+1))/2)/(hx/2)

               tab_ou(i, j) =                                              &  !
                     nq_i(xi, yi, 1, 1)*tab_in(ii   , jj   ) +             &  !  u1
                     nq_i(xi, yi, 2, 1)*tab_in(ii +1, jj   ) +             &  !  u2
                     nq_i(xi, yi, 3, 1)*tab_in(ii +1, jj +1) +             &  !  u3
                     nq_i(xi, yi, 4, 1)*tab_in(ii   , jj +1) +             &  !  u4

                     nq_i(xi, yi, 1, 2)*tab_in_dx(ii   , jj   )*hx/2 +     &  ! du1/dx
                     nq_i(xi, yi, 2, 2)*tab_in_dx(ii +1, jj   )*hx/2 +     &  ! du2/dx
                     nq_i(xi, yi, 3, 2)*tab_in_dx(ii +1, jj +1)*hx/2 +     &  ! du3/dx
                     nq_i(xi, yi, 4, 2)*tab_in_dx(ii   , jj +1)*hx/2 +     &  ! du4/dx

                     nq_i(xi, yi, 1, 3)*tab_in_dy(ii   , jj   )*hy/2 +     &  ! du1/dy
                     nq_i(xi, yi, 2, 3)*tab_in_dy(ii +1, jj   )*hy/2 +     &  ! du2/dy
                     nq_i(xi, yi, 3, 3)*tab_in_dy(ii +1, jj +1)*hy/2 +     &  ! du3/dy
                     nq_i(xi, yi, 4, 3)*tab_in_dy(ii   , jj +1)*hy/2 +     &  ! du4/dy

                     nq_i(xi, yi, 1, 4)*tab_in_xy(ii   , jj   )*hx*hy/4 +  &  ! du1/dxdy
                     nq_i(xi, yi, 2, 4)*tab_in_xy(ii +1, jj   )*hx*hy/4 +  &  ! du2/dxdy
                     nq_i(xi, yi, 3, 4)*tab_in_xy(ii +1, jj +1)*hx*hy/4 +  &  ! du3/dxdy
                     nq_i(xi, yi, 4, 4)*tab_in_xy(ii   , jj +1)*hx*hy/4       ! du4/dxdy

            enddo
         enddo

         call calcul_aire_hermite(   tab_in = tab_ou(1:long_new, 1:larg_new), &  !
                                       long = long_new,                       &  !
                                       larg = larg_new,                       &  !
                                         gw = gw(1:ng),                       &  !
                                    tab_dnq = tab_dnq(1:ng),                  &  !
                                         ng = ng,                             &  !
                                         hx = hhx,                            &  !
                                         hy = hhy,                            &  !
                                      width = width,                          &  !
                                     height = height,                         &  !
                                       aire = aire)                              !

         vec_s(it + 1) = log( aire )

!~          write(unit_out_her, *) vec_x(it+1), vec_s(it+1)

      enddo

      !$OMP END DO
      !$OMP END PARALLEL

      deallocate( x_new, y_new, x, y, tab_ou, tab_in_dx, tab_in_dy, tab_in_xy, tab_dnq, gx, gy, gw )

      if (out_lin) close(unit_out_lin)
      if (out_her) close(unit_out_her)
      !.............................................................

      k = 1
      do i = 1, npp
         if ( abs(vec_s(i)) > EPS_R8 ) then
            vec_y(k) = vec_s(i)
            vec_x(k) = vec_x(i)
            k = k + 1
         endif
      enddo
      nbpt = k - 1

      call interpolate_asfc( bt = beta(1:nb_beta),    &  !
                           n_bt = nb_beta,            &  !
                           n_pt = nbpt,               &  !
                            asf = asfc1,              &  !
                             r2 = asfc2 )                !

      asfc_res = [-1000*asfc1, asfc2]

   return

   contains

      subroutine interpolate_asfc(bt, n_bt, n_pt, asf, r2)
      !================================================================================================
      !< @note Function that fits the data points for the asfc determination.<br/>
      !        \(f\_boltz(x_i)=\beta_2 + \beta_3 \tanh \left( \dfrac{x_i -\beta_1}{\beta_4} \right) \)
      !
      !  @endnote
      !------------------------------------------------------------------------------------------------
      implicit none
      integer(kind=I4), intent(in )                    :: n_bt !! *number of parameters*
      integer(kind=I4), intent(in )                    :: n_pt !! *data vector length*
      real   (kind=R8), intent(out), dimension(1:n_bt) :: bt   !! *vector \(\beta\) of parameters*
      real   (kind=R8), intent(out)                    :: asf  !! *Asfc number*
      real   (kind=R8), intent(out)                    :: r2   !! *correlation number to assess validity of the Asfc calculus*

         real   (kind=R8), dimension(1:n_pt, 1:n_bt) :: jacob
         real   (kind=R8), dimension(1:n_pt)         :: v_x, v_y, res_y, pentes
         integer(kind=I4)                            :: i0, i, j, info
         real   (kind=R8)                            :: delta1, delta2, y_mean

         v_x(1:n_pt) = vec_x(1:n_pt)

         ! smoothing
         v_y(1) = vec_y(1)
         do i = 1 + 1, n_pt - 1

            v_y(i) = 0.25_R8 * ( vec_y(i - 1) + 2 * vec_y(i) + vec_y(i + 1) )

         enddo
         v_y(n_pt) = vec_y(n_pt)

         call init_beta_boltz(   bt = bt(1:n_bt),     & !
                               n_bt = n_bt,           & !
                                v_x = v_x(1:n_pt),    & !
                                v_y = v_y(1:n_pt),    & !
                               n_pt = n_pt)             !

         res_y(1:n_pt)         = 0._R8
         jacob(1:n_pt, 1:n_bt) = 0._R8

         call lmder1(   fcn = lmder1_f,               & !
                          m = n_pt,                   & !
                          n = n_bt,                   & !
                          x = bt(1:n_bt),             & !
                       fvec = res_y(1:n_pt),          & !
                       fjac = jacob(1:n_pt, 1:n_bt),  & !
                     ldfjac = n_pt,                   & !
                        tol = 1.0e-8_R8,              & !
                       info = info)                     !

         delta1 = 0._R8
         delta2 = 0._R8
         y_mean = sum( vec_y(1:n_pt) )/n_pt
         do i = 1, n_pt
            delta1 = delta1 +( vec_y(i) -f_boltz(xi = vec_x(i),   & !
                                               beta = bt(1:n_bt), & !
                                             n_beta = n_bt)       & !
                             )**2
            delta2 = delta2 +( vec_y(i) -y_mean )**2
         enddo
         r2 = 1._R8 -delta1/delta2

   !~       if (r2<0) then
   !~       do i=1,n_pt
   !~       write(99,*) v_x(i), v_y(i)
   !~       enddo
   !~       stop 'error'
   !~       endif

         i0 = locate(n = n_pt,        & !
                    xx = v_x(1:n_pt), & !
                     x = bt(1))         !

         j = 0
         do i = i0 -5, i0 +5
            if (i<1 .or. i>n_pt) cycle
            j = j +1
            pentes(j) = +df_boltz(  xi = v_x(i),     & !
                                  beta = bt(1:n_bt), & !
                                n_beta = n_bt,       & !
                                  ivar = 0)            !
         enddo

         asf = sum( pentes(1:j)/j )

      return
      endsubroutine interpolate_asfc

      subroutine lmder1_f(m, n, x, fvec, fjac, ldfjac, iflag)
      !================================================================================================
      !< @note Function called by [[lmder1]] as part of **minpack**, modified by
      !        [Burkardt](https://people.sc.fsu.edu/~jburkardt/f_src/minpack/minpack.html) <br/>
      !        According *iflag* value it calculates the function [[f_boltz]] at the data points
      !        or the jacobian.
      !
      !  @endnote
      !------------------------------------------------------------------------------------------------
      implicit none
      integer(kind=I4), intent(in   )                           :: m       !! *number of points*
      integer(kind=I4), intent(in   )                           :: n       !! *number of parameters*
      integer(kind=I4), intent(in   )                           :: ldfjac  !! *leading dimension of fjac, which must be at least n*
      integer(kind=I4), intent(in   )                           :: iflag   !! *which calculus to perform*
      real   (kind=R8), intent(out  ), dimension(1:m)           :: fvec    !! *vector of f_boltz(xi)*
      real   (kind=R8), intent(out  ), dimension(1:ldfjac, 1:n) :: fjac    !! *jacobian*
      real   (kind=R8), intent(inout), dimension(1:n)           :: x       !! *parameter values*

         integer(kind=I4) :: i, k

         select case (iflag)
            case(0)
               continue

            case(1)
               do i = 1, m
                  fvec(i) = f_boltz(     xi = vec_x(i),  &  !
                                       beta = x(1:n),    &  !
                                     n_beta = n)         &  !
                            -vec_y(i)
               enddo

            case(2)
               do i = 1, m
               do k = 1, n
                  fjac(i, k) = df_boltz(xi = vec_x(i), & !
                                      beta = x(1:n),   & !
                                      ivar = k,        & !
                                    n_beta = n)          !
               enddo
               enddo

            case default
               write ( *, '(a)' ) ' '
               write ( *, '(a)' ) 'LMDER1_F - Fatal error!'
               write ( *, '(a,i6)' ) '  Called with unexpected value of IFLAG = ', iflag
               stop

         endselect
      return
      endsubroutine lmder1_f

      subroutine init_aire_hermite(gx, gy, gw, tab_dnq, ng)
      !================================================================================================
      !< @note
      !
      !  @endnote
      !------------------------------------------------------------------------------------------------
      implicit none
      real   (kind=R8), intent(out), allocatable, dimension(:) :: gx
      real   (kind=R8), intent(out), allocatable, dimension(:) :: gy
      real   (kind=R8), intent(out), allocatable, dimension(:) :: gw
      real   (kind=R8), intent(out), allocatable, dimension(:) :: tab_dnq
      integer(kind=I4), intent(out)                            :: ng

         real   (kind=R8) :: x1, x2, y1, y2
         integer(kind=I4) :: i, k, nb_gauss_1d

         nb_gauss_1d = 2
         ng          = nb_gauss_1d**2

         allocate( gx(1:ng), gy(1:ng), gw(1:ng) )

         select case(nb_gauss_1d)
            case(1)
               x1 = 0._R8
               y1 = 2._R8

               gx(1:ng) = [   x1]
               gy(1:ng) = [   x1]
               gw(1:ng) = [y1*y1]
            case(2)
               x1 = sqrt(1._R8/3._R8)
               y1 = UN

               gx(1:ng) = [   -x1,    -x1,    +x1,    +x1]
               gy(1:ng) = [   -x1,    +x1,    -x1,    +x1]
               gw(1:ng) = [ y1*y1,  y1*y1,  y1*y1,  y1*y1]
            case(3)
               x1 = sqrt(3._R8/5.0_R8)
               x2 = 0._R8
               y1 = 5._R8/9._R8
               y2 = 8._R8/9._R8

               gx(1:ng) = [   -x1,   -x1,   -x1,    x2,    x2,    x2,   +x1,    +x1,   +x1]
               gy(1:ng) = [   -x1,    x2,   +x1,   -x1,    x2,   +x1,   -x1,     x2,   +x1]
               gw(1:ng) = [ y1*y1, y1*y2, y1*y1, y1*y2, y2*y2, y1*y2, y1*y1,  y1*y2, y1*y1]
         endselect

         allocate( tab_dnq(1:32*ng) )
         i = 1
         do k = 1, ng

            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 1, 1) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 2, 1) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 3, 1) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 4, 1) ; i = i +1

            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 1, 2) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 2, 2) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 3, 2) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 4, 2) ; i = i +1

            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 1, 3) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 2, 3) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 3, 3) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 4, 3) ; i = i +1

            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 1, 4) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 2, 4) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 3, 4) ; i = i +1
            tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 4, 4) ; i = i +1
            !!
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 1, 1) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 2, 1) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 3, 1) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 4, 1) ; i = i +1

            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 1, 2) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 2, 2) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 3, 2) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 4, 2) ; i = i +1

            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 1, 3) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 2, 3) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 3, 3) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 4, 3) ; i = i +1

            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 1, 4) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 2, 4) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 3, 4) ; i = i +1
            tab_dnq(i) = dnq_et_i(gx(k), gy(k), 4, 4) ; i = i +1

         enddo

      return
      endsubroutine init_aire_hermite

      subroutine calcul_tabd_hermite( tab_in, tab_dx, tab_dy, tab_xy, long, larg, hx, hy)
      !================================================================================================
      !< @note
      !
      !  @endnote
      !------------------------------------------------------------------------------------------------
      implicit none
      integer(kind=I4), intent(in )                            :: long
      integer(kind=I4), intent(in )                            :: larg
      real   (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab_in
      real   (kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_dx
      real   (kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_dy
      real   (kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_xy
      real   (kind=R8), intent(in )                            :: hx
      real   (kind=R8), intent(in )                            :: hy

         integer(kind=I4) :: i, im, ip, j, jm, jp
         real   (kind=R8) :: ui, uim, uip, ujm, ujp, upp, ump, upm, umm

         do j = 1, larg

            jm = max(j -1,    1)
            jp = min(j +1, larg)

            do i = 1, long
               im = max(i -1,    1)
               ip = min(i +1, long)

               ui  = tab_in(i , j )
               uim = tab_in(im, j )
               uip = tab_in(ip, j )
               ujm = tab_in(i , jm)
               ujp = tab_in(i , jp)
               upp = tab_in(ip, jp)
               ump = tab_in(im, jp)
               upm = tab_in(ip, jm)
               umm = tab_in(im, jm)

               tab_dx(i, j) = (uip -uim)/(2*hx)
               tab_dy(i, j) = (ujp -ujm)/(2*hy)
               tab_xy(i, j) = (upp -ump -upm +umm)/(4*hx*hy)
            enddo
         enddo

         tab_dx(   1, 1:larg) = ( tab_in(      2, 1:larg) -      &  !
                                  tab_in(      1, 1:larg) )/hx      !
         tab_dx(long, 1:larg) = ( tab_in(long   , 1:larg) -      &  !
                                  tab_in(long -1, 1:larg) )/hx      !

         tab_dy(1:long,    1) = ( tab_in(1:long,       2) -      &  !
                                  tab_in(1:long,       1) )/hy      !
         tab_dy(1:long, larg) = ( tab_in(1:long, larg   ) -      &  !
                                  tab_in(1:long, larg -1) )/hy      !

         tab_xy(   1, 1:larg) = ( tab_dy(      2, 1:larg) -      &  !
                                  tab_dy(      1, 1:larg) )/hx      !
         tab_xy(long, 1:larg) = ( tab_dy(long   , 1:larg) -      &  !
                                  tab_dy(long -1, 1:larg) )/hx      !

         tab_xy(1:long,    1) = ( tab_dx(1:long,       2) -      &  !
                                  tab_dx(1:long,       1) )/hy      !
         tab_xy(1:long, larg) = ( tab_dx(1:long, larg   ) -      &  !
                                  tab_dx(1:long, larg -1) )/hy      !

      return
      endsubroutine calcul_tabd_hermite

      subroutine calcul_aire_hermite(tab_in, long, larg, gw, tab_dnq, ng, hx, hy, width, height, aire)
      !================================================================================================
      !< @note
      !
      !  @endnote
      !------------------------------------------------------------------------------------------------
      implicit none
      integer(kind=I4), intent(in )                            :: long
      integer(kind=I4), intent(in )                            :: larg
      integer(kind=I4), intent(in )                            :: ng
      real   (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab_in
      real   (kind=R8), intent(in ), dimension(1:ng)           :: gw
      real   (kind=R8), intent(in ), dimension(1:32*ng)        :: tab_dnq
      real   (kind=R8), intent(in )                            :: hx
      real   (kind=R8), intent(in )                            :: hy
      real   (kind=R8), intent(in )                            :: width
      real   (kind=R8), intent(in )                            :: height
      real   (kind=R8), intent(out)                            :: aire

         integer(kind=I4) :: i, j, k, i1, i2, j1, j2
         real   (kind=R8) :: aire_tmp

         real(kind=R8), allocatable, dimension(:) :: dfx
         real(kind=R8), allocatable, dimension(:) :: dfy

         real(kind=R8), allocatable, dimension(:,:) :: tab_dx     ! new grid function evaluations
         real(kind=R8), allocatable, dimension(:,:) :: tab_dy     ! new grid function evaluations
         real(kind=R8), allocatable, dimension(:,:) :: tab_xy     ! new grid function evaluations

         allocate( dfx(1:ng),    &  !
                   dfy(1:ng) )      !

         allocate( tab_dx   (1:long, 1:larg),   &  !
                   tab_dy   (1:long, 1:larg),   &  !
                   tab_xy   (1:long, 1:larg) )     !

         call calcul_tabd_hermite( tab_in = tab_in(1:long, 1:larg),  &  !
                                   tab_dx = tab_dx(1:long, 1:larg),  &  !
                                   tab_dy = tab_dy(1:long, 1:larg),  &  !
                                   tab_xy = tab_xy(1:long, 1:larg),  &  !
                                     long = long,                    &  !
                                     larg = larg,                    &  !
                                       hx = hx,                      &  !
                                       hy = hy )                        !


         aire_tmp = 0._R8
         do j = 1, larg -1

            j1 = j ; j2 = j +1

            do i = 1, long -1

               i1 = i ; i2 = i +1

               do k = 1, ng

                  dfx(k) = (2._R8/hx)*(                              &  !
                  tab_dnq(32*(k -1) +01)*tab_in(i1, j1) +            &  !  u1
                  tab_dnq(32*(k -1) +02)*tab_in(i2, j1) +            &  !  u2
                  tab_dnq(32*(k -1) +03)*tab_in(i2, j2) +            &  !  u3
                  tab_dnq(32*(k -1) +04)*tab_in(i1, j2) +            &  !  u4

                  tab_dnq(32*(k -1) +05)*tab_dx(i1, j1)*hx/2 +       &  ! du1/dx
                  tab_dnq(32*(k -1) +06)*tab_dx(i2, j1)*hx/2 +       &  ! du2/dx
                  tab_dnq(32*(k -1) +07)*tab_dx(i2, j2)*hx/2 +       &  ! du3/dx
                  tab_dnq(32*(k -1) +08)*tab_dx(i1, j2)*hx/2 +       &  ! du4/dx

                  tab_dnq(32*(k -1) +09)*tab_dy(i1, j1)*hy/2 +       &  ! du1/dy
                  tab_dnq(32*(k -1) +10)*tab_dy(i2, j1)*hy/2 +       &  ! du2/dy
                  tab_dnq(32*(k -1) +11)*tab_dy(i2, j2)*hy/2 +       &  ! du3/dy
                  tab_dnq(32*(k -1) +12)*tab_dy(i1, j2)*hy/2 +       &  ! du4/dy

                  tab_dnq(32*(k -1) +13)*tab_xy(i1, j1)*hx*hy/4 +    &  ! du1/dxdy
                  tab_dnq(32*(k -1) +14)*tab_xy(i2, j1)*hx*hy/4 +    &  ! du2/dxdy
                  tab_dnq(32*(k -1) +15)*tab_xy(i2, j2)*hx*hy/4 +    &  ! du3/dxdy
                  tab_dnq(32*(k -1) +16)*tab_xy(i1, j2)*hx*hy/4 )       ! du4/dxdy

                  dfy(k) = (2._R8/hy)*(                              &  !
                  tab_dnq(32*(k -1) +17)*tab_in(i1, j1) +            &  !  u1
                  tab_dnq(32*(k -1) +18)*tab_in(i2, j1) +            &  !  u2
                  tab_dnq(32*(k -1) +19)*tab_in(i2, j2) +            &  !  u3
                  tab_dnq(32*(k -1) +20)*tab_in(i1, j2) +            &  !  u4

                  tab_dnq(32*(k -1) +21)*tab_dx(i1, j1)*hx/2 +       &  ! du1/dx
                  tab_dnq(32*(k -1) +22)*tab_dx(i2, j1)*hx/2 +       &  ! du2/dx
                  tab_dnq(32*(k -1) +23)*tab_dx(i2, j2)*hx/2 +       &  ! du3/dx
                  tab_dnq(32*(k -1) +24)*tab_dx(i1, j2)*hx/2 +       &  ! du4/dx

                  tab_dnq(32*(k -1) +25)*tab_dy(i1, j1)*hy/2 +       &  ! du1/dy
                  tab_dnq(32*(k -1) +26)*tab_dy(i2, j1)*hy/2 +       &  ! du2/dy
                  tab_dnq(32*(k -1) +27)*tab_dy(i2, j2)*hy/2 +       &  ! du3/dy
                  tab_dnq(32*(k -1) +28)*tab_dy(i1, j2)*hy/2 +       &  ! du4/dy

                  tab_dnq(32*(k -1) +29)*tab_xy(i1, j1)*hx*hy/4 +    &  ! du1/dxdy
                  tab_dnq(32*(k -1) +30)*tab_xy(i2, j1)*hx*hy/4 +    &  ! du2/dxdy
                  tab_dnq(32*(k -1) +31)*tab_xy(i2, j2)*hx*hy/4 +    &  ! du3/dxdy
                  tab_dnq(32*(k -1) +32)*tab_xy(i1, j2)*hx*hy/4 )       ! du4/dxdy

               enddo

               do k = 1, ng
                  aire_tmp = aire_tmp + gw(k) * sqrt( UN + dfx(k)**2 + dfy(k)**2 )
               enddo

            enddo

         enddo

         aire = aire_tmp*( hx/2 )*( hy/2 )
         aire = aire/(width*height)

         deallocate( tab_dx   ,  &  !
                     tab_dy   ,  &  !
                     tab_xy )       !

         deallocate( dfx, dfy )

      return
      endsubroutine calcul_aire_hermite

   endsubroutine calcul_asfc_hermite