MODULE wind_from_lonlat_mod
   USE omp_para
   USE icosa
   IMPLICIT NONE
   PRIVATE

   PUBLIC :: ulonlat2un, compute_wind2d_perp_from_lonlat_compound, &
             compute_wind_centered_from_lonlat_compound, compute_wind_perp_from_lonlat_compound

CONTAINS

   SUBROUTINE ulonlat2un(f_ulon, f_ulat, f_u)
      TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! IN : velocity reconstructed at hexagons
      TYPE(t_field), POINTER :: f_u(:)               ! OUT  : normal velocity components on edges

      REAL(rstd), POINTER :: u(:, :), ulon(:, :), ulat(:, :)
      INTEGER :: ind

      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_ulonlat2un(ulon, ulat, u)
      END DO

   END SUBROUTINE ulonlat2un

   SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u)
      USE grid_param
      REAL(rstd) :: u(3*iim*jjm, llm, 3)
      REAL(rstd) :: ulon(3*iim*jjm, llm)
      REAL(rstd) :: ulat(3*iim*jjm, llm)

      INTEGER :: i, j, ij, l

      DO l = ll_begin, ll_end
         DO j = jj_begin - 1, jj_end + 1
            DO i = ii_begin - 1, ii_end + 1
               ij = (j - 1)*iim + i
               u(ij + u_right, l, :) = ulon(ij + u_right, l)*elon_e(ij + u_right, :) + ulat(ij + u_right, l)*elat_e(ij + u_right, :)
               u(ij + u_lup, l, :) = ulon(ij + u_lup, l)*elon_e(ij + u_lup, :) + ulat(ij + u_lup, l)*elat_e(ij + u_lup, :)
               u(ij + u_ldown, l, :) = ulon(ij + u_ldown, l)*elon_e(ij + u_ldown, :) + ulat(ij + u_ldown, l)*elat_e(ij + u_ldown, :)
            END DO
         END DO
      END DO

   END SUBROUTINE compute_wind_from_lonlat_compound

   SUBROUTINE compute_wind_centered_from_lonlat_compound(ulon, ulat, u)
      USE grid_param
      REAL(rstd) :: u(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 - 1, jj_end + 1
            DO i = ii_begin - 1, ii_end + 1
               ij = (j - 1)*iim + i
               u(ij, l, :) = ulon(ij, l)*elon_i(ij, :) + ulat(ij, l)*elat_i(ij, :)
            END DO
         END DO
      END DO
   END SUBROUTINE compute_wind_centered_from_lonlat_compound

   SUBROUTINE compute_wind2D_from_lonlat_compound(ulon, ulat, u)
      USE grid_param
      REAL(rstd) :: u(3*iim*jjm, 3)
      REAL(rstd) :: ulon(3*iim*jjm)
      REAL(rstd) :: ulat(3*iim*jjm)

      INTEGER :: i, j, ij

      DO j = jj_begin - 1, jj_end + 1
         DO i = ii_begin - 1, ii_end + 1
            ij = (j - 1)*iim + i
            u(ij + u_right, :) = ulon(ij + u_right)*elon_e(ij + u_right, :) + ulat(ij + u_right)*elat_e(ij + u_right, :)
            u(ij + u_lup, :) = ulon(ij + u_lup)*elon_e(ij + u_lup, :) + ulat(ij + u_lup)*elat_e(ij + u_lup, :)
            u(ij + u_ldown, :) = ulon(ij + u_ldown)*elon_e(ij + u_ldown, :) + ulat(ij + u_ldown)*elat_e(ij + u_ldown, :)
         END DO
      END DO

   END SUBROUTINE compute_wind2D_from_lonlat_compound

   SUBROUTINE compute_wind_perp_from_lonlat_compound(ulon, ulat, up)
      USE grid_param
      REAL(rstd) :: up(3*iim*jjm, llm)
      REAL(rstd) :: ulon(3*iim*jjm, llm)
      REAL(rstd) :: ulat(3*iim*jjm, llm)
      REAL(rstd) :: u(3*iim*jjm, llm, 3)

      INTEGER :: i, j, ij, l

      CALL compute_wind_from_lonlat_compound(ulon, ulat, u)

      DO l = ll_begin, ll_end
         DO j = jj_begin - 1, jj_end + 1
            DO i = ii_begin - 1, ii_end + 1
               ij = (j - 1)*iim + i
               up(ij + u_right, l) = sum(u(ij + u_right, l, :)*ep_e(ij + u_right, :))
               up(ij + u_lup, l) = sum(u(ij + u_lup, l, :)*ep_e(ij + u_lup, :))
               up(ij + u_ldown, l) = sum(u(ij + u_ldown, l, :)*ep_e(ij + u_ldown, :))
            END DO
         END DO
      END DO

   END SUBROUTINE compute_wind_perp_from_lonlat_compound

   SUBROUTINE compute_wind2D_perp_from_lonlat_compound(ulon, ulat, up)
      USE grid_param
      REAL(rstd) :: up(3*iim*jjm)
      REAL(rstd) :: ulon(3*iim*jjm)
      REAL(rstd) :: ulat(3*iim*jjm)
      REAL(rstd) :: u(3*iim*jjm, 3)

      INTEGER :: i, j, ij

      CALL compute_wind2D_from_lonlat_compound(ulon, ulat, u)
      DO j = jj_begin - 1, jj_end + 1
         DO i = ii_begin - 1, ii_end + 1
            ij = (j - 1)*iim + i
            up(ij + u_right) = sum(u(ij + u_right, :)*ep_e(ij + u_right, :))
            up(ij + u_lup) = sum(u(ij + u_lup, :)*ep_e(ij + u_lup, :))
            up(ij + u_ldown) = sum(u(ij + u_ldown, :)*ep_e(ij + u_ldown, :))
         END DO
      END DO

   END SUBROUTINE compute_wind2D_perp_from_lonlat_compound

   SUBROUTINE compute_wind_centered_from_wind_lonlat_centered(ulon, ulat, uc)
      USE grid_param
      REAL(rstd) :: ulon(iim*jjm, llm)
      REAL(rstd) :: ulat(iim*jjm, llm)
      REAL(rstd) :: uc(iim*jjm, llm, 3)

      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
               uc(ij, l, :) = ulon(ij, l)*elon_i(ij, :) + ulat(ij, l)*elat_i(ij, :)
            END DO
         END DO
      END DO

   END SUBROUTINE compute_wind_centered_from_wind_lonlat_centered

   SUBROUTINE compute_wind_perp_from_wind_centered(uc, un)
      USE grid_param

      IMPLICIT NONE
      REAL(rstd), INTENT(IN)   :: uc(iim*jjm, llm, 3)
      REAL(rstd), INTENT(OUT)  :: un(3*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
               un(ij + u_right, l) = sum(0.5*(uc(ij, l, :) + uc(ij + t_right, l, :))*ep_e(ij + u_right, :))
               un(ij + u_lup, l) = sum(0.5*(uc(ij, l, :) + uc(ij + t_lup, l, :))*ep_e(ij + u_lup, :))
               un(ij + u_ldown, l) = sum(0.5*(uc(ij, l, :) + uc(ij + t_ldown, l, :))*ep_e(ij + u_ldown, :))
            END DO
         END DO
      END DO

   END SUBROUTINE compute_wind_perp_from_wind_centered

   SUBROUTINE compute_ulonlat2un(ulon, ulat, un)
      USE grid_param
      REAL(rstd), INTENT(IN)  :: ulon(iim*jjm, llm)
      REAL(rstd), INTENT(IN)  :: ulat(iim*jjm, llm)
      REAL(rstd), INTENT(OUT) :: un(3*iim*jjm, llm)
      REAL(rstd)             :: uc(iim*jjm, llm, 3)

      CALL compute_wind_centered_from_wind_lonlat_centered(ulon, ulat, uc)
      CALL compute_wind_perp_from_wind_centered(uc, un)

   END SUBROUTINE compute_ulonlat2un

END MODULE wind_from_lonlat_mod
