Index: LMDZ5/trunk/libf/dyn3d_common/caldyn0.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/caldyn0.F90	(revision 2335)
+++ LMDZ5/trunk/libf/dyn3d_common/caldyn0.F90	(revision 2336)
@@ -1,3 +1,3 @@
-SUBROUTINE caldyn0(itau,ucov,vcov,teta,ps,pk,phis,phi,w,time)
+SUBROUTINE caldyn0(itau,ucov,vcov,teta,ps,masse,pk,phis,phi,w,pbaru,pbarv,time)
 !
 !-------------------------------------------------------------------------------
@@ -15,19 +15,19 @@
 !===============================================================================
 ! Arguments:
-  INTEGER, INTENT(IN)  :: itau                      !--- 
+  INTEGER, INTENT(IN)  :: itau                      !--- TIME STEP INDEX
   REAL,    INTENT(IN)  :: vcov (ip1jm    ,llm)      !--- V COVARIANT WIND
   REAL,    INTENT(IN)  :: ucov (ip1jmp1  ,llm)      !--- U COVARIANT WIND
   REAL,    INTENT(IN)  :: teta (ip1jmp1  ,llm)      !--- POTENTIAL TEMPERATURE
   REAL,    INTENT(IN)  :: ps   (ip1jmp1)            !--- GROUND PRESSURE
+  REAL,    INTENT(OUT) :: masse(ip1jmp1  ,llm)      !--- MASS IN EACH CELL
   REAL,    INTENT(IN)  :: pk   (iip1,jjp1,llm)      !--- PRESSURE
   REAL,    INTENT(IN)  :: phis (ip1jmp1)            !--- GROUND GEOPOTENTIAL
   REAL,    INTENT(IN)  :: phi  (ip1jmp1  ,llm)      !--- 3D GEOPOTENTIAL
   REAL,    INTENT(OUT) :: w    (ip1jmp1  ,llm)      !--- VERTICAL WIND
+  REAL,    INTENT(OUT) :: pbaru(ip1jmp1  ,llm)      !--- U MASS FLUX
+  REAL,    INTENT(OUT) :: pbarv(ip1jm    ,llm)      !--- V MASS FLUX
   REAL,    INTENT(IN)  :: time                      !--- TIME
 !===============================================================================
 ! Local variables:
-  REAL masse(ip1jmp1  ,llm)      !--- MASS IN EACH CELL
-  REAL pbaru(ip1jmp1  ,llm)      !--- U MASS FLUX
-  REAL pbarv(ip1jm    ,llm)      !--- V MASS FLUX
   REAL, DIMENSION(ip1jmp1,llmp1) :: p
   REAL, DIMENSION(ip1jmp1,llm)   :: ucont, massebx, ang, ecin, convm, bern
Index: LMDZ5/trunk/libf/dyn3d_common/convmas.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/convmas.F	(revision 2335)
+++ 	(revision )
@@ -1,63 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE convmas (pbaru, pbarv, convm )
-c
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteurs:  P. Le Van , F. Hourdin  .
-c   -------
-c
-c   Objet:
-c   ------
-c
-c   ********************************************************************
-c   .... calcul de la convergence du flux de masse aux niveaux p ...
-c   ********************************************************************
-c
-c
-c     pbaru  et  pbarv  sont des arguments d'entree pour le s-pg  ....
-c      .....  convm      est  un argument de sortie pour le s-pg  ....
-c
-c    le calcul se fait de haut en bas, 
-c    la convergence de masse au niveau p(llm+1) est egale a 0. et
-c    n'est pas stockee dans le tableau convm .
-c
-c
-c=======================================================================
-c
-c   Declarations:
-c   -------------
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comvert.h"
-#include "logic.h"
-
-      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ),convm(  ip1jmp1,llm )
-      INTEGER   l,ij
-
-
-c-----------------------------------------------------------------------
-c    ....  calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ......
-
-      CALL  convflu( pbaru, pbarv, llm, convm )
-
-c-----------------------------------------------------------------------
-c   filtrage:
-c   ---------
-
-       CALL filtreg( convm, jjp1, llm, 2, 2, .true., 1 )
-
-c    integration de la convergence de masse de haut  en bas ......
-
-      DO      l      = llmm1, 1, -1
-        DO    ij     = 1, ip1jmp1
-         convm(ij,l) = convm(ij,l) + convm(ij,l+1)
-        ENDDO
-      ENDDO
-c
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3d_common/convmas.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/convmas.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3d_common/convmas.F90	(revision 2336)
@@ -0,0 +1,36 @@
+SUBROUTINE convmas (pbaru, pbarv, convm)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute mass flux convergence at p levels.
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+  include "logic.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: pbaru(ip1jmp1,llm)
+  REAL, INTENT(IN)  :: pbarv(ip1jm  ,llm)
+  REAL, INTENT(OUT) :: convm(ip1jmp1,llm)
+!===============================================================================
+! Method used:   Computation from top to bottom.
+!   Mass convergence at level llm is equal to zero and is not stored in convm.
+!===============================================================================
+! Local variables:
+  INTEGER :: l
+!===============================================================================
+
+!--- Computation of - (d(pbaru)/dx + d(pbarv)/dy )
+  CALL convflu( pbaru, pbarv, llm, convm )
+
+!--- Filter
+  CALL filtreg( convm, jjp1, llm, 2, 2, .TRUE., 1 )
+
+!--- Mass convergence is integrated from top to bottom
+  DO l=llmm1,1,-1
+    convm(:,l) = convm(:,l) + convm(:,l+1)
+  END DO
+
+END SUBROUTINE convmas
Index: LMDZ5/trunk/libf/dyn3d_common/enercin.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/enercin.F	(revision 2335)
+++ 	(revision )
@@ -1,98 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE enercin ( vcov, ucov, vcont, ucont, ecin )
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteur: P. Le Van
-c   -------
-c
-c   Objet:
-c   ------
-c
-c *********************************************************************
-c .. calcul de l'energie cinetique aux niveaux s  ......
-c *********************************************************************
-c  vcov, vcont, ucov et ucont sont des arguments d'entree pour le s-pg .
-c  ecin         est  un  argument de sortie pour le s-pg
-c
-c=======================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-
-      REAL vcov( ip1jm,llm ),vcont( ip1jm,llm ),
-     * ucov( ip1jmp1,llm ),ucont( ip1jmp1,llm ),ecin( ip1jmp1,llm )
-
-      REAL ecinni( iip1 ),ecinsi( iip1 )
-
-      REAL ecinpn, ecinps
-      INTEGER     l,ij,i
-
-      REAL        SSUM
-
-
-
-c                 . V
-c                i,j-1
-
-c      alpha4 .       . alpha1
-
-
-c        U .      . P     . U
-c       i-1,j    i,j      i,j
-
-c      alpha3 .       . alpha2
-
-
-c                 . V
-c                i,j
-
-c    
-c  L'energie cinetique au point scalaire P(i,j) ,autre que les poles, est :
-c       Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 )  +
-c              0.5 * U(i  ,j)**2 *( alpha1 + alpha2 )  +
-c              0.5 * V(i,j-1)**2 *( alpha1 + alpha4 )  +
-c              0.5 * V(i,  j)**2 *( alpha2 + alpha3 )
-
-
-      DO 5 l = 1,llm
-
-      DO 1  ij = iip2, ip1jm -1
-      ecin( ij+1, l )  =    0.5  *
-     * (   ucov( ij   ,l ) * ucont( ij   ,l ) * alpha3p4( ij +1 )   +
-     *     ucov( ij+1 ,l ) * ucont( ij+1 ,l ) * alpha1p2( ij +1 )   +
-     *     vcov(ij-iim,l ) * vcont(ij-iim,l ) * alpha1p4( ij +1 )   +
-     *     vcov( ij+ 1,l ) * vcont( ij+ 1,l ) * alpha2p3( ij +1 )   )
-   1  CONTINUE
-
-c    ... correction pour  ecin(1,j,l)  ....
-c    ...   ecin(1,j,l)= ecin(iip1,j,l) ...
-
-CDIR$ IVDEP
-      DO 2 ij = iip2, ip1jm, iip1
-      ecin( ij,l ) = ecin( ij + iim, l )
-   2  CONTINUE
-
-c     calcul aux poles  .......
-
-
-      DO 3 i = 1, iim
-      ecinni(i) = vcov(    i  ,  l) * vcont(    i    ,l) * aire(   i   )
-      ecinsi(i) = vcov(i+ip1jmi1,l) * vcont(i+ip1jmi1,l) * aire(i+ip1jm)
-   3  CONTINUE
-
-      ecinpn = 0.5 * SSUM( iim,ecinni,1 ) / apoln
-      ecinps = 0.5 * SSUM( iim,ecinsi,1 ) / apols
-
-      DO 4 ij = 1,iip1
-      ecin(   ij     , l ) = ecinpn
-      ecin( ij+ ip1jm, l ) = ecinps
-   4  CONTINUE
-
-   5  CONTINUE
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3d_common/enercin.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/enercin.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3d_common/enercin.F90	(revision 2336)
@@ -0,0 +1,71 @@
+SUBROUTINE enercin ( vcov, ucov, vcont, ucont, ecin )
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van.
+!-------------------------------------------------------------------------------
+! Purpose: Compute kinetic energy at sigma levels.
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: vcov    (ip1jm,  llm)
+  REAL, INTENT(IN)  :: ucov    (ip1jmp1,llm)
+  REAL, INTENT(IN)  :: vcont   (ip1jm,  llm)
+  REAL, INTENT(IN)  :: ucont   (ip1jmp1,llm)
+  REAL, INTENT(OUT) :: ecin    (ip1jmp1,llm)
+!===============================================================================
+! Notes:
+!                 . V
+!                i,j-1
+!
+!      alpha4 .       . alpha1
+!
+!
+!        U .      . P     . U
+!       i-1,j    i,j      i,j
+!
+!      alpha3 .       . alpha2
+!
+!
+!                 . V
+!                i,j
+!
+! Kinetic energy at scalar point P(i,j) (excluding poles) is:
+!       Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 )  +
+!              0.5 * U(i  ,j)**2 *( alpha1 + alpha2 )  +
+!              0.5 * V(i,j-1)**2 *( alpha1 + alpha4 )  +
+!              0.5 * V(i,  j)**2 *( alpha2 + alpha3 )
+!===============================================================================
+! Local variables:
+  INTEGER :: l, ij, i
+  REAL    :: ecinni(iip1), ecinsi(iip1), ecinpn, ecinps
+!===============================================================================
+  DO l=1,llm
+    DO ij = iip2, ip1jm -1
+      ecin(ij+1,l)=0.5*(ucov(ij    ,l)*ucont(ij    ,l)*alpha3p4(ij +1)          &
+                      + ucov(ij+1  ,l)*ucont(ij+1  ,l)*alpha1p2(ij +1)          &
+                      + vcov(ij-iim,l)*vcont(ij-iim,l)*alpha1p4(ij +1)          &
+                      + vcov(ij+1  ,l)*vcont(ij+1  ,l)*alpha2p3(ij +1) )
+    END DO
+    !--- Correction: ecin(1,j,l)= ecin(iip1,j,l)
+    DO ij=iip2,ip1jm,iip1; ecin(ij,l) = ecin(ij+iim,l); END DO
+
+    !--- North pole
+    DO i=1,iim
+      ecinni(i) = vcov(i,l)*vcont(i,l)*aire(i)
+    END DO
+    ecinpn = 0.5*SUM(ecinni(1:iim))/apoln
+    DO ij=1,iip1; ecin(ij,l)=ecinpn; END DO
+
+    !--- South pole
+    DO i=1,iim
+      ecinsi(i) = vcov(i+ip1jmi1,l)*vcont(i+ip1jmi1,l)*aire(i+ip1jm)
+    END DO
+    ecinps = 0.5*SUM(ecinsi(1:iim))/apols
+    DO ij=1,iip1; ecin(ij+ip1jm,l)=ecinps; END DO
+  END DO
+
+END SUBROUTINE enercin
+
Index: LMDZ5/trunk/libf/dyn3d_common/flumass.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/flumass.F	(revision 2335)
+++ 	(revision )
@@ -1,109 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE flumass (massebx,masseby, vcont, ucont, pbaru, pbarv )
-
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteurs:  P. Le Van, F. Hourdin  .
-c   -------
-c
-c   Objet:
-c   ------
-c
-c *********************************************************************
-c     .... calcul du flux de masse  aux niveaux s ......
-c *********************************************************************
-c   massebx,masseby,vcont et ucont sont des argum. d'entree pour le s-pg .
-c       pbaru  et pbarv            sont des argum.de sortie pour le s-pg .
-c
-c=======================================================================
-
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-
-      REAL massebx( ip1jmp1,llm ),masseby( ip1jm,llm ) ,
-     * vcont( ip1jm,llm ),ucont( ip1jmp1,llm ),pbaru( ip1jmp1,llm ),
-     * pbarv( ip1jm,llm )
-
-      REAL apbarun( iip1 ),apbarus( iip1 )
-
-      REAL sairen,saireun,saires,saireus,ctn,cts,ctn0,cts0
-      INTEGER  l,ij,i
-
-      REAL       SSUM
-
-
-      DO  5 l = 1,llm
-
-      DO  1 ij = iip2,ip1jm
-      pbaru( ij,l ) = massebx( ij,l ) * ucont( ij,l )
-   1  CONTINUE
-
-      DO 3 ij = 1,ip1jm
-      pbarv( ij,l ) = masseby( ij,l ) * vcont( ij,l )
-   3  CONTINUE
-
-   5  CONTINUE
-
-c    ................................................................
-c     calcul de la composante du flux de masse en x aux poles .......
-c    ................................................................
-c     par la resolution d'1 systeme de 2 equations .
-
-c     la premiere equat.decrivant le calcul de la divergence en 1 point i
-c     du pole,ce calcul etant itere de i=1 a i=im .
-c                 c.a.d   ,
-c     ( ( 0.5*pbaru(i)-0.5*pbaru(i-1) - pbarv(i))/aire(i)   =
-c                                           - somme de ( pbarv(n) )/aire pole
-
-c     l'autre equat.specifiant que la moyenne du flux de masse au pole est =0.
-c     c.a.d    somme de pbaru(n)*aire locale(n) = 0.
-
-c     on en revient ainsi a determiner la constante additive commune aux pbaru
-c     qui representait pbaru(0,j,l) dans l'equat.du calcul de la diverg.au pt
-c     i=1 .
-c     i variant de 1 a im
-c     n variant de 1 a im
-
-      sairen = SSUM( iim,  aire(   1     ), 1 )
-      saireun= SSUM( iim, aireu(   1     ), 1 )
-      saires = SSUM( iim,  aire( ip1jm+1 ), 1 )
-      saireus= SSUM( iim, aireu( ip1jm+1 ), 1 )
-
-      DO 20 l = 1,llm
-
-      ctn =  SSUM( iim, pbarv(    1     ,l),  1 )/ sairen
-      cts =  SSUM( iim, pbarv(ip1jmi1+ 1,l),  1 )/ saires
-
-      pbaru(    1   ,l )=   pbarv(    1     ,l ) - ctn * aire(    1    )
-      pbaru( ip1jm+1,l )= - pbarv( ip1jmi1+1,l ) + cts * aire( ip1jm+1 )
-
-      DO 11 i = 2,iim
-      pbaru(    i    ,l ) = pbaru(   i - 1   ,l )    +
-     *                      pbarv(    i      ,l ) - ctn * aire(   i    )
-
-      pbaru( i+ ip1jm,l ) = pbaru( i+ ip1jm-1,l )    -
-     *                      pbarv( i+ ip1jmi1,l ) + cts * aire(i+ ip1jm)
-  11  CONTINUE
-      DO 12 i = 1,iim
-      apbarun(i) = aireu(    i   ) * pbaru(   i    , l)
-      apbarus(i) = aireu(i +ip1jm) * pbaru(i +ip1jm, l)
-  12  CONTINUE
-      ctn0 = -SSUM( iim,apbarun,1 )/saireun
-      cts0 = -SSUM( iim,apbarus,1 )/saireus
-      DO 14 i = 1,iim
-      pbaru(   i    , l) = 2. * ( pbaru(   i    , l) + ctn0 )
-      pbaru(i+ ip1jm, l) = 2. * ( pbaru(i +ip1jm, l) + cts0 )
-  14  CONTINUE
-
-      pbaru(   iip1 ,l ) = pbaru(    1    ,l )
-      pbaru( ip1jmp1,l ) = pbaru( ip1jm +1,l )
-  20  CONTINUE
-
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3d_common/flumass.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/flumass.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3d_common/flumass.F90	(revision 2336)
@@ -0,0 +1,76 @@
+SUBROUTINE flumass (massebx,masseby, vcont, ucont, pbaru, pbarv )
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute mass flux at s levels.
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: massebx(ip1jmp1,llm)
+  REAL, INTENT(IN)  :: masseby(ip1jm  ,llm)
+  REAL, INTENT(IN)  :: vcont  (ip1jm  ,llm)
+  REAL, INTENT(IN)  :: ucont  (ip1jmp1,llm)
+  REAL, INTENT(OUT) :: pbaru  (ip1jmp1,llm)
+  REAL, INTENT(OUT) :: pbarv  (ip1jm  ,llm)
+!===============================================================================
+! Method used:   A 2 equations system is solved.
+!   * 1st one describes divergence computation at pole point nr. i (i=1 to im):
+!     (0.5*(pbaru(i)-pbaru(i-1))-pbarv(i))/aire(i) = - SUM(pbarv(n))/aire pole
+!   * 2nd one specifies that mean mass flux at pole is equal to 0:
+!     SUM(pbaru(n)*local_area(n))=0
+! This way, we determine additive constant common to pbary elements representing
+!   pbaru(0,j,l) in divergence computation equation for point i=1. (i=1 to im)
+!===============================================================================
+! Local variables:
+  REAL    :: sairen, saireun, ctn, ctn0, apbarun(iip1)
+  REAL    :: saires, saireus, cts, cts0, apbarus(iip1)
+  INTEGER :: l, i
+!===============================================================================
+  DO l=1,llm
+    pbaru(iip2:ip1jm,l)=massebx(iip2:ip1jm,l)*ucont(iip2:ip1jm,l)
+    pbarv(   1:ip1jm,l)=masseby(   1:ip1jm,l)*vcont(   1:ip1jm,l)
+  END DO
+
+  !--- NORTH POLE
+  sairen =SUM(aire (1:iim))
+  saireun=SUM(aireu(1:iim))
+  DO l = 1,llm
+    ctn=SUM(pbarv(1:iim,l))/sairen
+    pbaru(1,l)= pbarv(1,l)-ctn*aire(1)
+    DO i=2,iim
+      pbaru(i,l)=pbaru(i-1,l)+pbarv(i,l)-ctn*aire(i)
+    END DO
+    DO i=1,iim
+      apbarun(i)=aireu(i)*pbaru(i,l)
+    END DO
+    ctn0 = -SUM(apbarun(1:iim))/saireun
+    DO i = 1,iim
+      pbaru(i,l)=2.*(pbaru(i,l)+ctn0)
+    END DO
+    pbaru(iip1,l)=pbaru(1,l)
+  END DO
+
+  !--- SOUTH POLE
+  saires =SUM(aire (ip1jm+1:ip1jmp1-1))
+  saireus=SUM(aireu(ip1jm+1:ip1jmp1-1))
+  DO l = 1,llm
+    cts=SUM(pbarv(ip1jmi1+1:ip1jm-1,l))/saires
+    pbaru(1+ip1jm,l)=-pbarv(1+ip1jmi1,l)+cts*aire(1+ip1jm)
+    DO i=2,iim
+      pbaru(i+ip1jm,l)=pbaru(i-1+ip1jm,l)-pbarv(i+ip1jmi1,l)+cts*aire(i+ip1jm)
+    END DO
+    DO i=1,iim
+      apbarus(i)=aireu(i+ip1jm)*pbaru(i+ip1jm,l)
+    END DO
+    cts0 = -SUM(apbarus(1:iim))/saireus
+    DO i = 1,iim
+      pbaru(i+ip1jm,l)=2.*(pbaru(i+ip1jm,l)+cts0)
+    END DO
+    pbaru(ip1jmp1,l)=pbaru(1+ip1jm,l)
+  END DO
+
+END SUBROUTINE flumass
Index: LMDZ5/trunk/libf/dyn3d_common/massbar.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/massbar.F	(revision 2335)
+++ 	(revision )
@@ -1,102 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE massbar(  masse, massebx, masseby )
-      IMPLICIT NONE
-c
-c **********************************************************************
-c
-c  Calcule les moyennes en x et  y de la masse d'air dans chaque maille.
-c **********************************************************************
-c    Auteurs : P. Le Van , Fr. Hourdin  .
-c   ..........
-c
-c  ..  masse                 est  un argum. d'entree  pour le s-pg ...
-c  ..  massebx,masseby      sont des argum. de sortie pour le s-pg ...
-c     
-c
-c     IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-c
-      REAL    masse( ip1jmp1,llm ), massebx( ip1jmp1,llm )  ,
-     *      masseby(   ip1jm,llm )
-      INTEGER ij,l
-c
-c
-c   Methode pour calculer massebx et masseby .
-c   ----------------------------------------
-c
-c    A chaque point scalaire P (i,j) est affecte 4 coefficients d'aires
-c       alpha1(i,j)  calcule  au point ( i+1/4,j-1/4 )
-c       alpha2(i,j)  calcule  au point ( i+1/4,j+1/4 )
-c       alpha3(i,j)  calcule  au point ( i-1/4,j+1/4 )
-c       alpha4(i,j)  calcule  au point ( i-1/4,j-1/4 )
-c
-c    Avec  alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j)        
-c
-c    N.B .  Pour plus de details, voir s-pg  ...  iniconst ...
-c
-c
-c
-c   alpha4 .         . alpha1    . alpha4
-c    (i,j)             (i,j)       (i+1,j)
-c
-c             P .        U .          . P
-c           (i,j)       (i,j)         (i+1,j)
-c
-c   alpha3 .         . alpha2    .alpha3 
-c    (i,j)              (i,j)     (i+1,j)
-c
-c             V .        Z .          . V
-c           (i,j)
-c
-c   alpha4 .         . alpha1    .alpha4
-c   (i,j+1)            (i,j+1)   (i+1,j+1) 
-c
-c             P .        U .          . P
-c          (i,j+1)                    (i+1,j+1)
-c
-c
-c
-c                       On  a :
-c
-c    massebx(i,j) = masse(i  ,j) * ( alpha1(i  ,j) + alpha2(i,j))   +
-c                   masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) )
-c     localise  au point  ... U (i,j) ...
-c
-c    masseby(i,j) = masse(i,j  ) * ( alpha2(i,j  ) + alpha3(i,j  )  +
-c                   masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1)  
-c     localise  au point  ... V (i,j) ...
-c
-c
-c=======================================================================
-
-      DO   100    l = 1 , llm
-c
-        DO  ij = 1, ip1jmp1 - 1
-         massebx(ij,l) =  masse( ij, l) * alpha1p2( ij  )     + 
-     *                   masse(ij+1, l) * alpha3p4(ij+1 )
-        ENDDO
-
-c    .... correction pour massebx( iip1,j) .....
-c    ...    massebx(iip1,j)= massebx(1,j) ...
-c
-CDIR$ IVDEP
-        DO  ij = iip1, ip1jmp1, iip1
-         massebx( ij,l ) = massebx( ij - iim,l )
-        ENDDO
-
-
-         DO  ij = 1,ip1jm
-         masseby( ij,l ) = masse(  ij   , l ) * alpha2p3(   ij    )  +
-     *                     masse(ij+iip1, l ) * alpha1p4( ij+iip1 )
-         ENDDO
-
-100   CONTINUE
-c
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3d_common/massbar.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/massbar.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3d_common/massbar.F90	(revision 2336)
@@ -0,0 +1,66 @@
+SUBROUTINE massbar(masse,massebx,masseby)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute air mass mean along X and Y in each cell.
+! See iniconst for more details.
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: masse  (ip1jmp1,llm)
+  REAL, INTENT(OUT) :: massebx(ip1jmp1,llm)
+  REAL, INTENT(OUT) :: masseby(ip1jm  ,llm)
+!-------------------------------------------------------------------------------
+! Method used. Each scalar point is associated to 4 area coefficients:
+!    * alpha1(i,j) at point ( i+1/4,j-1/4 )
+!    * alpha2(i,j) at point ( i+1/4,j+1/4 )
+!    * alpha3(i,j) at point ( i-1/4,j+1/4 )
+!    * alpha4(i,j) at point ( i-1/4,j-1/4 )
+! where alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j)
+!
+!   alpha4 .         . alpha1    . alpha4
+!    (i,j)             (i,j)       (i+1,j)
+!
+!             P .        U .          . P
+!           (i,j)       (i,j)         (i+1,j)
+!
+!   alpha3 .         . alpha2    .alpha3 
+!    (i,j)              (i,j)     (i+1,j)
+!
+!             V .        Z .          . V
+!           (i,j)
+!
+!   alpha4 .         . alpha1    .alpha4
+!   (i,j+1)            (i,j+1)   (i+1,j+1) 
+!
+!             P .        U .          . P
+!          (i,j+1)                    (i+1,j+1)
+!
+!
+!    massebx(i,j) = masse(i  ,j) * ( alpha1(i  ,j) + alpha2(i,j))   +
+!                   masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) )
+!     localized at point  ... U (i,j) ...
+!
+!    masseby(i,j) = masse(i,j  ) * ( alpha2(i,j  ) + alpha3(i,j  )  +
+!                   masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1)  
+!     localized at point  ... V (i,j) ...
+!===============================================================================
+! Local variables:
+  INTEGER :: ij, l
+!===============================================================================
+  DO l=1,llm
+    DO ij=1,ip1jmp1-1
+      massebx(ij,l)=masse(ij,l)*alpha1p2(ij)+masse(ij+1   ,l)*alpha3p4(ij+1)
+    END DO
+    DO ij=iip1,ip1jmp1,iip1; massebx(ij,l)=massebx(ij-iim,l); END DO
+    DO ij=1,ip1jm
+      masseby(ij,l)=masse(ij,l)*alpha2p3(ij)+masse(ij+iip1,l)*alpha1p4(ij+iip1)
+    END DO
+  END DO
+
+END SUBROUTINE massbar
+
Index: LMDZ5/trunk/libf/dyn3d_common/massbarxy.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/massbarxy.F	(revision 2335)
+++ 	(revision )
@@ -1,49 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE massbarxy(  masse, massebxy )
-      IMPLICIT NONE
-c
-c **********************************************************************
-c
-c  Calcule les moyennes en x et  y de la masse d'air dans chaque maille.
-c **********************************************************************
-c    Auteurs : P. Le Van , Fr. Hourdin  .
-c   ..........
-c
-c  ..  masse          est  un  argum. d'entree  pour le s-pg ...
-c  ..  massebxy       est  un  argum. de sortie pour le s-pg ...
-c     
-c
-c     IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-c
-       REAL  masse( ip1jmp1,llm ), massebxy( ip1jm,llm )
-       INTEGER ij,l
-c
-
-      DO   100    l = 1 , llm
-c
-      DO 5 ij = 1, ip1jm - 1
-      massebxy( ij,l ) = masse(    ij  ,l ) * alpha2(   ij    )   +
-     +                   masse(   ij+1 ,l ) * alpha3(  ij+1   )   +
-     +                   masse( ij+iip1,l ) * alpha1( ij+iip1 )   +
-     +                   masse( ij+iip2,l ) * alpha4( ij+iip2 )
-   5  CONTINUE
-
-c    ....  correction pour     massebxy( iip1,j )  ........
-
-CDIR$ IVDEP
-
-      DO 7 ij = iip1, ip1jm, iip1
-      massebxy( ij,l ) = massebxy( ij - iim,l )
-   7  CONTINUE
-
-100   CONTINUE
-c
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3d_common/massbarxy.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/massbarxy.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3d_common/massbarxy.F90	(revision 2336)
@@ -0,0 +1,31 @@
+SUBROUTINE massbarxy(masse,massebxy)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute air mass mean along X and Y in each cell.
+! See iniconst for more details.
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comconst.h"
+  include "comgeom.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: masse   (ip1jmp1,llm)
+  REAL, INTENT(OUT) :: massebxy(ip1jm  ,llm)
+!===============================================================================
+! Local variables:
+  INTEGER :: ij, l
+!===============================================================================
+  DO l=1,llm
+    DO ij=1,ip1jm-1
+      massebxy(ij,l)=masse(ij     ,l)*alpha2(ij     ) + &
+     +               masse(ij+1   ,l)*alpha3(ij+1   ) + &
+     +               masse(ij+iip1,l)*alpha1(ij+iip1) + &
+     +               masse(ij+iip2,l)*alpha4(ij+iip2)
+    END DO
+    DO ij=iip1,ip1jm,iip1; massebxy(ij,l)=massebxy(ij-iim,l); END DO
+  END DO
+
+END SUBROUTINE massbarxy
Index: LMDZ5/trunk/libf/dyn3d_common/tourpot.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/tourpot.F	(revision 2335)
+++ 	(revision )
@@ -1,81 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE tourpot ( vcov, ucov, massebxy, vorpot )
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteur:  P. Le Van
-c   -------
-c
-c   Objet:
-c   ------
-c
-c    *******************************************************************
-c    .........      calcul du tourbillon potentiel             .........
-c    *******************************************************************
-c
-c     vcov,ucov,fext et pbarxyfl sont des argum. d'entree pour le s-pg .
-c             vorpot            est  un argum.de sortie pour le s-pg .
-c
-c=======================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "logic.h"
-
-      REAL  rot( ip1jm,llm )
-      REAL vcov( ip1jm,llm ),ucov( ip1jmp1,llm )
-      REAL massebxy( ip1jm,llm ),vorpot( ip1jm,llm )
-
-      INTEGER l, ij
-
-
-
-
-c  ... vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy ..
-
-
-
-c    ........  Calcul du rotationnel du vent V  puis filtrage  ........
-
-      DO 5 l = 1,llm
-
-      DO 2 ij = 1, ip1jm - 1
-      rot( ij,l ) = vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)
-   2  CONTINUE
-
-c    ....  correction pour  rot( iip1,j,l )  .....
-c    ....     rot(iip1,j,l) = rot(1,j,l)    .....
-
-CDIR$ IVDEP
-
-      DO 3 ij = iip1, ip1jm, iip1
-      rot( ij,l ) = rot( ij -iim, l )
-   3  CONTINUE
-
-   5  CONTINUE
-
-
-      CALL  filtreg( rot, jjm, llm, 2, 1, .FALSE., 1 )
-
-
-      DO 10 l = 1, llm
-
-      DO 6 ij = 1, ip1jm - 1
-      vorpot( ij,l ) = ( rot(ij,l) + fext(ij) ) / massebxy(ij,l)
-   6  CONTINUE
-
-c    ..... correction pour  vorpot( iip1,j,l)  .....
-c    ....   vorpot(iip1,j,l)= vorpot(1,j,l) ....
-CDIR$ IVDEP
-      DO 8 ij = iip1, ip1jm, iip1
-      vorpot( ij,l ) = vorpot( ij -iim,l )
-   8  CONTINUE
-
-  10  CONTINUE
-
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3d_common/tourpot.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/tourpot.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3d_common/tourpot.F90	(revision 2336)
@@ -0,0 +1,46 @@
+SUBROUTINE tourpot ( vcov, ucov, massebxy, vorpot )
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van.
+!-------------------------------------------------------------------------------
+! Purpose: Compute potential vorticity.
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+  include "logic.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: vcov    (ip1jm,  llm)
+  REAL, INTENT(IN)  :: ucov    (ip1jmp1,llm)
+  REAL, INTENT(IN)  :: massebxy(ip1jm,  llm)
+  REAL, INTENT(OUT) :: vorpot  (ip1jm,  llm)
+!===============================================================================
+! Method used:
+!   vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy
+!===============================================================================
+! Local variables:
+  INTEGER :: l, ij
+  REAL    :: rot(ip1jm,llm)
+!===============================================================================
+
+!--- Wind vorticity ; correction: rot(iip1,j,l) = rot(1,j,l)
+  DO l=1,llm
+    DO ij=1,ip1jm-1
+      rot(ij,l)=vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)
+    END DO
+    DO ij=iip1,ip1jm,iip1; rot(ij,l)=rot(ij-iim,l); END DO
+  END DO
+
+!--- Filter
+  CALL  filtreg(rot,jjm,llm,2,1,.FALSE.,1)
+
+!--- Potential vorticity ; correction: rot(iip1,j,l) = rot(1,j,l)
+  DO l=1,llm
+    DO ij=1,ip1jm-1
+      vorpot(ij,l)=(rot(ij,l)+fext(ij))/massebxy(ij,l)
+    END DO
+    DO ij=iip1,ip1jm,iip1; vorpot(ij,l)=vorpot(ij-iim,l); END DO
+  END DO
+
+END SUBROUTINE tourpot
Index: LMDZ5/trunk/libf/dyn3d_common/vitvert.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/vitvert.F	(revision 2335)
+++ 	(revision )
@@ -1,52 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE vitvert ( convm , w )
-c
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteurs:  P. Le Van , F. Hourdin .
-c   -------
-c
-c   Objet:
-c   ------
-c
-c    *******************************************************************
-c  .... calcul de la vitesse verticale aux niveaux sigma  ....
-c    *******************************************************************
-c     convm   est un argument  d'entree pour le s-pg  ......
-c       w     est un argument de sortie pour le s-pg  ......
-c
-c    la vitesse verticale est orientee de  haut en bas .
-c    au sol, au niveau sigma(1),   w(i,j,1) = 0.
-c    au sommet, au niveau sigma(llm+1) , la vit.verticale est aussi
-c    egale a 0. et n'est pas stockee dans le tableau w  .
-c
-c
-c=======================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comvert.h"
-
-      REAL w(ip1jmp1,llm),convm(ip1jmp1,llm)
-      INTEGER   l, ij
-
-
-
-      DO 2  l = 1,llmm1
-
-      DO 1 ij = 1,ip1jmp1
-      w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
-   1  CONTINUE
-
-   2  CONTINUE
-
-      DO 5 ij  = 1,ip1jmp1
-      w(ij,1)  = 0.
-5     CONTINUE
-
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3d_common/vitvert.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d_common/vitvert.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3d_common/vitvert.F90	(revision 2336)
@@ -0,0 +1,27 @@
+SUBROUTINE vitvert (convm, w)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute vertical speed at sigma levels.
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comvert.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: convm(ip1jmp1,llm)
+  REAL, INTENT(OUT) :: w    (ip1jmp1,llm)
+!===============================================================================
+! Notes: Vertical speed is oriented from bottom to top.
+!   * At ground - level sigma(1):     w(i,j,1) = 0.
+!   * At top    - level sigma(llm+1): w(i,j,l) = 0. (not stored in w)
+!===============================================================================
+! Local variables:
+  INTEGER :: l
+!===============================================================================
+  DO l=1,llmm1; w(:,l+1)=convm(:,l+1)-bp(l+1)*convm(:,1); END DO
+  w(:,1)=0.
+
+END SUBROUTINE vitvert
+
Index: LMDZ5/trunk/libf/dyn3dmem/ce0l.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/ce0l.F90	(revision 2335)
+++ 	(revision )
@@ -1,134 +1,0 @@
-!
-! $Id$
-!
-!-------------------------------------------------------------------------------
-!
-PROGRAM ce0l
-!-------------------------------------------------------------------------------
-! Purpose: Calls etat0, creates initial states and limit_netcdf
-!
-!     interbar=.T. for barycentric interpolation inter_barxy
-!     extrap  =.T. for data extrapolation, like for the SSTs when file does not
-!                  contain ocean points only.
-!     oldice  =.T. for old-style ice, obtained using grille_m (grid_atob).
-!     masque is created in etat0, passed to limit to ensure consistancy.
-!-------------------------------------------------------------------------------
-  USE control_mod
-#ifdef CPP_EARTH
-! This prog. is designed to work for Earth
-  USE dimphy
-  USE comgeomphy
-  USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
-  USE mod_const_mpi
-  USE infotrac
-  USE parallel_lmdz, ONLY: finalize_parallel  
-  USE indice_sol_mod
-
-#ifdef CPP_IOIPSL
-  USE ioipsl, ONLY: ioconf_calendar
-#endif
-
-#endif
-  IMPLICIT NONE
-#ifndef CPP_EARTH
-#include "iniprint.h"
-  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
-#else
-!-------------------------------------------------------------------------------
-! Local variables:
-  LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.
-#include "dimensions.h"
-#include "paramet.h"
-!#include "indicesol.h"
-#include "iniprint.h"
-#include "temps.h"
-#include "logic.h"
-#ifdef CPP_MPI
-      include 'mpif.h'
-#endif
-
-  INTEGER                       :: ierr
-  REAL,    DIMENSION(iip1,jjp1) :: masque
-  CHARACTER(LEN=15)             :: calnd
-  REAL,    DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol
-!-------------------------------------------------------------------------------
-  CALL conf_gcm( 99, .TRUE. )
-
-#ifdef CPP_MPI
-  CALL init_mpi
-#endif
-
-  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
-  WRITE(lunout,*)'---> klon=',klon
-  IF (mpi_size>1 .OR. omp_size>1) THEN
-       CALL abort_gcm('ce0l','In parallel mode,                         &
- &                 ce0l must be called only                             &
- &                 for 1 process and 1 task',1)
-  ENDIF
-
-  CALL InitComgeomphy
-
-#ifdef CPP_IOIPSL
-  SELECT CASE(calend)
-    CASE('earth_360d');CALL ioconf_calendar('360d');      calnd='a 360 jours/an'
-    CASE('earth_365d');CALL ioconf_calendar('noleap');    calnd='a 365 jours/an'
-    CASE('earth_366d');CALL ioconf_calendar('366d');      calnd='bissextile'
-    CASE('gregorian'); CALL ioconf_calendar('gregorian'); calnd='gregorien'
-    CASE('standard');  CALL ioconf_calendar('gregorian'); calnd='gregorien'
-    CASE('julian');    CALL ioconf_calendar('julian');    calnd='julien'
-    CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
-  !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
-    CASE DEFAULT
-      CALL abort_gcm('ce0l','Mauvais choix de calendrier',1)
-  END SELECT
-  WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre '//TRIM(calnd)
-#endif
-
-  IF (type_trac == 'inca') THEN
-#ifdef INCA
-      CALL init_const_lmdz( &
-         nbtr,anneeref,dayref,&
-         iphysiq,day_step,nday,& 
-         nbsrf, is_oce,is_sic,&
-         is_ter,is_lic)
-      
-#endif
-  END IF
-
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,*) '  ***  etat0_netcdf ***  '
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) ' interbar = ',interbar
-  CALL etat0_netcdf(interbar,masque,phis,ok_etat0)
-
-  IF(ok_limit) THEN
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,*) '  ***  Limit_netcdf ***  '
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,'(//)')
-  CALL limit_netcdf(interbar,extrap,oldice,masque)
-  END IF
-
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  ***************************  '
-  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
-  WRITE(lunout,*) '  ***************************  '
-  WRITE(lunout,'(//)')
-  CALL grilles_gcm_netcdf_sub(masque,phis)
-  
-#ifdef CPP_MPI
-!$OMP MASTER
-  CALL MPI_FINALIZE(ierr)
-  IF (ierr /= 0) CALL abort_gcm('ce0l','Error in MPI_FINALIZE',1)
-!$OMP END MASTER
-#endif
-
-#endif
-! of #ifndef CPP_EARTH #else
-
-END PROGRAM ce0l
-!
-!-------------------------------------------------------------------------------
Index: LMDZ5/trunk/libf/dyn3dmem/convmas1_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/convmas1_loc.F	(revision 2335)
+++ 	(revision )
@@ -1,64 +1,0 @@
-      SUBROUTINE convmas1_loc (pbaru, pbarv, convm )
-c
-      USE parallel_lmdz
-      USE mod_filtreg_p
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteurs:  P. Le Van , F. Hourdin  .
-c   -------
-c
-c   Objet:
-c   ------
-c
-c   ********************************************************************
-c   .... calcul de la convergence du flux de masse aux niveaux p ...
-c   ********************************************************************
-c
-c
-c     pbaru  et  pbarv  sont des arguments d'entree pour le s-pg  ....
-c      .....  convm      est  un argument de sortie pour le s-pg  ....
-c
-c    le calcul se fait de haut en bas, 
-c    la convergence de masse au niveau p(llm+1) est egale a 0. et
-c    n'est pas stockee dans le tableau convm .
-c
-c
-c=======================================================================
-c
-c   Declarations:
-c   -------------
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comvert.h"
-#include "logic.h"
-
-      REAL pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm )
-      REAL, target :: convm(  ijb_u:ije_u,llm )
-      INTEGER   l,ij
-
-      INTEGER ijb,ije,jjb,jje
- 
-      
-c-----------------------------------------------------------------------
-c    ....  calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ......
-
-      CALL  convflu_loc( pbaru, pbarv, llm, convm )
-
-c-----------------------------------------------------------------------
-c   filtrage:
-c   ---------
-       
-       jjb=jj_begin
-       jje=jj_end+1
-       if (pole_sud) jje=jj_end
- 
-       CALL filtreg_p( convm, jjb_u,jje_u,jjb, jje, jjp1, llm,
-     &                 2, 2, .true., 1 )
-
-c    integration de la convergence de masse de haut  en bas ......
-c
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/convmas1_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/convmas1_loc.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3dmem/convmas1_loc.F90	(revision 2336)
@@ -0,0 +1,38 @@
+SUBROUTINE convmas1_loc (pbaru, pbarv, convm)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute mass flux convergence at p levels.
+!          Equivalent to convmas_loc if convmas2_loc is called after.
+  USE parallel_lmdz
+  USE mod_filtreg_p
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+  include "logic.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: pbaru(ijb_u:ije_u,llm)
+  REAL, INTENT(IN)  :: pbarv(ijb_v:ije_v,llm)
+  REAL, TARGET, INTENT(OUT) :: convm(ijb_u:ije_u,llm)
+!===============================================================================
+! Method used:   Computation from top to bottom.
+!   Mass convergence at level llm is equal to zero and is not stored in convm.
+!===============================================================================
+! Local variables:
+  INTEGER :: l, jjb, jje
+!===============================================================================
+
+!--- Computation of - (d(pbaru)/dx + d(pbarv)/dy )
+  CALL convflu_loc( pbaru, pbarv, llm, convm )
+
+!--- Filter
+  jjb=jj_begin
+  jje=jj_end+1
+  IF(pole_sud) jje=jj_end
+  CALL filtreg_p(convm,jjb_u,jje_u,jjb,jje,jjp1,llm,2,2,.TRUE.,1)
+
+END SUBROUTINE convmas1_loc
+
Index: LMDZ5/trunk/libf/dyn3dmem/convmas2_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/convmas2_loc.F	(revision 2335)
+++ 	(revision )
@@ -1,56 +1,0 @@
-      SUBROUTINE convmas2_loc ( convm )
-c
-      USE parallel_lmdz
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteurs:  P. Le Van , F. Hourdin  .
-c   -------
-c
-c   Objet:
-c   ------
-c
-c   ********************************************************************
-c   .... calcul de la convergence du flux de masse aux niveaux p ...
-c   ********************************************************************
-c
-c
-c     pbaru  et  pbarv  sont des arguments d'entree pour le s-pg  ....
-c      .....  convm      est  un argument de sortie pour le s-pg  ....
-c
-c    le calcul se fait de haut en bas, 
-c    la convergence de masse au niveau p(llm+1) est egale a 0. et
-c    n'est pas stockee dans le tableau convm .
-c
-c
-c=======================================================================
-c
-c   Declarations:
-c   -------------
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comvert.h"
-#include "logic.h"
-
-      REAL pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm )
-      REAL :: convm(  ijb_u:ije_u,llm )
-      INTEGER   l,ij
-      INTEGER ijb,ije,jjb,jje
- 
-c$OMP MASTER
-c    integration de la convergence de masse de haut  en bas ......
-       ijb=ij_begin
-       ije=ij_end+iip1
-       if (pole_sud) ije=ij_end
-            
-      DO      l      = llmm1, 1, -1
-        DO    ij     = ijb, ije
-         convm(ij,l) = convm(ij,l) + convm(ij,l+1)
-        ENDDO
-      ENDDO
-c
-c$OMP END MASTER
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/convmas2_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/convmas2_loc.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3dmem/convmas2_loc.F90	(revision 2336)
@@ -0,0 +1,36 @@
+SUBROUTINE convmas2_loc (convm)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute mass flux convergence at p levels.
+!          Equivalent to convmas_loc if convmas1_loc is called before.
+  USE parallel_lmdz
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+  include "logic.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(INOUT) :: convm(ijb_u:ije_u,llm)
+!===============================================================================
+! Method used:   Computation from top to bottom.
+!   Mass convergence at level llm is equal to zero and is not stored in convm.
+!===============================================================================
+! Local variables:
+  INTEGER :: l, ijb, ije
+!===============================================================================
+
+!$OMP MASTER
+!--- Mass convergence is integrated from top to bottom
+  ijb=ij_begin
+  ije=ij_end+iip1
+  IF(pole_sud) ije=ij_end
+  DO l=llmm1,1,-1
+    convm(ijb:ije,l) = convm(ijb:ije,l) + convm(ijb:ije,l+1)
+  END DO
+!$OMP END MASTER
+
+END SUBROUTINE convmas2_loc
+
Index: LMDZ5/trunk/libf/dyn3dmem/convmas_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/convmas_loc.F	(revision 2335)
+++ 	(revision )
@@ -1,77 +1,0 @@
-      SUBROUTINE convmas_loc (pbaru, pbarv, convm )
-c
-      USE parallel_lmdz
-      USE mod_filtreg_p
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteurs:  P. Le Van , F. Hourdin  .
-c   -------
-c
-c   Objet:
-c   ------
-c
-c   ********************************************************************
-c   .... calcul de la convergence du flux de masse aux niveaux p ...
-c   ********************************************************************
-c
-c
-c     pbaru  et  pbarv  sont des arguments d'entree pour le s-pg  ....
-c      .....  convm      est  un argument de sortie pour le s-pg  ....
-c
-c    le calcul se fait de haut en bas, 
-c    la convergence de masse au niveau p(llm+1) est egale a 0. et
-c    n'est pas stockee dans le tableau convm .
-c
-c
-c=======================================================================
-c
-c   Declarations:
-c   -------------
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comvert.h"
-#include "logic.h"
-
-      REAL pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm )
-      REAL, target :: convm(  ijb_u:ije_u,llm )
-      INTEGER   l,ij
-
-      INTEGER ijb,ije,jjb,jje
- 
-      
-c-----------------------------------------------------------------------
-c    ....  calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ......
-
-      CALL  convflu_loc( pbaru, pbarv, llm, convm )
-
-c-----------------------------------------------------------------------
-c   filtrage:
-c   ---------
-       
-       jjb=jj_begin
-       jje=jj_end+1
-       if (pole_sud) jje=jj_end
- 
-       CALL filtreg_p(convm, jjb_u, jje_u,jjb, jje, jjp1, llm,
-     &                 2, 2, .true., 1 )
-
-c    integration de la convergence de masse de haut  en bas ......
-!$OMP BARRIER
-!$OMP MASTER
-       ijb=ij_begin
-       ije=ij_end+iip1
-       if (pole_sud) ije=ij_end
-            
-      DO      l      = llmm1, 1, -1
-        DO    ij     = ijb, ije
-         convm(ij,l) = convm(ij,l) + convm(ij,l+1)
-        ENDDO
-      ENDDO
-c
-!$OMP END MASTER
-!$OMP BARRIER
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/convmas_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/convmas_loc.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3dmem/convmas_loc.F90	(revision 2336)
@@ -0,0 +1,49 @@
+SUBROUTINE convmas_loc (pbaru, pbarv, convm)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute mass flux convergence at p levels.
+  USE parallel_lmdz
+  USE mod_filtreg_p
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+  include "logic.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: pbaru(ijb_u:ije_u,llm)
+  REAL, INTENT(IN)  :: pbarv(ijb_v:ije_v,llm)
+  REAL, INTENT(OUT) :: convm(ijb_u:ije_u,llm)
+!===============================================================================
+! Method used:   Computation from top to bottom.
+!   Mass convergence at level llm is equal to zero and is not stored in convm.
+!===============================================================================
+! Local variables:
+  INTEGER :: l, ijb, ije, jjb, jje
+!===============================================================================
+
+!--- Computation of - (d(pbaru)/dx + d(pbarv)/dy )
+  CALL convflu_loc( pbaru, pbarv, llm, convm )
+
+!--- Filter
+  jjb=jj_begin
+  jje=jj_end+1
+  IF(pole_sud) jje=jj_end
+  CALL filtreg_p(convm,jjb_u,jje_u,jjb,jje,jjp1,llm,2,2,.TRUE.,1)
+
+!--- Mass convergence is integrated from top to bottom
+!$OMP BARRIER
+!$OMP MASTER
+  ijb=ij_begin
+  ije=ij_end+iip1
+  IF(pole_sud) ije=ij_end
+  DO l=llmm1,1,-1
+    convm(ijb:ije,l) = convm(ijb:ije,l) + convm(ijb:ije,l+1)
+  END DO
+!$OMP END MASTER
+!$OMP BARRIER
+
+END SUBROUTINE convmas_loc
+
Index: LMDZ5/trunk/libf/dyn3dmem/enercin_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/enercin_loc.F	(revision 2335)
+++ 	(revision )
@@ -1,122 +1,0 @@
-      SUBROUTINE enercin_loc ( vcov, ucov, vcont, ucont, ecin )
-      USE parallel_lmdz
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteur: P. Le Van
-c   -------
-c
-c   Objet:
-c   ------
-c
-c *********************************************************************
-c .. calcul de l'energie cinetique aux niveaux s  ......
-c *********************************************************************
-c  vcov, vcont, ucov et ucont sont des arguments d'entree pour le s-pg .
-c  ecin         est  un  argument de sortie pour le s-pg
-c
-c=======================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-
-      REAL vcov( ijb_v:ije_v,llm ),vcont( ijb_v:ije_v,llm )
-      REAL ucov( ijb_u:ije_u,llm ),ucont( ijb_u:ije_u,llm )
-      REAL ecin( ijb_u:ije_u,llm )
-
-      REAL ecinni( iip1 ),ecinsi( iip1 )
-
-      REAL ecinpn, ecinps
-      INTEGER     l,ij,i,ijb,ije
-
-      EXTERNAL    SSUM
-      REAL        SSUM
-
-
-
-c                 . V
-c                i,j-1
-
-c      alpha4 .       . alpha1
-
-
-c        U .      . P     . U
-c       i-1,j    i,j      i,j
-
-c      alpha3 .       . alpha2
-
-
-c                 . V
-c                i,j
-
-c    
-c  L'energie cinetique au point scalaire P(i,j) ,autre que les poles, est :
-c       Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 )  +
-c              0.5 * U(i  ,j)**2 *( alpha1 + alpha2 )  +
-c              0.5 * V(i,j-1)**2 *( alpha1 + alpha4 )  +
-c              0.5 * V(i,  j)**2 *( alpha2 + alpha3 )
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO 5 l = 1,llm
-      
-      ijb=ij_begin
-      ije=ij_end+iip1
-      
-      IF (pole_nord) ijb=ij_begin+iip1
-      IF (pole_sud)  ije=ij_end-iip1
-      
-      DO 1  ij = ijb, ije -1
-      ecin( ij+1, l )  =    0.5  *
-     * (   ucov( ij   ,l ) * ucont( ij   ,l ) * alpha3p4( ij +1 )   +
-     *     ucov( ij+1 ,l ) * ucont( ij+1 ,l ) * alpha1p2( ij +1 )   +
-     *     vcov(ij-iim,l ) * vcont(ij-iim,l ) * alpha1p4( ij +1 )   +
-     *     vcov( ij+ 1,l ) * vcont( ij+ 1,l ) * alpha2p3( ij +1 )   )
-   1  CONTINUE
-
-c    ... correction pour  ecin(1,j,l)  ....
-c    ...   ecin(1,j,l)= ecin(iip1,j,l) ...
-
-CDIR$ IVDEP
-      DO 2 ij = ijb, ije, iip1
-      ecin( ij,l ) = ecin( ij + iim, l )
-   2  CONTINUE
-
-c     calcul aux poles  .......
-
-      IF (pole_nord) THEN
-    
-        DO  i = 1, iim
-         ecinni(i) = vcov(    i  ,  l) * 
-     *               vcont(    i    ,l) * aire(   i   )
-        ENDDO
-
-        ecinpn = 0.5 * SSUM( iim,ecinni,1 ) / apoln
-
-        DO ij = 1,iip1
-          ecin(   ij     , l ) = ecinpn
-        ENDDO
-   
-      ENDIF
-
-      IF (pole_sud) THEN
-    
-        DO  i = 1, iim
-         ecinsi(i) = vcov(i+ip1jmi1,l)* 
-     *               vcont(i+ip1jmi1,l) * aire(i+ip1jm)
-        ENDDO
-
-        ecinps = 0.5 * SSUM( iim,ecinsi,1 ) / apols
-
-        DO ij = 1,iip1
-          ecin( ij+ ip1jm, l ) = ecinps
-        ENDDO
-   
-      ENDIF
-
-      
-   5  CONTINUE
-c$OMP END DO NOWAIT
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/enercin_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/enercin_loc.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3dmem/enercin_loc.F90	(revision 2336)
@@ -0,0 +1,84 @@
+SUBROUTINE enercin_loc ( vcov, ucov, vcont, ucont, ecin )
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van.
+!-------------------------------------------------------------------------------
+! Purpose: Compute kinetic energy at sigma levels.
+  USE parallel_lmdz
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: vcov    (ijb_v:ije_v,llm)
+  REAL, INTENT(IN)  :: ucov    (ijb_u:ije_u,llm)
+  REAL, INTENT(IN)  :: vcont   (ijb_v:ije_v,llm)
+  REAL, INTENT(IN)  :: ucont   (ijb_u:ije_u,llm)
+  REAL, INTENT(OUT) :: ecin    (ijb_u:ije_u,llm)
+!===============================================================================
+! Notes:
+!                 . V
+!                i,j-1
+!
+!      alpha4 .       . alpha1
+!
+!
+!        U .      . P     . U
+!       i-1,j    i,j      i,j
+!
+!      alpha3 .       . alpha2
+!
+!
+!                 . V
+!                i,j
+!
+! Kinetic energy at scalar point P(i,j) (excluding poles) is:
+!       Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 )  +
+!              0.5 * U(i  ,j)**2 *( alpha1 + alpha2 )  +
+!              0.5 * V(i,j-1)**2 *( alpha1 + alpha4 )  +
+!              0.5 * V(i,  j)**2 *( alpha2 + alpha3 )
+!===============================================================================
+! Local variables:
+  INTEGER :: l, ij, i, ijb, ije
+  REAL    :: ecinni(iim), ecinsi(iim), ecinpn, ecinps
+!===============================================================================
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+  DO l=1,llm
+
+    ijb=ij_begin
+    ije=ij_end+iip1
+
+    IF(pole_nord) ijb=ij_begin+iip1
+    IF(pole_sud)  ije=ij_end-iip1
+
+    DO ij = ijb,ije-1
+      ecin(ij+1,l)=0.5*(ucov(ij    ,l)*ucont(ij    ,l)*alpha3p4(ij +1)          &
+                      + ucov(ij+1  ,l)*ucont(ij+1  ,l)*alpha1p2(ij +1)          &
+                      + vcov(ij-iim,l)*vcont(ij-iim,l)*alpha1p4(ij +1)          &
+                      + vcov(ij+1  ,l)*vcont(ij+1  ,l)*alpha2p3(ij +1) )
+    END DO
+
+    !--- Correction: ecin(1,j,l)= ecin(iip1,j,l)
+    DO ij=ijb,ije,iip1; ecin(ij,l) = ecin(ij+iim,l); END DO
+
+    !--- North pole
+    IF(pole_nord) THEN
+      ecinni(:) = vcov(1:iim,l)*vcont(1:iim,l)*aire(1:iim)
+      ecinpn = 0.5*SUM(ecinni)/apoln
+      ecin(1:iip1,l)=ecinpn
+    END IF
+
+    !--- South pole
+    IF(pole_sud) THEN
+      DO i=1,iim
+        ecinsi(i) = vcov(i+ip1jmi1,l)*vcont(i+ip1jmi1,l)*aire(i+ip1jm)
+      END DO
+      ecinps = 0.5*SUM(ecinsi)/apols
+      ecin(1+ip1jm:ip1jmp1,l)=ecinps
+    END IF
+  END DO
+!$OMP END DO NOWAIT
+
+END SUBROUTINE enercin_loc
+
Index: LMDZ5/trunk/libf/dyn3dmem/flumass_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/flumass_loc.F	(revision 2335)
+++ 	(revision )
@@ -1,152 +1,0 @@
-      SUBROUTINE flumass_loc(massebx,masseby,vcont,ucont,pbaru,pbarv)
-      USE parallel_lmdz
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteurs:  P. Le Van, F. Hourdin  .
-c   -------
-c
-c   Objet:
-c   ------
-c
-c *********************************************************************
-c     .... calcul du flux de masse  aux niveaux s ......
-c *********************************************************************
-c   massebx,masseby,vcont et ucont sont des argum. d'entree pour le s-pg .
-c       pbaru  et pbarv            sont des argum.de sortie pour le s-pg .
-c
-c=======================================================================
-
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-
-      REAL massebx( ijb_u:ije_u,llm ),masseby( ijb_v:ije_v,llm ) ,
-     * vcont( ijb_v:ije_v,llm ),ucont( ijb_u:ije_u,llm ),
-     * pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm )
-
-      REAL apbarun( iip1 ),apbarus( iip1 )
-
-      REAL sairen,saireun,saires,saireus,ctn,cts,ctn0,cts0
-      INTEGER  l,ij,i
-      INTEGER ijb,ije
-      
-      EXTERNAL   SSUM
-      REAL       SSUM
-      
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)    
-      DO  5 l = 1,llm
-
-        ijb=ij_begin
-        ije=ij_end+iip1
-      
-        if (pole_nord) ijb=ij_begin+iip1
-        if (pole_sud)  ije=ij_end-iip1
-        
-        DO  1 ij = ijb,ije
-          pbaru( ij,l ) = massebx( ij,l ) * ucont( ij,l )
-   1    CONTINUE
-
-        ijb=ij_begin-iip1
-        ije=ij_end+iip1
-      
-        if (pole_nord) ijb=ij_begin
-        if (pole_sud)  ije=ij_end-iip1
-        
-        DO 3 ij = ijb,ije
-          pbarv( ij,l ) = masseby( ij,l ) * vcont( ij,l )
-   3    CONTINUE
-
-   5  CONTINUE
-c$OMP END DO NOWAIT
-c    ................................................................
-c     calcul de la composante du flux de masse en x aux poles .......
-c    ................................................................
-c     par la resolution d'1 systeme de 2 equations .
-
-c     la premiere equat.decrivant le calcul de la divergence en 1 point i
-c     du pole,ce calcul etant itere de i=1 a i=im .
-c                 c.a.d   ,
-c     ( ( 0.5*pbaru(i)-0.5*pbaru(i-1) - pbarv(i))/aire(i)   =
-c                                           - somme de ( pbarv(n) )/aire pole
-
-c     l'autre equat.specifiant que la moyenne du flux de masse au pole est =0.
-c     c.a.d    somme de pbaru(n)*aire locale(n) = 0.
-
-c     on en revient ainsi a determiner la constante additive commune aux pbaru
-c     qui representait pbaru(0,j,l) dans l'equat.du calcul de la diverg.au pt
-c     i=1 .
-c     i variant de 1 a im
-c     n variant de 1 a im
-
-      IF (pole_nord) THEN
-     
-        sairen = SSUM( iim,  aire(   1     ), 1 )
-        saireun= SSUM( iim, aireu(   1     ), 1 )
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)    
-        DO l = 1,llm
- 
-          ctn =  SSUM( iim, pbarv(    1     ,l),  1 )/ sairen
-      
-          pbaru(1,l)=pbarv(1,l) - ctn * aire(1)
-        
-          DO i = 2,iim
-            pbaru(i,l) = pbaru(i- 1,l )    +
-     *                   pbarv(i,l) - ctn * aire(i )
-          ENDDO
-        
-          DO i = 1,iim
-            apbarun(i) = aireu(    i   ) * pbaru(   i    , l)
-          ENDDO
-      
-          ctn0 = -SSUM( iim,apbarun,1 )/saireun
-        
-          DO i = 1,iim
-            pbaru(   i    , l) = 2. * ( pbaru(   i    , l) + ctn0 )
-          ENDDO
-       
-          pbaru(   iip1 ,l ) = pbaru(    1    ,l )
-        
-        ENDDO
-c$OMP END DO NOWAIT              
-
-      ENDIF
-
-      
-      IF (pole_sud) THEN
-  
-        saires = SSUM( iim,  aire( ip1jm+1 ), 1 )
-        saireus= SSUM( iim, aireu( ip1jm+1 ), 1 )
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)    
-        DO  l = 1,llm
- 
-          cts =  SSUM( iim, pbarv(ip1jmi1+ 1,l),  1 )/ saires
-          pbaru(ip1jm+1,l)= - pbarv(ip1jmi1+1,l) + cts * aire(ip1jm+1)
-   
-          DO i = 2,iim
-            pbaru(i+ ip1jm,l) = pbaru(i+ip1jm-1,l)    -
-     *                          pbarv(i+ip1jmi1,l)+cts*aire(i+ip1jm)
-          ENDDO
-        
-          DO i = 1,iim
-            apbarus(i) = aireu(i +ip1jm) * pbaru(i +ip1jm, l)
-          ENDDO
-
-          cts0 = -SSUM( iim,apbarus,1 )/saireus
-
-          DO i = 1,iim
-            pbaru(i+ ip1jm, l) = 2. * ( pbaru(i +ip1jm, l) + cts0 )
-          ENDDO
-
-          pbaru( ip1jmp1,l ) = pbaru( ip1jm +1,l )
-       
-        ENDDO
-c$OMP END DO NOWAIT         
-      ENDIF
-      
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/flumass_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/flumass_loc.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3dmem/flumass_loc.F90	(revision 2336)
@@ -0,0 +1,91 @@
+SUBROUTINE flumass_loc(massebx,masseby, vcont, ucont, pbaru, pbarv )
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute mass flux at s levels.
+  USE parallel_lmdz
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: massebx(ijb_u:ije_u,llm)
+  REAL, INTENT(IN)  :: masseby(ijb_v:ije_v,llm)
+  REAL, INTENT(IN)  :: vcont  (ijb_v:ije_v,llm)
+  REAL, INTENT(IN)  :: ucont  (ijb_u:ije_u,llm)
+  REAL, INTENT(OUT) :: pbaru  (ijb_u:ije_u,llm)
+  REAL, INTENT(OUT) :: pbarv  (ijb_v:ije_v,llm)
+!===============================================================================
+! Method used:   A 2 equations system is solved.
+!   * 1st one describes divergence computation at pole point nr. i (i=1 to im):
+!     (0.5*(pbaru(i)-pbaru(i-1))-pbarv(i))/aire(i) = - SUM(pbarv(n))/aire pole
+!   * 2nd one specifies that mean mass flux at pole is equal to 0:
+!     SUM(pbaru(n)*local_area(n))=0
+! This way, we determine additive constant common to pbary elements representing
+!   pbaru(0,j,l) in divergence computation equation for point i=1. (i=1 to im)
+!===============================================================================
+! Local variables:
+  REAL    :: sairen, saireun, ctn, ctn0, apbarun(iim)
+  REAL    :: saires, saireus, cts, cts0, apbarus(iim)
+  INTEGER :: l, i, ij, ijb, ije
+!===============================================================================
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)    
+  DO l=1,llm
+
+    ijb=ij_begin
+    ije=ij_end+iip1
+    IF(pole_nord) ijb=ij_begin+iip1
+    IF(pole_sud)  ije=ij_end-iip1
+    pbaru(ijb:ije,l)=massebx(ijb:ije,l)*ucont(ijb:ije,l)
+
+    ijb=ij_begin-iip1
+    ije=ij_end+iip1
+    IF(pole_nord) ijb=ij_begin
+    IF(pole_sud)  ije=ij_end-iip1
+    pbarv(ijb:ije,l)=masseby(ijb:ije,l)*vcont(ijb:ije,l)
+
+  END DO
+!$OMP END DO NOWAIT
+
+  !--- North pole
+  IF(pole_nord) THEN
+    sairen =SUM(aire (1:iim))
+    saireun=SUM(aireu(1:iim))
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)    
+    DO l=1,llm
+      ctn=SUM(pbarv(1:iim,l))/sairen
+      pbaru(1,l)= pbarv(1,l)-ctn*aire(1)
+      DO i=2,iim
+        pbaru(i,l)=pbaru(i-1,l)+pbarv(i,l)-ctn*aire(i)
+      END DO
+      apbarun(:)=aireu(1:iim)*pbaru(1:iim,l)
+      ctn0 = -SUM(apbarun)/saireun
+      pbaru(1:iim,l)=2.*(pbaru(1:iim,l)+ctn0)
+      pbaru(iip1,l)=pbaru(1,l)
+    END DO
+!$OMP END DO NOWAIT              
+  END IF
+
+  !--- South pole
+  IF(pole_sud) THEN
+    saires =SUM(aire (ip1jm+1:ip1jmp1-1))
+    saireus=SUM(aireu(ip1jm+1:ip1jmp1-1))
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)    
+    DO l=1,llm
+      cts=SUM(pbarv(1+ip1jmi1:ip1jm-1,l))/saires
+      pbaru(1+ip1jm,l)=-pbarv(1+ip1jmi1,l)+cts*aire(1+ip1jm)
+      DO i=2,iim
+        pbaru(i+ip1jm,l)=pbaru(i-1+ip1jm,l)-pbarv(i+ip1jmi1,l)+cts*aire(i+ip1jm)
+      END DO
+      apbarus(:)=aireu(1+ip1jm:ip1jmp1-1)*pbaru(1+ip1jm:ip1jmp1-1,l)
+      cts0 = -SUM(apbarus)/saireus
+      pbaru(1+ip1jm:ip1jmp1-1,l)=2.*(pbaru(1+ip1jm:ip1jmp1-1,l)+cts0)
+      pbaru(ip1jmp1,l)=pbaru(1+ip1jm,l)
+    END DO
+!$OMP END DO NOWAIT         
+  END IF
+
+END SUBROUTINE flumass_loc
+
Index: LMDZ5/trunk/libf/dyn3dmem/massbar_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/massbar_loc.F	(revision 2335)
+++ 	(revision )
@@ -1,117 +1,0 @@
-      SUBROUTINE massbar_loc(  masse, massebx, masseby )
-     
-c
-c **********************************************************************
-c
-c  Calcule les moyennes en x et  y de la masse d'air dans chaque maille.
-c **********************************************************************
-c    Auteurs : P. Le Van , Fr. Hourdin  .
-c   ..........
-c
-c  ..  masse                 est  un argum. d'entree  pour le s-pg ...
-c  ..  massebx,masseby      sont des argum. de sortie pour le s-pg ...
-c     
-c
-      USE parallel_lmdz
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-c
-      REAL    masse( ijb_u:ije_u,llm ), massebx( ijb_u:ije_u,llm )  ,
-     *      masseby(   ijb_v:ije_v,llm )
-      INTEGER ij,l,ijb,ije
-c
-c
-c   Methode pour calculer massebx et masseby .
-c   ----------------------------------------
-c
-c    A chaque point scalaire P (i,j) est affecte 4 coefficients d'aires
-c       alpha1(i,j)  calcule  au point ( i+1/4,j-1/4 )
-c       alpha2(i,j)  calcule  au point ( i+1/4,j+1/4 )
-c       alpha3(i,j)  calcule  au point ( i-1/4,j+1/4 )
-c       alpha4(i,j)  calcule  au point ( i-1/4,j-1/4 )
-c
-c    Avec  alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j)        
-c
-c    N.B .  Pour plus de details, voir s-pg  ...  iniconst ...
-c
-c
-c
-c   alpha4 .         . alpha1    . alpha4
-c    (i,j)             (i,j)       (i+1,j)
-c
-c             P .        U .          . P
-c           (i,j)       (i,j)         (i+1,j)
-c
-c   alpha3 .         . alpha2    .alpha3 
-c    (i,j)              (i,j)     (i+1,j)
-c
-c             V .        Z .          . V
-c           (i,j)
-c
-c   alpha4 .         . alpha1    .alpha4
-c   (i,j+1)            (i,j+1)   (i+1,j+1) 
-c
-c             P .        U .          . P
-c          (i,j+1)                    (i+1,j+1)
-c
-c
-c
-c                       On  a :
-c
-c    massebx(i,j) = masse(i  ,j) * ( alpha1(i  ,j) + alpha2(i,j))   +
-c                   masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) )
-c     localise  au point  ... U (i,j) ...
-c
-c    masseby(i,j) = masse(i,j  ) * ( alpha2(i,j  ) + alpha3(i,j  )  +
-c                   masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1)  
-c     localise  au point  ... V (i,j) ...
-c
-c
-c=======================================================================
-      
-      
-      
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)  
-      DO   100    l = 1 , llm
-c
-        ijb=ij_begin
-        ije=ij_end+iip1
-        if (pole_sud) ije=ije-iip1
-        
-        DO  ij = ijb, ije - 1
-         massebx(ij,l) =  masse( ij, l) * alpha1p2( ij  )     + 
-     *                   masse(ij+1, l) * alpha3p4(ij+1 )
-        ENDDO
-
-c    .... correction pour massebx( iip1,j) .....
-c    ...    massebx(iip1,j)= massebx(1,j) ...
-c
-CDIR$ IVDEP
-
-        
-
-        DO  ij = ijb+iim, ije+iim, iip1
-         massebx( ij,l ) = massebx( ij - iim,l )
-        ENDDO
-
-
-      
-        ijb=ij_begin-iip1
-        ije=ij_end+iip1
-        if (pole_nord) ijb=ij_begin
-        if (pole_sud) ije=ij_end-iip1
-
-         DO  ij = ijb,ije
-         masseby( ij,l ) = masse(  ij   , l ) * alpha2p3(   ij    )  +
-     *                     masse(ij+iip1, l ) * alpha1p4( ij+iip1 )
-         ENDDO
-
-100   CONTINUE
-c$OMP END DO NOWAIT
-c
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/massbar_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/massbar_loc.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3dmem/massbar_loc.F90	(revision 2336)
@@ -0,0 +1,76 @@
+SUBROUTINE massbar_loc(masse,massebx,masseby)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute air mass mean along X and Y in each cell.
+! See iniconst for more details.
+  USE parallel_lmdz
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: masse  (ijb_u:ije_u,llm)
+  REAL, INTENT(OUT) :: massebx(ijb_u:ije_u,llm)
+  REAL, INTENT(OUT) :: masseby(ijb_v:ije_v,llm)
+!-------------------------------------------------------------------------------
+! Method used. Each scalar point is associated to 4 area coefficients:
+!    * alpha1(i,j) at point ( i+1/4,j-1/4 )
+!    * alpha2(i,j) at point ( i+1/4,j+1/4 )
+!    * alpha3(i,j) at point ( i-1/4,j+1/4 )
+!    * alpha4(i,j) at point ( i-1/4,j-1/4 )
+! where alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j)
+!
+!   alpha4 .         . alpha1    . alpha4
+!    (i,j)             (i,j)       (i+1,j)
+!
+!             P .        U .          . P
+!           (i,j)       (i,j)         (i+1,j)
+!
+!   alpha3 .         . alpha2    .alpha3 
+!    (i,j)              (i,j)     (i+1,j)
+!
+!             V .        Z .          . V
+!           (i,j)
+!
+!   alpha4 .         . alpha1    .alpha4
+!   (i,j+1)            (i,j+1)   (i+1,j+1) 
+!
+!             P .        U .          . P
+!          (i,j+1)                    (i+1,j+1)
+!
+!
+!    massebx(i,j) = masse(i  ,j) * ( alpha1(i  ,j) + alpha2(i,j))   +
+!                   masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) )
+!     localized at point  ... U (i,j) ...
+!
+!    masseby(i,j) = masse(i,j  ) * ( alpha2(i,j  ) + alpha3(i,j  )  +
+!                   masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1)  
+!     localized at point  ... V (i,j) ...
+!===============================================================================
+! Local variables:
+  INTEGER :: ij, l, ijb, ije
+!===============================================================================
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)  
+  DO l=1,llm
+    ijb=ij_begin
+    ije=ij_end+iip1
+    IF(pole_sud) ije=ije-iip1
+    DO ij=ijb,ije-1
+      massebx(ij,l)=masse(ij,l)*alpha1p2(ij)+masse(ij+1   ,l)*alpha3p4(ij+1)
+    END DO
+    DO ij=ijb+iim,ije+iim,iip1; massebx(ij,l)=massebx(ij-iim,l); END DO
+    ijb=ij_begin-iip1
+    ije=ij_end+iip1
+    IF(pole_nord) ijb=ij_begin
+    IF(pole_sud) ije=ij_end-iip1
+    DO ij=ijb,ije
+      masseby(ij,l)=masse(ij,l)*alpha2p3(ij)+masse(ij+iip1,l)*alpha1p4(ij+iip1)
+    END DO
+  END DO
+!$OMP END DO NOWAIT
+
+END SUBROUTINE massbar_loc
+
Index: LMDZ5/trunk/libf/dyn3dmem/massbarxy_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/massbarxy_loc.F	(revision 2335)
+++ 	(revision )
@@ -1,55 +1,0 @@
-      SUBROUTINE massbarxy_loc(  masse, massebxy )
-      USE parallel_lmdz
-      implicit none
-c **********************************************************************
-c
-c  Calcule les moyennes en x et  y de la masse d'air dans chaque maille.
-c **********************************************************************
-c    Auteurs : P. Le Van , Fr. Hourdin  .
-c   ..........
-c
-c  ..  masse          est  un  argum. d'entree  pour le s-pg ...
-c  ..  massebxy       est  un  argum. de sortie pour le s-pg ...
-c     
-c
-c     IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-c
-       REAL  masse( ijb_u:ije_u,llm ), massebxy( ijb_v:ije_v,llm )
-c
-      INTEGER ij,l,ijb,ije
-
-      
-      ijb=ij_begin-iip1
-      ije=ij_end
-      
-      if (pole_nord) ijb=ijb+iip1
-      if (pole_sud)  ije=ije-iip1
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)    
-      DO   100    l = 1 , llm
-c
-      DO 5 ij = ijb, ije - 1
-      massebxy( ij,l ) = masse(    ij  ,l ) * alpha2(   ij    )   +
-     +                   masse(   ij+1 ,l ) * alpha3(  ij+1   )   +
-     +                   masse( ij+iip1,l ) * alpha1( ij+iip1 )   +
-     +                   masse( ij+iip2,l ) * alpha4( ij+iip2 )
-   5  CONTINUE
-
-c    ....  correction pour     massebxy( iip1,j )  ........
-
-CDIR$ IVDEP
-
-      DO 7 ij = ijb+iip1-1, ije+iip1-1, iip1
-      massebxy( ij,l ) = massebxy( ij - iim,l )
-   7  CONTINUE
-
-100   CONTINUE
-c$OMP END DO NOWAIT
-c
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/massbarxy_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/massbarxy_loc.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3dmem/massbarxy_loc.F90	(revision 2336)
@@ -0,0 +1,39 @@
+SUBROUTINE massbarxy_loc(masse,massebxy)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute air mass mean along X and Y in each cell.
+! See iniconst for more details.
+  USE parallel_lmdz
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comconst.h"
+  include "comgeom.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: masse   (ijb_u:ije_u,llm)
+  REAL, INTENT(OUT) :: massebxy(ijb_v:ije_v,llm)
+!===============================================================================
+! Local variables:
+  INTEGER :: ij, l, ijb, ije
+!===============================================================================
+  ijb=ij_begin-iip1
+  ije=ij_end
+  IF(pole_nord) ijb=ijb+iip1
+  IF(pole_sud)  ije=ije-iip1
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)    
+  DO l=1,llm
+    DO ij=ijb,ije-1
+      massebxy(ij,l)=masse(ij     ,l)*alpha2(ij     ) + &
+     +               masse(ij+1   ,l)*alpha3(ij+1   ) + &
+     +               masse(ij+iip1,l)*alpha1(ij+iip1) + &
+     +               masse(ij+iip2,l)*alpha4(ij+iip2)
+    END DO
+    DO ij=ijb+iip1-1,ije+iip1-1,iip1; massebxy(ij,l)=massebxy(ij-iim,l); END DO
+  END DO
+!$OMP END DO NOWAIT
+
+END SUBROUTINE massbarxy_loc
+
Index: LMDZ5/trunk/libf/dyn3dmem/tourpot_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/tourpot_loc.F	(revision 2335)
+++ 	(revision )
@@ -1,95 +1,0 @@
-      SUBROUTINE tourpot_loc ( vcov, ucov, massebxy, vorpot )
-      USE parallel_lmdz
-      USE mod_filtreg_p
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteur:  P. Le Van
-c   -------
-c
-c   Objet:
-c   ------
-c
-c    *******************************************************************
-c    .........      calcul du tourbillon potentiel             .........
-c    *******************************************************************
-c
-c     vcov,ucov,fext et pbarxyfl sont des argum. d'entree pour le s-pg .
-c             vorpot            est  un argum.de sortie pour le s-pg .
-c
-c=======================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "logic.h"
-
-      REAL  rot( ijb_v:ije_v,llm )
-      REAL vcov( ijb_v:ije_v,llm ),ucov( ijb_u:ije_u,llm )
-      REAL massebxy( ijb_v:ije_v,llm ),vorpot( ijb_v:ije_v,llm )
-
-      INTEGER l, ij ,ije,ijb,jje,jjb
-
-
-      ijb=ij_begin-iip1
-      ije=ij_end
-      
-      if (pole_nord) ijb=ij_begin
-      
-      
-c  ... vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy ..
-
-
-
-c    ........  Calcul du rotationnel du vent V  puis filtrage  ........
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO 5 l = 1,llm
-
-      if (pole_sud)  ije=ij_end-iip1-1
-      DO 2 ij = ijb, ije 
-      rot( ij,l ) = vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)
-   2  CONTINUE
-
-c    ....  correction pour  rot( iip1,j,l )  .....
-c    ....     rot(iip1,j,l) = rot(1,j,l)    .....
-
-CDIR$ IVDEP
-
-      if (pole_sud)  ije=ij_end-iip1
-     
-      DO 3 ij = ijb+iip1-1, ije, iip1
-      rot( ij,l ) = rot( ij -iim, l )
-   3  CONTINUE
-
-   5  CONTINUE
-c$OMP END DO NOWAIT
-      jjb=jj_begin-1
-      jje=jj_end
-      
-      if (pole_nord) jjb=jjb+1
-      if (pole_sud)  jje=jje-1
-      CALL  filtreg_p( rot, jjb_v,jje_v,jjb,jje,jjm, llm, 
-     &                 2, 1, .FALSE., 1 )
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
-      DO 10 l = 1, llm
-      
-      if (pole_sud)  ije=ij_end-iip1-1  
-      
-      DO 6 ij = ijb, ije
-      vorpot( ij,l ) = ( rot(ij,l) + fext(ij) ) / massebxy(ij,l)
-   6  CONTINUE
-
-c    ..... correction pour  vorpot( iip1,j,l)  .....
-c    ....   vorpot(iip1,j,l)= vorpot(1,j,l) ....
-CDIR$ IVDEP
-      if (pole_sud)  ije=ij_end-iip1
-      DO 8 ij = ijb+iip1-1, ije, iip1
-      vorpot( ij,l ) = vorpot( ij -iim,l )
-   8  CONTINUE
-
-  10  CONTINUE
-c$OMP END DO NOWAIT
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/tourpot_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/tourpot_loc.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3dmem/tourpot_loc.F90	(revision 2336)
@@ -0,0 +1,65 @@
+SUBROUTINE tourpot_loc ( vcov, ucov, massebxy, vorpot )
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van.
+!-------------------------------------------------------------------------------
+! Purpose: Compute potential vorticity.
+  USE parallel_lmdz
+  USE mod_filtreg_p
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+  include "logic.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: vcov    (ijb_v:ije_v,llm)
+  REAL, INTENT(IN)  :: ucov    (ijb_u:ije_u,llm)
+  REAL, INTENT(IN)  :: massebxy(ijb_v:ije_v,llm)
+  REAL, INTENT(OUT) :: vorpot  (ijb_v:ije_v,llm)
+!===============================================================================
+! Method used:
+!   vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy
+!===============================================================================
+! Local variables:
+  INTEGER :: l, ij, ije, ijb, jje, jjb
+  REAL    :: rot(ijb_v:ije_v,llm)
+!===============================================================================
+
+  ijb=ij_begin-iip1
+  ije=ij_end
+  IF(pole_nord) ijb=ij_begin
+
+!--- Wind vorticity ; correction: rot(iip1,j,l) = rot(1,j,l)
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+  DO l=1,llm
+    IF(pole_sud) ije=ij_end-iip1-1
+    DO ij=ijb,ije
+      rot(ij,l)=vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)
+    END DO
+    IF(pole_sud) ije=ij_end-iip1
+    DO ij=ijb+iip1-1,ije,iip1; rot(ij,l)=rot(ij-iim,l); END DO
+  END DO
+!$OMP END DO NOWAIT
+
+!--- Filter
+  jjb=jj_begin-1
+  jje=jj_end
+  IF(pole_nord) jjb=jjb+1
+  IF(pole_sud)  jje=jje-1
+  CALL filtreg_p(rot,jjb_v,jje_v,jjb,jje,jjm,llm,2,1,.FALSE.,1)
+
+!--- Potential vorticity ; correction: rot(iip1,j,l) = rot(1,j,l)
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+  DO l=1,llm
+    IF(pole_sud) ije=ij_end-iip1-1
+    DO ij=ijb,ije
+      vorpot(ij,l)=(rot(ij,l)+fext(ij))/massebxy(ij,l)
+    END DO
+    IF(pole_sud) ije=ij_end-iip1
+    DO ij=ijb+iip1-1,ije,iip1; vorpot(ij,l)=vorpot(ij-iim,l); END DO
+  END DO
+!$OMP END DO NOWAIT
+
+END SUBROUTINE tourpot_loc
+
Index: LMDZ5/trunk/libf/dyn3dmem/vitvert_loc.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/vitvert_loc.F	(revision 2335)
+++ 	(revision )
@@ -1,56 +1,0 @@
-      SUBROUTINE vitvert_loc ( convm , w )
-c
-      USE parallel_lmdz
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Auteurs:  P. Le Van , F. Hourdin .
-c   -------
-c
-c   Objet:
-c   ------
-c
-c    *******************************************************************
-c  .... calcul de la vitesse verticale aux niveaux sigma  ....
-c    *******************************************************************
-c     convm   est un argument  d'entree pour le s-pg  ......
-c       w     est un argument de sortie pour le s-pg  ......
-c
-c    la vitesse verticale est orientee de  haut en bas .
-c    au sol, au niveau sigma(1),   w(i,j,1) = 0.
-c    au sommet, au niveau sigma(llm+1) , la vit.verticale est aussi
-c    egale a 0. et n'est pas stockee dans le tableau w  .
-c
-c
-c=======================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comvert.h"
-
-      REAL w(ijb_u:ije_u,llm),convm(ijb_u:ije_u,llm)
-      INTEGER   l, ij,ijb,ije
-
-
-      ijb=ij_begin
-      ije=ij_end+iip1
-      
-      if (pole_sud) ije=ij_end
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-      DO 2  l = 1,llmm1
-
-      DO 1 ij = ijb,ije
-      w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
-   1  CONTINUE
-
-   2  CONTINUE
-c$OMP END DO
-c$OMP MASTER
-      DO 5 ij  = ijb,ije
-      w(ij,1)  = 0.
-5     CONTINUE
-c$OMP END MASTER
-c$OMP BARRIER
-      RETURN
-      END
Index: LMDZ5/trunk/libf/dyn3dmem/vitvert_loc.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dmem/vitvert_loc.F90	(revision 2336)
+++ LMDZ5/trunk/libf/dyn3dmem/vitvert_loc.F90	(revision 2336)
@@ -0,0 +1,38 @@
+SUBROUTINE vitvert_loc(convm, w)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute vertical speed at sigma levels.
+  USE parallel_lmdz
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comvert.h"
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)  :: convm(ijb_u:ije_u,llm)
+  REAL, INTENT(OUT) :: w    (ijb_u:ije_u,llm)
+!===============================================================================
+! Notes: Vertical speed is oriented from bottom to top.
+!   * At ground - level sigma(1):     w(i,j,1) = 0.
+!   * At top    - level sigma(llm+1): w(i,j,l) = 0. (not stored in w)
+!===============================================================================
+! Local variables:
+  INTEGER :: l, ijb, ije
+!===============================================================================
+  ijb=ij_begin
+  ije=ij_end+iip1
+  IF(pole_sud) ije=ij_end
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+  DO l=1,llmm1
+    w(ijb:ije,l+1)=convm(ijb:ije,l+1)-bp(l+1)*convm(ijb:ije,1)
+  END DO
+!$OMP END DO
+!$OMP MASTER
+  w(ijb:ije,1)=0.
+!$OMP END MASTER
+!$OMP BARRIER
+
+END SUBROUTINE vitvert_loc
+
Index: LMDZ5/trunk/libf/dyn3dpar/ce0l.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/ce0l.F90	(revision 2335)
+++ 	(revision )
@@ -1,134 +1,0 @@
-!
-! $Id$
-!
-!-------------------------------------------------------------------------------
-!
-PROGRAM ce0l
-!-------------------------------------------------------------------------------
-! Purpose: Calls etat0, creates initial states and limit_netcdf
-!
-!     interbar=.T. for barycentric interpolation inter_barxy
-!     extrap  =.T. for data extrapolation, like for the SSTs when file does not
-!                  contain ocean points only.
-!     oldice  =.T. for old-style ice, obtained using grille_m (grid_atob).
-!     masque is created in etat0, passed to limit to ensure consistancy.
-!-------------------------------------------------------------------------------
-  USE control_mod
-#ifdef CPP_EARTH
-! This prog. is designed to work for Earth
-  USE dimphy
-  USE comgeomphy
-  USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
-  USE mod_const_mpi
-  USE infotrac
-  USE parallel_lmdz, ONLY: finalize_parallel
-  USE indice_sol_mod
-
-#ifdef CPP_IOIPSL
-  USE ioipsl, ONLY: ioconf_calendar
-#endif
-
-#endif
-  IMPLICIT NONE
-#ifndef CPP_EARTH
-#include "iniprint.h"
-  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
-#else
-!-------------------------------------------------------------------------------
-! Local variables:
-  LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.
-#include "dimensions.h"
-#include "paramet.h"
-!#include "indicesol.h"
-#include "iniprint.h"
-#include "temps.h"
-#include "logic.h"
-#ifdef CPP_MPI
-      include 'mpif.h'
-#endif
-
-  INTEGER                       :: ierr
-  REAL,    DIMENSION(iip1,jjp1) :: masque
-  CHARACTER(LEN=15)             :: calnd
-  REAL,    DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol
-!-------------------------------------------------------------------------------
-  CALL conf_gcm( 99, .TRUE. )
-
-#ifdef CPP_MPI
-  CALL init_mpi
-#endif
-
-  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
-  WRITE(lunout,*)'---> klon=',klon
-  IF (mpi_size>1 .OR. omp_size>1) THEN
-       CALL abort_gcm('ce0l','In parallel mode,                         &
- &                 ce0l must be called only                             &
- &                 for 1 process and 1 task',1)
-  ENDIF
-
-  CALL InitComgeomphy
-
-#ifdef CPP_IOIPSL
-  SELECT CASE(calend)
-    CASE('earth_360d');CALL ioconf_calendar('360d');      calnd='a 360 jours/an'
-    CASE('earth_365d');CALL ioconf_calendar('noleap');    calnd='a 365 jours/an'
-    CASE('earth_366d');CALL ioconf_calendar('366d');      calnd='bissextile'
-    CASE('gregorian'); CALL ioconf_calendar('gregorian'); calnd='gregorien'
-    CASE('standard');  CALL ioconf_calendar('gregorian'); calnd='gregorien'
-    CASE('julian');    CALL ioconf_calendar('julian');    calnd='julien'
-    CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
-  !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
-    CASE DEFAULT
-      CALL abort_gcm('ce0l','Mauvais choix de calendrier',1)
-  END SELECT
-  WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre '//TRIM(calnd)
-#endif
-
-  IF (type_trac == 'inca') THEN
-#ifdef INCA
-      CALL init_const_lmdz( &
-         nbtr,anneeref,dayref,&
-         iphysiq,day_step,nday,& 
-         nbsrf, is_oce,is_sic,&
-         is_ter,is_lic)
-      
-#endif
-  END IF
-
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,*) '  ***  etat0_netcdf ***  '
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) ' interbar = ',interbar
-  CALL etat0_netcdf(interbar,masque,phis,ok_etat0)
-
-  IF(ok_limit) THEN
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,*) '  ***  Limit_netcdf ***  '
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,'(//)')
-  CALL limit_netcdf(interbar,extrap,oldice,masque)
-  END IF
-
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  ***************************  '
-  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
-  WRITE(lunout,*) '  ***************************  '
-  WRITE(lunout,'(//)')
-  CALL grilles_gcm_netcdf_sub(masque,phis)
-  
-#ifdef CPP_MPI
-!$OMP MASTER
-  CALL MPI_FINALIZE(ierr)
-  IF (ierr /= 0) CALL abort_gcm('ce0l','Error in MPI_FINALIZE',1)
-!$OMP END MASTER
-#endif
-
-#endif
-! of #ifndef CPP_EARTH #else
-
-END PROGRAM ce0l
-!
-!-------------------------------------------------------------------------------
Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90	(revision 2335)
+++ LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90	(revision 2336)
@@ -1,28 +1,46 @@
 PROGRAM ce0l
-  !
-  ! Purpose: Calls etat0, creates initial states and limit_netcdf
-  !
-  ! interbar=.T. for barycentric interpolation inter_barxy
-  ! extrap  =.T. for data extrapolation, like for the SSTs when file does not
-  !                  contain ocean points only.
-  ! oldice  =.T. for old-style ice, obtained using grille_m (grid_atob).
-  ! masque is created in etat0, passed to limit to ensure consistancy.
-
-  USE control_mod, only: DAY_STEP, DAYREF, NSPLIT_PHYS
-  USE etat0dyn, only: etat0dyn_netcdf
-  USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
+!
+!-------------------------------------------------------------------------------
+! Purpose: Initial states and boundary conditions files creation:
+!     * start.nc    for dynamics    (using etat0dyn     routine)
+!     * startphy.nc for physics     (using etat0phys    routine)
+!     * limit.nc    for forced runs (using limit_netcdf routine)
+!-------------------------------------------------------------------------------
+! Notes:
+!     * extrap=.T. (default) for data extrapolation, like for the SSTs when file
+!                   does contain ocean points only.
+!     * "masque" can be:
+!       - read from file "o2a.nc"          (for coupled runs).
+!       - created in etat0phys or etat0dyn (for forced  runs).
+!     It is then passed to limit_netcdf to ensure consistancy.
+!-------------------------------------------------------------------------------
   USE ioipsl, ONLY: ioconf_calendar, getin, flininfo, flinopen, flinget, flinclo
-
-  USE etat0phys, only: etat0phys_netcdf
-  USE dimphy, only: KLON
-  USE infotrac, only: TYPE_TRAC, infotrac_init
+  USE control_mod,    ONLY: day_step, dayref, nsplit_phys
+  USE etat0dyn,       ONLY: etat0dyn_netcdf
+  USE etat0phys,      ONLY: etat0phys_netcdf
+  USE limit,          ONLY: limit_netcdf
+  USE netcdf,         ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
+  USE infotrac,       ONLY: type_trac, infotrac_init
+  USE dimphy,         ONLY: klon
   USE test_disvert_m, ONLY: test_disvert
+  USE filtreg_mod,    ONLY: inifilr
+#ifdef inca
+  USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic
+#endif
+#ifdef CPP_PARA
+  USE mod_const_mpi,  ONLY: init_const_mpi
+  USE parallel_lmdz,  ONLY: init_parallel, mpi_rank, omp_rank, mpi_size
+  USE bands,          ONLY: read_distrib, distrib_phys
+  USE mod_hallo,      ONLY: init_mod_hallo
+  USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys
+#endif
 
   IMPLICIT NONE
 
-  ! Local variables:
+!-------------------------------------------------------------------------------
+! Local variables:
   include "dimensions.h"
   include "paramet.h"
-  include "comgeom.h"
+  include "comgeom2.h"
   include "comconst.h"
   include "comvert.h"
@@ -30,21 +48,20 @@
   include "temps.h"
   include "logic.h"
-  REAL               :: masque(iip1, jjp1)             !--- CONTINENTAL MASK
-  REAL               :: phis  (iip1, jjp1)             !--- GROUND GEOPOTENTIAL
+  REAL               :: masque(iip1,jjp1)             !--- CONTINENTAL MASK
+  REAL               :: phis  (iip1,jjp1)             !--- GROUND GEOPOTENTIAL
   CHARACTER(LEN=256) :: modname, fmt, calnd           !--- CALENDAR TYPE
   LOGICAL            :: use_filtre_fft
-  LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.
-
-  !--- Local variables for ocean mask reading:
+  LOGICAL, PARAMETER :: extrap=.FALSE.
+
+!--- Local variables for ocean mask reading:
   INTEGER            :: nid_o2a, iml_omask, jml_omask, j
   INTEGER            :: fid, iret, llm_tmp, ttm_tmp, itaul(1)
-  REAL, ALLOCATABLE  :: lon_omask(:, :), dlon_omask(:), ocemask(:, :)
-  REAL, ALLOCATABLE  :: lat_omask(:, :), dlat_omask(:), ocetmp (:, :)
+  REAL, ALLOCATABLE  :: lon_omask(:,:), dlon_omask(:), ocemask(:,:)
+  REAL, ALLOCATABLE  :: lat_omask(:,:), dlat_omask(:), ocetmp (:,:)
   REAL               :: date, lev(1)
-
-  !----------------------------------------------------------------------
+!-------------------------------------------------------------------------------
   modname="ce0l"
 
-  !--- Constants
+!--- Constants
   pi     = 4. * ATAN(1.)
   rad    = 6371229.
@@ -59,145 +76,146 @@
 
   CALL conf_gcm( 99, .TRUE. )
-
   dtvr = daysec/REAL(day_step)
-  WRITE(lunout, *)'dtvr', dtvr
-
+  WRITE(lunout,*)'dtvr',dtvr
   CALL iniconst()
   CALL inigeom()
 
+!--- Calendar choice
 #ifdef CPP_IOIPSL
   calnd='gregorian'
   SELECT CASE(calend)
-  CASE('earth_360d')
-     CALL ioconf_calendar('360d')
-     calnd='with 360 days/year'
-  CASE('earth_365d')
-     CALL ioconf_calendar('noleap')
-     calnd='with no leap year'
-  CASE('earth_366d')
-     CALL ioconf_calendar('366d')
-     calnd='with leap years only'
-  CASE('gregorian')
-     CALL ioconf_calendar('gregorian')
-  CASE('standard')
-     CALL ioconf_calendar('gregorian')
-  CASE('julian')
-     CALL ioconf_calendar('julian')
-     calnd='julian'
-  CASE('proleptic_gregorian')
-     CALL ioconf_calendar('gregorian')
-     !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
-  CASE DEFAULT
-     CALL abort_gcm('ce0l', 'Bad choice for calendar', 1)
+    CASE('earth_360d');CALL ioconf_calendar('360d');   calnd='with 360 days/year'
+    CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd='with no leap year'
+    CASE('earth_366d');CALL ioconf_calendar('366d');   calnd='with leap years only'
+    CASE('gregorian'); CALL ioconf_calendar('gregorian')
+    CASE('standard');  CALL ioconf_calendar('gregorian')
+    CASE('julian');    CALL ioconf_calendar('julian'); calnd='julian'
+    CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
+  !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
+    CASE DEFAULT
+      CALL abort_gcm('ce0l','Bad choice for calendar',1)
   END SELECT
-  WRITE(lunout, *)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
-#endif
-
-  use_filtre_fft=.FALSE.
-  CALL getin('use_filtre_fft', use_filtre_fft)
-  IF(use_filtre_fft) THEN
-     WRITE(lunout, *)"FFT filter not available for sequential dynamics."
-     WRITE(lunout, *)"Your setting of variable use_filtre_fft is not used."
-  ENDIF
-
-  !--- LAND MASK. TWO CASES:
-  !   1) read from ocean model    file "o2a.nc"    (coupled runs)
-  !   2) computed from topography file "Relief.nc" (masque(:, :)=-99999.)
-  ! Coupled simulations (case 1) use the ocean model mask to compute the
-  ! weights to ensure ocean fractions are the same for atmosphere and ocean.
-
-  IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN
-     WRITE(lunout, *)'BEWARE !! No ocean mask "o2a.nc" file found'
-     WRITE(lunout, *)'Forced run.'
-     masque(:, :)=-99999.
-  ELSE
-     iret=NF90_CLOSE(nid_o2a)
-     WRITE(lunout, *)'BEWARE !! Ocean mask "o2a.nc" file found'
-     WRITE(lunout, *)'Coupled run.'
-     CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
-     IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
-        WRITE(lunout, *)'Mismatching dimensions for ocean mask'
-        WRITE(lunout, *)'iim  = ', iim , ' iml_omask = ', iml_omask
-        WRITE(lunout, *)'jjp1 = ', jjp1, ' jml_omask = ', jml_omask
-        CALL abort_gcm(modname, '', 1)
-     END IF
-     ALLOCATE(ocemask(iim, jjp1), lon_omask(iim, jjp1), dlon_omask(iim ))
-     ALLOCATE(ocetmp (iim, jjp1), lat_omask(iim, jjp1), dlat_omask(jjp1))
-     CALL flinopen("o2a.nc", .FALSE., iim, jjp1, llm_tmp, lon_omask, &
-          lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
-     CALL flinget(fid, "OceMask",    iim, jjp1, llm_tmp, ttm_tmp, 1, 1, ocetmp)
-     CALL flinclo(fid)
-     dlon_omask(1:iim ) = lon_omask(1:iim, 1)
-     dlat_omask(1:jjp1) = lat_omask(1, 1:jjp1)
-     ocemask = ocetmp
-     IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
-        DO j=1, jjp1
-           ocemask(:, j) = ocetmp(:, jjp1-j+1)
-        END DO
-     END IF
-     DEALLOCATE(ocetmp, lon_omask, lat_omask, dlon_omask, dlat_omask)
-     IF(prt_level>=1) THEN
-        WRITE(fmt, "(i4, 'i1)')")iim 
-        fmt='('//ADJUSTL(fmt)
-        WRITE(lunout, *)'OCEAN MASK :'
-        WRITE(lunout, fmt) NINT(ocemask)
-     END IF
-     masque(1:iim, :)=1.-ocemask(:, :)
-     masque(iip1 , :)=masque(1, :)
-     DEALLOCATE(ocemask)
-  END IF
-  phis(:, :)=-99999.
-
-  CALL Init_Phys_lmdz(iim, jjp1, llm, 1, (/(jjm-1)*iim+2/))
-  WRITE(lunout, *)'---> klon=', klon
-
-  call infotrac_init
-  CALL iniphysiq(iim, jjm, llm, daysec, dayref, dtphys / nsplit_phys, rlatu, &
-       rlonv, aire, cu, cv, rad, g, r, cpp, iflag_phys)
-
-  IF(pressure_exner) CALL test_disvert
-
+  WRITE(lunout,*)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
+#endif
+
+#ifdef CPP_PARA
+!--- Physical grid + parallel initializations
+  CALL init_const_mpi()
+  CALL init_parallel()
+  CALL read_distrib()
+  CALL init_mod_hallo()
+  CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
+  CALL init_interface_dyn_phys()
+#else
+  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
+#endif
+  WRITE(lunout,*)'---> klon=',klon
+
+!--- Tracers initializations
   IF (type_trac == 'inca') THEN
 #ifdef INCA
-     CALL init_const_lmdz(nbtr, anneeref, dayref, iphysiq, day_step, nday)
-     CALL init_inca_para(iim, jjm+1, klon, 1, klon_mpi_para_nb, 0)
-     WRITE(lunout, *)'nbtr =' , nbtr 
-#endif
-  END IF
+    CALL init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,&
+                         nbsrf,is_oce,is_sic,is_ter,is_lic,calend)
+    CALL init_inca_para(iim,jjp1,llm,klon_glo,mpi_size,distrib_phys,&
+                        COMM_LMDZ)
+    WRITE(lunout,*)'nbtr =' , nbtr 
+#endif
+  END IF
+  CALL infotrac_init()
+
+  CALL inifilr()
+  CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, &
+                 rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
+  IF(pressure_exner) CALL test_disvert
+
+#ifdef CPP_PARA
+  IF (mpi_rank==0.AND.omp_rank==0) THEN
+#endif
+  use_filtre_fft=.FALSE.
+  CALL getin('use_filtre_fft',use_filtre_fft)
+  IF(use_filtre_fft) THEN
+     WRITE(lunout,*)"FFT filter not available for sequential dynamics."
+     WRITE(lunout,*)"Your setting of variable use_filtre_fft is not used."
+  ENDIF
+
+!--- LAND MASK. TWO CASES:
+!   1) read from ocean model    file "o2a.nc"    (coupled runs)
+!   2) computed from topography file "Relief.nc" (masque(:,:)=-99999.)
+! Coupled simulations (case 1) use the ocean model mask to compute the
+! weights to ensure ocean fractions are the same for atmosphere and ocean.
+!*******************************************************************************
+  IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN
+    WRITE(lunout,*)'BEWARE !! No ocean mask "o2a.nc" file found'
+    WRITE(lunout,*)'Forced run.'
+    masque(:,:)=-99999.
+  ELSE
+    iret=NF90_CLOSE(nid_o2a)
+    WRITE(lunout,*)'BEWARE !! Ocean mask "o2a.nc" file found'
+    WRITE(lunout,*)'Coupled run.'
+    CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
+    IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
+      WRITE(lunout,*)'Mismatching dimensions for ocean mask'
+      WRITE(lunout,*)'iim  = ',iim ,' iml_omask = ',iml_omask
+      WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
+      CALL abort_gcm(modname,'',1)
+    END IF
+    ALLOCATE(ocemask(iim,jjp1),lon_omask(iim,jjp1),dlon_omask(iim ))
+    ALLOCATE(ocetmp (iim,jjp1),lat_omask(iim,jjp1),dlat_omask(jjp1))
+    CALL flinopen("o2a.nc", .FALSE.,iim,jjp1,llm_tmp,lon_omask,lat_omask,lev,  &
+     &                                       ttm_tmp,itaul,date,dt,fid)
+    CALL flinget(fid, "OceMask",    iim,jjp1,llm_tmp,ttm_tmp,1,1,ocetmp)
+    CALL flinclo(fid)
+    dlon_omask(1:iim ) = lon_omask(1:iim,1)
+    dlat_omask(1:jjp1) = lat_omask(1,1:jjp1)
+    ocemask = ocetmp
+    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
+      DO j=1,jjp1; ocemask(:,j) = ocetmp(:,jjp1-j+1); END DO
+    END IF
+    DEALLOCATE(ocetmp,lon_omask,lat_omask,dlon_omask,dlat_omask)
+    IF(prt_level>=1) THEN
+      WRITE(fmt,"(i4,'i1)')")iim ; fmt='('//ADJUSTL(fmt)
+      WRITE(lunout,*)'OCEAN MASK :'
+      WRITE(lunout,fmt) NINT(ocemask)
+    END IF
+    masque(1:iim,:)=1.-ocemask(:,:)
+    masque(iip1 ,:)=masque(1,:)
+    DEALLOCATE(ocemask)
+  END IF
+  phis(:,:)=-99999.
+
   IF(ok_etat0) THEN
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) '  ************************  '
-     WRITE(lunout, *) '  ***  etat0phy_netcdf ***  '
-     WRITE(lunout, *) '  ************************  '
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) ' interbar = ', interbar
-     CALL etat0phys_netcdf(interbar, masque, phis)
-  END IF
-
-  IF(ok_etat0) THEN
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) '  ************************  '
-     WRITE(lunout, *) '  ***  etat0dyn_netcdf ***  '
-     WRITE(lunout, *) '  ************************  '
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) ' interbar = ', interbar
-     CALL etat0dyn_netcdf(interbar, masque, phis)
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,*) '  ***  etat0phy_netcdf ***  '
+    WRITE(lunout,*) '  ************************  '
+    CALL etat0phys_netcdf(masque,phis)
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,*) '  ***  etat0dyn_netcdf ***  '
+    WRITE(lunout,*) '  ************************  '
+    CALL etat0dyn_netcdf(masque,phis)
   END IF
 
   IF(ok_limit) THEN
-     WRITE(lunout, '(//)')
-     WRITE(lunout, *) '  *********************  '
-     WRITE(lunout, *) '  ***  Limit_netcdf ***  '
-     WRITE(lunout, *) '  *********************  '
-     WRITE(lunout, '(//)')
-     CALL limit_netcdf(interbar, extrap, oldice, masque)
-  END IF
-
-  WRITE(lunout, '(//)')
-  WRITE(lunout, *) '  ***************************  '
-  WRITE(lunout, *) '  ***  grilles_gcm_netcdf ***  '
-  WRITE(lunout, *) '  ***************************  '
-  WRITE(lunout, '(//)')
-  CALL grilles_gcm_netcdf_sub(masque, phis)
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  *********************  '
+    WRITE(lunout,*) '  ***  Limit_netcdf ***  '
+    WRITE(lunout,*) '  *********************  '
+    WRITE(lunout,'(//)')
+    CALL limit_netcdf(masque,phis,extrap)
+  END IF
+
+  WRITE(lunout,'(//)')
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,'(//)')
+  CALL grilles_gcm_netcdf_sub(masque,phis)
+
+#ifdef CPP_PARA
+  END IF
+#endif
 
 END PROGRAM ce0l
+!
+!-------------------------------------------------------------------------------
Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90	(revision 2335)
+++ LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90	(revision 2336)
@@ -12,5 +12,5 @@
 !  routine (to be called after restget):
 !    CALL startget_dyn3d(varname, lon_in,  lat_in, pls, workvar,&
-!                          champ, lon_in2, lat_in2, ibar)
+!                          champ, lon_in2, lat_in2)
 !
 !    *  Variables should have the following names in the NetCDF files:
@@ -36,7 +36,5 @@
   USE ioipsl,         ONLY: flininfo, flinopen, flinget, flinclo, histclo
   USE assert_eq_m,    ONLY: assert_eq
-#ifdef CPP_PHYS
   USE indice_sol_mod, ONLY: epsfra
-#endif
   IMPLICIT NONE
 
@@ -54,11 +52,7 @@
   include "serre.h"
   REAL, SAVE :: deg2rad
-#ifndef CPP_PHYS
-  REAL, SAVE :: epsfra= 1.E-5
-#endif
   INTEGER,            SAVE      :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
   REAL, ALLOCATABLE,  SAVE      :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)
   CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc'
-  CHARACTER(LEN=120), PARAMETER :: orofname='Relief.nc'
 
 CONTAINS
@@ -66,5 +60,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE etat0dyn_netcdf(ib, masque, phis)
+SUBROUTINE etat0dyn_netcdf(masque, phis)
 !
 !-------------------------------------------------------------------------------
@@ -73,20 +67,17 @@
 ! Notes:  1) This routine is designed to work for Earth
 !         2) If masque(:,:)/=-99999., masque and phis are already known.
-!         Otherwise: read it from ocean model file (coupled run) or compute it.
+!         Otherwise: compute it.
 !-------------------------------------------------------------------------------
   USE control_mod
-#ifdef CPP_PHYS
   USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
   USE regr_pr_o3_m,   ONLY: regr_pr_o3
   USE press_coefoz_m, ONLY: press_coefoz
-#endif
   USE exner_hyb_m,    ONLY: exner_hyb
   USE exner_milieu_m, ONLY: exner_milieu
-  USE infotrac, only: NQTOT, TNAME
+  USE infotrac,       ONLY: nqtot, tname
   USE filtreg_mod
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
-  LOGICAL, INTENT(IN)    :: ib                  !--- Barycentric interpolation
   REAL,    INTENT(INOUT) :: masque(iip1,jjp1)   !--- Land-ocean mask
   REAL,    INTENT(INOUT) :: phis  (iip1,jjp1)   !--- Ground geopotential
@@ -111,12 +102,7 @@
   deg2rad = pi/180.0
 
-! Initializations for tracers and filter
-!*******************************************************************************
-  CALL inifilr()
-
 ! Compute ground geopotential and possibly the mask.
 !*******************************************************************************
   masque_tmp(:,:)=masque(:,:)
-  CALL start_init_orog0(rlonv, rlatu, phis, masque_tmp)
   WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
   IF(ALL(masque==-99999.)) THEN                         !--- KEEP NEW MASK
@@ -132,5 +118,5 @@
 ! Compute psol AND tsol, knowing phis.
 !*******************************************************************************
-  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, ib, phis, psol)
+  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol)
 
 ! Mid-levels pressure computation
@@ -147,10 +133,10 @@
 !*******************************************************************************
   uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
-  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv,ib)
+  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv)
   CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,      &
- &                           rlonu,rlatu(:jjm),ib)
-  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv,ib)
+ &                           rlonu,rlatu(:jjm))
+  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv)
   tpot(:,:,:)=t3d(:,:,:)
-  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv,ib)
+  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv)
 
   WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
@@ -165,9 +151,8 @@
 !  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
   qd (:,:,:) = 0.0
-  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv,ib)
+  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv)
   ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
   CALL flinclo(fid_dyn)
 
-#ifdef CPP_PHYS
 ! Parameterization of ozone chemistry:
 !*******************************************************************************
@@ -180,20 +165,5 @@
     q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29.                  !--- Mole->mass fraction         
   END IF
-
-#endif
   q3d(iip1,:,:,:)=q3d(1,:,:,:)
-
-! Intermediate computation
-!*******************************************************************************
-  CALL massdair(p3d,masse)
-  WRITE(lunout,*)' ALPHAX ',alphax
-  DO l=1,llm
-    xppn(:)=aire(1:iim,1   )*masse(1:iim,1   ,l)
-    xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
-    xpn=SUM(xppn)/apoln
-    xps=SUM(xpps)/apols
-    masse(:,1   ,l)=xpn
-    masse(:,jjp1,l)=xps
-  END DO
 
 ! Writing
@@ -215,10 +185,18 @@
   CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
   WRITE(lunout,*)'sortie geopot'
-  CALL caldyn0( itau, uvent, vvent, tpot, psol, pk, phis,               &
-                phi,  w, time+iday-dayref)
+  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
+                phi,  w, pbaru, pbarv, time+iday-dayref)
   WRITE(lunout,*)'sortie caldyn0'     
+#ifdef CPP_PARA
+  CALL dynredem0_loc( "start.nc", dayref, phis)
+#else
   CALL dynredem0( "start.nc", dayref, phis)
+#endif
   WRITE(lunout,*)'sortie dynredem0'
+#ifdef CPP_PARA
+  CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
+#else
   CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
+#endif
   WRITE(lunout,*)'sortie dynredem1' 
   CALL histclo()
@@ -232,5 +210,5 @@
 !
 SUBROUTINE startget_dyn3d(var, lon_in,  lat_in,  pls,  workvar,&
-                        champ, lon_in2, lat_in2, ibar)
+                        champ, lon_in2, lat_in2)
 !-------------------------------------------------------------------------------
   IMPLICIT NONE
@@ -252,5 +230,4 @@
   REAL,             INTENT(IN)    :: lon_in2(:)       ! dim (iml)
   REAL,             INTENT(IN)    :: lat_in2(:)       ! dim (jml2)
-  LOGICAL,          INTENT(IN)    :: ibar
 !-------------------------------------------------------------------------------
 ! Local variables:
@@ -287,5 +264,5 @@
 !--- INTERPOLATE 3D FIELD IF NEEDED
   IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,      &
-                                                  lat_in2,pls,champ,ibar)
+                                                  lat_in2,pls,champ)
 
 !--- COMPUTE THE REQUIRED FILED
@@ -317,212 +294,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
-!
-!-------------------------------------------------------------------------------
-  USE conf_dat_m, ONLY: conf_dat2d
-  IMPLICIT NONE
-!===============================================================================
-! Purpose:  Compute "phis" just like it would be in start_init_orog.
-!===============================================================================
-! Arguments:
-  REAL, INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
-  REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=256) :: modname="start_init_orog0"
-  CHARACTER(LEN=256) :: title="RELIEF"
-  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
-  REAL               :: lev(1), date, dt
-  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
-  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
-!-------------------------------------------------------------------------------
-  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
-  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
-  IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
-  IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
-  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
-  IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
-
-!--- HIGH RESOLUTION OROGRAPHY
-  CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
-
-  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
-  CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
-                lev, ttm_tmp, itau, date, dt, fid)
-  ALLOCATE(relief_hi(iml_rel,jml_rel))
-  CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
-  CALL flinclo(fid)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
-  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
-
-!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
-  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
-  DEALLOCATE(lon_ini,lat_ini)
-
-!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
-  WRITE(lunout,*)
-  WRITE(lunout,*)'*** Compute surface geopotential ***'
-
-!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
-  CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
-  phis = phis * 9.81
-  phis(iml,:) = phis(1,:)
-  DEALLOCATE(relief_hi,lon_rad,lat_rad)
-
-END SUBROUTINE start_init_orog0
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
-!
-!===============================================================================
-! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
-!          without any call to physics subroutines.
-!===============================================================================
-  IMPLICIT NONE 
-!-------------------------------------------------------------------------------
-! Arguments:
-  REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
-  REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
-  REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
-  REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
-  REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=256) :: modname="grid_noro0"
-  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
-  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
-  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
-  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
-  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
-  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
-  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
-  LOGICAL :: masque_lu
-  INTEGER :: i, ii, imdp, imar, iext
-  INTEGER :: j, jj, jmdp, jmar, nn
-  REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
-  REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
-!-------------------------------------------------------------------------------
-  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
-  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
-  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
-  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
-  IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
-  IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
-  iext=imdp/10
-  xpi = ACOS(-1.)
-  rad = 6371229.
-
-!--- ARE WE USING A READ MASK ?
-  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
-  WRITE(lunout,*)'Masque lu: ',masque_lu
-
-!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
-  ALLOCATE(xusn(imdp+2*iext))
-  xusn(1     +iext:imdp  +iext)=xd(:)
-  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
-  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
-
-  ALLOCATE(yusn(jmdp+2))
-  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
-  yusn(2:jmdp+1)=yd(:)
-  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
-
-  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
-  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
-  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
-  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
-  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
-  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
-  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
-  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
-
-!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
-  ALLOCATE(a(imar+1),b(imar+1))
-  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
-  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
-  a(1)=x(1)-(x(2)-x(1))/2.0
-  a(2:imar+1)= b(1:imar)
-
-  ALLOCATE(c(jmar),d(jmar))
-  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
-  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
-  c(1)=y(1)-(y(2)-y(1))/2.0
-  c(2:jmar)=d(1:jmar-1)
-
-!--- INITIALIZATIONS:
-  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
-  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
-
-!--- SUMMATION OVER GRIDPOINT AREA
-  zleny=xpi/REAL(jmdp)*rad
-  xincr=xpi/REAL(jmdp)/2.
-  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
-  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
-  DO ii = 1, imar+1
-    DO jj = 1, jmar
-      DO j = 2,jmdp+1 
-        zlenx  =zleny  *COS(yusn(j))
-        zbordnor=(xincr+c(jj)-yusn(j))*rad
-        zbordsud=(xincr-d(jj)+yusn(j))*rad
-        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
-        IF(weighy/=0) THEN
-          DO i = 2, imdp+2*iext-1
-            zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
-            zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
-            weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
-            IF(weighx/=0)THEN
-              num_tot(ii,jj)=num_tot(ii,jj)+1.0
-              IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
-              weight(ii,jj)=weight(ii,jj)+weighx*weighy
-              zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
-            END IF
-          END DO
-        END IF
-      END DO
-    END DO
-  END DO
-
-!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
-  IF(.NOT.masque_lu) THEN
-    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
-  END IF
-  nn=COUNT(weight(:,1:jmar-1)==0.0)
-  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
-  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
-
-!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
-  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
-  WHERE(mask>=0.1) mask_tmp = 1.
-  WHERE(weight(:,:)/=0.0)
-    zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
-    zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
-  END WHERE
-  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
-
-!--- Values at poles
-  zphi(imar+1,:)=zphi(1,:)
-
-  zweinor=SUM(weight(1:imar,   1),DIM=1)
-  zweisud=SUM(weight(1:imar,jmar),DIM=1)
-  zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
-  zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
-  zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
-
-END SUBROUTINE grid_noro0
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,ibar,zs,psol)
+SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol)
 !
 !-------------------------------------------------------------------------------
@@ -534,5 +304,4 @@
   REAL,    INTENT(IN)  :: lon_in (:),  lat_in (:)    ! dim (iml) (jml)
   REAL,    INTENT(IN)  :: lon_in2(:),  lat_in2(:)    ! dim (iml) (jml2)
-  LOGICAL, INTENT(IN)  :: ibar
   REAL,    INTENT(IN)  :: zs  (:,:)                  ! dim (iml,jml)
   REAL,    INTENT(OUT) :: psol(:,:)                  ! dim (iml,jml)
@@ -606,7 +375,7 @@
     ALLOCATE(field(iml,jml))
     CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
-    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)
-    CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,       &
-                              lon_in, lat_in, lon_in2, lat_in2, field)
+    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
+    CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana,              &
+                                        lon_in, lat_in, lon_in2, lat_in2, field)
   ELSE IF(SIZE(field)/=SIZE(z)) THEN
     msg='The '//TRIM(msg)//' field we have does not have the right size'
@@ -625,5 +394,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d,ibar)
+SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d)
 !
 !-------------------------------------------------------------------------------
@@ -639,5 +408,4 @@
   REAL,    INTENT(IN)  :: pls_in(:,:,:)           ! dim (iml,jml,lml)
   REAL,    INTENT(OUT) :: var3d (:,:,:)           ! dim (iml,jml,lml)
-  LOGICAL, INTENT(IN)  :: ibar
 !-------------------------------------------------------------------------------
 ! Local variables:
@@ -667,5 +435,5 @@
   ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn))
   CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini,                           &
-                       lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
+                       lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.)
   DEALLOCATE(lon_ini, lat_ini)
 
@@ -673,5 +441,5 @@
   ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
   DO il = 1,llm_dyn
-    CALL interp_startvar(var,ibar,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),     &
+    CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),          &
                           lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
   END DO
@@ -708,14 +476,13 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
+SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
 !
 !-------------------------------------------------------------------------------
   USE inter_barxy_m, ONLY: inter_barxy
-  USE grid_atob_m,   ONLY: grille_m
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
   CHARACTER(LEN=*), INTENT(IN)  :: nam
-  LOGICAL,          INTENT(IN)  :: ibar, ibeg
+  LOGICAL,          INTENT(IN)  :: ibeg
   REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
   REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
@@ -735,14 +502,10 @@
   j2=SIZE(lat2)
   ALLOCATE(vtmp(i1-1,j1))
-  IF(ibar) THEN
-    IF(ibeg.AND.prt_level>1) THEN
-      WRITE(lunout,*)"---------------------------------------------------------"
-      WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
-      WRITE(lunout,*)"---------------------------------------------------------"
-    END IF
-    CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
-  ELSE
-    CALL grille_m   (lon, lat,        vari, lon1,        lat1, vtmp)
-  END IF
+  IF(ibeg.AND.prt_level>1) THEN
+    WRITE(lunout,*)"---------------------------------------------------------"
+    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
+    WRITE(lunout,*)"---------------------------------------------------------"
+  END IF
+  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
   CALL gr_int_dyn(vtmp, varo, i1-1, j1)
 
Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90	(revision 2335)
+++ LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90	(revision 2336)
@@ -36,5 +36,6 @@
   USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
   USE assert_eq_m,        ONLY: assert_eq
-  USE dimphy
+  USE dimphy,             ONLY: klon
+  USE conf_dat_m,         ONLY: conf_dat2d
   USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
     rlon, solsw, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
@@ -51,18 +52,15 @@
   include "paramet.h"
   include "comgeom2.h"
-  include "comvert.h"
   include "comconst.h"
   include "dimsoil.h"
   include "temps.h"
-  include "comdissnew.h"
-  include "serre.h"
   include "clesphys.h"
   REAL, SAVE :: deg2rad
   REAL, SAVE, ALLOCATABLE :: tsol(:)
-!  REAL, SAVE, ALLOCATABLE :: rugo(:,:)  ! ??? COMPUTED BUT NOT USED ???
   INTEGER,            SAVE      :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
   REAL, ALLOCATABLE,  SAVE      :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:)
-  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc"
-  CHARACTER(LEN=256), PARAMETER :: title="RELIEF"
+  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF"
+  CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc",  psrfvar="SP"
+  CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW",       tsrfvar="ST"
 
 
@@ -72,10 +70,12 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE etat0phys_netcdf(ib, masque, phis)
+SUBROUTINE etat0phys_netcdf(masque, phis)
 !
 !-------------------------------------------------------------------------------
 ! Purpose: Creates initial states
 !-------------------------------------------------------------------------------
-! Note: This routine is designed to work for Earth
+! Notes:  1) This routine is designed to work for Earth
+!         2) If masque(:,:)/=-99999., masque and phis are already known.
+!         Otherwise: compute it.
 !-------------------------------------------------------------------------------
   USE control_mod
@@ -84,13 +84,9 @@
   USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
   USE indice_sol_mod
-  USE conf_phys_m,    ONLY: conf_phys
-  USE exner_hyb_m,    ONLY: exner_hyb
-  USE exner_milieu_m, ONLY: exner_milieu
-  USE test_disvert_m, ONLY: test_disvert
-  USE grid_atob_m,    ONLY: grille_m
+  USE conf_phys_m, ONLY: conf_phys
+  USE init_ssrf_m, ONLY: start_init_subsurf
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
-  LOGICAL, INTENT(IN)    :: ib          !--- Barycentric interpolation
   REAL,    INTENT(INOUT) :: masque(:,:) !--- Land mask           dim(iip1,jjp1)
   REAL,    INTENT(INOUT) :: phis  (:,:) !--- Ground geopotential dim(iip1,jjp1)
@@ -99,18 +95,10 @@
   CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt
   INTEGER            :: i, j, l, ji, iml, jml
-  REAL               :: phystep
+  LOGICAL            :: read_mask
+  REAL               :: phystep, dummy
   REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp
   REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
   REAL, DIMENSION(klon,nbsrf)         :: qsolsrf, snsrf
   REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
-
-!--- Local variables for sea-ice reading:
-  LOGICAL           :: read_mask
-  INTEGER           :: iml_lic, jml_lic, isst(klon-2)
-  INTEGER           :: fid, llm_tmp, ttm_tmp, itaul(1)
-  REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic(:,:)
-  REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:)
-  REAL              :: date, lev(1), dummy
-  REAL              :: flic_tmp(SIZE(masque,1),SIZE(masque,2))
 
 !--- Arguments for conf_phys
@@ -132,5 +120,5 @@
   jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml")
 
-! Grid construction and miscellanous initializations.
+! Physics configuration
 !*******************************************************************************
   CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
@@ -145,18 +133,8 @@
                    read_climoz,                                         &
                    alp_offset)
-
   CALL phys_state_var_init(read_climoz)
 
-  co2_ppm0 = co2_ppm  !--- Initial atmospheric CO2 conc. from .def file
-
-  rlat(1) = pi/2.
-  DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);     END DO
-  rlat(klon) = - pi/2.
-  rlat(:)=rlat(:)*(180.0/pi)
-
-  rlon(1) = 0.0
-  DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO
-  rlon(klon) = 0.0
-  rlon(:)=rlon(:)*(180.0/pi)
+!--- Initial atmospheric CO2 conc. from .def file
+  co2_ppm0 = co2_ppm
 
 ! Compute ground geopotential, sub-cells quantities and possibly the mask.
@@ -174,9 +152,9 @@
     WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
   END IF
- CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
+  CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
 
 ! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
 !*******************************************************************************
-  CALL start_init_phys(rlonv, rlatu, rlonu, rlatv, ib, phis)
+  CALL start_init_phys(rlonu, rlatv, phis)
 
 ! Some initializations.
@@ -188,56 +166,7 @@
     CALL regr_lat_time_climoz(read_climoz)
 
-! Sub-surfaces initialization
-!*******************************************************************************
-!--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
-  CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
-  ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
-  ALLOCATE(dlat_lic(jml_lic),       dlon_lic(iml_lic))
-  ALLOCATE( fraclic(iml_lic,jml_lic))
-  CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp,  &
- &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
-  CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
-  CALL flinclo(fid)
-  WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic
-  IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. !--- Conversion to degrees
-  IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
-  dlon_lic(:)=lon_lic(:,1)
-  dlat_lic(:)=lat_lic(1,:)
-  CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim,:) )
-  flic_tmp(iml,:)=flic_tmp(1,:)
-
-!--- To the physical grid
-  pctsrf(:,:) = 0.
-  CALL gr_dyn_fi(1, iml, jml, klon, flic_tmp, pctsrf(:,is_lic))
-
-!--- Adequation with soil/sea mask
-  WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 
-  WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
-  pctsrf(:,is_ter)=zmasq(:)
-  DO ji=1,klon
-    IF(zmasq(ji)>EPSFRA) THEN 
-      IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
-        pctsrf(ji,is_lic)=zmasq(ji)
-        pctsrf(ji,is_ter)=0.
-      ELSE
-        pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
-        IF(pctsrf(ji,is_ter)<EPSFRA) THEN
-          pctsrf(ji,is_ter)=0.
-          pctsrf(ji,is_lic)=zmasq(ji)
-        END IF 
-      END IF 
-    END IF 
-  END DO 
-
-!--- Sub-surface ocean and sea ice (sea ice set to zero for start).
-  pctsrf(:,is_oce)=(1.-zmasq(:))
-  WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
-  IF(read_mask) pctsrf(:,is_oce)=1-zmasq(:)
-  isst=0
-  WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
-
-!--- It is checked that the sub-surfaces sum is equal to 1.
-  ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
-  IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points'
+! Sub-surfaces initialization.
+!*******************************************************************************
+  CALL start_init_subsurf(read_mask)
 
 ! Write physical initial state
@@ -321,6 +250,4 @@
 !   described in LOTT & MILLER (1997) and LOTT(1999).
 !===============================================================================
-  USE conf_dat_m,  ONLY: conf_dat2d
-!  USE grid_atob_m, ONLY: rugsoro
   USE grid_noro_m, ONLY: grid_noro
   IMPLICIT NONE
@@ -331,6 +258,5 @@
 !-------------------------------------------------------------------------------
 ! Local variables:
-  CHARACTER(LEN=256) :: modname="start_init_orog"
-  CHARACTER(LEN=256) :: title="RELIEF"
+  CHARACTER(LEN=256) :: modname
   INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
   REAL               :: lev(1), date, dt
@@ -340,4 +266,5 @@
   REAL, ALLOCATABLE  :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:)
 !-------------------------------------------------------------------------------
+  modname="start_init_orog"
   iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
   jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
@@ -350,5 +277,5 @@
                 lev, ttm_tmp, itau, date, dt, fid)
   ALLOCATE(relief_hi(iml_rel,jml_rel))
-  CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
+  CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi)
   CALL flinclo(fid)
 
@@ -360,5 +287,5 @@
 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
   ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
-  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
+  CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.)
   DEALLOCATE(lon_ini,lat_ini)
 
@@ -378,12 +305,4 @@
   phis = phis * 9.81
   phis(iml,:) = phis(1,:)
-
-!--- COMPUTE SURFACE ROUGHNESS
-!  WRITE(lunout,*)
-!  WRITE(lunout,*)'*** Compute surface roughness induced by the orography ***'
-!  ALLOCATE(tmp_var(iml-1,jml))
-!  CALL rugsoro(lon_rad, lat_rad, relief_hi, lon_in(1:iml-1), lat_in, tmp_var)
-!  ALLOCATE(rugo(iml,jml)); rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
-!  DEALLOCATE(tmp_var)
   DEALLOCATE(relief_hi,lon_rad,lat_rad)
 
@@ -405,5 +324,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_init_phys(lon_in,lat_in,lon_in2,lat_in2,ibar,phis)
+SUBROUTINE start_init_phys(lon_in,lat_in,phis)
 !
 !===============================================================================
@@ -413,11 +332,9 @@
 !-------------------------------------------------------------------------------
 ! Arguments:
-  REAL,    INTENT(IN) :: lon_in (:),  lat_in (:)     ! dim (iml) (jml)
-  REAL,    INTENT(IN) :: lon_in2(:),  lat_in2(:)     ! dim (iml) (jml2)
-  LOGICAL, INTENT(IN) :: ibar
+  REAL,    INTENT(IN) :: lon_in(:),  lat_in(:)       ! dim (iml) (jml2)
   REAL,    INTENT(IN) :: phis(:,:)                   ! dim (iml,jml)
 !-------------------------------------------------------------------------------
 ! Local variables:
-  CHARACTER(LEN=256) :: modname="start_init_phys", physfname="ECPHY.nc"
+  CHARACTER(LEN=256) :: modname
   REAL               :: date, dt
   INTEGER            :: iml, jml, jml2, itau(1)
@@ -426,15 +343,15 @@
   REAL, ALLOCATABLE  :: ts(:,:), qs(:,:)
 !-------------------------------------------------------------------------------
-  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(lon_in2),TRIM(modname)//" iml")
-  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),              TRIM(modname)//" jml")
-  jml2=SIZE(lat_in2)
+  modname="start_init_phys"
+  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml")
+  jml=SIZE(phis,2); jml2=SIZE(lat_in)
 
   WRITE(lunout,*)'Opening the surface analysis'
-  CALL flininfo(physfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
-  WRITE(lunout,*) 'Values read: ',   iml_phys, jml_phys, llm_phys, ttm_phys
+  CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
+  WRITE(lunout,*) 'Values read: ',  iml_phys, jml_phys, llm_phys, ttm_phys
 
   ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys))
   ALLOCATE(levphys_ini(llm_phys))
-  CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
+  CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
                 lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys)
 
@@ -445,6 +362,6 @@
 
   ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
-  CALL get_var_phys('ST'  ,ts)                   !--- SURFACE TEMPERATURE
-  CALL get_var_phys('CDSW',qs)                   !--- SOIL MOISTURE
+  CALL get_var_phys(tsrfvar,ts)                   !--- SURFACE TEMPERATURE
+  CALL get_var_phys(qsolvar,qs)                   !--- SOIL MOISTURE
   CALL flinclo(fid_phys)
   DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
@@ -463,5 +380,4 @@
 !
 !-------------------------------------------------------------------------------
-  USE conf_dat_m, ONLY: conf_dat2d
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
@@ -474,13 +390,13 @@
 !-------------------------------------------------------------------------------
   SELECT CASE(title)
-    CASE('SP');        tllm=0
-    CASE('ST','CDSW'); tllm=llm_phys
+    CASE(psrfvar);         tllm=0
+    CASE(tsrfvar,qsolvar); tllm=llm_phys
   END SELECT
   IF(ALLOCATED(field)) RETURN
   ALLOCATE(field(iml,jml)); field(:,:)=0.
   CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana)
-  CALL conf_dat2d(title, lon_ini, lat_ini,  lon_rad, lat_rad, var_ana, ibar)
-  CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,         &
-                            lon_in, lat_in, lon_in2, lat_in2, field)
+  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
+  CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana,               &
+                                      lon_in,  lat_in,  field)
 
 END SUBROUTINE get_var_phys
@@ -495,41 +411,35 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
+SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo)
 !
 !-------------------------------------------------------------------------------
   USE inter_barxy_m, ONLY: inter_barxy
-  USE grid_atob_m,   ONLY: grille_m
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
   CHARACTER(LEN=*), INTENT(IN)  :: nam
-  LOGICAL,          INTENT(IN)  :: ibar, ibeg
+  LOGICAL,          INTENT(IN)  :: ibeg
   REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
   REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
-  REAL,             INTENT(IN)  :: lon1(:), lat1(:) ! dim (i1) (j1)
   REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
   REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
 !-------------------------------------------------------------------------------
 ! Local variables:
-  CHARACTER(LEN=256) :: modname="interp_startvar"
+  CHARACTER(LEN=256) :: modname
   INTEGER            :: ii, jj, i1, j1, j2
   REAL, ALLOCATABLE  :: vtmp(:,:)
 !-------------------------------------------------------------------------------
-  ii=assert_eq(SIZE(lon),            SIZE(vari,1),TRIM(modname)//" ii")
-  jj=assert_eq(SIZE(lat),            SIZE(vari,2),TRIM(modname)//" jj")
-  i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
-  j1=assert_eq(SIZE(lat1),           SIZE(varo,2),TRIM(modname)//" j1")
-  j2=SIZE(lat2)
+  modname="interp_startvar"
+  ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")
+  jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")
+  i1=assert_eq(SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
+  j1=SIZE(varo,2); j2=SIZE(lat2)
   ALLOCATE(vtmp(i1-1,j1))
-  IF(ibar) THEN
-    IF(ibeg.AND.prt_level>1) THEN
-      WRITE(lunout,*)"--------------------------------------------------------"
-      WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
-      WRITE(lunout,*)"--------------------------------------------------------"
-    END IF
-    CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
-  ELSE
-    CALL grille_m   (lon, lat,        vari, lon1,        lat1, vtmp)
+  IF(ibeg.AND.prt_level>1) THEN
+    WRITE(lunout,*)"--------------------------------------------------------"
+    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
+    WRITE(lunout,*)"--------------------------------------------------------"
   END IF
+  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
   CALL gr_int_dyn(vtmp, varo, i1-1, j1)
 
Index: LMDZ5/trunk/libf/phylmd/limit_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/limit_netcdf.F90	(revision 2335)
+++ LMDZ5/trunk/libf/phylmd/limit_netcdf.F90	(revision 2336)
@@ -1,5 +1,5 @@
-SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
-!
-!-------------------------------------------------------------------------------
+MODULE limit
+!
+!*******************************************************************************
 ! Author : L. Fairhead, 27/01/94
 !-------------------------------------------------------------------------------
@@ -16,35 +16,14 @@
 !  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
 !-------------------------------------------------------------------------------
-#ifndef CPP_1D
+
+  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo,          &
+   ioconf_calendar, ioget_calendar, lock_calendar, ioget_mon_len, ioget_year_len
+  USE assert_eq_m,        ONLY: assert_eq
+  USE conf_dat_m,         ONLY: conf_dat2d, conf_dat3d
+  USE dimphy,             ONLY: klon, zmasq
+  USE phys_state_var_mod, ONLY: pctsrf, rlon, rlat
   USE control_mod
-  USE indice_sol_mod
-  USE dimphy
-  USE ioipsl,             ONLY : ioget_year_len
-  USE phys_state_var_mod, ONLY : pctsrf, rlon, rlat
-  USE netcdf,             ONLY : NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,       &
-                   NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
-                   NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
-		   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
-  USE inter_barxy_m, only: inter_barxy
-  USE netcdf95,    ONLY: nf95_def_var, nf95_put_att, nf95_put_var
-  USE grid_atob_m, ONLY: grille_m, rugosite, sea_ice
-  USE print_control_mod, ONLY: prt_level,lunout
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  include "dimensions.h"
-  include "paramet.h"
-  LOGICAL,                    INTENT(IN) :: interbar ! barycentric interpolation
-  LOGICAL,                    INTENT(IN) :: extrap   ! SST extrapolation flag
-  LOGICAL,                    INTENT(IN) :: oldice   ! old way ice computation
-  REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
-!-------------------------------------------------------------------------------
-! Local variables:
-  include "logic.h"
-  include "comgeom2.h"
-  include "comconst.h"
-
-!--- INPUT NETCDF FILES NAMES --------------------------------------------------
-  CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam
+  USE init_ssrf_m,        ONLY: start_init_subsurf
+
   CHARACTER(LEN=20), PARAMETER :: &
   fsst(4)=['amipbc_sst_1x1.nc   ','cpl_atm_sst.nc      ','histmth_sst.nc      '&
@@ -56,8 +35,56 @@
   vsst(4)=['tosbcs    ','SISUTESW  ','tsol_oce  ','sstk      '], &
   vsic(4)=['sicbcs    ','SIICECOV  ','pourc_sic ','ci        ']
-  CHARACTER(LEN=20), PARAMETER :: frugo='Rugos.nc            ',                &
-                                  falbe='Albedo.nc           '
-  CHARACTER(LEN=10), PARAMETER :: vrug='RUGOS     ', valb='ALBEDO    '
+  CHARACTER(LEN=10), PARAMETER :: &
+  frugo='Rugos.nc  ', falbe='Albedo.nc ', frelf='Relief.nc ',    &
+   vrug='RUGOS     ',  valb='ALBEDO    ',  vrel='RELIEF    '
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE limit_netcdf(masque, phis, extrap)
+!
+!-------------------------------------------------------------------------------
+! Author : L. Fairhead, 27/01/94
+!-------------------------------------------------------------------------------
+! Purpose: Boundary conditions files building for new model using climatologies.
+!          Both grids have to be regular.
+!-------------------------------------------------------------------------------
+! Note: This routine is designed to work for Earth
+!-------------------------------------------------------------------------------
+! Modification history:
+!  * 23/03/1994: Z. X. Li
+!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
+!  *    07/2001: P. Le Van
+!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
+!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
+!-------------------------------------------------------------------------------
+#ifndef CPP_1D
+  USE indice_sol_mod
+  USE netcdf,             ONLY: NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,        &
+                  NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,      &
+                  NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,       &
+		  NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
+  USE inter_barxy_m,      ONLY: inter_barxy
+  USE netcdf95,           ONLY: nf95_def_var, nf95_put_att, nf95_put_var
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
+  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: phis   ! ground geopotential
+  LOGICAL,                    INTENT(IN)    :: extrap ! SST extrapolation flag
+!-------------------------------------------------------------------------------
+! Local variables:
+  include "logic.h"
+  include "comgeom2.h"
+  include "comconst.h"
+
+!--- INPUT NETCDF FILES NAMES --------------------------------------------------
+  CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam
   CHARACTER(LEN=10) :: varname
+
 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
   REAL               :: fi_ice(klon), verif(klon)
@@ -82,4 +109,12 @@
   CALL inigeom
 
+!--- MASK, GROUND GEOPOT. & SUBSURFACES COMPUTATION (IN CASE ok_etat0==.FALSE.)
+   IF(ALL(masque==-99999.)) THEN
+    CALL start_init_orog0(rlonv,rlatu,phis,masque)
+    CALL gr_dyn_fi(1,iip1,jjp1,klon,masque,zmasq)          !--- To physical grid
+    ALLOCATE(rlon(klon),rlat(klon),pctsrf(klon,nbsrf))
+    CALL start_init_subsurf(.FALSE.)
+  END IF
+
 !--- Beware: anneeref (from gcm.def) is used to determine output time sampling
   ndays=ioget_year_len(anneeref)
@@ -87,5 +122,5 @@
 !--- RUGOSITY TREATMENT --------------------------------------------------------
   CALL msg(1,'Traitement de la rugosite')
-  CALL get_2Dfield(frugo,vrug,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
+  CALL get_2Dfield(frugo,vrug,'RUG',ndays,phy_rug,mask=masque(1:iim,:))
 
 !--- OCEAN TREATMENT -----------------------------------------------------------
@@ -108,5 +143,5 @@
   CALL msg(-1,'Fichier choisi pour la glace de mer:'//TRIM(icefile))
 
-  CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
+  CALL get_2Dfield(icefile,varname, 'SIC',ndays,phy_ice)
 
   ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
@@ -167,9 +202,9 @@
   CALL msg(-1,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile))
 
-  CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
+  CALL get_2Dfield(sstfile,varname,'SST',ndays,phy_sst,flag=extrap)
 
 !--- ALBEDO TREATMENT ----------------------------------------------------------
   CALL msg(1,'Traitement de l albedo')
-  CALL get_2Dfield(falbe,valb,'ALB',interbar,ndays,phy_alb)
+  CALL get_2Dfield(falbe,valb,'ALB',ndays,phy_alb)
 
 !--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
@@ -251,5 +286,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask)
+SUBROUTINE get_2Dfield(fnam, varname, mode, ndays, champo, flag, mask)
 !
 !-----------------------------------------------------------------------------
@@ -263,8 +298,4 @@
        NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
        NF90_GET_ATT
-  USE dimphy, ONLY : klon
-  USE phys_state_var_mod, ONLY : pctsrf
-  USE conf_dat_m, ONLY: conf_dat2d
-  USE control_mod
   USE pchsp_95_m, only: pchsp_95
   USE pchfe_95_m, only: pchfe_95
@@ -281,5 +312,4 @@
   CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
   CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
-  LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
   INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
   REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
@@ -312,5 +342,4 @@
   CHARACTER(LEN=25) :: title              ! for messages
   LOGICAL           :: extrp              ! flag for extrapolation
-  LOGICAL           :: oldice             ! flag for old way ice computation 
   REAL              :: chmin, chmax
   INTEGER ierr
@@ -328,12 +357,5 @@
   CASE('ALB'); title='Albedo'
   END SELECT
-  
-
-  extrp=.FALSE. 
-  oldice=.FALSE.
-  IF ( PRESENT(flag) ) THEN 
-    IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 
-    IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. 
-  END IF
+  extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
 
 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
@@ -400,25 +422,16 @@
   DO l=1, lmdep
     CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam)
-    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, ibar)
+    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
     IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
-
-    IF(ibar .AND. .NOT.oldice) THEN
-      IF(l==1) THEN
+    IF(l==1) THEN
       CALL msg(5,"----------------------------------------------------------")
       CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$")
       CALL msg(5,"----------------------------------------------------------")
-      END IF
-      IF(mode=='RUG') champ=LOG(champ)
-      CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
-      IF(mode=='RUG') THEN
-        champint=EXP(champint)
-        WHERE(NINT(mask)/=1) champint=0.001
-      END IF
-    ELSE
-      SELECT CASE(mode)
-        CASE('RUG');       CALL rugosite(dlon,dlat,champ,rlonv,rlatu,champint,mask)
-        CASE('SIC');       CALL sea_ice (dlon,dlat,champ,rlonv,rlatu,champint)
-        CASE('SST','ALB'); CALL grille_m(dlon,dlat,champ,rlonv,rlatu,champint)
-      END SELECT
+    END IF
+    IF(mode=='RUG') champ=LOG(champ)
+    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
+    IF(mode=='RUG') THEN
+      champint=EXP(champint)
+      WHERE(NINT(mask)/=1) champint=0.001
     END IF
     champtime(:, :, l)=champint
@@ -501,8 +514,212 @@
 !-------------------------------------------------------------------------------
 !
+SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!===============================================================================
+! Purpose:  Compute "phis" just like it would be in start_init_orog.
+!===============================================================================
+! Arguments:
+  REAL,             INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
+  REAL,             INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="start_init_orog0"
+  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
+  REAL               :: lev(1), date, dt, deg2rad
+  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
+  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
+!-------------------------------------------------------------------------------
+  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
+  IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
+  IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
+  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
+  IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
+
+!--- HIGH RESOLUTION OROGRAPHY
+  CALL flininfo(frelf, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
+
+  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
+  CALL flinopen(frelf, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,   &
+                lev, ttm_tmp, itau, date, dt, fid)
+  ALLOCATE(relief_hi(iml_rel,jml_rel))
+  CALL flinget(fid, vrel, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
+  CALL flinclo(fid)
+
+!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
+  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
+  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
+  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
+
+!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
+  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
+  CALL conf_dat2d(vrel, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
+  DEALLOCATE(lon_ini,lat_ini)
+
+!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
+  WRITE(lunout,*)
+  WRITE(lunout,*)'*** Compute surface geopotential ***'
+
+!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
+  CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
+  phis = phis * 9.81
+  phis(iml,:) = phis(1,:)
+  DEALLOCATE(relief_hi,lon_rad,lat_rad)
+
+END SUBROUTINE start_init_orog0
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
+!
+!===============================================================================
+! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
+!          without any call to physics subroutines.
+!===============================================================================
+  IMPLICIT NONE 
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
+  REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
+  REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
+  REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
+  REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="grid_noro0"
+  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
+  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
+  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
+  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
+  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
+  LOGICAL :: masque_lu
+  INTEGER :: i, ii, imdp, imar, iext
+  INTEGER :: j, jj, jmdp, jmar, nn
+  REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
+  REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
+!-------------------------------------------------------------------------------
+  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
+  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
+  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
+  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
+  IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
+  IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
+  iext=imdp/10
+  xpi = ACOS(-1.)
+  rad = 6371229.
+
+!--- ARE WE USING A READ MASK ?
+  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
+  WRITE(lunout,*)'Masque lu: ',masque_lu
+
+!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
+  ALLOCATE(xusn(imdp+2*iext))
+  xusn(1     +iext:imdp  +iext)=xd(:)
+  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
+  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
+
+  ALLOCATE(yusn(jmdp+2))
+  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
+  yusn(2:jmdp+1)=yd(:)
+  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
+
+  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
+  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
+  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
+  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
+  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
+  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
+  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
+  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
+
+!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
+  ALLOCATE(a(imar+1),b(imar+1))
+  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
+  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
+  a(1)=x(1)-(x(2)-x(1))/2.0
+  a(2:imar+1)= b(1:imar)
+
+  ALLOCATE(c(jmar),d(jmar))
+  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
+  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
+  c(1)=y(1)-(y(2)-y(1))/2.0
+  c(2:jmar)=d(1:jmar-1)
+
+!--- INITIALIZATIONS:
+  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
+  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
+
+!--- SUMMATION OVER GRIDPOINT AREA
+  zleny=xpi/REAL(jmdp)*rad
+  xincr=xpi/REAL(jmdp)/2.
+  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
+  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
+  DO ii = 1, imar+1
+    DO jj = 1, jmar
+      DO j = 2,jmdp+1 
+        zlenx  =zleny  *COS(yusn(j))
+        zbordnor=(xincr+c(jj)-yusn(j))*rad
+        zbordsud=(xincr-d(jj)+yusn(j))*rad
+        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
+        IF(weighy/=0) THEN
+          DO i = 2, imdp+2*iext-1
+            zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
+            zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
+            weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
+            IF(weighx/=0)THEN
+              num_tot(ii,jj)=num_tot(ii,jj)+1.0
+              IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
+              weight(ii,jj)=weight(ii,jj)+weighx*weighy
+              zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
+            END IF
+          END DO
+        END IF
+      END DO
+    END DO
+  END DO
+
+!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
+  IF(.NOT.masque_lu) THEN
+    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
+  END IF
+  nn=COUNT(weight(:,1:jmar-1)==0.0)
+  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
+  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
+
+!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
+  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
+  WHERE(mask>=0.1) mask_tmp = 1.
+  WHERE(weight(:,:)/=0.0)
+    zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
+    zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
+  END WHERE
+  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
+
+!--- Values at poles
+  zphi(imar+1,:)=zphi(1,:)
+
+  zweinor=SUM(weight(1:imar,   1),DIM=1)
+  zweisud=SUM(weight(1:imar,jmar),DIM=1)
+  zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
+  zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
+  zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
+
+END SUBROUTINE grid_noro0
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
 FUNCTION year_len(y,cal_in)
 !
 !-------------------------------------------------------------------------------
-  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
@@ -537,5 +754,4 @@
 !
 !-------------------------------------------------------------------------------
-  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
@@ -624,5 +840,4 @@
 !-------------------------------------------------------------------------------
   USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
-  USE print_control_mod, ONLY: lunout
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
@@ -643,2 +858,6 @@
 ! of #ifndef CPP_1D
 END SUBROUTINE limit_netcdf
+
+END MODULE limit
+!
+!*******************************************************************************
Index: LMDZ5/trunk/makelmdz_fcm
===================================================================
--- LMDZ5/trunk/makelmdz_fcm	(revision 2335)
+++ LMDZ5/trunk/makelmdz_fcm	(revision 2336)
@@ -184,12 +184,4 @@
 done
 
-if [[ $code = ce0l && $paramem = mem ]]
-then
-    echo "There is no parallel version of ce0l at the moment."
-    echo "Please compile the sequential version of the code to produce the" \
-	"executable ce0l."
-    exit 1
-fi
-
 ###############################################################
 # path to fcm 
