MODULE diags_unstructured_mod
   USE ISO_C_BINDING
   USE OMP_LIB
   USE data_unstructured_mod
   USE transfer_unstructured_mod
   IMPLICIT NONE
   PRIVATE
   SAVE

CONTAINS

   SUBROUTINE wind_centered(un, uxyz) BIND(C, name='dynamico_wind_centered')
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn)  :: uxyz ! OUT
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)      :: un    ! IN, normal components
      !$OMP PARALLEL NUM_THREADS(nb_threads)
      CALL update_halo(transfer_edge, un)
      CALL compute_wind_centered(un, uxyz)
      !$OMP END PARALLEL
   END SUBROUTINE wind_centered

   SUBROUTINE compute_wind_centered(ue, ucenter)
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn)  :: ucenter ! OUT
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)      :: ue    ! IN, normal components
      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
      REAL(C_DOUBLE) :: cx, cy, cz, ux, uy, uz, ue_le, fac
      REAL(C_DOUBLE), parameter :: scale = 1.

      !--------------------------------------------------------------------------
      !---------------------------- wind_centered ----------------------------------
      ! Perot reconstruction based on Gauss theorem
      ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
      !$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)
            sign1 = primal_ne(1, ij)
            sign2 = primal_ne(2, ij)
            sign3 = primal_ne(3, ij)
            sign4 = primal_ne(4, ij)
            ij_up1 = up(edge1)
            ij_up2 = up(edge2)
            ij_up3 = up(edge3)
            ij_up4 = up(edge4)
            ij_down1 = down(edge1)
            ij_down2 = down(edge2)
            ij_down3 = down(edge3)
            ij_down4 = down(edge4)
            !DIR$ SIMD
            DO l = 1, llm
               cx = centroid(ij, 1)
               cy = centroid(ij, 2)
               cz = centroid(ij, 3)
               ux = 0.; uy = 0.; uz = 0.
               ue_le = sign1*ue(l, edge1)*le(edge1)
               ux = ux + ue_le*(.5*(xyz_v(ij_up1, 1) + xyz_v(ij_down1, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up1, 2) + xyz_v(ij_down1, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up1, 3) + xyz_v(ij_down1, 3)) - cz)
               ue_le = sign2*ue(l, edge2)*le(edge2)
               ux = ux + ue_le*(.5*(xyz_v(ij_up2, 1) + xyz_v(ij_down2, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up2, 2) + xyz_v(ij_down2, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up2, 3) + xyz_v(ij_down2, 3)) - cz)
               ue_le = sign3*ue(l, edge3)*le(edge3)
               ux = ux + ue_le*(.5*(xyz_v(ij_up3, 1) + xyz_v(ij_down3, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up3, 2) + xyz_v(ij_down3, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up3, 3) + xyz_v(ij_down3, 3)) - cz)
               ue_le = sign4*ue(l, edge4)*le(edge4)
               ux = ux + ue_le*(.5*(xyz_v(ij_up4, 1) + xyz_v(ij_down4, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up4, 2) + xyz_v(ij_down4, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up4, 3) + xyz_v(ij_down4, 3)) - cz)
               fac = scale*(1./Ai(ij))
               ucenter(l, ij, 1) = ux*fac
               ucenter(l, ij, 2) = uy*fac
               ucenter(l, ij, 3) = uz*fac
            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)
            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)
            ij_up1 = up(edge1)
            ij_up2 = up(edge2)
            ij_up3 = up(edge3)
            ij_up4 = up(edge4)
            ij_up5 = up(edge5)
            ij_down1 = down(edge1)
            ij_down2 = down(edge2)
            ij_down3 = down(edge3)
            ij_down4 = down(edge4)
            ij_down5 = down(edge5)
            !DIR$ SIMD
            DO l = 1, llm
               cx = centroid(ij, 1)
               cy = centroid(ij, 2)
               cz = centroid(ij, 3)
               ux = 0.; uy = 0.; uz = 0.
               ue_le = sign1*ue(l, edge1)*le(edge1)
               ux = ux + ue_le*(.5*(xyz_v(ij_up1, 1) + xyz_v(ij_down1, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up1, 2) + xyz_v(ij_down1, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up1, 3) + xyz_v(ij_down1, 3)) - cz)
               ue_le = sign2*ue(l, edge2)*le(edge2)
               ux = ux + ue_le*(.5*(xyz_v(ij_up2, 1) + xyz_v(ij_down2, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up2, 2) + xyz_v(ij_down2, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up2, 3) + xyz_v(ij_down2, 3)) - cz)
               ue_le = sign3*ue(l, edge3)*le(edge3)
               ux = ux + ue_le*(.5*(xyz_v(ij_up3, 1) + xyz_v(ij_down3, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up3, 2) + xyz_v(ij_down3, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up3, 3) + xyz_v(ij_down3, 3)) - cz)
               ue_le = sign4*ue(l, edge4)*le(edge4)
               ux = ux + ue_le*(.5*(xyz_v(ij_up4, 1) + xyz_v(ij_down4, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up4, 2) + xyz_v(ij_down4, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up4, 3) + xyz_v(ij_down4, 3)) - cz)
               ue_le = sign5*ue(l, edge5)*le(edge5)
               ux = ux + ue_le*(.5*(xyz_v(ij_up5, 1) + xyz_v(ij_down5, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up5, 2) + xyz_v(ij_down5, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up5, 3) + xyz_v(ij_down5, 3)) - cz)
               fac = scale*(1./Ai(ij))
               ucenter(l, ij, 1) = ux*fac
               ucenter(l, ij, 2) = uy*fac
               ucenter(l, ij, 3) = uz*fac
            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)
            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)
            ij_up1 = up(edge1)
            ij_up2 = up(edge2)
            ij_up3 = up(edge3)
            ij_up4 = up(edge4)
            ij_up5 = up(edge5)
            ij_up6 = up(edge6)
            ij_down1 = down(edge1)
            ij_down2 = down(edge2)
            ij_down3 = down(edge3)
            ij_down4 = down(edge4)
            ij_down5 = down(edge5)
            ij_down6 = down(edge6)
            !DIR$ SIMD
            DO l = 1, llm
               cx = centroid(ij, 1)
               cy = centroid(ij, 2)
               cz = centroid(ij, 3)
               ux = 0.; uy = 0.; uz = 0.
               ue_le = sign1*ue(l, edge1)*le(edge1)
               ux = ux + ue_le*(.5*(xyz_v(ij_up1, 1) + xyz_v(ij_down1, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up1, 2) + xyz_v(ij_down1, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up1, 3) + xyz_v(ij_down1, 3)) - cz)
               ue_le = sign2*ue(l, edge2)*le(edge2)
               ux = ux + ue_le*(.5*(xyz_v(ij_up2, 1) + xyz_v(ij_down2, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up2, 2) + xyz_v(ij_down2, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up2, 3) + xyz_v(ij_down2, 3)) - cz)
               ue_le = sign3*ue(l, edge3)*le(edge3)
               ux = ux + ue_le*(.5*(xyz_v(ij_up3, 1) + xyz_v(ij_down3, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up3, 2) + xyz_v(ij_down3, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up3, 3) + xyz_v(ij_down3, 3)) - cz)
               ue_le = sign4*ue(l, edge4)*le(edge4)
               ux = ux + ue_le*(.5*(xyz_v(ij_up4, 1) + xyz_v(ij_down4, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up4, 2) + xyz_v(ij_down4, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up4, 3) + xyz_v(ij_down4, 3)) - cz)
               ue_le = sign5*ue(l, edge5)*le(edge5)
               ux = ux + ue_le*(.5*(xyz_v(ij_up5, 1) + xyz_v(ij_down5, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up5, 2) + xyz_v(ij_down5, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up5, 3) + xyz_v(ij_down5, 3)) - cz)
               ue_le = sign6*ue(l, edge6)*le(edge6)
               ux = ux + ue_le*(.5*(xyz_v(ij_up6, 1) + xyz_v(ij_down6, 1)) - cx)
               uy = uy + ue_le*(.5*(xyz_v(ij_up6, 2) + xyz_v(ij_down6, 2)) - cy)
               uz = uz + ue_le*(.5*(xyz_v(ij_up6, 3) + xyz_v(ij_down6, 3)) - cz)
               fac = scale*(1./Ai(ij))
               ucenter(l, ij, 1) = ux*fac
               ucenter(l, ij, 2) = uy*fac
               ucenter(l, ij, 3) = uz*fac
            END DO
         CASE DEFAULT
            !DIR$ SIMD
            DO l = 1, llm
               cx = centroid(ij, 1)
               cy = centroid(ij, 2)
               cz = centroid(ij, 3)
               ux = 0.; uy = 0.; uz = 0.
               DO iedge = 1, primal_deg(ij)
                  edge = primal_edge(iedge, ij)
                  ij_up = up(edge)
                  ij_down = down(edge)
                  ue_le = primal_ne(iedge, ij)*ue(l, edge)*le(edge)
                  ux = ux + ue_le*(.5*(xyz_v(ij_up, 1) + xyz_v(ij_down, 1)) - cx)
                  uy = uy + ue_le*(.5*(xyz_v(ij_up, 2) + xyz_v(ij_down, 2)) - cy)
                  uz = uz + ue_le*(.5*(xyz_v(ij_up, 3) + xyz_v(ij_down, 3)) - cz)
               END DO
               fac = scale*(1./Ai(ij))
               ucenter(l, ij, 1) = ux*fac
               ucenter(l, ij, 2) = uy*fac
               ucenter(l, ij, 3) = uz*fac
            END DO
         END SELECT
      END DO
      !$OMP END DO
      !---------------------------- wind_centered ----------------------------------
      !--------------------------------------------------------------------------
   END SUBROUTINE compute_wind_centered

END MODULE diags_unstructured_mod
