MODULE geopotential_mod

CONTAINS

  SUBROUTINE geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
  USE transfert_mod
  USE field_mod
  USE dimensions
  USE geometry
  USE domain_mod
  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
      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 dimensions
  USE geometry
  USE metric
  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

  ! for first layer
   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
   
  END SUBROUTINE compute_geopotential

END MODULE geopotential_mod