SUBMODULE(compute_caldyn_mod) compute_pvort_only_mod

   IMPLICIT NONE

CONTAINS

   MODULE SUBROUTINE compute_pvort_only(u, rhodz, qu, qv)
      USE prec, ONLY: rstd
      USE grid_param, ONLY: grid_type, grid_ico
      USE geometry, ONLY: Av, Riv2, fv
      REAL(rstd), INTENT(IN)  :: u(:, :)
      REAL(rstd), INTENT(INOUT) :: rhodz(:, :)
      REAL(rstd), INTENT(OUT) :: qu(:, :), qv(:, :)
      IF (grid_type == grid_ico) THEN
         CALL compute_pvort_only_hex_(u, rhodz, qu, qv, Av, Riv2, fv)
      ELSE
         CALL compute_pvort_only_unst_(u, rhodz, qu, qv)
      END IF
   END SUBROUTINE compute_pvort_only

   SUBROUTINE compute_pvort_only_unst_(u, rhodz, qu, qv)
      USE ISO_C_BINDING, only: C_DOUBLE, C_FLOAT
      USE geometry, ONLY: R_iv2, Av, fv
      USE grid_param, ONLY: llm, dual_num, edge_num, primal_num
      USE data_unstructured_mod, ONLY: enter_trace, exit_trace, id_pvort_only, &
                                       dual_deg, dual_edge, dual_ne, dual_vertex, up, down
      REAL(C_DOUBLE), DIMENSION(llm, primal_num) :: rhodz
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)    :: u, qu
      REAL(C_DOUBLE), DIMENSION(llm, dual_num)    :: qv
      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, hv
      CALL enter_trace(id_pvort_only, 8*llm*((1)*primal_num + (1)*dual_num + (2)*edge_num)) ! primal, dual, edge

      !--------------------------------------------------------------------------
      !---------------------------- pvort_only ----------------------------------
      !$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)
            sign1 = dual_ne(1, ij)
            sign2 = dual_ne(2, ij)
            sign3 = dual_ne(3, ij)
            vertex1 = dual_vertex(1, ij)
            vertex2 = dual_vertex(2, ij)
            vertex3 = dual_vertex(3, ij)
            !DIR$ SIMD
            DO l = 1, llm
               etav = 0.d0
               etav = etav + sign1*u(l, edge1)
               etav = etav + sign2*u(l, edge2)
               etav = etav + sign3*u(l, edge3)
               hv = 0.
               hv = hv + R_iv2(1, ij)*rhodz(l, vertex1)
               hv = hv + R_iv2(2, ij)*rhodz(l, vertex2)
               hv = hv + R_iv2(3, ij)*rhodz(l, vertex3)
               qv(l, ij) = (etav + fv(ij)*Av(ij))/(hv*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)
            sign1 = dual_ne(1, ij)
            sign2 = dual_ne(2, ij)
            sign3 = dual_ne(3, ij)
            sign4 = dual_ne(4, ij)
            vertex1 = dual_vertex(1, ij)
            vertex2 = dual_vertex(2, ij)
            vertex3 = dual_vertex(3, ij)
            vertex4 = dual_vertex(4, ij)
            !DIR$ SIMD
            DO l = 1, llm
               etav = 0.d0
               etav = etav + sign1*u(l, edge1)
               etav = etav + sign2*u(l, edge2)
               etav = etav + sign3*u(l, edge3)
               etav = etav + sign4*u(l, edge4)
               hv = 0.
               hv = hv + R_iv2(1, ij)*rhodz(l, vertex1)
               hv = hv + R_iv2(2, ij)*rhodz(l, vertex2)
               hv = hv + R_iv2(3, ij)*rhodz(l, vertex3)
               hv = hv + R_iv2(4, ij)*rhodz(l, vertex4)
               qv(l, ij) = (etav + fv(ij)*Av(ij))/(hv*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)*u(l, edge)
               END DO
               hv = 0.
               DO ivertex = 1, dual_deg(ij)
                  vertex = dual_vertex(ivertex, ij)
                  hv = hv + R_iv2(ivertex, ij)*rhodz(l, vertex)
               END DO
               qv(l, ij) = (etav + fv(ij)*Av(ij))/(hv*Av(ij))
            END DO
         END SELECT
      END DO
      !$OMP END DO
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         ij_up = up(edge)
         ij_down = down(edge)
         !DIR$ SIMD
         DO l = 1, llm
            qu(l, edge) = 0.5d0*(qv(l, ij_down) + qv(l, ij_up))
         END DO
      END DO
      !$OMP END DO
      !---------------------------- pvort_only ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_pvort_only_unst_

   SUBROUTINE compute_pvort_only_hex_(u, rhodz, qu, qv, Av_, Riv2_, fv_)
      USE icosa
      USE trace, ONLY: trace_start, trace_end
      USE caldyn_vars_mod, ONLY: dysl_pvort_only
      USE omp_para, ONLY: ll_begin, ll_end
      REAL(rstd), INTENT(IN)  :: u(iim*3*jjm, llm)
      REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm, llm)
      REAL(rstd), INTENT(OUT) :: qu(iim*3*jjm, llm)
      REAL(rstd), INTENT(OUT) :: qv(iim*2*jjm, llm)
      REAL(rstd), INTENT(IN)  :: Av_(2*iim*jjm)
      REAL(rstd), INTENT(IN)  :: Riv2_(iim*jjm, 6)
      REAL(rstd), INTENT(IN)  :: fv_(2*iim*jjm)

      INTEGER :: ij, l
      REAL(rstd) :: etav, hv, radius_m2

      CALL trace_start("compute_pvort_only")

      !!! Compute shallow-water potential vorticity

      !$acc data async &
      !$acc present(rhodz(:,:), u(:,:), qu(:,:), qv(:,:), Av_(:), Riv2_(:,:), fv_(:))

      radius_m2 = radius**(-2)
      !$acc parallel loop collapse(2) async
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij = ij_begin_ext, ij_end_ext
            etav = 1./Av_(ij + z_up)*(ne_rup*u(ij + u_rup, l) &
                                      + ne_left*u(ij + t_rup + u_left, l) &
                                      - ne_lup*u(ij + u_lup, l))
            hv = Riv2_(ij, vup)*rhodz(ij, l) &
                 + Riv2_(ij + t_rup, vldown)*rhodz(ij + t_rup, l) &
                 + Riv2_(ij + t_lup, vrdown)*rhodz(ij + t_lup, l)
            qv(ij + z_up, l) = (etav + fv_(ij + z_up))/hv

            etav = 1./Av_(ij + z_down)*(ne_ldown*u(ij + u_ldown, l) &
                                        + ne_right*u(ij + t_ldown + u_right, l) &
                                        - ne_rdown*u(ij + u_rdown, l))
            hv = Riv2_(ij, vdown)*rhodz(ij, l) &
                 + Riv2_(ij + t_ldown, vrup)*rhodz(ij + t_ldown, l) &
                 + Riv2_(ij + t_rdown, vlup)*rhodz(ij + t_rdown, l)
            qv(ij + z_down, l) = (etav + fv_(ij + z_down))/hv
         END DO
      END DO

      !$acc parallel loop collapse(2) async
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij = ij_begin, ij_end
            qu(ij + u_right, l) = 0.5*(qv(ij + z_rdown, l) + qv(ij + z_rup, l))
            qu(ij + u_lup, l) = 0.5*(qv(ij + z_up, l) + qv(ij + z_lup, l))
            qu(ij + u_ldown, l) = 0.5*(qv(ij + z_ldown, l) + qv(ij + z_down, l))
         END DO
      END DO

      !$acc end data

      CALL trace_end("compute_pvort_only")

   END SUBROUTINE compute_pvort_only_hex_

END SUBMODULE compute_pvort_only_mod
