MODULE compute_divergence_mod
   USE grid_param
   USE prec, ONLY: rstd
   IMPLICIT NONE
   PRIVATE

   PUBLIC :: divergence, compute_divergence_hex, compute_divergence_unst

CONTAINS

   SUBROUTINE divergence(f_ue, f_div)
      USE icosa
      USE compute_diagnostics_mod, ONLY: compute_divergence
      TYPE(t_field), POINTER :: f_ue(:)
      TYPE(t_field), POINTER :: f_div(:)

      REAL(rstd), POINTER :: ue(:, :)
      REAL(rstd), POINTER :: div(:, :)
      INTEGER :: ind

      CALL transfert_request(f_ue, req_e1_vect)

      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
         ue = f_ue(ind)
         div = f_div(ind)
         CALL compute_divergence(ue, div)
      END DO

   END SUBROUTINE divergence

!-------------- Wrappers for F2008 conformity -----------------

   SUBROUTINE compute_divergence_unst(ue, div)
      REAL(rstd), INTENT(IN)  :: ue(:, :)
      REAL(rstd), INTENT(OUT) :: div(:, :)
      CALL compute_divergence_unst_(ue, div)
   END SUBROUTINE compute_divergence_unst

   SUBROUTINE compute_divergence_hex(ue, div)
      REAL(rstd), INTENT(IN)  :: ue(:, :)
      REAL(rstd), INTENT(OUT) :: div(:, :)
      CALL compute_divergence_hex_(ue, div)
   END SUBROUTINE compute_divergence_hex

!--------------------------------------------------------------

   SUBROUTINE compute_divergence_unst_(ue, div)
      USE ISO_C_BINDING, only: C_DOUBLE, C_FLOAT
      USE geometry, ONLY: le, Ai
      USE data_unstructured_mod, ONLY: primal_deg, primal_edge, primal_ne
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)    :: ue
      REAL(C_DOUBLE), DIMENSION(llm, primal_num) :: div
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: sign1, sign2, sign3, sign4, sign5, sign6, le_de1, le_de2, le_de3, le_de4, le_de5, le_de6, &
                        de1, de2, de3, de4, de5, de6, le1, le2, le3, le4, le5, le6
      INTEGER :: edge1, edge2, edge3, edge4, edge5, edge6, ij_up1, ij_up2, ij_up3, ij_up4, ij_up5, ij_up6, &
                 ij_down1, ij_down2, ij_down3, ij_down4, ij_down5, ij_down6
      INTEGER ::vertex1, vertex2, vertex3, vertex4, vertex5, vertex6
      REAL(C_DOUBLE) :: div_ij

      !--------------------------------------------------------------------------
      !---------------------------- compute_divergence ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         ! this VLOOP iterates over primal cell edges
         SELECT CASE (primal_deg(ij))
         CASE (4)
            edge1 = primal_edge(1, ij)
            edge2 = primal_edge(2, ij)
            edge3 = primal_edge(3, ij)
            edge4 = primal_edge(4, ij)
            le1 = le(edge1)
            le2 = le(edge2)
            le3 = le(edge3)
            le4 = le(edge4)
            sign1 = primal_ne(1, ij)
            sign2 = primal_ne(2, ij)
            sign3 = primal_ne(3, ij)
            sign4 = primal_ne(4, ij)
            !DIR$ SIMD
            DO l = 1, llm
               div_ij = 0.d0
               div_ij = div_ij + sign1*ue(l, edge1)*le1
               div_ij = div_ij + sign2*ue(l, edge2)*le2
               div_ij = div_ij + sign3*ue(l, edge3)*le3
               div_ij = div_ij + sign4*ue(l, edge4)*le4
               div(l, ij) = div_ij/Ai(ij)
            END DO
         CASE (5)
            edge1 = primal_edge(1, ij)
            edge2 = primal_edge(2, ij)
            edge3 = primal_edge(3, ij)
            edge4 = primal_edge(4, ij)
            edge5 = primal_edge(5, ij)
            le1 = le(edge1)
            le2 = le(edge2)
            le3 = le(edge3)
            le4 = le(edge4)
            le5 = le(edge5)
            sign1 = primal_ne(1, ij)
            sign2 = primal_ne(2, ij)
            sign3 = primal_ne(3, ij)
            sign4 = primal_ne(4, ij)
            sign5 = primal_ne(5, ij)
            !DIR$ SIMD
            DO l = 1, llm
               div_ij = 0.d0
               div_ij = div_ij + sign1*ue(l, edge1)*le1
               div_ij = div_ij + sign2*ue(l, edge2)*le2
               div_ij = div_ij + sign3*ue(l, edge3)*le3
               div_ij = div_ij + sign4*ue(l, edge4)*le4
               div_ij = div_ij + sign5*ue(l, edge5)*le5
               div(l, ij) = div_ij/Ai(ij)
            END DO
         CASE (6)
            edge1 = primal_edge(1, ij)
            edge2 = primal_edge(2, ij)
            edge3 = primal_edge(3, ij)
            edge4 = primal_edge(4, ij)
            edge5 = primal_edge(5, ij)
            edge6 = primal_edge(6, ij)
            le1 = le(edge1)
            le2 = le(edge2)
            le3 = le(edge3)
            le4 = le(edge4)
            le5 = le(edge5)
            le6 = le(edge6)
            sign1 = primal_ne(1, ij)
            sign2 = primal_ne(2, ij)
            sign3 = primal_ne(3, ij)
            sign4 = primal_ne(4, ij)
            sign5 = primal_ne(5, ij)
            sign6 = primal_ne(6, ij)
            !DIR$ SIMD
            DO l = 1, llm
               div_ij = 0.d0
               div_ij = div_ij + sign1*ue(l, edge1)*le1
               div_ij = div_ij + sign2*ue(l, edge2)*le2
               div_ij = div_ij + sign3*ue(l, edge3)*le3
               div_ij = div_ij + sign4*ue(l, edge4)*le4
               div_ij = div_ij + sign5*ue(l, edge5)*le5
               div_ij = div_ij + sign6*ue(l, edge6)*le6
               div(l, ij) = div_ij/Ai(ij)
            END DO
         CASE DEFAULT
            !DIR$ SIMD
            DO l = 1, llm
               div_ij = 0.d0
               DO iedge = 1, primal_deg(ij)
                  edge = primal_edge(iedge, ij)
                  div_ij = div_ij + primal_ne(iedge, ij)*ue(l, edge)*le(edge)
               END DO
               div(l, ij) = div_ij/Ai(ij)
            END DO
         END SELECT
      END DO
      !$OMP END DO
      !---------------------------- compute_divergence ----------------------------------
      !--------------------------------------------------------------------------
   END SUBROUTINE compute_divergence_unst_

   SUBROUTINE compute_divergence_hex_(ue, div)
      USE icosa
      USE omp_para
      REAL(rstd), INTENT(IN)  :: ue(3*iim*jjm, llm)
      REAL(rstd), INTENT(OUT) :: div(iim*jjm, llm)
      REAL(rstd) :: div_ij
      INTEGER :: ij, l

      !--------------------------------------------------------------------------
      !---------------------------- compute_divergence ----------------------------------
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij = ij_begin_ext, ij_end_ext
            div_ij = 0.d0
            div_ij = div_ij + ne_rup*ue(ij + u_rup, l)*le(ij + u_rup)
            div_ij = div_ij + ne_lup*ue(ij + u_lup, l)*le(ij + u_lup)
            div_ij = div_ij + ne_left*ue(ij + u_left, l)*le(ij + u_left)
            div_ij = div_ij + ne_ldown*ue(ij + u_ldown, l)*le(ij + u_ldown)
            div_ij = div_ij + ne_rdown*ue(ij + u_rdown, l)*le(ij + u_rdown)
            div_ij = div_ij + ne_right*ue(ij + u_right, l)*le(ij + u_right)
            div(ij, l) = div_ij/Ai(ij)
         END DO
      END DO
      !---------------------------- compute_divergence ----------------------------------
      !--------------------------------------------------------------------------
   END SUBROUTINE compute_divergence_hex_

END MODULE compute_divergence_mod
