module exner_milieu_loc_m

  IMPLICIT NONE

CONTAINS

  SUBROUTINE exner_milieu_loc( ngrid, ps, p, pks, pk, pkf )

    !     Auteurs :  F. Forget , Y. Wanherdrick
    ! P.Le Van  , Fr. Hourdin  .
    !    ..........

    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...

    !   ************************************************************************
    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    !   ************************************************************************
    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    !    la pression et la fonction d'Exner  au  sol  .

    !     WARNING : CECI est une version speciale de exner_hyb originale
    !               Utilise dans la version martienne pour pouvoir 
    !               tourner avec des coordonnees verticales complexe
    !              => Il ne verifie PAS la condition la proportionalite en 
    !              energie totale/ interne / potentielle (F.Forget 2001)
    !    ( voir note de Fr.Hourdin )  ,

    USE parallel_lmdz
    USE lmdz_filtreg_p
    USE comconst_mod, ONLY: cpp, kappa, r, jmp1
    USE comvert_mod, ONLY: preff

    
    IMPLICIT NONE

    include "dimensions.h"
    include "paramet.h"
    include "comgeom.h"

    INTEGER  ngrid
    REAL p(ijb_u:ije_u,llmp1),pk(ijb_u:ije_u,llm)
    REAL, optional:: pkf(ijb_u:ije_u,llm)
    REAL ps(ijb_u:ije_u),pks(ijb_u:ije_u)

    !    .... variables locales   ...

    INTEGER l, ij
    REAL dum1

    INTEGER ije,ijb,jje,jjb
    logical,save :: firstcall=.TRUE.
    !$OMP THREADPRIVATE(firstcall) 
    CHARACTER(LEN=*),parameter :: modname="exner_milieu_loc"

    ! Sanity check
    IF (firstcall) THEN
       ! sanity checks for Shallow Water case (1 vertical layer)
       IF (llm==1) THEN
          IF (kappa/=1) THEN
             CALL abort_gcm(modname, &
                  "kappa!=1 , but running in Shallow Water mode!!",42)
          endif
          IF (cpp/=r) THEN
             CALL abort_gcm(modname, &
                  "cpp!=r , but running in Shallow Water mode!!",42)
          endif
       endif ! of if (llm.EQ.1)

       firstcall=.FALSE.
    endif ! of if (firstcall)

    !$OMP BARRIER

    ! Specific behaviour for Shallow Water (1 vertical layer) case:
    IF (llm==1) THEN
       ! Compute pks(:),pk(:),pkf(:)
       ijb=ij_begin
       ije=ij_end
       !$OMP DO SCHEDULE(STATIC)
       DO ij=ijb, ije
          pks(ij) = (cpp/preff) * ps(ij)
          pk(ij,1) = .5*pks(ij)
          IF (present(pkf)) pkf(ij,1)=pk(ij,1)
       ENDDO
       !$OMP ENDDO

       !$OMP BARRIER
       IF (present(pkf)) THEN
          jjb=jj_begin
          jje=jj_end
          CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, &
               2, 1, .TRUE., 1 )
       end if

       ! our work is done, exit routine
       RETURN
    endif ! of if (llm.EQ.1)

    ! General case:

    !     -------------
    !     Calcul de pks
    !     -------------

    ijb=ij_begin
    ije=ij_end

    !$OMP DO SCHEDULE(STATIC)
    DO   ij  = ijb, ije
       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    ENDDO
    !$OMP ENDDO
    ! Synchro OPENMP ici

    !$OMP BARRIER


    !    .... Calcul de pk  pour la couche l 
    !    --------------------------------------------

    dum1 = cpp * (2*preff)**(-kappa) 
    DO l = 1, llm-1
       !$OMP DO SCHEDULE(STATIC)
       DO   ij   = ijb, ije
          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
       ENDDO
       !$OMP ENDDO NOWAIT
    ENDDO

    !    .... Calcul de pk  pour la couche l = llm ..
    !    (on met la meme distance (en log pression)  entre Pk(llm)
    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)

    !$OMP DO SCHEDULE(STATIC)
    DO   ij   = ijb, ije
       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
    ENDDO
    !$OMP ENDDO NOWAIT        

    IF (present(pkf)) THEN
       !    calcul de pkf

       DO l = 1, llm
          !$OMP DO SCHEDULE(STATIC)
          DO   ij   = ijb, ije
             pkf(ij,l)=pk(ij,l)
          ENDDO
          !$OMP ENDDO NOWAIT
       ENDDO

       !$OMP BARRIER

       jjb=jj_begin
       jje=jj_end
       CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, &
            2, 1, .TRUE., 1 )
    end if

  END SUBROUTINE exner_milieu_loc

END MODULE exner_milieu_loc_m