SUBMODULE(compute_transport_mod) compute_backward_traj_mod

   IMPLICIT NONE

CONTAINS

   MODULE SUBROUTINE compute_backward_traj(normal, tangent, ue, tau, cc, &
   &                                       xyz_e, de, wee, le) ! metrics terms
      USE grid_param, ONLY: grid_type, grid_ico

      REAL(rstd), INTENT(IN)    :: normal(:, :)
      REAL(rstd), INTENT(IN)    :: tangent(:, :)
      REAL(rstd), INTENT(IN)    :: ue(:, :)
      REAL(rstd), INTENT(OUT)   :: cc(:, :, :) ! start of backward trajectory
      REAL(rstd), INTENT(IN)    :: tau
      ! metrics terms
      REAL(rstd), INTENT(IN)    :: xyz_e(:, :)
      REAL(rstd), INTENT(IN)    :: de(:)
      REAL(rstd), INTENT(IN)    :: wee(:, :)
      REAL(rstd), INTENT(IN)    :: le(:)

      IF (grid_type == grid_ico) THEN
         CALL compute_backward_traj_hex(normal, tangent, ue, tau, cc, &
      &                              xyz_e, de, wee, le)
      END IF
   END SUBROUTINE

! Backward trajectories, for use with Miura approach
   SUBROUTINE compute_backward_traj_hex(normal, tangent, ue, tau, cc, &
                                        xyz_e, de, wee, le) ! metrics terms
      USE trace
      USE omp_para
      USE dimensions
      USE grid_param, ONLY: llm
      IMPLICIT NONE
      REAL(rstd), INTENT(IN)    :: normal(3*iim*jjm, 3)
      REAL(rstd), INTENT(IN)    :: tangent(3*iim*jjm, 3)
      REAL(rstd), INTENT(IN)    :: ue(iim*3*jjm, llm)
      REAL(rstd), INTENT(OUT)   :: cc(3*iim*jjm, llm, 3) ! start of backward trajectory
      REAL(rstd), INTENT(IN)    :: tau
! metrics terms
      REAL(rstd), INTENT(IN)    :: xyz_e(iim*3*jjm, 3)
      REAL(rstd), INTENT(IN)    :: de(iim*3*jjm)
      REAL(rstd), INTENT(IN)    :: wee(iim*3*jjm, 5, 2)
      REAL(rstd), INTENT(IN)    :: le(iim*3*jjm)

      REAL(rstd) :: v_e(3), up_e
      INTEGER :: ij, l

      CALL trace_start("compute_backward_traj")

      ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement

      !$acc data present(ue(:,:), cc(:,:,:), normal(:,:), tangent(:,:), xyz_e(:,:), de(:), wee(:,:,:), le(:)) async

      ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
      !$acc parallel loop private(up_e, v_e) collapse(2) gang vector async
      DO l = ll_begin, ll_end
!DIR$ IVDEP
         DO ij = ij_begin, ij_end
            up_e = 1/de(ij + u_right)*( &
                   wee(ij + u_right, 1, 1)*le(ij + u_rup)*ue(ij + u_rup, l) + &
                   wee(ij + u_right, 2, 1)*le(ij + u_lup)*ue(ij + u_lup, l) + &
                   wee(ij + u_right, 3, 1)*le(ij + u_left)*ue(ij + u_left, l) + &
                   wee(ij + u_right, 4, 1)*le(ij + u_ldown)*ue(ij + u_ldown, l) + &
                   wee(ij + u_right, 5, 1)*le(ij + u_rdown)*ue(ij + u_rdown, l) + &
                   wee(ij + u_right, 1, 2)*le(ij + t_right + u_ldown)*ue(ij + t_right + u_ldown, l) + &
                   wee(ij + u_right, 2, 2)*le(ij + t_right + u_rdown)*ue(ij + t_right + u_rdown, l) + &
                   wee(ij + u_right, 3, 2)*le(ij + t_right + u_right)*ue(ij + t_right + u_right, l) + &
                   wee(ij + u_right, 4, 2)*le(ij + t_right + u_rup)*ue(ij + t_right + u_rup, l) + &
                   wee(ij + u_right, 5, 2)*le(ij + t_right + u_lup)*ue(ij + t_right + u_lup, l) &
                   )
            v_e = ue(ij + u_right, l)*normal(ij + u_right, :) + up_e*tangent(ij + u_right, :)
            cc(ij + u_right, l, :) = xyz_e(ij + u_right, :) - v_e*tau

            up_e = 1/de(ij + u_lup)*( &
                   wee(ij + u_lup, 1, 1)*le(ij + u_left)*ue(ij + u_left, l) + &
                   wee(ij + u_lup, 2, 1)*le(ij + u_ldown)*ue(ij + u_ldown, l) + &
                   wee(ij + u_lup, 3, 1)*le(ij + u_rdown)*ue(ij + u_rdown, l) + &
                   wee(ij + u_lup, 4, 1)*le(ij + u_right)*ue(ij + u_right, l) + &
                   wee(ij + u_lup, 5, 1)*le(ij + u_rup)*ue(ij + u_rup, l) + &
                   wee(ij + u_lup, 1, 2)*le(ij + t_lup + u_right)*ue(ij + t_lup + u_right, l) + &
                   wee(ij + u_lup, 2, 2)*le(ij + t_lup + u_rup)*ue(ij + t_lup + u_rup, l) + &
                   wee(ij + u_lup, 3, 2)*le(ij + t_lup + u_lup)*ue(ij + t_lup + u_lup, l) + &
                   wee(ij + u_lup, 4, 2)*le(ij + t_lup + u_left)*ue(ij + t_lup + u_left, l) + &
                   wee(ij + u_lup, 5, 2)*le(ij + t_lup + u_ldown)*ue(ij + t_lup + u_ldown, l) &
                   )
            v_e = ue(ij + u_lup, l)*normal(ij + u_lup, :) + up_e*tangent(ij + u_lup, :)
            cc(ij + u_lup, l, :) = xyz_e(ij + u_lup, :) - v_e*tau

            up_e = 1/de(ij + u_ldown)*( &
                   wee(ij + u_ldown, 1, 1)*le(ij + u_rdown)*ue(ij + u_rdown, l) + &
                   wee(ij + u_ldown, 2, 1)*le(ij + u_right)*ue(ij + u_right, l) + &
                   wee(ij + u_ldown, 3, 1)*le(ij + u_rup)*ue(ij + u_rup, l) + &
                   wee(ij + u_ldown, 4, 1)*le(ij + u_lup)*ue(ij + u_lup, l) + &
                   wee(ij + u_ldown, 5, 1)*le(ij + u_left)*ue(ij + u_left, l) + &
                   wee(ij + u_ldown, 1, 2)*le(ij + t_ldown + u_lup)*ue(ij + t_ldown + u_lup, l) + &
                   wee(ij + u_ldown, 2, 2)*le(ij + t_ldown + u_left)*ue(ij + t_ldown + u_left, l) + &
                   wee(ij + u_ldown, 3, 2)*le(ij + t_ldown + u_ldown)*ue(ij + t_ldown + u_ldown, l) + &
                   wee(ij + u_ldown, 4, 2)*le(ij + t_ldown + u_rdown)*ue(ij + t_ldown + u_rdown, l) + &
                   wee(ij + u_ldown, 5, 2)*le(ij + t_ldown + u_right)*ue(ij + t_ldown + u_right, l) &
                   )
            v_e = ue(ij + u_ldown, l)*normal(ij + u_ldown, :) + up_e*tangent(ij + u_ldown, :)
            cc(ij + u_ldown, l, :) = xyz_e(ij + u_ldown, :) - v_e*tau
         END DO
      END DO
      !$acc end data
      CALL trace_end("compute_backward_traj")

   END SUBROUTINE compute_backward_traj_hex

END SUBMODULE
