c
c $Id$
c
      SUBROUTINE  exner_hyb_loc(ngrid, ps, p,alpha,beta, pks,pk,pkf)
c
c     Auteurs :  P.Le Van  , Fr. Hourdin  .
c    ..........
c
c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
c
c   ************************************************************************
c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des 
c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
c   ************************************************************************
c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
c    la pression et la fonction d'Exner  au  sol  .
c
c                                 -------- z                                   
c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
c    ( voir note de Fr.Hourdin )  ,
c
c    on determine successivement , du haut vers le bas des couches, les 
c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
c
c
      USE parallel_lmdz
      USE mod_filtreg_p
      USE write_field_loc
      IMPLICIT NONE
c
#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
#include "comgeom.h"
#include "comvert.h"
#include "serre.h"

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

c    .... variables locales   ...

      INTEGER l, ij
      REAL unpl2k,dellta

      REAL ppn(iim),pps(iim)
      REAL xpn, xps
      REAL SSUM
      EXTERNAL SSUM
      INTEGER ije,ijb,jje,jjb
      logical,save :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall) 
      character(len=*),parameter :: modname="exner_hyb_loc"
c
c$OMP BARRIER           

      ! Sanity check
      if (firstcall) then
        ! sanity checks for Shallow Water case (1 vertical layer)
        if (llm.eq.1) then
          if (kappa.ne.1) then
            call abort_gcm(modname,
     &      "kappa!=1 , but running in Shallow Water mode!!",42)
          endif
          if (cpp.ne.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)

c$OMP BARRIER

! Specific behaviour for Shallow Water (1 vertical layer) case
      if (llm.eq.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)
          pkf(ij,1)=pk(ij,1)
        ENDDO
!$OMP ENDDO

!$OMP MASTER
      if (pole_nord) then
        DO  ij   = 1, iim
          ppn(ij) = aire(   ij   ) * pks(  ij     )
        ENDDO
        xpn      = SSUM(iim,ppn,1) /apoln
  
        DO ij   = 1, iip1
          pks(   ij     )  =  xpn
          pk(ij,1) = .5*pks(ij)
          pkf(ij,1)=pk(ij,1)
        ENDDO
      endif
      
      if (pole_sud) then
        DO  ij   = 1, iim
          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
        ENDDO
        xps      = SSUM(iim,pps,1) /apols 
  
        DO ij   = 1, iip1
          pks( ij+ip1jm )  =  xps
          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
        ENDDO
      endif
!$OMP END MASTER
!$OMP BARRIER
        jjb=jj_begin
        jje=jj_end
        CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,
     &                 2, 1, .TRUE., 1 )

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


      unpl2k    = 1.+ 2.* kappa
c
      ijb=ij_begin
      ije=ij_end

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

c$OMP MASTER
      if (pole_nord) then
        DO  ij   = 1, iim
          ppn(ij) = aire(   ij   ) * pks(  ij     )
        ENDDO
        xpn      = SSUM(iim,ppn,1) /apoln
  
        DO ij   = 1, iip1
          pks(   ij     )  =  xpn
        ENDDO
      endif
      
      if (pole_sud) then
        DO  ij   = 1, iim
          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
        ENDDO
        xps      = SSUM(iim,pps,1) /apols 
  
        DO ij   = 1, iip1
          pks( ij+ip1jm )  =  xps
        ENDDO
      endif
c$OMP END MASTER
c$OMP BARRIER
c
c
c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
c
c$OMP DO SCHEDULE(STATIC)
      DO     ij      = ijb,ije
       alpha(ij,llm) = 0.
       beta (ij,llm) = 1./ unpl2k
      ENDDO
c$OMP ENDDO NOWAIT
c
c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
c
      DO l = llm -1 , 2 , -1
c
c$OMP DO SCHEDULE(STATIC)
        DO ij = ijb, ije
        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
        beta (ij,l)  =   p(ij,l  ) / dellta   
        ENDDO
c$OMP ENDDO NOWAIT
c
      ENDDO

c
c  ***********************************************************************
c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
c
c$OMP DO SCHEDULE(STATIC)
      DO   ij   = ijb, ije
       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
      ENDDO
c$OMP ENDDO NOWAIT
c
c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
c
      DO l = 2, llm
c$OMP DO SCHEDULE(STATIC)
        DO   ij   = ijb, ije
         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
        ENDDO
c$OMP ENDDO NOWAIT        
      ENDDO
c
c
c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
      DO l = 1, llm
c$OMP DO SCHEDULE(STATIC)
         DO   ij   = ijb, ije
           pkf(ij,l)=pk(ij,l)
         ENDDO
c$OMP ENDDO NOWAIT             
      ENDDO

c$OMP BARRIER
      
      jjb=jj_begin
      jje=jj_end
#ifdef DEBUG_IO    
      call WriteField_u('pkf',pkf)
#endif
      CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,
     &                 2, 1, .TRUE., 1 )
#ifdef DEBUG_IO    
      call WriteField_u('pkf',pkf)
#endif      

      RETURN
      END