SUBMODULE(compute_caldyn_mod) compute_caldyn_slow_hydro_mod
   USE prec, ONLY: rstd
   USE grid_param
   USE earth_const
   USE disvert_mod
   USE omp_para
   USE trace
   IMPLICIT NONE

CONTAINS

   MODULE SUBROUTINE compute_caldyn_slow_hydro(zero, u, rhodz, berni, hflux, du)
      USE grid_param, ONLY: grid_type, grid_ico
      USE geometry, ONLY: Ai, le_de
      LOGICAL, INTENT(IN)      :: zero
      REAL(rstd), INTENT(IN)    :: u(:, :), rhodz(:, :)
      REAL(rstd), INTENT(OUT)  :: berni(:, :), hflux(:, :)
      REAL(rstd), INTENT(INOUT) :: du(:, :)
      IF (grid_type == grid_ico) THEN
         CALL compute_caldyn_slow_hydro_hex_(zero, u, rhodz, berni, hflux, du, Ai, le_de)
      ELSE
         CALL compute_caldyn_slow_hydro_unst_(zero, u, rhodz, berni, hflux, du)
      END IF
   END SUBROUTINE compute_caldyn_slow_hydro

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

   SUBROUTINE compute_caldyn_slow_hydro_unst_(zero, u, rhodz, berni, hflux, du)
      USE ISO_C_BINDING, only: C_DOUBLE, C_FLOAT
      USE data_unstructured_mod, ONLY: enter_trace, exit_trace, &
                                       id_slow_hydro, left, right, primal_deg, primal_edge
      LOGICAL, INTENT(IN) :: zero
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)  :: rhodz, berni  ! IN, IN, IN, BUF
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)     :: u, hflux, du ! IN, OUT, OUT
      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) :: ke, uu
      CALL enter_trace(id_slow_hydro, 8*llm*((3)*primal_num + (0)*dual_num + (3)*edge_num))

      !--------------------------------------------------------------------------
      !---------------------------- caldyn_slow_hydro ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         ij_left = left(edge)
         ij_right = right(edge)
         !DIR$ SIMD
         DO l = 1, llm
            uu = .5*(rhodz(l, ij_left) + rhodz(l, ij_right))*u(l, edge)
            hflux(l, edge) = uu*le_de(edge)
         END DO
      END DO
      !$OMP END DO
      !$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)
            le_de1 = le_de(edge1)
            le_de2 = le_de(edge2)
            le_de3 = le_de(edge3)
            le_de4 = le_de(edge4)
            !DIR$ SIMD
            DO l = 1, llm
               ke = 0.d0
               ke = ke + le_de1*u(l, edge1)**2
               ke = ke + le_de2*u(l, edge2)**2
               ke = ke + le_de3*u(l, edge3)**2
               ke = ke + le_de4*u(l, edge4)**2
               BERNI(l, ij) = ke*(.25/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)
            le_de1 = le_de(edge1)
            le_de2 = le_de(edge2)
            le_de3 = le_de(edge3)
            le_de4 = le_de(edge4)
            le_de5 = le_de(edge5)
            !DIR$ SIMD
            DO l = 1, llm
               ke = 0.d0
               ke = ke + le_de1*u(l, edge1)**2
               ke = ke + le_de2*u(l, edge2)**2
               ke = ke + le_de3*u(l, edge3)**2
               ke = ke + le_de4*u(l, edge4)**2
               ke = ke + le_de5*u(l, edge5)**2
               BERNI(l, ij) = ke*(.25/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)
            le_de1 = le_de(edge1)
            le_de2 = le_de(edge2)
            le_de3 = le_de(edge3)
            le_de4 = le_de(edge4)
            le_de5 = le_de(edge5)
            le_de6 = le_de(edge6)
            !DIR$ SIMD
            DO l = 1, llm
               ke = 0.d0
               ke = ke + le_de1*u(l, edge1)**2
               ke = ke + le_de2*u(l, edge2)**2
               ke = ke + le_de3*u(l, edge3)**2
               ke = ke + le_de4*u(l, edge4)**2
               ke = ke + le_de5*u(l, edge5)**2
               ke = ke + le_de6*u(l, edge6)**2
               BERNI(l, ij) = ke*(.25/Ai(ij))
            END DO
         CASE DEFAULT
            !DIR$ SIMD
            DO l = 1, llm
               ke = 0.d0
               DO iedge = 1, primal_deg(ij)
                  edge = primal_edge(iedge, ij)
                  ke = ke + le_de(edge)*u(l, edge)**2
               END DO
               BERNI(l, ij) = ke*(.25/Ai(ij))
            END DO
         END SELECT
      END DO
      !$OMP END DO
      IF (zero) THEN
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            ij_left = left(edge)
            ij_right = right(edge)
            !DIR$ SIMD
            DO l = 1, llm
               du(l, edge) = 1.*(berni(l, ij_left) - berni(l, ij_right)) ! minus gradient
            END DO
         END DO
         !$OMP END DO
      ELSE
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            ij_left = left(edge)
            ij_right = right(edge)
            !DIR$ SIMD
            DO l = 1, llm
               du(l, edge) = du(l, edge) + 1.*(berni(l, ij_left) - berni(l, ij_right)) ! minus gradient
            END DO
         END DO
         !$OMP END DO
      END IF
      !---------------------------- caldyn_slow_hydro ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_caldyn_slow_hydro_unst_

   SUBROUTINE compute_caldyn_slow_hydro_hex_(zero, u, rhodz, berni, hflux, du, Ai, le_de)
      USE icosa, hide_Ai => Ai, hide_lede => le_de
      USE caldyn_vars_mod
      LOGICAL, INTENT(IN) :: zero
      REAL(rstd), INTENT(IN)    :: u(3*iim*jjm, llm)    ! prognostic "velocity"
      REAL(rstd), INTENT(IN)    :: rhodz(iim*jjm, llm)
      REAL(rstd), INTENT(OUT)   :: berni(iim*jjm, llm)  ! Bernoulli function
      REAL(rstd), INTENT(OUT)   :: hflux(3*iim*jjm, llm) ! hflux in kg/s
      REAL(rstd), INTENT(INOUT) :: du(3*iim*jjm, llm)
      REAL(rstd), INTENT(IN)    :: Ai(iim*jjm)
      REAL(rstd), INTENT(IN)    :: le_de(3*iim*jjm)

      REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
      INTEGER :: ij, l

      CALL trace_start("compute_caldyn_slow_hydro")

      !$acc data async &
      !$acc present(u(:,ll_begin:ll_end), rhodz(:,ll_begin:ll_end), berni(:,ll_begin:ll_end), hflux(:,:), du(:,ll_begin:ll_end), Ai(:), le_de(:))

      !$acc parallel loop collapse(2) async
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij = ij_begin_ext, ij_end_ext
            uu = .5*(rhodz(ij, l) + rhodz(ij + t_right, l))*u(ij + u_right, l)
            hflux(ij + u_right, l) = uu*le_de(ij + u_right)
            uu = .5*(rhodz(ij, l) + rhodz(ij + t_lup, l))*u(ij + u_lup, l)
            hflux(ij + u_lup, l) = uu*le_de(ij + u_lup)
            uu = .5*(rhodz(ij, l) + rhodz(ij + t_ldown, l))*u(ij + u_ldown, l)
            hflux(ij + u_ldown, l) = uu*le_de(ij + u_ldown)
         END DO
      END DO
      !$acc parallel loop collapse(2) async
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij = ij_begin, ij_end
            ke = 0.d0
            ke = ke + le_de(ij + u_rup)*u(ij + u_rup, l)**2
            ke = ke + le_de(ij + u_lup)*u(ij + u_lup, l)**2
            ke = ke + le_de(ij + u_left)*u(ij + u_left, l)**2
            ke = ke + le_de(ij + u_ldown)*u(ij + u_ldown, l)**2
            ke = ke + le_de(ij + u_rdown)*u(ij + u_rdown, l)**2
            ke = ke + le_de(ij + u_right)*u(ij + u_right, l)**2
            berni(ij, l) = ke*(.25/Ai(ij))
         END DO
      END DO
      IF (zero) THEN
         !$acc parallel loop collapse(2) async
         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               du(ij + u_right, l) = ne_right*(berni(ij, l) - berni(ij + t_right, l)) ! minus gradient
               du(ij + u_lup, l) = ne_lup*(berni(ij, l) - berni(ij + t_lup, l)) ! minus gradient
               du(ij + u_ldown, l) = ne_ldown*(berni(ij, l) - berni(ij + t_ldown, l)) ! minus gradient
            END DO
         END DO
      ELSE
         !$acc parallel loop collapse(2) async
         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               du(ij + u_right, l) = du(ij + u_right, l) + ne_right*(berni(ij, l) - berni(ij + t_right, l)) ! minus gradient
               du(ij + u_lup, l) = du(ij + u_lup, l) + ne_lup*(berni(ij, l) - berni(ij + t_lup, l)) ! minus gradient
               du(ij + u_ldown, l) = du(ij + u_ldown, l) + ne_ldown*(berni(ij, l) - berni(ij + t_ldown, l)) ! minus gradient
            END DO
         END DO
      END IF

      !$acc end data

      CALL trace_end("compute_caldyn_slow_hydro")

   END SUBROUTINE compute_caldyn_slow_hydro_hex_

   SUBROUTINE compute_caldyn_slow_hydro_manual(zero, u, rhodz, hv, Kv, berni, hflux, du)
      USE icosa
      USE caldyn_vars_mod
      LOGICAL, INTENT(IN) :: zero
      REAL(rstd), INTENT(IN)  :: u(3*iim*jjm, llm)    ! prognostic "velocity"
      REAL(rstd), INTENT(IN)  :: Kv(2*iim*jjm, llm)   ! kinetic energy at vertices
      REAL(rstd), INTENT(IN)  :: hv(2*iim*jjm, llm)   ! height/mass averaged to vertices
      REAL(rstd), INTENT(IN)  :: rhodz(iim*jjm, llm)
      REAL(rstd), INTENT(OUT) :: hflux(3*iim*jjm, llm) ! hflux in kg/s
      REAL(rstd), INTENT(INOUT) :: du(3*iim*jjm, llm)

      REAL(rstd) :: berni(iim*jjm, llm)  ! Bernoulli function
      REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
      REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
      INTEGER :: ij, l

      CALL trace_start("compute_caldyn_slow_hydro")

      DO l = ll_begin, ll_end
         !  Compute mass fluxes
         IF (caldyn_conserv == conserv_energy) CALL test_message(req_qu)

         IF (caldyn_kinetic == kinetic_trisk) THEN
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               uu_right = 0.5*(rhodz(ij, l) + rhodz(ij + t_right, l))*u(ij + u_right, l)
               uu_lup = 0.5*(rhodz(ij, l) + rhodz(ij + t_lup, l))*u(ij + u_lup, l)
               uu_ldown = 0.5*(rhodz(ij, l) + rhodz(ij + t_ldown, l))*u(ij + u_ldown, l)
               uu_right = uu_right*le_de(ij + u_right)
               uu_lup = uu_lup*le_de(ij + u_lup)
               uu_ldown = uu_ldown*le_de(ij + u_ldown)
               hflux(ij + u_right, l) = uu_right
               hflux(ij + u_lup, l) = uu_lup
               hflux(ij + u_ldown, l) = uu_ldown
            END DO
         ELSE ! mass flux deriving from consistent kinetic energy
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               uu_right = 0.5*(hv(ij + z_rup, l) + hv(ij + z_rdown, l))*u(ij + u_right, l)
               uu_lup = 0.5*(hv(ij + z_up, l) + hv(ij + z_lup, l))*u(ij + u_lup, l)
               uu_ldown = 0.5*(hv(ij + z_ldown, l) + hv(ij + z_down, l))*u(ij + u_ldown, l)
               uu_right = uu_right*le_de(ij + u_right)
               uu_lup = uu_lup*le_de(ij + u_lup)
               uu_ldown = uu_ldown*le_de(ij + u_ldown)
               hflux(ij + u_right, l) = uu_right
               hflux(ij + u_lup, l) = uu_lup
               hflux(ij + u_ldown, l) = uu_ldown
            END DO
         END IF

         ! Compute Bernoulli=kinetic energy
         IF (caldyn_kinetic == kinetic_trisk) THEN
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               berni1(ij) = &
                  1/(4*Ai(ij))*(le_de(ij + u_right)*u(ij + u_right, l)**2 + &
                                le_de(ij + u_rup)*u(ij + u_rup, l)**2 + &
                                le_de(ij + u_lup)*u(ij + u_lup, l)**2 + &
                                le_de(ij + u_left)*u(ij + u_left, l)**2 + &
                                le_de(ij + u_ldown)*u(ij + u_ldown, l)**2 + &
                                le_de(ij + u_rdown)*u(ij + u_rdown, l)**2)
            END DO
         ELSE
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               berni1(ij) = Riv(ij, vup)*Kv(ij + z_up, l) + &
                            Riv(ij, vlup)*Kv(ij + z_lup, l) + &
                            Riv(ij, vldown)*Kv(ij + z_ldown, l) + &
                            Riv(ij, vdown)*Kv(ij + z_down, l) + &
                            Riv(ij, vrdown)*Kv(ij + z_rdown, l) + &
                            Riv(ij, vrup)*Kv(ij + z_rup, l)
            END DO
         END IF
         ! Compute du=-grad(Bernoulli)
         IF (zero) THEN
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               du(ij + u_right, l) = ne_right*(berni1(ij) - berni1(ij + t_right))
               du(ij + u_lup, l) = ne_lup*(berni1(ij) - berni1(ij + t_lup))
               du(ij + u_ldown, l) = ne_ldown*(berni1(ij) - berni1(ij + t_ldown))
            END DO
         ELSE
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               du(ij + u_right, l) = du(ij + u_right, l) + &
                                     ne_right*(berni1(ij) - berni1(ij + t_right))
               du(ij + u_lup, l) = du(ij + u_lup, l) + &
                                   ne_lup*(berni1(ij) - berni1(ij + t_lup))
               du(ij + u_ldown, l) = du(ij + u_ldown, l) + &
                                     ne_ldown*(berni1(ij) - berni1(ij + t_ldown))
            END DO
         END IF
      END DO

      CALL trace_end("compute_caldyn_slow_hydro")
   END SUBROUTINE compute_caldyn_slow_hydro_manual

END SUBMODULE compute_caldyn_slow_hydro_mod
