MODULE geopotential_mod
   IMPLICIT NONE
   PRIVATE

   PUBLIC :: compute_geopotential
CONTAINS

   SUBROUTINE thetarhodz2geopot(f_ps, f_phis, f_theta_rhodz, f_p, f_theta, f_phi) ! FIXME : never called, dry only
      USE icosa
      USE omp_para
      USE pression_mod
      USE theta2theta_rhodz_mod
      TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
                                f_p(:), f_theta(:), f_phi(:)               ! OUT
      REAL(rstd), POINTER :: p(:, :), theta(:, :, :), theta_rhodz(:, :, :), &
                             phi(:, :), phis(:), ps(:)
      INTEGER :: ind

      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
         ps = f_ps(ind)
         p = f_p(ind)
!$OMP BARRIER
         CALL compute_pression(ps, p, 0)
!$OMP BARRIER
         theta_rhodz = f_theta_rhodz(ind)
         theta = f_theta(ind)
         CALL compute_theta_rhodz2theta(ps, theta_rhodz, theta, 0, ondevice=.FALSE.)
         phis = f_phis(ind)
         phi = f_phi(ind)
         CALL compute_geopotential(phis, p, theta, phi, 0)
      END DO

   END SUBROUTINE thetarhodz2geopot

   SUBROUTINE compute_geopotential(phis, p, theta, phi, offset)
      USE icosa
      USE omp_para
      REAL(rstd), INTENT(IN) :: p(iim*jjm, llm + 1)
      REAL(rstd), INTENT(IN) :: phis(iim*jjm)
      REAL(rstd), INTENT(IN) :: theta(iim*jjm, llm, nqdyn)
      REAL(rstd), INTENT(OUT) :: phi(iim*jjm, llm)
      INTEGER, INTENT(IN) :: offset
      INTEGER :: i, j, ij, l
      REAL(rstd) :: mg_ij, p_ij, exner_ij

!!! Compute geopotential
      ! LMD-Z : dh=pi.dtheta-vdp  => (1/rho) dp = pi.dtheta - d(h=pi.theta) = -dtheta.gradpi
      ! g.dz =  -(1/rho)dp =  theta.gradpi
      ! direct (variational)
      ! dz = -(1/rho.g)dp, rho = RT/p = (R/C_p).theta.pi/p

      ! for first layer

! flush pks, pk thetha, phis
!$OMP BARRIER
      IF (is_omp_level_master) THEN
         DO j = jj_begin - offset, jj_end + offset
            DO i = ii_begin - offset, ii_end + offset
               ij = (j - 1)*iim + i
               phi(ij, 1) = phis(ij)
            END DO
         END DO

         ! for other layers
         DO l = 1, llm - 1
            DO j = jj_begin - offset, jj_end + offset
               DO i = ii_begin - offset, ii_end + offset
                  ij = (j - 1)*iim + i
                  mg_ij = p(ij, l + 1) - p(ij, l)
                  ! v=RT/p=(kappa.cpp).Theta.(p/preff)**kappa /p
                  p_ij = .5*(p(ij, l + 1) + p(ij, l))
                  exner_ij = (p_ij/preff)**kappa ! NB : no cpp factor
                  phi(ij, l + 1) = phi(ij, l) + (kappa*cpp)*mg_ij*theta(ij, l, 1)*exner_ij/p_ij
               END DO
            END DO
         END DO
      END IF
! flush phi
!$OMP BARRIER

   END SUBROUTINE compute_geopotential

END MODULE geopotential_mod
