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

   PUBLIC :: vorticity, compute_vorticity_hex, compute_vorticity_unst

CONTAINS

   SUBROUTINE vorticity(f_ue, f_vort)
      USE icosa
      USE compute_diagnostics_mod, ONLY: compute_vorticity
      TYPE(t_field), POINTER :: f_ue(:)
      TYPE(t_field), POINTER :: f_vort(:)

      REAL(rstd), POINTER :: ue(:, :)
      REAL(rstd), POINTER :: vort(:, :)
      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)
         vort = f_vort(ind)
         CALL compute_vorticity(ue, vort)
      END DO

   END SUBROUTINE vorticity

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

   SUBROUTINE compute_vorticity_unst(ue, vort)
      REAL(rstd), INTENT(IN)  :: ue(:, :)
      REAL(rstd), INTENT(OUT) :: vort(:, :)
      CALL compute_vorticity_unst_(ue, vort)
   END SUBROUTINE compute_vorticity_unst

   SUBROUTINE compute_vorticity_hex(ue, vort)
      REAL(rstd), INTENT(IN)  :: ue(:, :)
      REAL(rstd), INTENT(OUT) :: vort(:, :)
      CALL compute_vorticity_hex_(ue, vort)
   END SUBROUTINE compute_vorticity_hex

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

   SUBROUTINE compute_vorticity_unst_(ue, vort)
      USE ISO_C_BINDING, only: C_DOUBLE, C_FLOAT
      USE geometry, ONLY: de, Av
      USE data_unstructured_mod, ONLY: dual_deg, dual_edge, dual_ne
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)    :: ue
      REAL(C_DOUBLE), DIMENSION(llm, dual_num)    :: vort
      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) :: etav

      !--------------------------------------------------------------------------
      !---------------------------- compute_vorticity ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, dual_num
         ! this VLOOP iterates over dual cell edges
         SELECT CASE (dual_deg(ij))
         CASE (3)
            edge1 = dual_edge(1, ij)
            edge2 = dual_edge(2, ij)
            edge3 = dual_edge(3, ij)
            de1 = de(edge1)
            de2 = de(edge2)
            de3 = de(edge3)
            sign1 = dual_ne(1, ij)
            sign2 = dual_ne(2, ij)
            sign3 = dual_ne(3, ij)
            !DIR$ SIMD
            DO l = 1, llm
               etav = 0.d0
               etav = etav + sign1*ue(l, edge1)*de1
               etav = etav + sign2*ue(l, edge2)*de2
               etav = etav + sign3*ue(l, edge3)*de3
               vort(l, ij) = etav/Av(ij)
            END DO
         CASE (4)
            edge1 = dual_edge(1, ij)
            edge2 = dual_edge(2, ij)
            edge3 = dual_edge(3, ij)
            edge4 = dual_edge(4, ij)
            de1 = de(edge1)
            de2 = de(edge2)
            de3 = de(edge3)
            de4 = de(edge4)
            sign1 = dual_ne(1, ij)
            sign2 = dual_ne(2, ij)
            sign3 = dual_ne(3, ij)
            sign4 = dual_ne(4, ij)
            !DIR$ SIMD
            DO l = 1, llm
               etav = 0.d0
               etav = etav + sign1*ue(l, edge1)*de1
               etav = etav + sign2*ue(l, edge2)*de2
               etav = etav + sign3*ue(l, edge3)*de3
               etav = etav + sign4*ue(l, edge4)*de4
               vort(l, ij) = etav/Av(ij)
            END DO
         CASE DEFAULT
            !DIR$ SIMD
            DO l = 1, llm
               etav = 0.d0
               DO iedge = 1, dual_deg(ij)
                  edge = dual_edge(iedge, ij)
                  etav = etav + dual_ne(iedge, ij)*ue(l, edge)*de(edge)
               END DO
               vort(l, ij) = etav/Av(ij)
            END DO
         END SELECT
      END DO
      !$OMP END DO
      !---------------------------- compute_vorticity ----------------------------------
      !--------------------------------------------------------------------------
   END SUBROUTINE compute_vorticity_unst_

   SUBROUTINE compute_vorticity_hex_(ue, vort)
      USE icosa
      USE omp_para
      REAL(rstd), INTENT(IN)  :: ue(3*iim*jjm, llm)
      REAL(rstd), INTENT(OUT) :: vort(2*iim*jjm, llm)
      REAL(rstd) :: etav
      INTEGER :: ij, l

      !--------------------------------------------------------------------------
      !---------------------------- compute_vorticity ----------------------------------
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij = ij_begin_ext, ij_end_ext
            etav = 0.d0
            etav = etav + ne_rup*ue(ij + u_rup, l)*de(ij + u_rup)
            etav = etav + ne_left*ue(ij + t_rup + u_left, l)*de(ij + t_rup + u_left)
            etav = etav + (-ne_lup)*ue(ij + u_lup, l)*de(ij + u_lup)
            vort(ij + z_up, l) = etav/Av(ij + z_up)
            etav = 0.d0
            etav = etav + (-ne_rdown)*ue(ij + u_rdown, l)*de(ij + u_rdown)
            etav = etav + ne_right*ue(ij + t_ldown + u_right, l)*de(ij + t_ldown + u_right)
            etav = etav + ne_ldown*ue(ij + u_ldown, l)*de(ij + u_ldown)
            vort(ij + z_down, l) = etav/Av(ij + z_down)
         END DO
      END DO
      !---------------------------- compute_vorticity ----------------------------------
      !--------------------------------------------------------------------------
   END SUBROUTINE compute_vorticity_hex_

END MODULE compute_vorticity_mod
