MODULE transport_unstructured_mod
   USE ISO_C_BINDING
   USE data_unstructured_mod
   IMPLICIT NONE
   PRIVATE
   SAVE

CONTAINS

   SUBROUTINE remap_unstructured(rhodz, thetarhodz, u) BIND(C, name='dynamico_remap')
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)   :: rhodz, thetarhodz    ! INOUT
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)      :: u                    ! INOUT
      ! DIAG
      REAL(C_DOUBLE), DIMENSION(primal_num)     :: new_rhodz_cum, mass_col
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)   :: old_rhodz
      REAL(C_DOUBLE), DIMENSION(llm + 1, primal_num) :: eta, rhodz_cum, thetarhodz_cum, new_thetarhodz_cum
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)      :: urhodz
      REAL(C_DOUBLE), DIMENSION(llm + 1, edge_num)     :: urhodz_cum, new_urhodz_cum
      INTEGER, DIMENSION(primal_num) :: cur_lev

      REAL(C_DOUBLE)          :: rhodz_cum_target, rhodz_cum_lev, rhodz_cum_levp1, X
      REAL(C_DOUBLE)         :: time1, time2
      INTEGER :: ij, l, level, edge, ij_left, ij_right

      time1 = OMP_GET_WTIME()

      !--------------------------------------------------------------------------
      !---------------------------- remap_eta ----------------------------------
      ! IN : rhodz, mass_dak, mass_dbk
      ! TMP : mass_col, cur_lev, new_rhodz_cum
      ! OUT : rhodz, old_rhodz, eta
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         rhodz_cum(1, ij) = 0.
         cur_lev(ij) = 1
         eta(1, ij) = 1.
         new_rhodz_cum(ij) = 0.
         DO l = 1, llm
            rhodz_cum(l + 1, ij) = rhodz_cum(l, ij) + rhodz(l, ij)
         END DO
         mass_col(ij) = rhodz_cum(llm + 1, ij)
         DO l = 1, llm
            old_rhodz(l, ij) = rhodz(l, ij)
            rhodz(l, ij) = mass_dak(l, ij) + mass_col(ij)*mass_dbk(l, ij)
            rhodz_cum_target = new_rhodz_cum(ij) + rhodz(l, ij)
            DO level = cur_lev(ij), llm
               rhodz_cum_levp1 = rhodz_cum(level + 1, ij)
               IF (rhodz_cum_target <= rhodz_cum_levp1) EXIT
            END DO
            IF (level > llm) level = llm
            rhodz_cum_lev = rhodz_cum(level, ij)
            ! now rhodz_cum_lev <= rhodz_cum_target <= rhodz_cum_levp1
            cur_lev(ij) = level
            new_rhodz_cum(ij) = rhodz_cum_target
            eta(l + 1, ij) = level + (rhodz_cum_target - rhodz_cum_lev)/(rhodz_cum_levp1 - rhodz_cum_lev)
         END DO
      END DO
      !$OMP END DO
      !---------------------------- remap_eta ----------------------------------
      !--------------------------------------------------------------------------

      !--------------------------------------------------------------------------
      !---------------------------- remap_theta ----------------------------------
      ! IN : thetarhodz, eta
      ! TMP : thetarhodz_cum, new_thetarhodz_cum
      ! OUT : thetarhodz
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         thetarhodz_cum(1, ij) = 0.
         DO l = 1, llm
            thetarhodz_cum(l + 1, ij) = thetarhodz_cum(l, ij) + thetarhodz(l, ij)
         END DO
         DO l = 1, llm + 1
            X = eta(l, ij)
            level = MIN(llm, FLOOR(X)) ! eta=llm+1 => level=llm, X=1
            X = X - level
            new_thetarhodz_cum(l, ij) = thetarhodz_cum(level, ij) + X*thetarhodz(level, ij)
         END DO
         DO l = 1, llm
            thetarhodz(l, ij) = new_thetarhodz_cum(l + 1, ij) - new_thetarhodz_cum(l, ij)
         END DO
      END DO
      !$OMP END DO
      !---------------------------- remap_theta ----------------------------------
      !--------------------------------------------------------------------------

      !--------------------------------------------------------------------------
      !---------------------------- remap_u ----------------------------------
      ! IN : u, old_rhodz, rhodz, eta
      ! TMP : urhodz_cum, new_urhodz_cum
      ! OUT : u
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         ij_left = left(edge)
         ij_right = right(edge)
         urhodz_cum(1, edge) = 0.
         DO l = 1, llm
            urhodz(l, edge) = u(l, edge)*(old_rhodz(l, ij_left) + old_rhodz(l, ij_right))
            urhodz_cum(l + 1, edge) = urhodz_cum(l, edge) + urhodz(l, edge)
         END DO
         DO l = 1, llm + 1
            X = .5*(eta(l, ij_left) + eta(l, ij_right))
            level = MIN(llm, FLOOR(X))
            X = X - level
            new_urhodz_cum(l, edge) = urhodz_cum(level, edge) + X*urhodz(level, edge)
         END DO
         DO l = 1, llm
            u(l, edge) = (new_urhodz_cum(l + 1, edge) - new_urhodz_cum(l, edge))/(rhodz(l, ij_left) + rhodz(l, ij_right))
         END DO
      END DO
      !$OMP END DO
      !---------------------------- remap_u ----------------------------------
      !--------------------------------------------------------------------------

      time2 = OMP_GET_WTIME()
      IF (time2 > time1) elapsed = elapsed + time2 - time1

      CALL print_trace()
   END SUBROUTINE remap_unstructured

END MODULE transport_unstructured_mod
