MODULE geopotential_mod

CONTAINS

  SUBROUTINE geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field), POINTER :: f_phis(:)
    TYPE(t_field), POINTER :: f_pks(:)
    TYPE(t_field), POINTER :: f_pk(:)
    TYPE(t_field), POINTER :: f_theta(:)
    TYPE(t_field), POINTER :: f_phi(:)
  
    REAL(rstd), POINTER :: phis(:)
    REAL(rstd), POINTER :: pks(:)
    REAL(rstd), POINTER :: pk(:,:)
    REAL(rstd), POINTER :: theta(:,:)
    REAL(rstd), POINTER :: phi(:,:)
    INTEGER :: ind

    DO ind=1,ndomain
      IF (.NOT. assigned_domain(ind)) CYCLE
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      phis=f_phis(ind)
      pks=f_pks(ind)
      pk=f_pk(ind)
      theta=f_theta(ind)
      phi=f_phi(ind)
      CALL compute_geopotential(phis,pks,pk,theta,phi,0)
    ENDDO
  
  END SUBROUTINE geopotential
  
  SUBROUTINE compute_geopotential(phis,pks,pk,theta,phi,offset)
  USE icosa
  USE omp_para
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
    REAL(rstd),INTENT(IN) :: pks(iim*jjm)
    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
    REAL(rstd),INTENT(OUT) :: phi(iim*jjm,llm)
    INTEGER,INTENT(IN) :: offset
    INTEGER :: i,j,ij,l

!!! 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 ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
       ENDDO
     ENDDO
          
    ! for other layers
     DO l = 2, llm
       DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
           ij=(j-1)*iim+i
           phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
                                         * (  pk(ij,l-1) -  pk(ij,l)    )
         ENDDO
       ENDDO
     ENDDO
   ENDIF
! flush phi
!$OMP BARRIER
   
  END SUBROUTINE compute_geopotential

END MODULE geopotential_mod