SUBMODULE(compute_transport_mod) compute_gradq3d_mod

   IMPLICIT NONE

CONTAINS

   MODULE SUBROUTINE compute_gradq3d(qi, sqrt_leng, gradq3d, xyz_i, xyz_v)
      USE grid_param, ONLY: grid_type, grid_ico
      REAL(rstd), INTENT(IN)  :: qi(:, :)
      REAL(rstd), INTENT(IN)  :: sqrt_leng(:)
      REAL(rstd), INTENT(IN)  :: xyz_i(:, :)
      REAL(rstd), INTENT(IN)  :: xyz_v(:, :)
      REAL(rstd), INTENT(OUT) :: gradq3d(:, :, :)

      IF (grid_type == grid_ico) THEN
         CALL compute_gradq3d_hex(qi, sqrt_leng, gradq3d, xyz_i, xyz_v)
      END IF

   END SUBROUTINE compute_gradq3d

!=======================================================================================

   SUBROUTINE compute_gradq3d_hex(qi, sqrt_leng, gradq3d, xyz_i, xyz_v)
      USE trace
      USE omp_para
      USE dimensions
      USE grid_param, ONLY: llm
      IMPLICIT NONE
      REAL(rstd), INTENT(IN)  :: qi(iim*jjm, llm)
      REAL(rstd), INTENT(IN)  :: sqrt_leng(iim*jjm)
      REAL(rstd), INTENT(IN)  :: xyz_i(iim*jjm, 3)
      REAL(rstd), INTENT(IN)  :: xyz_v(2*iim*jjm, 3)
      REAL(rstd), INTENT(OUT) :: gradq3d(iim*jjm, llm, 3)
      REAL(rstd) :: maxq, minq, minq_c, maxq_c
      REAL(rstd) :: alphamx, alphami, alpha, maggrd
      REAL(rstd) :: arr(2*iim*jjm)
      REAL(rstd) :: ar(iim*jjm)
      REAL(rstd) :: gradtri(2*iim*jjm, llm, 3)
      INTEGER :: ij, k, l
      REAL(rstd) :: detx, dety, detz, det
      REAL(rstd) :: A(3, 3), a11, a12, a13, a21, a22, a23, a31, a32, a33
      REAL(rstd) :: x1, x2, x3
      REAL(rstd) :: dq(3)

      CALL trace_start("compute_gradq3d1")

      ! TODO : precompute ar, drop arr as output argument of gradq ?

      !========================================================================================== GRADIENT
      ! Compute gradient at triangles solving a linear system
      ! arr = area of triangle joining centroids of hexagons
!     DO l = ll_begin,ll_end
!!DIR$ SIMD
!      DO ij=ij_begin_ext,ij_end_ext
!!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up))
!!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down))
!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
!       END DO
!    END DO

      !$acc data create(gradtri(:,:,:), arr(:), ar(:)) present(sqrt_leng(:), xyz_i(:,:), xyz_v(:,:), qi(:,:), gradq3d(:,:,:)) async

      !$acc parallel loop collapse(2) async private(A, dq)
      DO l = ll_begin, ll_end
!DIR$ SIMD
         DO ij = ij_begin_ext, ij_end_ext
!       CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))

            A(1, 1) = xyz_i(ij + t_rup, 1) - xyz_i(ij, 1); A(1, 2) = xyz_i(ij + t_rup, 2) - xyz_i(ij, 2); A(1, 3) = xyz_i(ij + t_rup, 3) - xyz_i(ij, 3)
            A(2, 1) = xyz_i(ij + t_lup, 1) - xyz_i(ij, 1); A(2, 2) = xyz_i(ij + t_lup, 2) - xyz_i(ij, 2); A(2, 3) = xyz_i(ij + t_lup, 3) - xyz_i(ij, 3)
            A(3, 1) = xyz_v(ij + z_up, 1); A(3, 2) = xyz_v(ij + z_up, 2); A(3, 3) = xyz_v(ij + z_up, 3)

            dq(1) = qi(ij + t_rup, l) - qi(ij, l)
            dq(2) = qi(ij + t_lup, l) - qi(ij, l)
            dq(3) = 0.0

!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)

            a11 = A(1, 1); a12 = A(2, 1); a13 = A(3, 1)
            a21 = A(1, 2); a22 = A(2, 2); a23 = A(3, 2)
            a31 = A(1, 3); a32 = A(2, 3); a33 = A(3, 3)

            x1 = a11*(a22*a33 - a23*a32)
            x2 = a12*(a21*a33 - a23*a31)
            x3 = a13*(a21*a32 - a22*a31)
            det = x1 - x2 + x3

!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)

            a11 = dq(1); a12 = dq(2); a13 = dq(3)
            a21 = A(1, 2); a22 = A(2, 2); a23 = A(3, 2)
            a31 = A(1, 3); a32 = A(2, 3); a33 = A(3, 3)

            x1 = a11*(a22*a33 - a23*a32)
            x2 = a12*(a21*a33 - a23*a31)
            x3 = a13*(a21*a32 - a22*a31)
            detx = x1 - x2 + x3

!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)

            a11 = A(1, 1); a12 = A(2, 1); a13 = A(3, 1)
            a21 = dq(1); a22 = dq(2); a23 = dq(3)
            a31 = A(1, 3); a32 = A(2, 3); a33 = A(3, 3)

            x1 = a11*(a22*a33 - a23*a32)
            x2 = a12*(a21*a33 - a23*a31)
            x3 = a13*(a21*a32 - a22*a31)
            dety = x1 - x2 + x3

!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)

            a11 = A(1, 1); a12 = A(2, 1); a13 = A(3, 1)
            a21 = A(1, 2); a22 = A(2, 2); a23 = A(3, 2)
            a31 = dq(1); a32 = dq(2); a33 = dq(3)

            x1 = a11*(a22*a33 - a23*a32)
            x2 = a12*(a21*a33 - a23*a31)
            x3 = a13*(a21*a32 - a22*a31)
            detz = x1 - x2 + x3

            gradtri(ij + z_up, l, 1) = detx
            gradtri(ij + z_up, l, 2) = dety
            gradtri(ij + z_up, l, 3) = detz
            arr(ij + z_up) = det

         END DO
      END DO

      !$acc parallel loop collapse(2) async private(A, dq)
      DO l = ll_begin, ll_end
         DO ij = ij_begin_ext, ij_end_ext

!        CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))

            A(1, 1) = xyz_i(ij + t_ldown, 1) - xyz_i(ij, 1); A(1, 2) = xyz_i(ij + t_ldown, 2) - xyz_i(ij, 2); A(1, 3) = xyz_i(ij + t_ldown, 3) - xyz_i(ij, 3)
            A(2, 1) = xyz_i(ij + t_rdown, 1) - xyz_i(ij, 1); A(2, 2) = xyz_i(ij + t_rdown, 2) - xyz_i(ij, 2); A(2, 3) = xyz_i(ij + t_rdown, 3) - xyz_i(ij, 3)
            A(3, 1) = xyz_v(ij + z_down, 1); A(3, 2) = xyz_v(ij + z_down, 2); A(3, 3) = xyz_v(ij + z_down, 3)

            dq(1) = qi(ij + t_ldown, l) - qi(ij, l)
            dq(2) = qi(ij + t_rdown, l) - qi(ij, l)
            dq(3) = 0.0

!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)

            a11 = A(1, 1); a12 = A(2, 1); a13 = A(3, 1)
            a21 = A(1, 2); a22 = A(2, 2); a23 = A(3, 2)
            a31 = A(1, 3); a32 = A(2, 3); a33 = A(3, 3)

            x1 = a11*(a22*a33 - a23*a32)
            x2 = a12*(a21*a33 - a23*a31)
            x3 = a13*(a21*a32 - a22*a31)
            det = x1 - x2 + x3

!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)

            a11 = dq(1); a12 = dq(2); a13 = dq(3)
            a21 = A(1, 2); a22 = A(2, 2); a23 = A(3, 2)
            a31 = A(1, 3); a32 = A(2, 3); a33 = A(3, 3)

            x1 = a11*(a22*a33 - a23*a32)
            x2 = a12*(a21*a33 - a23*a31)
            x3 = a13*(a21*a32 - a22*a31)
            detx = x1 - x2 + x3

!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)

            a11 = A(1, 1); a12 = A(2, 1); a13 = A(3, 1)
            a21 = dq(1); a22 = dq(2); a23 = dq(3)
            a31 = A(1, 3); a32 = A(2, 3); a33 = A(3, 3)

            x1 = a11*(a22*a33 - a23*a32)
            x2 = a12*(a21*a33 - a23*a31)
            x3 = a13*(a21*a32 - a22*a31)
            dety = x1 - x2 + x3

!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)

            a11 = A(1, 1); a12 = A(2, 1); a13 = A(3, 1)
            a21 = A(1, 2); a22 = A(2, 2); a23 = A(3, 2)
            a31 = dq(1); a32 = dq(2); a33 = dq(3)

            x1 = a11*(a22*a33 - a23*a32)
            x2 = a12*(a21*a33 - a23*a31)
            x3 = a13*(a21*a32 - a22*a31)
            detz = x1 - x2 + x3

            gradtri(ij + z_down, l, 1) = detx
            gradtri(ij + z_down, l, 2) = dety
            gradtri(ij + z_down, l, 3) = detz
            arr(ij + z_down) = det

         END DO
      END DO

!DIR$ SIMD
      !$acc parallel loop async
      DO ij = ij_begin, ij_end
     ar(ij) = arr(ij + z_up) + arr(ij + z_lup) + arr(ij + z_ldown) + arr(ij + z_down) + arr(ij + z_rdown) + arr(ij + z_rup) + 1.e-50
      END DO
      CALL trace_end("compute_gradq3d1")
      CALL trace_start2("compute_gradq3d2")

      !$acc parallel loop collapse(3) async
      DO k = 1, 3
         DO l = ll_begin, ll_end
!DIR$ SIMD
            DO ij = ij_begin, ij_end
               gradq3d(ij, l, k) = (gradtri(ij + z_up, l, k) + gradtri(ij + z_down, l, k) + &
                                    gradtri(ij + z_rup, l, k) + gradtri(ij + z_ldown, l, k) + &
                                    gradtri(ij + z_lup, l, k) + gradtri(ij + z_rdown, l, k))/ar(ij)
            END DO
         END DO
      END DO
      CALL trace_end2("compute_gradq3d2")
      CALL trace_start("compute_gradq3d3")

      !============================================================================================= LIMITING
      !$acc parallel loop collapse(2) async
      DO l = ll_begin, ll_end
!DIR$ SIMD
         DO ij = ij_begin, ij_end
!             maggrd =  dot_product_3d(gradq3d(ij,l,:),gradq3d(ij,l,:))
            maggrd = gradq3d(ij, l, 1)*gradq3d(ij, l, 1) + gradq3d(ij, l, 2)*gradq3d(ij, l, 2) + gradq3d(ij, l, 3)*gradq3d(ij, l, 3)
            maggrd = sqrt(maggrd)
            maxq_c = qi(ij, l) + maggrd*sqrt_leng(ij)
            minq_c = qi(ij, l) - maggrd*sqrt_leng(ij)
            maxq = max(qi(ij, l), qi(ij + t_right, l), qi(ij + t_lup, l), qi(ij + t_rup, l), qi(ij + t_left, l), &
                       qi(ij + t_rdown, l), qi(ij + t_ldown, l))
            minq = min(qi(ij, l), qi(ij + t_right, l), qi(ij + t_lup, l), qi(ij + t_rup, l), qi(ij + t_left, l), &
                       qi(ij + t_rdown, l), qi(ij + t_ldown, l))
            IF ((maxq_c - qi(ij, l)) /= 0.0) THEN
               alphamx = (maxq - qi(ij, l)); alphamx = alphamx/(maxq_c - qi(ij, l))
               alphamx = max(alphamx, 0.0)
            ELSE
               alphamx = 0.0
            END IF
            IF ((minq_c - qi(ij, l)) /= 0.0) THEN
               alphami = (minq - qi(ij, l)); alphami = alphami/(minq_c - qi(ij, l))
               alphami = max(alphami, 0.0)
            ELSE
               alphami = 0.0
            END IF
            alpha = min(alphamx, alphami, 1.0)
!             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
            gradq3d(ij, l, 1) = alpha*gradq3d(ij, l, 1)
            gradq3d(ij, l, 2) = alpha*gradq3d(ij, l, 2)
            gradq3d(ij, l, 3) = alpha*gradq3d(ij, l, 3)
         END DO
      END DO

      CALL trace_end("compute_gradq3d3")

      !$acc end data

   CONTAINS

      SUBROUTINE gradq(n0, l, n1, n2, n3, q, dq1, dq2, dq3, det)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: n0, l, n1, n2, n3
         REAL(rstd), INTENT(IN)     :: q(iim*jjm, llm)
!    REAL(rstd), INTENT(OUT)    :: dq(3), det
         REAL(rstd), INTENT(OUT)    :: dq1, dq2, dq3, det
         REAL(rstd)    :: dq(3)

         REAL(rstd) :: A(3, 3)

         ! TODO : replace A by A1,A2,A3
         A(1, 1) = xyz_i(n1, 1) - xyz_i(n0, 1); A(1, 2) = xyz_i(n1, 2) - xyz_i(n0, 2); A(1, 3) = xyz_i(n1, 3) - xyz_i(n0, 3)
         A(2, 1) = xyz_i(n2, 1) - xyz_i(n0, 1); A(2, 2) = xyz_i(n2, 2) - xyz_i(n0, 2); A(2, 3) = xyz_i(n2, 3) - xyz_i(n0, 3)
         A(3, 1) = xyz_v(n3, 1); A(3, 2) = xyz_v(n3, 2); A(3, 3) = xyz_v(n3, 3)

         dq(1) = q(n1, l) - q(n0, l)
         dq(2) = q(n2, l) - q(n0, l)
         dq(3) = 0.0

         !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
!    CALL determinant(A(:,1),A(:,2),A(:,3),det)
!    CALL determinant(dq,A(:,2),A(:,3),detx)
!    CALL determinant(A(:,1),dq,A(:,3),dety)
!    CALL determinant(A(:,1),A(:,2),dq,detz)
!    dq(1) = detx
!    dq(2) = dety
!    dq(3) = detz

         CALL determinant(A(1, 1), A(2, 1), A(3, 1), A(1, 2), A(2, 2), A(3, 2), A(1, 3), A(2, 3), A(3, 3), det)
         CALL determinant(dq(1), dq(2), dq(3), A(1, 2), A(2, 2), A(3, 2), A(1, 3), A(2, 3), A(3, 3), dq1)
         CALL determinant(A(1, 1), A(2, 1), A(3, 1), dq(1), dq(2), dq(3), A(1, 3), A(2, 3), A(3, 3), dq2)
         CALL determinant(A(1, 1), A(2, 1), A(3, 1), A(1, 2), A(2, 2), A(3, 2), dq(1), dq(2), dq(3), dq3)

      END SUBROUTINE gradq

      !==========================================================================
!  PURE SUBROUTINE determinant(a1,a2,a3,det)
!    IMPLICIT NONE
!    REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3
!    REAL(rstd), INTENT(OUT) :: det
!    REAL(rstd) ::  x1,x2,x3
!    x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
!    x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
!    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
!    det =  x1 - x2 + x3
!  END SUBROUTINE determinant

      SUBROUTINE determinant(a11, a12, a13, a21, a22, a23, a31, a32, a33, det)
         IMPLICIT NONE
         REAL(rstd), INTENT(IN) :: a11, a12, a13, a21, a22, a23, a31, a32, a33
         REAL(rstd), INTENT(OUT) :: det
         REAL(rstd) ::  x1, x2, x3
         x1 = a11*(a22*a33 - a23*a32)
         x2 = a12*(a21*a33 - a23*a31)
         x3 = a13*(a21*a32 - a22*a31)
         det = x1 - x2 + x3
      END SUBROUTINE determinant

   END SUBROUTINE compute_gradq3d_hex

END SUBMODULE compute_gradq3d_mod
