Index: LMDZ5/trunk/libf/dyn3d/gcm.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/gcm.F	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3d/gcm.F	(revision 2021)
@@ -105,7 +105,4 @@
       REAL ps(ip1jmp1)                       ! pression  au sol
       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-      REAL pks(ip1jmp1)                      ! exner au  sol
-      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL masse(ip1jmp1,llm)                ! masse d'air
       REAL phis(ip1jmp1)                     ! geopotentiel au sol
@@ -131,5 +128,4 @@
       data call_iniphys/.true./
 
-      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
 c+jld variables test conservation energie
 c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
Index: LMDZ5/trunk/libf/dyn3d/guide_mod.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d/guide_mod.F90	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3d/guide_mod.F90	(revision 2021)
@@ -589,4 +589,6 @@
   SUBROUTINE guide_interp(psi,teta)
   
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
   IMPLICIT NONE
 
@@ -610,6 +612,5 @@
   REAL, DIMENSION (iip1,jjm,llm)     :: pbary 
   ! Variables pour fonction Exner (P milieu couche)
-  REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
-  REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
+  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
   REAL, DIMENSION (iip1,jjp1)        :: pks    
   REAL                               :: prefkap,unskap
@@ -676,7 +677,7 @@
     CALL pression( ip1jmp1, ap, bp, psi, p )
     if (pressure_exner) then
-      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+      CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
     else
-      CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
+      CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
     endif
 !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
Index: LMDZ5/trunk/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d/iniacademic.F90	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3d/iniacademic.F90	(revision 2021)
@@ -14,4 +14,6 @@
 #endif
   USE Write_Field
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -54,5 +56,4 @@
   REAL pks(ip1jmp1)                      ! exner au  sol
   REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
   REAL phi(ip1jmp1,llm)                  ! geopotentiel
   REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
@@ -70,5 +71,5 @@
   integer idum
 
-  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+  REAL zdtvr
   
   character(len=*),parameter :: modname="iniacademic"
@@ -223,7 +224,7 @@
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
         if (pressure_exner) then
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        else
-          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+          CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
+        else
+          call exner_milieu(ip1jmp1,ps,p,pks,pk)
         endif
         CALL massdair(p,masse)
Index: LMDZ5/trunk/libf/dyn3d/leapfrog.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/leapfrog.F	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3d/leapfrog.F	(revision 2021)
@@ -19,4 +19,7 @@
      &                       iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
      &                       periodav, ok_dyn_ave, output_grads_dyn
+      use exner_hyb_m, only: exner_hyb
+      use exner_milieu_m, only: exner_milieu
+
       IMPLICIT NONE
 
@@ -158,5 +161,4 @@
       character*10 string10
 
-      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
       REAL :: flxw(ip1jmp1,llm)  ! flux de masse verticale
 
@@ -217,7 +219,7 @@
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       if (pressure_exner) then
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
       else
-        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
       endif
 
@@ -373,7 +375,7 @@
          CALL pression (  ip1jmp1, ap, bp, ps,  p      )
          if (pressure_exner) then
-           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+           CALL exner_hyb(  ip1jmp1, ps, p,pks, pk, pkf )
          else
-           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+           CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
          endif
 
@@ -448,7 +450,7 @@
           CALL massdair(p,masse)
           if (pressure_exner) then
-            CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+            CALL exner_hyb(ip1jmp1,ps,p,pks,pk,pkf)
           else
-            CALL exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+            CALL exner_milieu(ip1jmp1,ps,p,pks,pk,pkf)
           endif
 
@@ -506,7 +508,7 @@
         CALL pression ( ip1jmp1, ap, bp, ps, p                  )
         if (pressure_exner) then
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+          CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
         else
-          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+          CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
         endif
         CALL massdair(p,masse)
Index: LMDZ5/trunk/libf/dyn3d_common/exner_hyb.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/exner_hyb.F	(revision 2020)
+++ 	(revision )
@@ -1,152 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_hyb ( 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
-      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(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
-
-c    .... variables locales   ...
-
-      INTEGER l, ij
-      REAL unpl2k,dellta
-
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-c
-      logical,save :: firstcall=.true.
-      character(len=*),parameter :: modname="exner_hyb"
-      
-      ! 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)
-
-      if (llm.eq.1) then
-        
-        ! Compute pks(:),pk(:),pkf(:)
-        
-        DO   ij  = 1, ngrid
-          pks(ij) = (cpp/preff) * ps(ij) 
-          pk(ij,1) = .5*pks(ij)
-        ENDDO
-        
-        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
-        
-        ! our work is done, exit routine
-        return
-
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-     
-      unpl2k    = 1.+ 2.* kappa
-c
-      DO   ij  = 1, ngrid
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-
-      DO  ij   = 1, iim
-        ppn(ij) = aire(   ij   ) * pks(  ij     )
-        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-      ENDDO
-      xpn      = SSUM(iim,ppn,1) /apoln
-      xps      = SSUM(iim,pps,1) /apols
-
-      DO ij   = 1, iip1
-        pks(   ij     )  =  xpn
-        pks( ij+ip1jm )  =  xps
-      ENDDO
-c
-c
-c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
-c
-      DO     ij      = 1, ngrid
-       alpha(ij,llm) = 0.
-       beta (ij,llm) = 1./ unpl2k
-      ENDDO
-c
-c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
-c
-      DO l = llm -1 , 2 , -1
-c
-        DO ij = 1, ngrid
-        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
-      ENDDO
-c
-c  ***********************************************************************
-c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
-c
-      DO   ij   = 1, ngrid
-       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
-c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
-c
-      DO l = 2, llm
-        DO   ij   = 1, ngrid
-         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-        ENDDO
-      ENDDO
-c
-c
-      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3d_common/exner_hyb_m.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/exner_hyb_m.F90	(revision 2021)
+++ LMDZ5/trunk/libf/dyn3d_common/exner_hyb_m.F90	(revision 2021)
@@ -0,0 +1,145 @@
+module exner_hyb_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_hyb ( ngrid, ps, p, pks, pk, pkf )
+
+    !     Auteurs :  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  .
+    !
+    !                                 -------- z
+    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
+    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !    on determine successivement , du haut vers le bas des couches, les 
+    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
+    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
+    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
+    !
+    !
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    real, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL unpl2k,dellta
+
+    logical,save :: firstcall=.true.
+    character(len=*),parameter :: modname="exner_hyb"
+
+    ! 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)
+
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+
+       DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij)
+          pk(ij,1) = .5*pks(ij)
+       ENDDO
+
+       if (present(pkf)) then
+          pkf = pk
+          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+       end if
+
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
+
+    ! General case:
+
+    unpl2k    = 1.+ 2.* kappa
+
+    !     -------------
+    !     Calcul de pks
+    !     -------------
+
+    DO   ij  = 1, ngrid
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+
+    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
+    !
+    DO     ij      = 1, ngrid
+       alpha(ij,llm) = 0.
+       beta (ij,llm) = 1./ unpl2k
+    ENDDO
+    !
+    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
+    !
+    DO l = llm -1 , 2 , -1
+       !
+       DO ij = 1, ngrid
+          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
+    ENDDO
+
+    !  ***********************************************************************
+    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
+    !
+    DO   ij   = 1, ngrid
+       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
+    !
+    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
+    !
+    DO l = 2, llm
+       DO   ij   = 1, ngrid
+          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+       ENDDO
+    ENDDO
+
+    if (present(pkf)) then
+       !    calcul de pkf
+       pkf = pk
+       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_hyb
+
+end module exner_hyb_m
Index: LMDZ5/trunk/libf/dyn3d_common/exner_milieu.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/exner_milieu.F	(revision 2020)
+++ 	(revision )
@@ -1,144 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
-c
-c     Auteurs :  F. Forget , Y. Wanherdrick
-c P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** 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     WARNING : CECI est une version speciale de exner_hyb originale
-c               Utilise dans la version martienne pour pouvoir 
-c               tourner avec des coordonnees verticales complexe
-c              => Il ne verifie PAS la condition la proportionalite en 
-c              energie totale/ interne / potentielle (F.Forget 2001)
-c    ( voir note de Fr.Hourdin )  ,
-c
-      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(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
-
-c    .... variables locales   ...
-
-      INTEGER l, ij
-      REAL dum1
-
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-      EXTERNAL SSUM
-      logical,save :: firstcall=.true.
-      character(len=*),parameter :: modname="exner_milieu"
-
-      ! 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)
-
-!!!! Specific behaviour for Shallow Water (1 vertical layer) case:
-      if (llm.eq.1) then
-      
-        ! Compute pks(:),pk(:),pkf(:)
-        
-        DO   ij  = 1, ngrid
-          pks(ij) = (cpp/preff) * ps(ij) 
-          pk(ij,1) = .5*pks(ij)
-        ENDDO
-        
-        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
-        
-        ! our work is done, exit routine
-        return
-
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-
-c     -------------
-c     Calcul de pks
-c     -------------
-   
-      DO   ij  = 1, ngrid
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-
-      DO  ij   = 1, iim
-        ppn(ij) = aire(   ij   ) * pks(  ij     )
-        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-      ENDDO
-      xpn      = SSUM(iim,ppn,1) /apoln
-      xps      = SSUM(iim,pps,1) /apols
-
-      DO ij   = 1, iip1
-        pks(   ij     )  =  xpn
-        pks( ij+ip1jm )  =  xps
-      ENDDO
-c
-c
-c    .... Calcul de pk  pour la couche l 
-c    --------------------------------------------
-c
-      dum1 = cpp * (2*preff)**(-kappa) 
-      DO l = 1, llm-1
-        DO   ij   = 1, ngrid
-         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
-        ENDDO
-      ENDDO
-
-c    .... Calcul de pk  pour la couche l = llm ..
-c    (on met la meme distance (en log pression)  entre Pk(llm)
-c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
-
-      DO   ij   = 1, ngrid
-         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
-      ENDDO
-
-
-c    calcul de pkf
-c    -------------
-      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-c    EST-CE UTILE ?? : calcul de beta
-c    --------------------------------
-      DO l = 2, llm
-        DO   ij   = 1, ngrid
-          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
-        ENDDO
-      ENDDO
-
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3d_common/exner_milieu_m.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/exner_milieu_m.F90	(revision 2021)
+++ LMDZ5/trunk/libf/dyn3d_common/exner_milieu_m.F90	(revision 2021)
@@ -0,0 +1,124 @@
+module exner_milieu_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_milieu ( 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 )  ,
+    !
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    real, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL dum1
+
+    logical,save :: firstcall=.true.
+    character(len=*),parameter :: modname="exner_milieu"
+
+    ! 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)
+
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+
+       DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij)
+          pk(ij,1) = .5*pks(ij)
+       ENDDO
+
+       if (present(pkf)) then
+          pkf = pk
+          CALL filtreg ( pkf, 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
+    !     -------------
+
+    DO   ij  = 1, ngrid
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+
+    !    .... Calcul de pk  pour la couche l 
+    !    --------------------------------------------
+    !
+    dum1 = cpp * (2*preff)**(-kappa) 
+    DO l = 1, llm-1
+       DO   ij   = 1, ngrid
+          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
+       ENDDO
+    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)
+
+    DO   ij   = 1, ngrid
+       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
+    ENDDO
+
+    if (present(pkf)) then
+       !    calcul de pkf
+       pkf = pk
+       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_milieu
+
+end module exner_milieu_m
Index: LMDZ5/trunk/libf/dyn3dmem/call_calfis_mod.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/call_calfis_mod.F90	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dmem/call_calfis_mod.F90	(revision 2021)
@@ -12,6 +12,4 @@
 
     REAL,POINTER,SAVE :: p(:,:) 
-    REAL,POINTER,SAVE :: alpha(:,:) 
-    REAL,POINTER,SAVE :: beta(:,:) 
     REAL,POINTER,SAVE :: pks(:) 
     REAL,POINTER,SAVE :: pk(:,:) 
@@ -53,6 +51,4 @@
     CALL allocate_u(flxw,llm,d)
     CALL allocate_u(p,llmp1,d)
-    CALL allocate_u(alpha,llm,d)
-    CALL allocate_u(beta,llm,d)
     CALL allocate_u(pks,d)
     CALL allocate_u(pk,llm,d)
@@ -75,4 +71,6 @@
                          phis_dyn,q_dyn,flxw_dyn)
   USE dimensions_mod
+  use exner_hyb_loc_m, only: exner_hyb_loc
+  use exner_milieu_loc_m, only: exner_milieu_loc
   USE parallel_lmdz
   USE times
@@ -201,5 +199,5 @@
 
   !$OMP BARRIER
-    CALL exner_hyb_loc(  ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+    CALL exner_hyb_loc(  ip1jmp1, ps, p, pks, pk, pkf )
   !$OMP BARRIER
     CALL geopot_loc  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
@@ -343,7 +341,7 @@
 !$OMP BARRIER
     if (pressure_exner) then
-      CALL exner_hyb_loc(ijnb_u,ps,p,alpha,beta,pks,pk,pkf)
+      CALL exner_hyb_loc(ijnb_u,ps,p,pks,pk,pkf)
     else 
-      CALL exner_milieu_loc(ijnb_u,ps,p,beta,pks,pk,pkf)
+      CALL exner_milieu_loc(ijnb_u,ps,p,pks,pk,pkf)
     endif
 !$OMP BARRIER
Index: LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc.F	(revision 2020)
+++ 	(revision )
@@ -1,242 +1,0 @@
-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
Index: LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc_m.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc_m.F90	(revision 2021)
+++ LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc_m.F90	(revision 2021)
@@ -0,0 +1,198 @@
+module exner_hyb_loc_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_hyb_loc(ngrid, ps, p, pks,pk,pkf)
+
+    !     Auteurs :  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  .
+    !
+    !                                 -------- z
+    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
+    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !    on determine successivement , du haut vers le bas des couches, les 
+    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
+    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
+    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
+    !
+    !
+    USE parallel_lmdz
+    USE mod_filtreg_p
+    USE write_field_loc
+    !
+    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, optional:: 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)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL unpl2k,dellta
+
+    INTEGER ije,ijb,jje,jjb
+    logical,save :: firstcall=.true.
+    !$OMP THREADPRIVATE(firstcall) 
+    character(len=*),parameter :: modname="exner_hyb_loc"
+    !
+    !$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)
+
+    !$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)
+          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:
+
+    unpl2k    = 1.+ 2.* kappa
+
+    !     -------------
+    !     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 des coeff. alpha et beta  pour la couche l = llm ..
+    !
+    !$OMP DO SCHEDULE(STATIC)
+    DO     ij      = ijb,ije
+       alpha(ij,llm) = 0.
+       beta (ij,llm) = 1./ unpl2k
+    ENDDO
+    !$OMP ENDDO NOWAIT
+    !
+    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
+    !
+    DO l = llm -1 , 2 , -1
+       !
+       !$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
+       !$OMP ENDDO NOWAIT
+    ENDDO
+
+    !  ***********************************************************************
+    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
+    !
+    !$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
+    !$OMP ENDDO NOWAIT
+    !
+    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
+    !
+    DO l = 2, llm
+       !$OMP DO SCHEDULE(STATIC)
+       DO   ij   = ijb, ije
+          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+       ENDDO
+       !$OMP ENDDO NOWAIT        
+    ENDDO
+
+    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
+#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     
+    end if
+
+  END SUBROUTINE exner_hyb_loc
+
+end module exner_hyb_loc_m
Index: LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc.F	(revision 2020)
+++ 	(revision )
@@ -1,224 +1,0 @@
-!
-! $Id$
-!
-      SUBROUTINE  exner_milieu_loc ( ngrid, ps, p,beta, pks, pk, pkf )
-c
-c     Auteurs :  F. Forget , Y. Wanherdrick
-c P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** 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     WARNING : CECI est une version speciale de exner_hyb originale
-c               Utilise dans la version martienne pour pouvoir 
-c               tourner avec des coordonnees verticales complexe
-c              => Il ne verifie PAS la condition la proportionalite en 
-c              energie totale/ interne / potentielle (F.Forget 2001)
-c    ( voir note de Fr.Hourdin )  ,
-c
-      USE parallel_lmdz
-      USE mod_filtreg_p
-      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 dum1
-
-      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_milieu_loc"
-
-      ! 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)
-
-!!!! General case:
-
-c     -------------
-c     Calcul de pks
-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 de pk  pour la couche l 
-c    --------------------------------------------
-c
-      dum1 = cpp * (2*preff)**(-kappa) 
-      DO l = 1, llm-1
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
-        ENDDO
-c$OMP ENDDO NOWAIT        
-      ENDDO
-
-c    .... Calcul de pk  pour la couche l = llm ..
-c    (on met la meme distance (en log pression)  entre Pk(llm)
-c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
-
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij   = ijb, ije
-         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
-      ENDDO
-c$OMP ENDDO NOWAIT        
-
-
-c    calcul de pkf
-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
-      CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,
-     &                 2, 1, .TRUE., 1 )
-      
-c    EST-CE UTILE ?? : calcul de beta
-c    --------------------------------
-      DO l = 2, llm
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
-        ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
-
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc_m.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc_m.F90	(revision 2021)
+++ LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc_m.F90	(revision 2021)
@@ -0,0 +1,163 @@
+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 mod_filtreg_p
+    !
+    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, 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.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)
+
+    !$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)
+          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
Index: LMDZ5/trunk/libf/dyn3dmem/gcm.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/gcm.F	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dmem/gcm.F	(revision 2021)
@@ -98,7 +98,4 @@
       REAL,ALLOCATABLE,SAVE  :: ps(:)         ! pression  au sol
 c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-c      REAL pks(ip1jmp1)                      ! exner au  sol
-c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-c      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL,ALLOCATABLE,SAVE  :: masse(:,:)    ! masse d'air
       REAL,ALLOCATABLE,SAVE  :: phis(:)       ! geopotentiel au sol
@@ -124,5 +121,4 @@
       data call_iniphys/.true./
 
-c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
 c+jld variables test conservation energie
 c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
Index: LMDZ5/trunk/libf/dyn3dmem/guide_loc_mod.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/guide_loc_mod.F90	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dmem/guide_loc_mod.F90	(revision 2021)
@@ -329,4 +329,6 @@
 !=======================================================================
   SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
+    use exner_hyb_loc_m, only: exner_hyb_loc
+    use exner_milieu_loc_m, only: exner_milieu_loc
     USE parallel_lmdz
     USE control_mod
@@ -353,6 +355,5 @@
     REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: f_addv ! var aux: champ de guidage
     ! Variables pour fonction Exner (P milieu couche)
-    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)    :: pk, pkf
-    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)    :: alpha, beta
+    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)    :: pk
     REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)        :: pks    
     REAL                               :: unskap
@@ -399,7 +400,4 @@
         ALLOCATE(f_addv(ijb_v:ije_v,llm) ) 
         ALLOCATE(pk(iip1,jjb_u:jje_u,llm)  ) 
-        ALLOCATE(pkf(iip1,jjb_u:jje_u,llm)  ) 
-        ALLOCATE(alpha(iip1,jjb_u:jje_u,llm)  ) 
-        ALLOCATE(beta(iip1,jjb_u:jje_u,llm)  ) 
         ALLOCATE(pks(iip1,jjb_u:jje_u)  ) 
         ALLOCATE(p(ijb_u:ije_u,llmp1) ) 
@@ -539,7 +537,7 @@
 	CALL pression_loc( ijnb_u, ap, bp, ps, p )
 	if (pressure_exner) then
-          CALL exner_hyb_loc(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+          CALL exner_hyb_loc(ip1jmp1,ps,p,pks,pk)
 	else
-          CALL exner_milieu_loc(ip1jmp1,ps,p,beta,pks,pk,pkf)
+          CALL exner_milieu_loc(ip1jmp1,ps,p,pks,pk)
         endif
 !$OMP BARRIER        
@@ -894,4 +892,6 @@
 !=======================================================================
   SUBROUTINE guide_interp(psi,teta)
+    use exner_hyb_loc_m, only: exner_hyb_loc
+    use exner_milieu_loc_m, only: exner_milieu_loc
   USE parallel_lmdz
   USE mod_hallo
@@ -919,6 +919,5 @@
   REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)     :: pbary 
   ! Variables pour fonction Exner (P milieu couche)
-  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)    :: pk, pkf
-  REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)    :: alpha, beta
+  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)    :: pk
   REAL ,ALLOCATABLE, SAVE, DIMENSION (:,:)        :: pks    
   REAL                               :: unskap
@@ -949,7 +948,4 @@
       ALLOCATE(pbary(iip1,jjb_v:jje_v,llm) )    
       ALLOCATE(pk(iip1,jjb_u:jje_u,llm) )    
-      ALLOCATE(pkf(iip1,jjb_u:jje_u,llm)  )    
-      ALLOCATE(alpha(iip1,jjb_u:jje_u,llm) )    
-      ALLOCATE(beta(iip1,jjb_u:jje_u,llm) )    
       ALLOCATE(pks (iip1,jjb_u:jje_u) )    
       ALLOCATE(qsat(ijb_u:ije_u,llm) )    
@@ -1032,7 +1028,7 @@
 	CALL pression_loc( ijnb_u, ap, bp, psi, p )
 	if (disvert_type==1) then
-          CALL exner_hyb_loc(ijnb_u,psi,p,alpha,beta,pks,pk,pkf)
+          CALL exner_hyb_loc(ijnb_u,psi,p,pks,pk)
         else ! we assume that we are in the disvert_type==2 case
-          CALL exner_milieu_loc(ijnb_u,psi,p,beta,pks,pk,pkf)
+          CALL exner_milieu_loc(ijnb_u,psi,p,pks,pk)
         endif
 	unskap=1./kappa
Index: LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 2021)
@@ -4,4 +4,6 @@
 SUBROUTINE iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
 
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
   USE filtreg_mod
   USE infotrac, ONLY : nqtot
@@ -58,5 +60,4 @@
   REAL pks(ip1jmp1)                      ! exner au  sol
   REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
   REAL phi(ip1jmp1,llm)                  ! geopotentiel
   REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
@@ -75,5 +76,4 @@
 
   REAL zdtvr
-  real,allocatable :: alpha(:,:),beta(:,:)
   
   character(len=*),parameter :: modname="iniacademic"
@@ -219,6 +219,4 @@
        allocate(masse_glo(ip1jmp1,llm))
        allocate(phis_glo(ip1jmp1))
-       allocate(alpha(ip1jmp1,llm))
-       allocate(beta(ip1jmp1,llm))
 
         ! surface pressure
@@ -238,7 +236,7 @@
         CALL pression ( ip1jmp1, ap, bp, ps_glo, p       )
         if (pressure_exner) then
-          CALL exner_hyb( ip1jmp1, ps_glo, p,alpha,beta, pks, pk, pkf )
+          CALL exner_hyb( ip1jmp1, ps_glo, p, pks, pk )
         else
-          call exner_milieu(ip1jmp1,ps_glo,p,beta,pks,pk,pkf)
+          call exner_milieu(ip1jmp1,ps_glo,p,pks,pk)
         endif
         CALL massdair(p,masse_glo)
@@ -301,6 +299,4 @@
         deallocate(ps_glo)
         deallocate(phis_glo)
-        deallocate(alpha)
-        deallocate(beta)
      ENDIF ! of IF (.NOT. read_start)
   endif academic_case
Index: LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F	(revision 2021)
@@ -31,4 +31,6 @@
        USE call_calfis_mod, ONLY : call_calfis
        USE leapfrog_mod
+       use exner_hyb_loc_m, only: exner_hyb_loc
+       use exner_milieu_loc_m, only: exner_milieu_loc
       IMPLICIT NONE
 
@@ -156,5 +158,4 @@
       character*10 string10
 
-!      REAL,SAVE,ALLOCATABLE :: alpha(:,:),beta(:,:)
 !      REAL,SAVE,ALLOCATABLE :: flxw(:,:) ! flux de masse verticale
 
@@ -261,5 +262,4 @@
 !      ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
 !      ALLOCATE(finvmaold(ijb_u:ije_u,llm))
-!      ALLOCATE(alpha(ijb_u:ije_u,llm),beta(ijb_u:ije_u,llm))
 !      ALLOCATE(flxw(ijb_u:ije_u,llm))
 !      ALLOCATE(ecin(ijb_u:ije_u,llm),ecin0(ijb_u:ije_u,llm))
@@ -284,7 +284,7 @@
 c$OMP END MASTER
       if (pressure_exner) then
-      CALL exner_hyb_loc( ijnb_u, ps, p,alpha,beta, pks, pk, pkf)
+      CALL exner_hyb_loc( ijnb_u, ps, p, pks, pk, pkf)
       else 
-        CALL exner_milieu_loc( ijnb_u, ps, p, beta, pks, pk, pkf )
+        CALL exner_milieu_loc( ijnb_u, ps, p, pks, pk, pkf )
       endif
 c-----------------------------------------------------------------------
@@ -780,5 +780,5 @@
 
 ! c$OMP BARRIER
-!          CALL exner_hyb_loc(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+!          CALL exner_hyb_loc(  ip1jmp1, ps, p,pks, pk, pkf )
 ! c$OMP BARRIER
 !            jD_cur = jD_ref + day_ini - day_ref
@@ -1135,7 +1135,7 @@
 c$OMP BARRIER
         if (pressure_exner) then
-        CALL exner_hyb_loc( ijnb_u, ps, p,alpha,beta, pks, pk, pkf )
+        CALL exner_hyb_loc( ijnb_u, ps, p, pks, pk, pkf )
         else 
-          CALL exner_milieu_loc( ijnb_u, ps, p, beta, pks, pk, pkf )
+          CALL exner_milieu_loc( ijnb_u, ps, p, pks, pk, pkf )
         endif
 c$OMP BARRIER
Index: LMDZ5/trunk/libf/dyn3dmem/leapfrog_mod.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/leapfrog_mod.F90	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dmem/leapfrog_mod.F90	(revision 2021)
@@ -27,6 +27,4 @@
   REAL,POINTER,SAVE :: dq(:,:,:)
   REAL,POINTER,SAVE :: finvmaold(:,:)
-  REAL,POINTER,SAVE :: alpha(:,:)
-  REAL,POINTER,SAVE :: beta(:,:)
   REAL,POINTER,SAVE :: flxw(:,:)
   REAL,POINTER,SAVE :: unat(:,:)
@@ -79,6 +77,4 @@
     CALL allocate_u(dq,llm,nqtot,d)
     CALL allocate_u(finvmaold,llm,d)
-    CALL allocate_u(alpha,llm,d)
-    CALL allocate_u(beta,llm,d)
     CALL allocate_u(flxw,llm,d)
     CALL allocate_u(unat,llm,d)
@@ -129,6 +125,4 @@
     CALL switch_u(dq,distrib_caldyn,dist)
     CALL switch_u(finvmaold,distrib_caldyn,dist)
-    CALL switch_u(alpha,distrib_caldyn,dist)
-    CALL switch_u(beta,distrib_caldyn,dist)
     CALL switch_u(flxw,distrib_caldyn,dist)
     CALL switch_u(unat,distrib_caldyn,dist)
Index: LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p.F	(revision 2020)
+++ 	(revision )
@@ -1,231 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_hyb_p ( 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
-      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(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,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_p"
-c
-
-      ! 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,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-
-        ! our work is done, exit routine
-        return
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-
-      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
-      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p_m.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p_m.F90	(revision 2021)
+++ LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p_m.F90	(revision 2021)
@@ -0,0 +1,186 @@
+module exner_hyb_p_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_hyb_p ( ngrid, ps, p, pks, pk, pkf )
+
+    !     Auteurs :  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  .
+    !
+    !                                 -------- z
+    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
+    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !    on determine successivement , du haut vers le bas des couches, les 
+    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
+    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
+    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
+    !
+    !
+    USE parallel_lmdz
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    REAL, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
+    REAL alpha(ngrid,llm),beta(ngrid,llm)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL unpl2k,dellta
+
+    INTEGER ije,ijb,jje,jjb
+    logical,save :: firstcall=.true.
+    !$OMP THREADPRIVATE(firstcall) 
+    character(len=*),parameter :: modname="exner_hyb_p"
+
+    ! 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)
+
+    !$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)
+          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,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+       end if
+
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
+
+    ! General case:
+
+    unpl2k    = 1.+ 2.* kappa
+
+    !     -------------
+    !     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 des coeff. alpha et beta  pour la couche l = llm ..
+    !
+    !$OMP DO SCHEDULE(STATIC)
+    DO     ij      = ijb,ije
+       alpha(ij,llm) = 0.
+       beta (ij,llm) = 1./ unpl2k
+    ENDDO
+    !$OMP ENDDO NOWAIT
+    !
+    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
+    !
+    DO l = llm -1 , 2 , -1
+       !
+       !$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
+       !$OMP ENDDO NOWAIT
+    ENDDO
+
+    !  ***********************************************************************
+    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
+    !
+    !$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
+    !$OMP ENDDO NOWAIT
+    !
+    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
+    !
+    DO l = 2, llm
+       !$OMP DO SCHEDULE(STATIC)
+       DO   ij   = ijb, ije
+          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+       ENDDO
+       !$OMP ENDDO NOWAIT        
+    ENDDO
+
+    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,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_hyb_p
+
+end module exner_hyb_p_m
Index: LMDZ5/trunk/libf/dyn3dpar/exner_milieu_p.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/exner_milieu_p.F	(revision 2020)
+++ 	(revision )
@@ -1,218 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_milieu_p ( ngrid, ps, p,beta, pks, pk, pkf )
-c
-c     Auteurs :  F. Forget , Y. Wanherdrick
-c P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** 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     WARNING : CECI est une version speciale de exner_hyb originale
-c               Utilise dans la version martienne pour pouvoir 
-c               tourner avec des coordonnees verticales complexe
-c              => Il ne verifie PAS la condition la proportionalite en 
-c              energie totale/ interne / potentielle (F.Forget 2001)
-c    ( voir note de Fr.Hourdin )  ,
-c
-      USE parallel_lmdz
-      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(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
-
-c    .... variables locales   ...
-
-      INTEGER l, ij
-      REAL dum1
-
-      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_milieu_p"
-
-      ! 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,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-
-        ! our work is done, exit routine
-        return
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-
-c     -------------
-c     Calcul de pks
-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 de pk  pour la couche l 
-c    --------------------------------------------
-c
-      dum1 = cpp * (2*preff)**(-kappa) 
-      DO l = 1, llm-1
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
-        ENDDO
-c$OMP ENDDO NOWAIT        
-      ENDDO
-
-c    .... Calcul de pk  pour la couche l = llm ..
-c    (on met la meme distance (en log pression)  entre Pk(llm)
-c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
-
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij   = ijb, ije
-         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
-      ENDDO
-c$OMP ENDDO NOWAIT        
-
-
-c    calcul de pkf
-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
-      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-c    EST-CE UTILE ?? : calcul de beta
-c    --------------------------------
-      DO l = 2, llm
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
-        ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
-
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dpar/exner_milieu_p_m.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/exner_milieu_p_m.F90	(revision 2021)
+++ LMDZ5/trunk/libf/dyn3dpar/exner_milieu_p_m.F90	(revision 2021)
@@ -0,0 +1,159 @@
+module exner_milieu_p_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_milieu_p ( 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
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    REAL, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL dum1
+
+    logical,save :: firstcall=.true.
+    !$OMP THREADPRIVATE(firstcall) 
+    character(len=*),parameter :: modname="exner_milieu_p"
+
+    ! 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)
+
+    !$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)
+          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,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,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_milieu_p
+
+end module exner_milieu_p_m
Index: LMDZ5/trunk/libf/dyn3dpar/gcm.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/gcm.F	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dpar/gcm.F	(revision 2021)
@@ -99,7 +99,4 @@
       REAL ps(ip1jmp1)                       ! pression  au sol
 c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-c      REAL pks(ip1jmp1)                      ! exner au  sol
-c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-c      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL masse(ip1jmp1,llm)                ! masse d'air
       REAL phis(ip1jmp1)                     ! geopotentiel au sol
@@ -125,5 +122,4 @@
       data call_iniphys/.true./
 
-c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
 c+jld variables test conservation energie
 c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
Index: LMDZ5/trunk/libf/dyn3dpar/guide_p_mod.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/guide_p_mod.F90	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dpar/guide_p_mod.F90	(revision 2021)
@@ -328,4 +328,6 @@
 !=======================================================================
   SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
+    use exner_hyb_p_m, only: exner_hyb_p
+    use exner_milieu_p_m, only: exner_milieu_p
     USE parallel_lmdz
     USE control_mod
@@ -349,6 +351,5 @@
     REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
     ! Variables pour fonction Exner (P milieu couche)
-    REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
-    REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
+    REAL, DIMENSION (iip1,jjp1,llm)    :: pk
     REAL, DIMENSION (iip1,jjp1)        :: pks    
     REAL                               :: unskap
@@ -493,7 +494,7 @@
 	CALL pression_p( ip1jmp1, ap, bp, ps, p )
 	if (pressure_exner) then
-          CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+          CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
 	else
-          CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
+          CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
         endif
         unskap=1./kappa
@@ -689,4 +690,6 @@
 !=======================================================================
   SUBROUTINE guide_interp(psi,teta)
+    use exner_hyb_p_m, only: exner_hyb_p
+    use exner_milieu_p_m, only: exner_milieu_p
   USE parallel_lmdz
   USE mod_hallo
@@ -713,6 +716,5 @@
   REAL, DIMENSION (iip1,jjm,llm)     :: pbary 
   ! Variables pour fonction Exner (P milieu couche)
-  REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
-  REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
+  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
   REAL, DIMENSION (iip1,jjp1)        :: pks    
   REAL                               :: unskap
@@ -793,7 +795,7 @@
 	CALL pression_p( ip1jmp1, ap, bp, psi, p )
 	if (pressure_exner) then
-          CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+          CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
         else
-          CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
+          CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
         endif
 	unskap=1./kappa
Index: LMDZ5/trunk/libf/dyn3dpar/iniacademic.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/iniacademic.F90	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dpar/iniacademic.F90	(revision 2021)
@@ -4,4 +4,6 @@
 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
 
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
   USE filtreg_mod
   USE infotrac, ONLY : nqtot
@@ -54,5 +56,4 @@
   REAL pks(ip1jmp1)                      ! exner au  sol
   REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
   REAL phi(ip1jmp1,llm)                  ! geopotentiel
   REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
@@ -70,5 +71,5 @@
   integer idum
 
-  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+  REAL zdtvr
   
   character(len=*),parameter :: modname="iniacademic"
@@ -223,7 +224,7 @@
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
         if (pressure_exner) then
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        else
-          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+          CALL exner_hyb( ip1jmp1, ps, p, pks, pk )
+        else
+          call exner_milieu(ip1jmp1,ps,p,pks,pk)
         endif
         CALL massdair(p,masse)
Index: LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 2020)
+++ LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 2021)
@@ -8,4 +8,8 @@
      &                    time_0)
 
+      use exner_hyb_m, only: exner_hyb
+      use exner_milieu_m, only: exner_milieu
+      use exner_hyb_p_m, only: exner_hyb_p
+      use exner_milieu_p_m, only: exner_milieu_p
        USE misc_mod
        USE parallel_lmdz
@@ -149,5 +153,4 @@
       character*10 string10
 
-      REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
       REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
 
@@ -241,7 +244,7 @@
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       if (pressure_exner) then
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
       else
-        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
       endif
 c$OMP END MASTER
@@ -705,7 +708,7 @@
 c$OMP BARRIER
          if (pressure_exner) then
-           CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+           CALL exner_hyb_p(  ip1jmp1, ps, p,pks, pk, pkf )
          else
-           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+           CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
          endif
 c$OMP BARRIER
@@ -918,7 +921,7 @@
 c$OMP BARRIER
           if (pressure_exner) then
-            CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+            CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk,pkf)
           else
-            CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
+            CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk,pkf)
           endif
 c$OMP BARRIER
@@ -1059,7 +1062,7 @@
 c$OMP BARRIER
         if (pressure_exner) then
-          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+          CALL exner_hyb_p( ip1jmp1, ps, p, pks, pk, pkf )
         else
-          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+          CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
         endif
 c$OMP BARRIER
Index: LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90	(revision 2020)
+++ LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90	(revision 2021)
@@ -29,4 +29,7 @@
   USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
   USE indice_sol_mod
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
+  use test_disvert_m, only: test_disvert
 #endif
   IMPLICIT NONE
@@ -74,5 +77,5 @@
   CHARACTER(LEN=80)                        :: x, fmt
   INTEGER                                  :: i, j, l, ji
-  REAL,    DIMENSION(iip1,jjp1,llm)        :: alpha, beta, pk, pls, y
+  REAL,    DIMENSION(iip1,jjp1,llm)        :: pk, pls, y
   REAL,    DIMENSION(ip1jmp1)              :: pks
 
@@ -150,4 +153,5 @@
 
   CALL iniconst()
+  call test_disvert
   CALL inigeom()
 
@@ -253,7 +257,7 @@
   CALL pression(ip1jmp1, ap, bp, psol, p3d)
   if (pressure_exner) then
-    CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
+    CALL exner_hyb(ip1jmp1, psol, p3d, pks, pk)
   else
-    CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
+    CALL exner_milieu(ip1jmp1,psol,p3d, pks,pk)
   endif
   pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)
Index: LMDZ5/trunk/libf/phylmd/test_disvert_m.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/test_disvert_m.F90	(revision 2021)
+++ LMDZ5/trunk/libf/phylmd/test_disvert_m.F90	(revision 2021)
@@ -0,0 +1,70 @@
+module test_disvert_m
+
+  implicit none
+
+contains
+
+  subroutine test_disvert
+
+    ! Author: Lionel GUEZ
+
+    ! This procedure tests the order of pressure values at half-levels
+    ! and full levels. We arbitrarily choose to test ngrid values of
+    ! the surface pressure, which sample possible values on Earth.
+
+    use exner_hyb_m, only: exner_hyb
+
+    ! For llm:
+    include "dimensions.h"
+
+    ! For ap, bp, preff:
+    include "comvert.h"
+
+    ! For kappa, cpp:
+    include "comconst.h"
+
+    ! Local:
+    integer l, i
+    integer, parameter:: ngrid = 7
+    real p(ngrid, llm + 1) ! pressure at half-level, in Pa
+    real pks(ngrid) ! exner function at the surface, in J K-1 kg-1 
+    real pk(ngrid, llm) ! exner function at full level, in J K-1 kg-1 
+    real ps(ngrid) ! surface pressure, in Pa
+    real p_lay(ngrid, llm) ! pressure at full level, in Pa
+    real delta_ps ! in Pa
+
+    !---------------------
+
+    print *, "Call sequence information: test_disvert"
+
+    delta_ps = 6e4 / (ngrid - 1)
+    ps = (/(5e4 + delta_ps * i, i = 0, ngrid - 1)/)
+    forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
+    call exner_hyb(ngrid, ps, p, pks, pk)
+    p_lay = preff * (pk / cpp)**(1. / kappa)
+
+    ! Are pressure values in the right order?
+    if (any(p(:, :llm) <= p_lay .or. p_lay <= p(:, 2:))) then
+       ! List details and stop:
+       do l = 1, llm
+          do i = 1, ngrid
+             if (p(i, l) <= p_lay(i, l)) then
+                print 1000, "p(", i, ", ", l, ") = ", p(i, l) / 100., &
+                     " hPa <= p_lay(", i, ", ", l, ") = ", &
+                     p_lay(i, l) / 100., " hPa"
+             end if
+             if (p_lay(i, l) <= p(i, l + 1)) then
+                print 1000, "p_lay(", i, ", ", l, ") = ", p_lay(i, l) / 100., &
+                     " hPa <= p(", i, ", ", l + 1, ") = ", &
+                     p(i, l + 1) / 100., " hPa"
+             end if
+          end do
+       end do
+       call abort_gcm("test_disvert", "bad order of pressure values", 1)
+    end if
+
+1000 format (2(a, 2(i0, a), g13.6), a)
+
+  end subroutine test_disvert
+
+end module test_disvert_m
