MODULE wind_to_lonlat_mod
   USE omp_para
   USE icosa
   IMPLICIT NONE
   PRIVATE

   PUBLIC :: un2ulonlat, compute_wind_centered, compute_flux_centered, compute_wind_centered_lonlat, &
             compute_wind_on_edge, compute_un2ulonlat_unst, compute_un2ulonlat_hex

CONTAINS

   SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat, scale_)
      USE compute_diagnostics_mod, ONLY: compute_un2ulonlat
      TYPE(t_field), POINTER :: f_u(:) ! IN  : normal velocity components on edges
      TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
      REAL(rstd), POINTER    :: u(:, :), ulon(:, :), ulat(:, :)
      REAL(rstd), OPTIONAL   :: scale_
      REAL(rstd)             :: scale
      INTEGER :: ind
      scale = 1.
      IF (present(scale_)) scale = scale_
      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
         u = f_u(ind)
         ulon = f_ulon(ind)
         ulat = f_ulat(ind)
         CALL compute_un2ulonlat(u, ulon, ulat, scale)
      END DO

   END SUBROUTINE un2ulonlat

   SUBROUTINE compute_wind_centered(ue, ucenter, scale_)
      USE grid_param
      REAL(rstd) :: ue(3*iim*jjm, llm)
      REAL(rstd) :: ucenter(iim*jjm, llm, 3)
      REAL(rstd), INTENT(IN), OPTIONAL :: scale_
      INTEGER :: ij, l
      REAL(rstd) :: scale, fac, ue_le, cx, cy, cz, ux, uy, uz
      scale = 1.
      IF (present(scale_)) scale = scale_
      !--------------------------------------------------------------------------
      !---------------------------- wind_centered ----------------------------------
      ! Perot reconstruction based on Gauss theorem
      ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij = ij_begin, ij_end
            cx = centroid(ij, 1)
            cy = centroid(ij, 2)
            cz = centroid(ij, 3)
            ux = 0.; uy = 0.; uz = 0.
            ue_le = ne_rup*ue(ij + u_rup, l)*le(ij + u_rup)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_rup, 1) + xyz_v(ij + z_up, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_rup, 2) + xyz_v(ij + z_up, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_rup, 3) + xyz_v(ij + z_up, 3)) - cz)
            ue_le = ne_lup*ue(ij + u_lup, l)*le(ij + u_lup)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_lup, 1) + xyz_v(ij + z_up, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_lup, 2) + xyz_v(ij + z_up, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_lup, 3) + xyz_v(ij + z_up, 3)) - cz)
            ue_le = ne_left*ue(ij + u_left, l)*le(ij + u_left)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_lup, 1) + xyz_v(ij + z_ldown, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_lup, 2) + xyz_v(ij + z_ldown, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_lup, 3) + xyz_v(ij + z_ldown, 3)) - cz)
            ue_le = ne_ldown*ue(ij + u_ldown, l)*le(ij + u_ldown)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_ldown, 1) + xyz_v(ij + z_down, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_ldown, 2) + xyz_v(ij + z_down, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_ldown, 3) + xyz_v(ij + z_down, 3)) - cz)
            ue_le = ne_rdown*ue(ij + u_rdown, l)*le(ij + u_rdown)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_rdown, 1) + xyz_v(ij + z_down, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_rdown, 2) + xyz_v(ij + z_down, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_rdown, 3) + xyz_v(ij + z_down, 3)) - cz)
            ue_le = ne_right*ue(ij + u_right, l)*le(ij + u_right)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_rup, 1) + xyz_v(ij + z_rdown, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_rup, 2) + xyz_v(ij + z_rdown, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_rup, 3) + xyz_v(ij + z_rdown, 3)) - cz)
            fac = scale*(1./Ai(ij))
            ucenter(ij, l, 1) = ux*fac
            ucenter(ij, l, 2) = uy*fac
            ucenter(ij, l, 3) = uz*fac
         END DO
      END DO
      !---------------------------- wind_centered ----------------------------------
      !--------------------------------------------------------------------------
   END SUBROUTINE compute_wind_centered

   SUBROUTINE compute_flux_centered(scale, ue, ucenter)
      USE grid_param
      REAL(rstd), INTENT(IN) :: scale
      REAL(rstd) :: ue(3*iim*jjm, llm)
      REAL(rstd) :: ucenter(iim*jjm, llm, 3)
      INTEGER :: ij, l
      REAL(rstd) :: fac, ue_le, cx, cy, cz, ux, uy, uz

      !--------------------------------------------------------------------------
      !---------------------------- flux_centered ----------------------------------
      ! NB : here the input data is a flux and has already the factor l_e in it
      ! Input data is rescaled by factor scale
      ! Perot reconstruction based on Gauss theorem
      ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij = ij_begin, ij_end
            cx = centroid(ij, 1)
            cy = centroid(ij, 2)
            cz = centroid(ij, 3)
            ux = 0.; uy = 0.; uz = 0.
            ue_le = ne_rup*ue(ij + u_rup, l)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_rup, 1) + xyz_v(ij + z_up, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_rup, 2) + xyz_v(ij + z_up, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_rup, 3) + xyz_v(ij + z_up, 3)) - cz)
            ue_le = ne_lup*ue(ij + u_lup, l)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_lup, 1) + xyz_v(ij + z_up, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_lup, 2) + xyz_v(ij + z_up, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_lup, 3) + xyz_v(ij + z_up, 3)) - cz)
            ue_le = ne_left*ue(ij + u_left, l)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_lup, 1) + xyz_v(ij + z_ldown, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_lup, 2) + xyz_v(ij + z_ldown, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_lup, 3) + xyz_v(ij + z_ldown, 3)) - cz)
            ue_le = ne_ldown*ue(ij + u_ldown, l)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_ldown, 1) + xyz_v(ij + z_down, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_ldown, 2) + xyz_v(ij + z_down, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_ldown, 3) + xyz_v(ij + z_down, 3)) - cz)
            ue_le = ne_rdown*ue(ij + u_rdown, l)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_rdown, 1) + xyz_v(ij + z_down, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_rdown, 2) + xyz_v(ij + z_down, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_rdown, 3) + xyz_v(ij + z_down, 3)) - cz)
            ue_le = ne_right*ue(ij + u_right, l)
            ux = ux + ue_le*(.5*(xyz_v(ij + z_rup, 1) + xyz_v(ij + z_rdown, 1)) - cx)
            uy = uy + ue_le*(.5*(xyz_v(ij + z_rup, 2) + xyz_v(ij + z_rdown, 2)) - cy)
            uz = uz + ue_le*(.5*(xyz_v(ij + z_rup, 3) + xyz_v(ij + z_rdown, 3)) - cz)
            fac = scale*(1./Ai(ij))
            ucenter(ij, l, 1) = ux*fac
            ucenter(ij, l, 2) = uy*fac
            ucenter(ij, l, 3) = uz*fac
         END DO
      END DO
      !---------------------------- flux_centered ----------------------------------
      !--------------------------------------------------------------------------
   END SUBROUTINE compute_flux_centered

   SUBROUTINE compute_wind_on_edge(ue, uedge)
      USE grid_param
      REAL(rstd) :: ue(3*iim*jjm, llm)
      REAL(rstd) :: uedge(3*iim*jjm, llm, 3)

      REAL(rstd) :: ut(3*iim*jjm, llm)
      INTEGER :: i, j, ij, l

      CALL compute_tangential_compound(ue, ut)

      DO l = ll_begin, ll_end
         DO j = jj_begin, jj_end
            DO i = ii_begin, ii_end
               ij = (j - 1)*iim + i
               uedge(ij + u_right, l, :) = ue(ij + u_right, l)*ep_e(ij + u_right, :)*ne(ij, right) &
                                           + ut(ij + u_right, l)*et_e(ij + u_right, :)*ne(ij, right)
               uedge(ij + u_lup, l, :) = ue(ij + u_lup, l)*ep_e(ij + u_lup, :)*ne(ij, lup) &
                                         + ut(ij + u_lup, l)*et_e(ij + u_lup, :)*ne(ij, lup)
               uedge(ij + u_ldown, l, :) = ue(ij + u_ldown, l)*ep_e(ij + u_ldown, :)*ne(ij, ldown) &
                                           + ut(ij + u_ldown, l)*et_e(ij + u_ldown, :)*ne(ij, ldown)
            END DO
         END DO
      END DO

   END SUBROUTINE compute_wind_on_edge

   SUBROUTINE compute_tangential_compound(ue, ut)
      USE grid_param
      REAL(rstd) :: ue(3*iim*jjm, llm)
      REAL(rstd) :: ut(3*iim*jjm, llm)
      INTEGER :: i, j, l, ij

      DO l = ll_begin, ll_end
         DO j = jj_begin, jj_end
            DO i = ii_begin, ii_end
               ij = (j - 1)*iim + i

               ut(ij + u_right, l) = 1/de(ij + u_right)* &
                                     (wee(ij + u_right, 1)*ue(ij + u_rup, l) + &
                                      wee(ij + u_right, 2)*ue(ij + u_lup, l) + &
                                      wee(ij + u_right, 3)*ue(ij + u_left, l) + &
                                      wee(ij + u_right, 4)*ue(ij + u_ldown, l) + &
                                      wee(ij + u_right, 5)*ue(ij + u_rdown, l) + &
                                      wee(ij + u_right, 6)*ue(ij + t_right + u_ldown, l) + &
                                      wee(ij + u_right, 7)*ue(ij + t_right + u_rdown, l) + &
                                      wee(ij + u_right, 8)*ue(ij + t_right + u_right, l) + &
                                      wee(ij + u_right, 9)*ue(ij + t_right + u_rup, l) + &
                                      wee(ij + u_right, 10)*ue(ij + t_right + u_lup, l))

               ut(ij + u_lup, l) = 1/de(ij + u_lup)* &
                                   (wee(ij + u_lup, 1)*ue(ij + u_left, l) + &
                                    wee(ij + u_lup, 2)*ue(ij + u_ldown, l) + &
                                    wee(ij + u_lup, 3)*ue(ij + u_rdown, l) + &
                                    wee(ij + u_lup, 4)*ue(ij + u_right, l) + &
                                    wee(ij + u_lup, 5)*ue(ij + u_rup, l) + &
                                    wee(ij + u_lup, 6)*ue(ij + t_lup + u_right, l) + &
                                    wee(ij + u_lup, 7)*ue(ij + t_lup + u_rup, l) + &
                                    wee(ij + u_lup, 8)*ue(ij + t_lup + u_lup, l) + &
                                    wee(ij + u_lup, 9)*ue(ij + t_lup + u_left, l) + &
                                    wee(ij + u_lup, 10)*ue(ij + t_lup + u_ldown, l))

               ut(ij + u_ldown, l) = 1/de(ij + u_ldown)* &
                                     (wee(ij + u_ldown, 1)*ue(ij + u_rdown, l) + &
                                      wee(ij + u_ldown, 2)*ue(ij + u_right, l) + &
                                      wee(ij + u_ldown, 3)*ue(ij + u_rup, l) + &
                                      wee(ij + u_ldown, 4)*ue(ij + u_lup, l) + &
                                      wee(ij + u_ldown, 5)*ue(ij + u_left, l) + &
                                      wee(ij + u_ldown, 6)*ue(ij + t_ldown + u_lup, l) + &
                                      wee(ij + u_ldown, 7)*ue(ij + t_ldown + u_left, l) + &
                                      wee(ij + u_ldown, 8)*ue(ij + t_ldown + u_ldown, l) + &
                                      wee(ij + u_ldown, 9)*ue(ij + t_ldown + u_rdown, l) + &
                                      wee(ij + u_ldown, 10)*ue(ij + t_ldown + u_right, l))

            END DO
         END DO
      END DO

   END SUBROUTINE compute_tangential_compound

   SUBROUTINE compute_wind_centered_lonlat(uc, ulon, ulat)
      USE grid_param
      REAL(rstd) :: uc(iim*jjm, llm, 3)
      REAL(rstd) :: ulon(iim*jjm, llm)
      REAL(rstd) :: ulat(iim*jjm, llm)
      INTEGER :: i, j, ij, l

      DO l = ll_begin, ll_end
         DO j = jj_begin, jj_end
            DO i = ii_begin, ii_end
               ij = (j - 1)*iim + i
               ulon(ij, l) = sum(uc(ij, l, :)*elon_i(ij, :))
               ulat(ij, l) = sum(uc(ij, l, :)*elat_i(ij, :))
            END DO
         END DO
      END DO
   END SUBROUTINE compute_wind_centered_lonlat

!-------------- Wrappers for F2008 conformity -----------------

   SUBROUTINE compute_un2ulonlat_unst(un, ulon, ulat, scale)
      REAL(rstd), INTENT(IN)  :: un(:, :)
      REAL(rstd), INTENT(OUT) :: ulon(:, :), ulat(:, :)
      REAL(rstd), INTENT(IN)  :: scale
      CALL compute_un2ulonlat_unst_(un, ulon, ulat, scale)
   END SUBROUTINE compute_un2ulonlat_unst

   SUBROUTINE compute_un2ulonlat_hex(un, ulon, ulat, scale)
      REAL(rstd), INTENT(IN)  :: un(:, :)
      REAL(rstd), INTENT(OUT) :: ulon(:, :), ulat(:, :)
      REAL(rstd), INTENT(IN)  :: scale
      CALL compute_un2ulonlat_hex_(un, ulon, ulat, scale)
   END SUBROUTINE compute_un2ulonlat_hex

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

   SUBROUTINE compute_un2ulonlat_unst_(un, ulon, ulat, scale)
      USE grid_param, ONLY: primal_num, llm
      USE data_unstructured_mod, ONLY: primal_deg, primal_edge
      REAL(C_DOUBLE), DIMENSION(llm, edge_num), INTENT(IN)  :: un
      REAL(C_DOUBLE), DIMENSION(llm, primal_num), INTENT(OUT) :: ulon, ulat
      REAL(rstd)             :: scale
      REAL(rstd) :: u3d(3), perot(SIZE(perot_ie, 1), 3), coslon, coslat, sinlon, sinlat
      INTEGER :: cell, deg, l, iedge, edge
      DO cell = 1, primal_num
         deg = primal_deg(cell)
         coslon = COS(lon_i(cell))*scale
         coslat = COS(lat_i(cell))*scale
         sinlon = SIN(lon_i(cell))*scale
         sinlat = SIN(lat_i(cell))
         perot(1:deg, 1:3) = perot_ie(1:deg, cell, 1:3)
         DO l = 1, llm
            u3d(:) = 0.
            DO iedge = 1, deg
               edge = primal_edge(iedge, cell)
               ! The coefficients grad_ie are computed assuming
               ! that un is contravariant (= normal component * cell edge length).
               ! Input un is the normal component of the vector field.
               ! => we multiply input un by cell edge length le(edge)
               u3d(:) = u3d(:) + (un(l, edge)*le(edge))*perot(iedge, :)
            END DO
            ! At this point we have the 3D components of the Perot-reconstructed vector
            ! We project it onto the local basis :
            !    elon = (-sinlon, coslon, 0),
            !    elat = (-coslon*sinlat, -sinlon*sinlat, coslat)
            ulon(l, cell) = coslon*u3d(2) - sinlon*u3d(1)
            ulat(l, cell) = coslat*u3d(3) - sinlat*(coslon*u3d(1) + sinlon*u3d(2))
         END DO
      END DO
   END SUBROUTINE compute_un2ulonlat_unst_

   SUBROUTINE compute_un2ulonlat_hex_(un, ulon, ulat, scale)
      USE grid_param
      REAL(rstd), INTENT(IN)  :: un(3*iim*jjm, llm)
      REAL(rstd), INTENT(OUT) :: ulon(iim*jjm, llm)
      REAL(rstd), INTENT(OUT) :: ulat(iim*jjm, llm)
      REAL(rstd)             :: uc(iim*jjm, llm, 3)
      REAL(rstd)             :: scale
      CALL compute_wind_centered(un, uc, scale)
      CALL compute_wind_centered_lonlat(uc, ulon, ulat)
   END SUBROUTINE compute_un2ulonlat_hex_
END MODULE wind_to_lonlat_mod
