MODULE vertical_interp_mod
   USE prec, ONLY: rstd
   IMPLICIT NONE
   PRIVATE
   SAVE

   PUBLIC init_vertical_interp, vertical_interp, vertical_interp_4d

CONTAINS

   SUBROUTINE init_vertical_interp
   END SUBROUTINE init_vertical_interp

   SUBROUTINE vertical_interp(f_pmid, f_in, f_out, pval)
      USE icosa
      TYPE(t_field), POINTER :: f_pmid(:)
      TYPE(t_field), POINTER :: f_in(:)
      TYPE(t_field), POINTER :: f_out(:)
      REAL(rstd), INTENT(IN) :: pval

      REAL(rstd), POINTER :: in(:, :)
      REAL(rstd), POINTER :: out(:)
      REAL(rstd), POINTER :: pmid(:, :)

      INTEGER :: ind

      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
         pmid = f_pmid(ind)
         in = f_in(ind)
         out = f_out(ind)
         CALL compute_vertical_interp(pmid, in, out, pval)
      END DO

   END SUBROUTINE vertical_interp

   SUBROUTINE vertical_interp_4d(f_pmid, f_in, iq, f_out, pval)
      USE icosa
      TYPE(t_field), POINTER :: f_pmid(:)
      TYPE(t_field), POINTER :: f_in(:)
      TYPE(t_field), POINTER :: f_out(:)
      INTEGER,    INTENT(IN) :: iq
      REAL(rstd), INTENT(IN) :: pval

      REAL(rstd), POINTER :: in(:, :, :)
      REAL(rstd), POINTER :: out(:)
      REAL(rstd), POINTER :: pmid(:, :)

      INTEGER :: ind

      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
         pmid = f_pmid(ind)
         in = f_in(ind)
         out = f_out(ind)
         CALL compute_vertical_interp(pmid, in(:, :, iq), out, pval)
      END DO

   END SUBROUTINE vertical_interp_4d

   SUBROUTINE compute_vertical_interp(pmid, in, out, pval)
      USE grid_param, ONLY: grid_type, grid_ico, grid_unst
      REAL(rstd), INTENT(IN) :: pmid(:, :), in(:, :), pval
      REAL(rstd), INTENT(OUT):: out(:)
      IF (grid_type == grid_ico) CALL compute_vertical_interp_hex_(pmid, in, out, pval)
      IF (grid_type == grid_unst) CALL compute_vertical_interp_unst_(pmid, in, out, pval)
   END SUBROUTINE compute_vertical_interp

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

   SUBROUTINE compute_vertical_interp_hex_(pmid, in, out, pval)
      USE icosa
      USE omp_para, ONLY: is_omp_level_master
      REAL(rstd), INTENT(IN) :: pmid(iim*jjm, llm)
      REAL(rstd), INTENT(IN) :: in(iim*jjm, llm)
      REAL(rstd), INTENT(OUT):: out(iim*jjm)
      REAL(rstd), INTENT(IN) :: pval
      INTEGER :: i, j, ij

!$OMP BARRIER
      IF (is_omp_level_master) THEN

         DO j = jj_begin - 1, jj_end + 1
            DO i = ii_begin - 1, ii_end + 1
               ij = (j - 1)*iim + i
               CALL interp_1d(pmid(ij, :), in(ij, :), out(ij), pval)
            END DO
         END DO
      END IF
!$OMP BARRIER

   END SUBROUTINE compute_vertical_interp_hex_

   SUBROUTINE compute_vertical_interp_unst_(pmid, in, out, pval)
      USE icosa
      REAL(rstd), INTENT(IN) :: pmid(llm, primal_num)
      REAL(rstd), INTENT(IN) :: in(llm, primal_num)
      REAL(rstd), INTENT(OUT) :: out(primal_num)
      REAL(rstd) :: pval
      INTEGER :: ij
!$OMP BARRIER
!$OMP MASTER
      DO ij = 1, primal_num
         CALL interp_1d(pmid(:, ij), in(:, ij), out(ij), pval)
      END DO
!$OMP END MASTER
!$OMP BARRIER
   END SUBROUTINE compute_vertical_interp_unst_

   PURE SUBROUTINE interp_1d(pmid, in, out, pval)
      USE icosa
      REAL(rstd), INTENT(IN)  :: pmid(:)
      REAL(rstd), INTENT(IN)  :: in(:)
      REAL(rstd), INTENT(OUT) :: out
      REAL(rstd), INTENT(IN)  :: pval
      REAL(rstd) :: coeff
      INTEGER :: l
      l = llm - 1
      DO WHILE (pmid(l) < pval .AND. l > 1)
         l = l - 1
      END DO
      coeff = (pval - pmid(l))/(pmid(l) - pmid(l + 1))
      out = in(l) + coeff*(in(l) - in(l + 1))
   END SUBROUTINE interp_1d

END MODULE vertical_interp_mod
