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)
       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
         ENDDO
       ENDDO
     ENDDO
   ENDIF
! flush phi
!$OMP BARRIER
   
  END SUBROUTINE compute_geopotential

END MODULE geopotential_mod