Index: /trunk/DOC/chantiers/commit_importants.log
===================================================================
--- /trunk/DOC/chantiers/commit_importants.log	(revision 1522)
+++ /trunk/DOC/chantiers/commit_importants.log	(revision 1523)
@@ -1618,2 +1618,29 @@
 * added a "print_control_mod.F90 (as in Earth GCM) in phy_common
 * made similar changes in LMDZ.GENERIC and LMDZ.MARS
+
+**********************
+**** commit_v1523 ****
+**********************
+Ehouarn: More updates to make planetary codes (+Earth) setups converge.
+
+* in dyn3d_common:
+- convmas.F => convmas.F90
+- enercin.F => enercin.F90
+- flumass.F => flumass.F90
+- massbar.F => massbar.F90
+- tourpot.F => tourpot.F90
+- vitvert.F => vitvert.F90
+
+* in misc:
+- move "q_sat" from "dyn3d_common" to "misc" (in Earth model, it is also called
+  by the physics)
+- move "write_field" from "dyn3d_common" to "misc"(may be called from physics
+  or dynamics and depends on neither).
+
+* in phy_common:
+- move "write_field_phy" here since it may be called from any physics package)
+- add module "regular_lonlat_mod" to store global information on lon-lat grid
+
+* in dynlonlat_phylonlat/phy*:
+- turn "iniphysiq.F90" into module "iniphysiq_mod.F90"
+  (and of course adapt gcm.F[90] and 1D models accordingly)
Index: /trunk/LMDZ.COMMON/libf/dyn3d/gcm.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/gcm.F90	(revision 1522)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/gcm.F90	(revision 1523)
@@ -41,6 +41,5 @@
 ! Ehouarn: the following are needed with (parallel) physics:
 #ifdef CPP_PHYS
-!      USE dimphy
-!      USE comgeomphy, ONLY: initcomgeomphy
+  USE iniphysiq_mod, ONLY: iniphysiq
 #endif
 
@@ -447,5 +446,5 @@
 #ifdef CPP_PHYS
          CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, &
-                      rlatu,rlonv,aire,cu,cv,rad,g,r,cpp, &
+                      rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
                       iflag_phys)
 #endif
Index: unk/LMDZ.COMMON/libf/dyn3d_common/convmas.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/convmas.F	(revision 1522)
+++ 	(revision )
@@ -1,61 +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"
-
-      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: /trunk/LMDZ.COMMON/libf/dyn3d_common/convmas.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/convmas.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/convmas.F90	(revision 1523)
@@ -0,0 +1,35 @@
+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"
+!===============================================================================
+! 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: unk/LMDZ.COMMON/libf/dyn3d_common/enercin.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/enercin.F	(revision 1522)
+++ 	(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: /trunk/LMDZ.COMMON/libf/dyn3d_common/enercin.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/enercin.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/enercin.F90	(revision 1523)
@@ -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: unk/LMDZ.COMMON/libf/dyn3d_common/flumass.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/flumass.F	(revision 1522)
+++ 	(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: /trunk/LMDZ.COMMON/libf/dyn3d_common/flumass.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/flumass.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/flumass.F90	(revision 1523)
@@ -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: unk/LMDZ.COMMON/libf/dyn3d_common/massbar.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/massbar.F	(revision 1522)
+++ 	(revision )
@@ -1,101 +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 "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: /trunk/LMDZ.COMMON/libf/dyn3d_common/massbar.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/massbar.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/massbar.F90	(revision 1523)
@@ -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: unk/LMDZ.COMMON/libf/dyn3d_common/massbarxy.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/massbarxy.F	(revision 1522)
+++ 	(revision )
@@ -1,48 +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 "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: /trunk/LMDZ.COMMON/libf/dyn3d_common/massbarxy.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/massbarxy.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/massbarxy.F90	(revision 1523)
@@ -0,0 +1,30 @@
+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 "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: unk/LMDZ.COMMON/libf/dyn3d_common/q_sat.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/q_sat.F	(revision 1522)
+++ 	(revision )
@@ -1,73 +1,0 @@
-!
-! $Header$
-!
-!
-!
-
-      subroutine q_sat(np,temp,pres,qsat)
-!
-      IMPLICIT none
-!======================================================================
-! Autheur(s): Z.X. Li (LMD/CNRS)
-!  reecriture vectorisee par F. Hourdin.
-! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
-!======================================================================
-! Arguments:
-! kelvin---input-R: temperature en Kelvin
-! millibar--input-R: pression en mb
-!
-! q_sat----output-R: vapeur d'eau saturante en kg/kg
-!======================================================================
-!
-      integer np
-      REAL temp(np),pres(np),qsat(np)
-!
-      REAL r2es
-      PARAMETER (r2es=611.14 *18.0153/28.9644)
-!
-      REAL r3les, r3ies, r3es
-      PARAMETER (R3LES=17.269)
-      PARAMETER (R3IES=21.875)
-!
-      REAL r4les, r4ies, r4es
-      PARAMETER (R4LES=35.86)
-      PARAMETER (R4IES=7.66)
-!
-      REAL rtt
-      PARAMETER (rtt=273.16)
-!
-      REAL retv
-      PARAMETER (retv=28.9644/18.0153 - 1.0)
-
-      real zqsat
-      integer ip
-!
-!     ------------------------------------------------------------------
-!
-!
-
-      do ip=1,np
-
-!      write(*,*)'kelvin,millibar=',kelvin,millibar
-!       write(*,*)'temp,pres=',temp(ip),pres(ip)
-!
-         IF (temp(ip) .LE. rtt) THEN
-            r3es = r3ies
-            r4es = r4ies
-         ELSE
-            r3es = r3les
-            r4es = r4les
-         ENDIF
-!
-         zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
-         zqsat=MIN(0.5,ZQSAT)
-         zqsat=zqsat/(1.-retv *zqsat)
-!
-         qsat(ip)= zqsat
-!      write(*,*)'qsat=',qsat(ip)
-
-      enddo
-!
-      RETURN
-      END
-
Index: unk/LMDZ.COMMON/libf/dyn3d_common/tourpot.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/tourpot.F	(revision 1522)
+++ 	(revision )
@@ -1,80 +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"
-
-      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: /trunk/LMDZ.COMMON/libf/dyn3d_common/tourpot.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/tourpot.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/tourpot.F90	(revision 1523)
@@ -0,0 +1,45 @@
+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"
+!===============================================================================
+! 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: unk/LMDZ.COMMON/libf/dyn3d_common/vitvert.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/vitvert.F	(revision 1522)
+++ 	(revision )
@@ -1,53 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE vitvert ( convm , w )
-c
-      USE comvert_mod, ONLY: bp
-
-      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"
-
-      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: /trunk/LMDZ.COMMON/libf/dyn3d_common/vitvert.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/vitvert.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/vitvert.F90	(revision 1523)
@@ -0,0 +1,30 @@
+!
+! $Header$
+!
+SUBROUTINE vitvert (convm, w)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , Fr. Hourdin.
+!-------------------------------------------------------------------------------
+! Purpose: Compute vertical speed at sigma levels.
+  USE comvert_mod, ONLY: bp
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.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: unk/LMDZ.COMMON/libf/dyn3d_common/write_field.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/write_field.F90	(revision 1522)
+++ 	(revision )
@@ -1,326 +1,0 @@
-!
-! $Id: write_field.F90 1279 2009-12-10 09:02:56Z fairhead $
-!
-module write_field
-implicit none
-
-  integer, parameter :: MaxWriteField = 100
-  integer, dimension(MaxWriteField),save :: FieldId
-  integer, dimension(MaxWriteField),save :: FieldVarId
-  integer, dimension(MaxWriteField),save :: FieldIndex
-  character(len=255), dimension(MaxWriteField) ::  FieldName 
-   
-  integer,save :: NbField = 0
-  
-  interface WriteField
-    module procedure WriteField3d,WriteField2d,WriteField1d
-  end interface WriteField
-  contains
-  
-    function GetFieldIndex(name)
-    implicit none
-      integer          :: GetFieldindex
-      character(len=*) :: name
-    
-      character(len=255) :: TrueName
-      integer            :: i
-       
-      
-      TrueName=TRIM(ADJUSTL(name))
-    
-      GetFieldIndex=-1
-      do i=1,NbField
-        if (TrueName==FieldName(i)) then
-          GetFieldIndex=i
-          exit
-        endif
-      enddo
-    end function GetFieldIndex
- 
-    subroutine WriteField3d(name,Field)
-    implicit none
-      character(len=*) :: name
-      real, dimension(:,:,:) :: Field 
-      integer, dimension(3) :: Dim
-      
-      Dim=shape(Field)
-      call WriteField_gen(name,Field,Dim(1),Dim(2),Dim(3))  
-  
-    end subroutine WriteField3d
-    
-    subroutine WriteField2d(name,Field)
-    implicit none
-      character(len=*) :: name
-      real, dimension(:,:) :: Field 
-      integer, dimension(2) :: Dim
-      
-      Dim=shape(Field)
-      call WriteField_gen(name,Field,Dim(1),Dim(2),1)  
-  
-    end subroutine WriteField2d
-    
-    subroutine WriteField1d(name,Field)
-    implicit none
-      character(len=*) :: name
-      real, dimension(:) :: Field 
-      integer, dimension(1) :: Dim
-      
-      Dim=shape(Field)
-      call WriteField_gen(name,Field,Dim(1),1,1)  
-  
-    end subroutine WriteField1d
-        
-    subroutine WriteField_gen(name,Field,dimx,dimy,dimz)
-    implicit none
-    include 'netcdf.inc'
-      character(len=*) :: name
-      integer :: dimx,dimy,dimz
-      real,dimension(dimx,dimy,dimz) :: Field
-      integer,dimension(dimx*dimy*dimz) :: ndex
-      integer :: status
-      integer :: index
-      integer :: start(4)
-      integer :: count(4)
-      
-           
-      Index=GetFieldIndex(name)
-      if (Index==-1) then
-        call CreateNewField(name,dimx,dimy,dimz)
-	Index=GetFieldIndex(name)
-      else
-        FieldIndex(Index)=FieldIndex(Index)+1.
-      endif
-      
-      start(1)=1
-      start(2)=1
-      start(3)=1
-      start(4)=FieldIndex(Index)
-
-      count(1)=dimx
-      count(2)=dimy
-      count(3)=dimz
-      count(4)=1
-
-      status = NF_PUT_VARA_DOUBLE(FieldId(Index),FieldVarId(Index),start,count,Field)
-      status = NF_SYNC(FieldId(Index))
-      
-    end subroutine WriteField_gen
-       
-    subroutine CreateNewField(name,dimx,dimy,dimz)
-    implicit none
-    include 'netcdf.inc'  
-      character(len=*) :: name
-      integer :: dimx,dimy,dimz
-      integer :: TabDim(4)
-      integer :: status
-      
-      
-      NbField=NbField+1
-      FieldName(NbField)=TRIM(ADJUSTL(name))
-      FieldIndex(NbField)=1
-      
-      
-      status = NF_CREATE(TRIM(ADJUSTL(name))//'.nc', NF_CLOBBER, FieldId(NbField))
-      status = NF_DEF_DIM(FieldId(NbField),'X',dimx,TabDim(1))
-      status = NF_DEF_DIM(FieldId(NbField),'Y',dimy,TabDim(2))
-      status = NF_DEF_DIM(FieldId(NbField),'Z',dimz,TabDim(3))
-      status = NF_DEF_DIM(FieldId(NbField),'iter',NF_UNLIMITED,TabDim(4))
-      status = NF_DEF_VAR(FieldId(NbField),FieldName(NbField),NF_DOUBLE,4,TabDim,FieldVarId(NbField))
-      status = NF_ENDDEF(FieldId(NbField))
-
-    end subroutine CreateNewField
-    
-    
-    
-  subroutine write_field1D(name,Field)
-    implicit none
-  
-    integer, parameter :: MaxDim=1
-    character(len=*)   :: name
-    real, dimension(:) :: Field
-    real, dimension(:),allocatable :: New_Field
-    character(len=20) :: str
-    integer, dimension(MaxDim) :: Dim
-    integer :: i,nb
-    integer, parameter :: id=10
-    integer, parameter :: NbCol=4
-    integer :: ColumnSize 
-    integer :: pos
-    character(len=255) :: form
-    character(len=255) :: MaxLen
-    
-    
-    open(unit=id,file=name//'.field',form='formatted',status='replace')
-    write (id,'("----- Field '//name//'",//)')
-    Dim=shape(Field)
-    MaxLen=int2str(len(trim(int2str(Dim(1)))))
-    ColumnSize=20+6+3+len(trim(int2str(Dim(1))))
-    Nb=0
-    Pos=2
-    do i=1,Dim(1)
-      nb=nb+1
-      
-      if (MOD(nb,NbCol)==0) then
-        form='(t'//trim(int2str(pos))// ',i'//trim(MaxLen) //'," ---> ",g22.16,/)'
-        Pos=2
-      else
-        form='(t'//trim(int2str(pos))// ',i'//trim(MaxLen) //'," ---> ",g22.16," | ",)'
-        Pos=Pos+ColumnSize
-      endif
-      write (id,form,advance='no') i,Field(i)
-    enddo
-     
-    close(id)
-
-  end subroutine write_field1D
-
-  subroutine write_field2D(name,Field)
-    implicit none
-  
-    integer, parameter :: MaxDim=2
-    character(len=*)   :: name
-    real, dimension(:,:) :: Field
-    real, dimension(:,:),allocatable :: New_Field
-    character(len=20) :: str
-    integer, dimension(MaxDim) :: Dim
-    integer :: i,j,nb
-    integer, parameter :: id=10
-    integer, parameter :: NbCol=4
-    integer :: ColumnSize 
-    integer :: pos,offset
-    character(len=255) :: form
-    character(len=255) :: spacing
-    
-    open(unit=id,file=name//'.field',form='formatted',status='replace')
-    write (id,'("----- Field '//name//'",//)')
-    
-    Dim=shape(Field)
-    offset=len(trim(int2str(Dim(1))))+len(trim(int2str(Dim(2))))+3
-    ColumnSize=20+6+3+offset
-
-    spacing='(t2,"'//repeat('-',ColumnSize*NbCol)//'")'
-    
-    do i=1,Dim(2)
-      nb=0
-      Pos=2
-      do j=1,Dim(1)
-        nb=nb+1
-      
-        if (MOD(nb,NbCol)==0) then
-          form='(t'//trim(int2str(pos))//            &
-               ',"('//trim(int2str(j))//','          &
-                    //trim(int2str(i))//')",t'       & 
-                    //trim(int2str(pos+offset))     &    
-                    //'," ---> ",g22.16,/)'
-          Pos=2
-        else
-          form='(t'//trim(int2str(pos))//            &
-               ',"('//trim(int2str(j))//','          &
-                    //trim(int2str(i))//')",t'       & 
-                    //trim(int2str(pos+offset))     &    
-                    //'," ---> ",g22.16," | ")'
-          Pos=Pos+ColumnSize
-        endif
-        write (id,form,advance='no') Field(j,i)
-      enddo
-      if (MOD(nb,NbCol)==0) then
-        write (id,spacing)
-      else
-        write (id,'("")')
-        write (id,spacing)
-      endif
-    enddo
-     
-  end subroutine write_field2D
-  
-  subroutine write_field3D(name,Field)
-    implicit none
-  
-    integer, parameter :: MaxDim=3
-    character(len=*)   :: name
-    real, dimension(:,:,:) :: Field
-    real, dimension(:,:,:),allocatable :: New_Field
-    integer, dimension(MaxDim) :: Dim
-    integer :: i,j,k,nb
-    integer, parameter :: id=10
-    integer, parameter :: NbCol=4
-    integer :: ColumnSize 
-    integer :: pos,offset
-    character(len=255) :: form
-    character(len=255) :: spacing
-
-    open(unit=id,file=name//'.field',form='formatted',status='replace')
-    write (id,'("----- Field '//name//'"//)')
-    
-    Dim=shape(Field)
-    offset=len(trim(int2str(Dim(1))))+len(trim(int2str(Dim(2))))+len(trim(int2str(Dim(3))))+4
-    ColumnSize=22+6+3+offset
-
-!    open(unit=id,file=name,form=formatted
-   
-    spacing='(t2,"'//repeat('-',ColumnSize*NbCol)//'")'
-    
-    do i=1,Dim(3)
-    
-      do j=1,Dim(2)
-        nb=0
-        Pos=2
-        
-        do k=1,Dim(1)
-        nb=nb+1
-      
-          if (MOD(nb,NbCol)==0) then
-            form='(t'//trim(int2str(pos))//            &
-                 ',"('//trim(int2str(k))//','          &
-                      //trim(int2str(j))//','          &
-                      //trim(int2str(i))//')",t'       & 
-                      //trim(int2str(pos+offset))      &    
-                      //'," ---> ",g22.16,/)'
-           Pos=2
-          else
-            form='(t'//trim(int2str(pos))//            &
-                 ',"('//trim(int2str(k))//','          &
-                      //trim(int2str(j))//','          &
-                      //trim(int2str(i))//')",t'       & 
-                      //trim(int2str(pos+offset))      &    
-                      //'," ---> ",g22.16," | ")'
-! dépent de l'implémention, sur compaq, c'est necessaire
-!            Pos=Pos+ColumnSize
-          endif
-          write (id,form,advance='no') Field(k,j,i)
-        enddo
-        if (MOD(nb,NbCol)==0) then
-          write (id,spacing)
-        else
-          write (id,'("")')
-          write (id,spacing)
-        endif
-      enddo
-      write (id,spacing)
-    enddo
-    
-    close(id)
-  
-  end subroutine write_field3D  
-  
-  function int2str(int)
-    implicit none
-    integer, parameter :: MaxLen=10
-    integer,intent(in) :: int
-    character(len=MaxLen) :: int2str
-    logical :: flag
-    integer :: i
-    flag=.true.
-    
-    i=int
-    
-    int2str=''
-    do while (flag)
-      int2str=CHAR(MOD(i,10)+48)//int2str
-      i=i/10
-      if (i==0) flag=.false.
-    enddo
-  end function int2str
-
-end module write_field
-  
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1522)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1523)
@@ -14,5 +14,7 @@
       USE parallel_lmdz
       USE infotrac
+#ifdef CPP_PHYS
       USE mod_interface_dyn_phys
+#endif
       USE mod_hallo
       USE Bands
@@ -33,4 +35,5 @@
 ! Ehouarn: the following are needed with (parallel) physics:
 #ifdef CPP_PHYS
+      USE iniphysiq_mod, ONLY: iniphysiq
       USE mod_grid_phy_lmdz
 !      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
@@ -495,5 +498,5 @@
 !     &                iflag_phys)
          CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
-     &                rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,
+     &                rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,
      &                iflag_phys)
 #endif
Index: unk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/iniphysiq.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/iniphysiq.F90	(revision 1522)
+++ 	(revision )
@@ -1,195 +1,0 @@
-
-! $Id: iniphysiq.F90 2225 2015-03-11 14:55:23Z emillour $
-
-
-SUBROUTINE iniphysiq(iim,jjm,nlayer,punjours, pdayref,ptimestep,         &
-                     rlatu,rlonv,aire,cu,cv,                             &
-                     prad,pg,pr,pcpp,iflag_phys)
-  USE dimphy, ONLY: klev ! number of atmospheric levels
-  USE mod_grid_phy_lmdz, ONLY: klon_glo ! number of atmospheric columns
-                                        ! (on full grid)
-  USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
-                                klon_omp_begin, & ! start index of local omp subgrid
-                                klon_omp_end, & ! end index of local omp subgrid
-                                klon_mpi_begin ! start indes of columns (on local mpi grid)
-  USE comgeomphy, ONLY: initcomgeomphy, &
-                        airephy, & ! physics grid area (m2)
-                        cuphy, & ! cu coeff. (u_covariant = cu * u)
-                        cvphy, & ! cv coeff. (v_covariant = cv * v)
-                        rlond, & ! longitudes
-                        rlatd ! latitudes
-  IMPLICIT NONE
-
-  ! =======================================================================
-  ! Initialisation of the physical constants and some positional and
-  ! geometrical arrays for the physics
-  ! =======================================================================
-
-  include "YOMCST.h"
-  include "iniprint.h"
-
-  REAL, INTENT (IN) :: prad ! radius of the planet (m)
-  REAL, INTENT (IN) :: pg ! gravitational acceleration (m/s2)
-  REAL, INTENT (IN) :: pr ! ! reduced gas constant R/mu
-  REAL, INTENT (IN) :: pcpp ! specific heat Cp
-  REAL, INTENT (IN) :: punjours ! length (in s) of a standard day
-  INTEGER, INTENT (IN) :: nlayer ! number of atmospheric layers
-  INTEGER, INTENT (IN) :: iim ! number of atmospheric columns along longitudes
-  INTEGER, INTENT (IN) :: jjm ! number of atompsheric columns along latitudes
-  REAL, INTENT (IN) :: rlatu(jjm+1) ! latitudes of the physics grid
-  REAL, INTENT (IN) :: rlonv(iim+1) ! longitudes of the physics grid
-  REAL, INTENT (IN) :: aire(iim+1,jjm+1) ! area of the dynamics grid (m2)
-  REAL, INTENT (IN) :: cu((iim+1)*(jjm+1)) ! cu coeff. (u_covariant = cu * u)
-  REAL, INTENT (IN) :: cv((iim+1)*jjm) ! cv coeff. (v_covariant = cv * v)
-  INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
-  REAL, INTENT (IN) :: ptimestep !physics time step (s)
-  INTEGER, INTENT (IN) :: iflag_phys ! type of physics to be called
-
-  INTEGER :: ibegin, iend, offset
-  INTEGER :: i,j
-  CHARACTER (LEN=20) :: modname = 'iniphysiq'
-  CHARACTER (LEN=80) :: abort_message
-  REAL :: total_area_phy, total_area_dyn
-
-
-  ! global array, on full physics grid:
-  REAL,ALLOCATABLE :: latfi(:)
-  REAL,ALLOCATABLE :: lonfi(:)
-  REAL,ALLOCATABLE :: cufi(:)
-  REAL,ALLOCATABLE :: cvfi(:)
-  REAL,ALLOCATABLE :: airefi(:)
-
-  IF (nlayer/=klev) THEN
-    WRITE (lunout, *) 'STOP in ', trim(modname)
-    WRITE (lunout, *) 'Problem with dimensions :'
-    WRITE (lunout, *) 'nlayer     = ', nlayer
-    WRITE (lunout, *) 'klev   = ', klev
-    abort_message = ''
-    CALL abort_gcm(modname, abort_message, 1)
-  END IF
-
-  !call init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
-  
-  ! Generate global arrays on full physics grid
-  ALLOCATE(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
-  ALLOCATE(airefi(klon_glo))
-
-  IF (klon_glo>1) THEN ! general case
-    ! North pole
-    latfi(1)=rlatu(1)
-    lonfi(1)=0.
-    cufi(1) = cu(1)
-    cvfi(1) = cv(1)
-    DO j=2,jjm
-      DO i=1,iim
-        latfi((j-2)*iim+1+i)= rlatu(j)
-        lonfi((j-2)*iim+1+i)= rlonv(i)
-        cufi((j-2)*iim+1+i) = cu((j-1)*(iim+1)+i)
-        cvfi((j-2)*iim+1+i) = cv((j-1)*(iim+1)+i)
-      ENDDO
-    ENDDO
-    ! South pole
-    latfi(klon_glo)= rlatu(jjm+1)
-    lonfi(klon_glo)= 0.
-    cufi(klon_glo) = cu((iim+1)*jjm+1)
-    cvfi(klon_glo) = cv((iim+1)*jjm-iim)
-
-    ! build airefi(), mesh area on physics grid
-    CALL gr_dyn_fi(1,iim+1,jjm+1,klon_glo,aire,airefi)
-    ! Poles are single points on physics grid
-    airefi(1)=sum(aire(1:iim,1))
-    airefi(klon_glo)=sum(aire(1:iim,jjm+1))
-
-    ! Sanity check: do total planet area match between physics and dynamics?
-    total_area_dyn=sum(aire(1:iim,1:jjm+1))
-    total_area_phy=sum(airefi(1:klon_glo))
-    IF (total_area_dyn/=total_area_phy) THEN
-      WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
-      WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
-      WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
-      IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
-        ! stop here if the relative difference is more than 0.001%
-        abort_message = 'planet total surface discrepancy'
-        CALL abort_gcm(modname, abort_message, 1)
-      ENDIF
-    ENDIF
-  ELSE ! klon_glo==1, running the 1D model
-    ! just copy over input values
-    latfi(1)=rlatu(1)
-    lonfi(1)=rlonv(1)
-    cufi(1)=cu(1)
-    cvfi(1)=cv(1)
-    airefi(1)=aire(1,1)
-  ENDIF ! of IF (klon_glo>1)
-
-!$OMP PARALLEL 
-  ! Now generate local lon/lat/cu/cv/area arrays
-  CALL initcomgeomphy
-
-  offset = klon_mpi_begin - 1
-  airephy(1:klon_omp) = airefi(offset+klon_omp_begin:offset+klon_omp_end)
-  cuphy(1:klon_omp) = cufi(offset+klon_omp_begin:offset+klon_omp_end)
-  cvphy(1:klon_omp) = cvfi(offset+klon_omp_begin:offset+klon_omp_end)
-  rlond(1:klon_omp) = lonfi(offset+klon_omp_begin:offset+klon_omp_end)
-  rlatd(1:klon_omp) = latfi(offset+klon_omp_begin:offset+klon_omp_end)
-
-  ! Initialize some physical constants
-  call suphec
-
-!$OMP END PARALLEL
-
-  ! check that physical constants set in 'suphec' are coherent
-  ! with values set in the dynamics:
-  IF (rday/=punjours) THEN
-    WRITE (lunout, *) 'iniphysiq: length of day discrepancy!!!'
-    WRITE (lunout, *) '  in the dynamics punjours=', punjours
-    WRITE (lunout, *) '   but in the physics RDAY=', rday
-    IF (abs(rday-punjours)>0.01*punjours) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'length of day discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-
-  IF (rg/=pg) THEN
-    WRITE (lunout, *) 'iniphysiq: gravity discrepancy !!!'
-    WRITE (lunout, *) '     in the dynamics pg=', pg
-    WRITE (lunout, *) '  but in the physics RG=', rg
-    IF (abs(rg-pg)>0.01*pg) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'gravity discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-  IF (ra/=prad) THEN
-    WRITE (lunout, *) 'iniphysiq: planet radius discrepancy !!!'
-    WRITE (lunout, *) '   in the dynamics prad=', prad
-    WRITE (lunout, *) '  but in the physics RA=', ra
-    IF (abs(ra-prad)>0.01*prad) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'planet radius discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-  IF (rd/=pr) THEN
-    WRITE (lunout, *) 'iniphysiq: reduced gas constant discrepancy !!!'
-    WRITE (lunout, *) '     in the dynamics pr=', pr
-    WRITE (lunout, *) '  but in the physics RD=', rd
-    IF (abs(rd-pr)>0.01*pr) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'reduced gas constant discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-  IF (rcpd/=pcpp) THEN
-    WRITE (lunout, *) 'iniphysiq: specific heat discrepancy !!!'
-    WRITE (lunout, *) '     in the dynamics pcpp=', pcpp
-    WRITE (lunout, *) '  but in the physics RCPD=', rcpd
-    IF (abs(rcpd-pcpp)>0.01*pcpp) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'specific heat discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-
-END SUBROUTINE iniphysiq
Index: /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/iniphysiq_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/iniphysiq_mod.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/iniphysiq_mod.F90	(revision 1523)
@@ -0,0 +1,232 @@
+!
+! $Id: iniphysiq.F90 2225 2015-03-11 14:55:23Z emillour $
+!
+MODULE iniphysiq_mod
+
+CONTAINS
+
+SUBROUTINE iniphysiq(iim,jjm,nlayer,punjours, pdayref,ptimestep,         &
+                     rlatu,rlatv,rlonu,rlonv,aire,cu,cv,                 &
+                     prad,pg,pr,pcpp,iflag_phys)
+  USE dimphy, ONLY: klev ! number of atmospheric levels
+  USE mod_grid_phy_lmdz, ONLY: klon_glo ! number of atmospheric columns
+                                        ! (on full grid)
+  USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
+                                klon_omp_begin, & ! start index of local omp subgrid
+                                klon_omp_end, & ! end index of local omp subgrid
+                                klon_mpi_begin ! start indes of columns (on local mpi grid)
+  USE comgeomphy, ONLY: initcomgeomphy, &
+                        airephy, & ! physics grid area (m2)
+                        cuphy, & ! cu coeff. (u_covariant = cu * u)
+                        cvphy, & ! cv coeff. (v_covariant = cv * v)
+                        rlond, & ! longitudes
+                        rlatd ! latitudes
+  USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
+                                 east, west, north, south, &
+                                 north_east, north_west, &
+                                 south_west, south_east
+  USE nrtype, ONLY: pi
+  IMPLICIT NONE
+
+  ! =======================================================================
+  ! Initialisation of the physical constants and some positional and
+  ! geometrical arrays for the physics
+  ! =======================================================================
+
+  include "YOMCST.h"
+  include "iniprint.h"
+
+  REAL, INTENT (IN) :: prad ! radius of the planet (m)
+  REAL, INTENT (IN) :: pg ! gravitational acceleration (m/s2)
+  REAL, INTENT (IN) :: pr ! ! reduced gas constant R/mu
+  REAL, INTENT (IN) :: pcpp ! specific heat Cp
+  REAL, INTENT (IN) :: punjours ! length (in s) of a standard day
+  INTEGER, INTENT (IN) :: nlayer ! number of atmospheric layers
+  INTEGER, INTENT (IN) :: iim ! number of atmospheric columns along longitudes
+  INTEGER, INTENT (IN) :: jjm ! number of atompsheric columns along latitudes
+  REAL, INTENT (IN) :: rlatu(jjm+1) ! latitudes of the physics grid
+  REAL, INTENT (IN) :: rlatv(jjm) ! latitude boundaries of the physics grid
+  REAL, INTENT (IN) :: rlonv(iim+1) ! longitudes of the physics grid
+  REAL, INTENT (IN) :: rlonu(iim+1) ! longitude boundaries of the physics grid
+  REAL, INTENT (IN) :: aire(iim+1,jjm+1) ! area of the dynamics grid (m2)
+  REAL, INTENT (IN) :: cu((iim+1)*(jjm+1)) ! cu coeff. (u_covariant = cu * u)
+  REAL, INTENT (IN) :: cv((iim+1)*jjm) ! cv coeff. (v_covariant = cv * v)
+  INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
+  REAL, INTENT (IN) :: ptimestep !physics time step (s)
+  INTEGER, INTENT (IN) :: iflag_phys ! type of physics to be called
+
+  INTEGER :: ibegin, iend, offset
+  INTEGER :: i,j
+  CHARACTER (LEN=20) :: modname = 'iniphysiq'
+  CHARACTER (LEN=80) :: abort_message
+  REAL :: total_area_phy, total_area_dyn
+
+  ! boundaries, on global grid
+  REAL,ALLOCATABLE :: boundslon_reg(:,:)
+  REAL,ALLOCATABLE :: boundslat_reg(:,:)
+
+  ! global array, on full physics grid:
+  REAL,ALLOCATABLE :: latfi(:)
+  REAL,ALLOCATABLE :: lonfi(:)
+  REAL,ALLOCATABLE :: cufi(:)
+  REAL,ALLOCATABLE :: cvfi(:)
+  REAL,ALLOCATABLE :: airefi(:)
+
+  IF (nlayer/=klev) THEN
+    WRITE (lunout, *) 'STOP in ', trim(modname)
+    WRITE (lunout, *) 'Problem with dimensions :'
+    WRITE (lunout, *) 'nlayer     = ', nlayer
+    WRITE (lunout, *) 'klev   = ', klev
+    abort_message = ''
+    CALL abort_gcm(modname, 'Problem with dimensions', 1)
+  END IF
+
+  !call init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
+  
+  ! init regular global longitude-latitude grid points and boundaries
+  ALLOCATE(boundslon_reg(iim,2))
+  ALLOCATE(boundslat_reg(jjm+1,2))
+  
+  DO i=1,iim
+   boundslon_reg(i,east)=rlonu(i) 
+   boundslon_reg(i,west)=rlonu(i+1) 
+  ENDDO
+
+  boundslat_reg(1,north)= PI/2 
+  boundslat_reg(1,south)= rlatv(1)
+  DO j=2,jjm
+   boundslat_reg(j,north)=rlatv(j-1) 
+   boundslat_reg(j,south)=rlatv(j) 
+  ENDDO
+  boundslat_reg(jjm+1,north)= rlatv(jjm) 
+  boundslat_reg(jjm+1,south)= -PI/2
+
+  ! Write values in module regular_lonlat_mod
+  CALL init_regular_lonlat(iim,jjm+1, rlonv(1:iim), rlatu, &
+                           boundslon_reg, boundslat_reg)
+
+  ! Generate global arrays on full physics grid
+  ALLOCATE(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
+  ALLOCATE(airefi(klon_glo))
+
+  IF (klon_glo>1) THEN ! general case
+    ! North pole
+    latfi(1)=rlatu(1)
+    lonfi(1)=0.
+    cufi(1) = cu(1)
+    cvfi(1) = cv(1)
+    DO j=2,jjm
+      DO i=1,iim
+        latfi((j-2)*iim+1+i)= rlatu(j)
+        lonfi((j-2)*iim+1+i)= rlonv(i)
+        cufi((j-2)*iim+1+i) = cu((j-1)*(iim+1)+i)
+        cvfi((j-2)*iim+1+i) = cv((j-1)*(iim+1)+i)
+      ENDDO
+    ENDDO
+    ! South pole
+    latfi(klon_glo)= rlatu(jjm+1)
+    lonfi(klon_glo)= 0.
+    cufi(klon_glo) = cu((iim+1)*jjm+1)
+    cvfi(klon_glo) = cv((iim+1)*jjm-iim)
+
+    ! build airefi(), mesh area on physics grid
+    CALL gr_dyn_fi(1,iim+1,jjm+1,klon_glo,aire,airefi)
+    ! Poles are single points on physics grid
+    airefi(1)=sum(aire(1:iim,1))
+    airefi(klon_glo)=sum(aire(1:iim,jjm+1))
+
+    ! Sanity check: do total planet area match between physics and dynamics?
+    total_area_dyn=sum(aire(1:iim,1:jjm+1))
+    total_area_phy=sum(airefi(1:klon_glo))
+    IF (total_area_dyn/=total_area_phy) THEN
+      WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
+      WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
+      WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
+      IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
+        ! stop here if the relative difference is more than 0.001%
+        abort_message = 'planet total surface discrepancy'
+        CALL abort_gcm(modname, abort_message, 1)
+      ENDIF
+    ENDIF
+  ELSE ! klon_glo==1, running the 1D model
+    ! just copy over input values
+    latfi(1)=rlatu(1)
+    lonfi(1)=rlonv(1)
+    cufi(1)=cu(1)
+    cvfi(1)=cv(1)
+    airefi(1)=aire(1,1)
+  ENDIF ! of IF (klon_glo>1)
+
+!$OMP PARALLEL 
+  ! Now generate local lon/lat/cu/cv/area arrays
+  CALL initcomgeomphy
+
+  offset = klon_mpi_begin - 1
+  airephy(1:klon_omp) = airefi(offset+klon_omp_begin:offset+klon_omp_end)
+  cuphy(1:klon_omp) = cufi(offset+klon_omp_begin:offset+klon_omp_end)
+  cvphy(1:klon_omp) = cvfi(offset+klon_omp_begin:offset+klon_omp_end)
+  rlond(1:klon_omp) = lonfi(offset+klon_omp_begin:offset+klon_omp_end)
+  rlatd(1:klon_omp) = latfi(offset+klon_omp_begin:offset+klon_omp_end)
+
+  ! Initialize some physical constants
+  call suphec
+
+!$OMP END PARALLEL
+
+  ! check that physical constants set in 'suphec' are coherent
+  ! with values set in the dynamics:
+  IF (rday/=punjours) THEN
+    WRITE (lunout, *) 'iniphysiq: length of day discrepancy!!!'
+    WRITE (lunout, *) '  in the dynamics punjours=', punjours
+    WRITE (lunout, *) '   but in the physics RDAY=', rday
+    IF (abs(rday-punjours)>0.01*punjours) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'length of day discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+
+  IF (rg/=pg) THEN
+    WRITE (lunout, *) 'iniphysiq: gravity discrepancy !!!'
+    WRITE (lunout, *) '     in the dynamics pg=', pg
+    WRITE (lunout, *) '  but in the physics RG=', rg
+    IF (abs(rg-pg)>0.01*pg) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'gravity discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+  IF (ra/=prad) THEN
+    WRITE (lunout, *) 'iniphysiq: planet radius discrepancy !!!'
+    WRITE (lunout, *) '   in the dynamics prad=', prad
+    WRITE (lunout, *) '  but in the physics RA=', ra
+    IF (abs(ra-prad)>0.01*prad) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'planet radius discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+  IF (rd/=pr) THEN
+    WRITE (lunout, *) 'iniphysiq: reduced gas constant discrepancy !!!'
+    WRITE (lunout, *) '     in the dynamics pr=', pr
+    WRITE (lunout, *) '  but in the physics RD=', rd
+    IF (abs(rd-pr)>0.01*pr) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'reduced gas constant discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+  IF (rcpd/=pcpp) THEN
+    WRITE (lunout, *) 'iniphysiq: specific heat discrepancy !!!'
+    WRITE (lunout, *) '     in the dynamics pcpp=', pcpp
+    WRITE (lunout, *) '  but in the physics RCPD=', rcpd
+    IF (abs(rcpd-pcpp)>0.01*pcpp) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'specific heat discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+
+END SUBROUTINE iniphysiq
+
+END MODULE iniphysiq_mod
Index: /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/newstart.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/newstart.F	(revision 1522)
+++ /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/newstart.F	(revision 1523)
@@ -286,5 +286,9 @@
       ELSE
         CALL getin('ysinus',ysinus)
-        tab_cntrl(27) = ysinus
+        IF (ysinus) THEN
+          tab_cntrl(27) = 1
+        ELSE
+          tab_cntrl(27) = 0
+        ENDIF
       ENDIF
 
Index: /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/startvar.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/startvar.F90	(revision 1522)
+++ /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/startvar.F90	(revision 1523)
@@ -363,4 +363,5 @@
 SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in)
       USE comconst_mod
+      USE grid_atob_m, ONLY: rugsoro
 
 !-------------------------------------------------------------------------------
@@ -427,5 +428,5 @@
   ALLOCATE(rugo   (iml  ,jml))
   ALLOCATE(tmp_var(iml-1,jml))
-  CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml,      &
+  CALL rugsoro(lon_rad, lat_rad, relief_hi,      &
        lon_in, lat_in, tmp_var)
   rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
@@ -727,4 +728,5 @@
 
   USE inter_barxy_m, only: inter_barxy
+  USE grid_atob_m, only: grille_m
 
 ! Arguments:
@@ -757,5 +759,5 @@
     CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2(:j2), vtmp)
   ELSE
-    CALL grille_m   (ii, jj,   lon, lat, vari, i1-1, j1, lon1, lat1,     vtmp)
+    CALL grille_m   (lon, lat, vari, lon1, lat1, vtmp)
   END IF
   CALL gr_int_dyn(vtmp, varo, i1-1, j1)
Index: unk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/iniphysiq.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/iniphysiq.F90	(revision 1522)
+++ 	(revision )
@@ -1,195 +1,0 @@
-
-! $Id: iniphysiq.F90 2225 2015-03-11 14:55:23Z emillour $
-
-
-SUBROUTINE iniphysiq(iim,jjm,nlayer,punjours, pdayref,ptimestep,         &
-                     rlatu,rlonv,aire,cu,cv,                             &
-                     prad,pg,pr,pcpp,iflag_phys)
-  USE dimphy, ONLY: klev ! number of atmospheric levels
-  USE mod_grid_phy_lmdz, ONLY: klon_glo ! number of atmospheric columns
-                                        ! (on full grid)
-  USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
-                                klon_omp_begin, & ! start index of local omp subgrid
-                                klon_omp_end, & ! end index of local omp subgrid
-                                klon_mpi_begin ! start indes of columns (on local mpi grid)
-  USE comgeomphy, ONLY: initcomgeomphy, &
-                        airephy, & ! physics grid area (m2)
-                        cuphy, & ! cu coeff. (u_covariant = cu * u)
-                        cvphy, & ! cv coeff. (v_covariant = cv * v)
-                        rlond, & ! longitudes
-                        rlatd ! latitudes
-  IMPLICIT NONE
-
-  ! =======================================================================
-  ! Initialisation of the physical constants and some positional and
-  ! geometrical arrays for the physics
-  ! =======================================================================
-
-  include "YOMCST.h"
-  include "iniprint.h"
-
-  REAL, INTENT (IN) :: prad ! radius of the planet (m)
-  REAL, INTENT (IN) :: pg ! gravitational acceleration (m/s2)
-  REAL, INTENT (IN) :: pr ! ! reduced gas constant R/mu
-  REAL, INTENT (IN) :: pcpp ! specific heat Cp
-  REAL, INTENT (IN) :: punjours ! length (in s) of a standard day
-  INTEGER, INTENT (IN) :: nlayer ! number of atmospheric layers
-  INTEGER, INTENT (IN) :: iim ! number of atmospheric columns along longitudes
-  INTEGER, INTENT (IN) :: jjm ! number of atompsheric columns along latitudes
-  REAL, INTENT (IN) :: rlatu(jjm+1) ! latitudes of the physics grid
-  REAL, INTENT (IN) :: rlonv(iim+1) ! longitudes of the physics grid
-  REAL, INTENT (IN) :: aire(iim+1,jjm+1) ! area of the dynamics grid (m2)
-  REAL, INTENT (IN) :: cu((iim+1)*(jjm+1)) ! cu coeff. (u_covariant = cu * u)
-  REAL, INTENT (IN) :: cv((iim+1)*jjm) ! cv coeff. (v_covariant = cv * v)
-  INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
-  REAL, INTENT (IN) :: ptimestep !physics time step (s)
-  INTEGER, INTENT (IN) :: iflag_phys ! type of physics to be called
-
-  INTEGER :: ibegin, iend, offset
-  INTEGER :: i,j
-  CHARACTER (LEN=20) :: modname = 'iniphysiq'
-  CHARACTER (LEN=80) :: abort_message
-  REAL :: total_area_phy, total_area_dyn
-
-
-  ! global array, on full physics grid:
-  REAL,ALLOCATABLE :: latfi(:)
-  REAL,ALLOCATABLE :: lonfi(:)
-  REAL,ALLOCATABLE :: cufi(:)
-  REAL,ALLOCATABLE :: cvfi(:)
-  REAL,ALLOCATABLE :: airefi(:)
-
-  IF (nlayer/=klev) THEN
-    WRITE (lunout, *) 'STOP in ', trim(modname)
-    WRITE (lunout, *) 'Problem with dimensions :'
-    WRITE (lunout, *) 'nlayer     = ', nlayer
-    WRITE (lunout, *) 'klev   = ', klev
-    abort_message = ''
-    CALL abort_gcm(modname, abort_message, 1)
-  END IF
-
-  !call init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
-  
-  ! Generate global arrays on full physics grid
-  ALLOCATE(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
-  ALLOCATE(airefi(klon_glo))
-
-  IF (klon_glo>1) THEN ! general case
-    ! North pole
-    latfi(1)=rlatu(1)
-    lonfi(1)=0.
-    cufi(1) = cu(1)
-    cvfi(1) = cv(1)
-    DO j=2,jjm
-      DO i=1,iim
-        latfi((j-2)*iim+1+i)= rlatu(j)
-        lonfi((j-2)*iim+1+i)= rlonv(i)
-        cufi((j-2)*iim+1+i) = cu((j-1)*(iim+1)+i)
-        cvfi((j-2)*iim+1+i) = cv((j-1)*(iim+1)+i)
-      ENDDO
-    ENDDO
-    ! South pole
-    latfi(klon_glo)= rlatu(jjm+1)
-    lonfi(klon_glo)= 0.
-    cufi(klon_glo) = cu((iim+1)*jjm+1)
-    cvfi(klon_glo) = cv((iim+1)*jjm-iim)
-
-    ! build airefi(), mesh area on physics grid
-    CALL gr_dyn_fi(1,iim+1,jjm+1,klon_glo,aire,airefi)
-    ! Poles are single points on physics grid
-    airefi(1)=sum(aire(1:iim,1))
-    airefi(klon_glo)=sum(aire(1:iim,jjm+1))
-
-    ! Sanity check: do total planet area match between physics and dynamics?
-    total_area_dyn=sum(aire(1:iim,1:jjm+1))
-    total_area_phy=sum(airefi(1:klon_glo))
-    IF (total_area_dyn/=total_area_phy) THEN
-      WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
-      WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
-      WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
-      IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
-        ! stop here if the relative difference is more than 0.001%
-        abort_message = 'planet total surface discrepancy'
-        CALL abort_gcm(modname, abort_message, 1)
-      ENDIF
-    ENDIF
-  ELSE ! klon_glo==1, running the 1D model
-    ! just copy over input values
-    latfi(1)=rlatu(1)
-    lonfi(1)=rlonv(1)
-    cufi(1)=cu(1)
-    cvfi(1)=cv(1)
-    airefi(1)=aire(1,1)
-  ENDIF ! of IF (klon_glo>1)
-
-!$OMP PARALLEL 
-  ! Now generate local lon/lat/cu/cv/area arrays
-  CALL initcomgeomphy
-
-  offset = klon_mpi_begin - 1
-  airephy(1:klon_omp) = airefi(offset+klon_omp_begin:offset+klon_omp_end)
-  cuphy(1:klon_omp) = cufi(offset+klon_omp_begin:offset+klon_omp_end)
-  cvphy(1:klon_omp) = cvfi(offset+klon_omp_begin:offset+klon_omp_end)
-  rlond(1:klon_omp) = lonfi(offset+klon_omp_begin:offset+klon_omp_end)
-  rlatd(1:klon_omp) = latfi(offset+klon_omp_begin:offset+klon_omp_end)
-
-  ! Initialize some physical constants
-  call suphec
-
-!$OMP END PARALLEL
-
-  ! check that physical constants set in 'suphec' are coherent
-  ! with values set in the dynamics:
-  IF (rday/=punjours) THEN
-    WRITE (lunout, *) 'iniphysiq: length of day discrepancy!!!'
-    WRITE (lunout, *) '  in the dynamics punjours=', punjours
-    WRITE (lunout, *) '   but in the physics RDAY=', rday
-    IF (abs(rday-punjours)>0.01*punjours) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'length of day discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-
-  IF (rg/=pg) THEN
-    WRITE (lunout, *) 'iniphysiq: gravity discrepancy !!!'
-    WRITE (lunout, *) '     in the dynamics pg=', pg
-    WRITE (lunout, *) '  but in the physics RG=', rg
-    IF (abs(rg-pg)>0.01*pg) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'gravity discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-  IF (ra/=prad) THEN
-    WRITE (lunout, *) 'iniphysiq: planet radius discrepancy !!!'
-    WRITE (lunout, *) '   in the dynamics prad=', prad
-    WRITE (lunout, *) '  but in the physics RA=', ra
-    IF (abs(ra-prad)>0.01*prad) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'planet radius discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-  IF (rd/=pr) THEN
-    WRITE (lunout, *) 'iniphysiq: reduced gas constant discrepancy !!!'
-    WRITE (lunout, *) '     in the dynamics pr=', pr
-    WRITE (lunout, *) '  but in the physics RD=', rd
-    IF (abs(rd-pr)>0.01*pr) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'reduced gas constant discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-  IF (rcpd/=pcpp) THEN
-    WRITE (lunout, *) 'iniphysiq: specific heat discrepancy !!!'
-    WRITE (lunout, *) '     in the dynamics pcpp=', pcpp
-    WRITE (lunout, *) '  but in the physics RCPD=', rcpd
-    IF (abs(rcpd-pcpp)>0.01*pcpp) THEN
-        ! stop here if the relative difference is more than 1%
-      abort_message = 'specific heat discrepancy'
-      CALL abort_gcm(modname, abort_message, 1)
-    END IF
-  END IF
-
-END SUBROUTINE iniphysiq
Index: /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/iniphysiq_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/iniphysiq_mod.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/iniphysiq_mod.F90	(revision 1523)
@@ -0,0 +1,232 @@
+!
+! $Id: iniphysiq.F90 2225 2015-03-11 14:55:23Z emillour $
+!
+MODULE iniphysiq_mod
+
+CONTAINS
+
+SUBROUTINE iniphysiq(iim,jjm,nlayer,punjours, pdayref,ptimestep,         &
+                     rlatu,rlatv,rlonu,rlonv,aire,cu,cv,                 &
+                     prad,pg,pr,pcpp,iflag_phys)
+  USE dimphy, ONLY: klev ! number of atmospheric levels
+  USE mod_grid_phy_lmdz, ONLY: klon_glo ! number of atmospheric columns
+                                        ! (on full grid)
+  USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
+                                klon_omp_begin, & ! start index of local omp subgrid
+                                klon_omp_end, & ! end index of local omp subgrid
+                                klon_mpi_begin ! start indes of columns (on local mpi grid)
+  USE comgeomphy, ONLY: initcomgeomphy, &
+                        airephy, & ! physics grid area (m2)
+                        cuphy, & ! cu coeff. (u_covariant = cu * u)
+                        cvphy, & ! cv coeff. (v_covariant = cv * v)
+                        rlond, & ! longitudes
+                        rlatd ! latitudes
+  USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
+                                 east, west, north, south, &
+                                 north_east, north_west, &
+                                 south_west, south_east
+  USE nrtype, ONLY: pi
+  IMPLICIT NONE
+
+  ! =======================================================================
+  ! Initialisation of the physical constants and some positional and
+  ! geometrical arrays for the physics
+  ! =======================================================================
+
+  include "YOMCST.h"
+  include "iniprint.h"
+
+  REAL, INTENT (IN) :: prad ! radius of the planet (m)
+  REAL, INTENT (IN) :: pg ! gravitational acceleration (m/s2)
+  REAL, INTENT (IN) :: pr ! ! reduced gas constant R/mu
+  REAL, INTENT (IN) :: pcpp ! specific heat Cp
+  REAL, INTENT (IN) :: punjours ! length (in s) of a standard day
+  INTEGER, INTENT (IN) :: nlayer ! number of atmospheric layers
+  INTEGER, INTENT (IN) :: iim ! number of atmospheric columns along longitudes
+  INTEGER, INTENT (IN) :: jjm ! number of atompsheric columns along latitudes
+  REAL, INTENT (IN) :: rlatu(jjm+1) ! latitudes of the physics grid
+  REAL, INTENT (IN) :: rlatv(jjm) ! latitude boundaries of the physics grid
+  REAL, INTENT (IN) :: rlonv(iim+1) ! longitudes of the physics grid
+  REAL, INTENT (IN) :: rlonu(iim+1) ! longitude boundaries of the physics grid
+  REAL, INTENT (IN) :: aire(iim+1,jjm+1) ! area of the dynamics grid (m2)
+  REAL, INTENT (IN) :: cu((iim+1)*(jjm+1)) ! cu coeff. (u_covariant = cu * u)
+  REAL, INTENT (IN) :: cv((iim+1)*jjm) ! cv coeff. (v_covariant = cv * v)
+  INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
+  REAL, INTENT (IN) :: ptimestep !physics time step (s)
+  INTEGER, INTENT (IN) :: iflag_phys ! type of physics to be called
+
+  INTEGER :: ibegin, iend, offset
+  INTEGER :: i,j
+  CHARACTER (LEN=20) :: modname = 'iniphysiq'
+  CHARACTER (LEN=80) :: abort_message
+  REAL :: total_area_phy, total_area_dyn
+
+  ! boundaries, on global grid
+  REAL,ALLOCATABLE :: boundslon_reg(:,:)
+  REAL,ALLOCATABLE :: boundslat_reg(:,:)
+
+  ! global array, on full physics grid:
+  REAL,ALLOCATABLE :: latfi(:)
+  REAL,ALLOCATABLE :: lonfi(:)
+  REAL,ALLOCATABLE :: cufi(:)
+  REAL,ALLOCATABLE :: cvfi(:)
+  REAL,ALLOCATABLE :: airefi(:)
+
+  IF (nlayer/=klev) THEN
+    WRITE (lunout, *) 'STOP in ', trim(modname)
+    WRITE (lunout, *) 'Problem with dimensions :'
+    WRITE (lunout, *) 'nlayer     = ', nlayer
+    WRITE (lunout, *) 'klev   = ', klev
+    abort_message = ''
+    CALL abort_gcm(modname, 'Problem with dimensions', 1)
+  END IF
+
+  !call init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
+  
+  ! init regular global longitude-latitude grid points and boundaries
+  ALLOCATE(boundslon_reg(iim,2))
+  ALLOCATE(boundslat_reg(jjm+1,2))
+  
+  DO i=1,iim
+   boundslon_reg(i,east)=rlonu(i) 
+   boundslon_reg(i,west)=rlonu(i+1) 
+  ENDDO
+
+  boundslat_reg(1,north)= PI/2 
+  boundslat_reg(1,south)= rlatv(1)
+  DO j=2,jjm
+   boundslat_reg(j,north)=rlatv(j-1) 
+   boundslat_reg(j,south)=rlatv(j) 
+  ENDDO
+  boundslat_reg(jjm+1,north)= rlatv(jjm) 
+  boundslat_reg(jjm+1,south)= -PI/2
+
+  ! Write values in module regular_lonlat_mod
+  CALL init_regular_lonlat(iim,jjm+1, rlonv(1:iim), rlatu, &
+                           boundslon_reg, boundslat_reg)
+
+  ! Generate global arrays on full physics grid
+  ALLOCATE(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
+  ALLOCATE(airefi(klon_glo))
+
+  IF (klon_glo>1) THEN ! general case
+    ! North pole
+    latfi(1)=rlatu(1)
+    lonfi(1)=0.
+    cufi(1) = cu(1)
+    cvfi(1) = cv(1)
+    DO j=2,jjm
+      DO i=1,iim
+        latfi((j-2)*iim+1+i)= rlatu(j)
+        lonfi((j-2)*iim+1+i)= rlonv(i)
+        cufi((j-2)*iim+1+i) = cu((j-1)*(iim+1)+i)
+        cvfi((j-2)*iim+1+i) = cv((j-1)*(iim+1)+i)
+      ENDDO
+    ENDDO
+    ! South pole
+    latfi(klon_glo)= rlatu(jjm+1)
+    lonfi(klon_glo)= 0.
+    cufi(klon_glo) = cu((iim+1)*jjm+1)
+    cvfi(klon_glo) = cv((iim+1)*jjm-iim)
+
+    ! build airefi(), mesh area on physics grid
+    CALL gr_dyn_fi(1,iim+1,jjm+1,klon_glo,aire,airefi)
+    ! Poles are single points on physics grid
+    airefi(1)=sum(aire(1:iim,1))
+    airefi(klon_glo)=sum(aire(1:iim,jjm+1))
+
+    ! Sanity check: do total planet area match between physics and dynamics?
+    total_area_dyn=sum(aire(1:iim,1:jjm+1))
+    total_area_phy=sum(airefi(1:klon_glo))
+    IF (total_area_dyn/=total_area_phy) THEN
+      WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
+      WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
+      WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
+      IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
+        ! stop here if the relative difference is more than 0.001%
+        abort_message = 'planet total surface discrepancy'
+        CALL abort_gcm(modname, abort_message, 1)
+      ENDIF
+    ENDIF
+  ELSE ! klon_glo==1, running the 1D model
+    ! just copy over input values
+    latfi(1)=rlatu(1)
+    lonfi(1)=rlonv(1)
+    cufi(1)=cu(1)
+    cvfi(1)=cv(1)
+    airefi(1)=aire(1,1)
+  ENDIF ! of IF (klon_glo>1)
+
+!$OMP PARALLEL 
+  ! Now generate local lon/lat/cu/cv/area arrays
+  CALL initcomgeomphy
+
+  offset = klon_mpi_begin - 1
+  airephy(1:klon_omp) = airefi(offset+klon_omp_begin:offset+klon_omp_end)
+  cuphy(1:klon_omp) = cufi(offset+klon_omp_begin:offset+klon_omp_end)
+  cvphy(1:klon_omp) = cvfi(offset+klon_omp_begin:offset+klon_omp_end)
+  rlond(1:klon_omp) = lonfi(offset+klon_omp_begin:offset+klon_omp_end)
+  rlatd(1:klon_omp) = latfi(offset+klon_omp_begin:offset+klon_omp_end)
+
+  ! Initialize some physical constants
+  call suphec
+
+!$OMP END PARALLEL
+
+  ! check that physical constants set in 'suphec' are coherent
+  ! with values set in the dynamics:
+  IF (rday/=punjours) THEN
+    WRITE (lunout, *) 'iniphysiq: length of day discrepancy!!!'
+    WRITE (lunout, *) '  in the dynamics punjours=', punjours
+    WRITE (lunout, *) '   but in the physics RDAY=', rday
+    IF (abs(rday-punjours)>0.01*punjours) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'length of day discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+
+  IF (rg/=pg) THEN
+    WRITE (lunout, *) 'iniphysiq: gravity discrepancy !!!'
+    WRITE (lunout, *) '     in the dynamics pg=', pg
+    WRITE (lunout, *) '  but in the physics RG=', rg
+    IF (abs(rg-pg)>0.01*pg) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'gravity discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+  IF (ra/=prad) THEN
+    WRITE (lunout, *) 'iniphysiq: planet radius discrepancy !!!'
+    WRITE (lunout, *) '   in the dynamics prad=', prad
+    WRITE (lunout, *) '  but in the physics RA=', ra
+    IF (abs(ra-prad)>0.01*prad) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'planet radius discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+  IF (rd/=pr) THEN
+    WRITE (lunout, *) 'iniphysiq: reduced gas constant discrepancy !!!'
+    WRITE (lunout, *) '     in the dynamics pr=', pr
+    WRITE (lunout, *) '  but in the physics RD=', rd
+    IF (abs(rd-pr)>0.01*pr) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'reduced gas constant discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+  IF (rcpd/=pcpp) THEN
+    WRITE (lunout, *) 'iniphysiq: specific heat discrepancy !!!'
+    WRITE (lunout, *) '     in the dynamics pcpp=', pcpp
+    WRITE (lunout, *) '  but in the physics RCPD=', rcpd
+    IF (abs(rcpd-pcpp)>0.01*pcpp) THEN
+        ! stop here if the relative difference is more than 1%
+      abort_message = 'specific heat discrepancy'
+      CALL abort_gcm(modname, abort_message, 1)
+    END IF
+  END IF
+
+END SUBROUTINE iniphysiq
+
+END MODULE iniphysiq_mod
Index: /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/newstart.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/newstart.F	(revision 1522)
+++ /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/newstart.F	(revision 1523)
@@ -291,5 +291,9 @@
       ELSE
         CALL getin('ysinus',ysinus)
-        tab_cntrl(27) = ysinus
+        IF (ysinus) THEN
+          tab_cntrl(27) = 1
+        ELSE
+          tab_cntrl(27) = 0
+        ENDIF
       ENDIF
 
Index: /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/startvar.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/startvar.F90	(revision 1522)
+++ /trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/startvar.F90	(revision 1523)
@@ -364,4 +364,5 @@
 
       USE comconst_mod
+      USE grid_atob_m, ONLY: rugsoro
 
 !-------------------------------------------------------------------------------
@@ -428,5 +429,5 @@
   ALLOCATE(rugo   (iml  ,jml))
   ALLOCATE(tmp_var(iml-1,jml))
-  CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml,      &
+  CALL rugsoro(lon_rad, lat_rad, relief_hi,      &
        lon_in, lat_in, tmp_var)
   rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
@@ -728,4 +729,5 @@
 
   USE inter_barxy_m, only: inter_barxy
+  USE grid_atob_m, only: grille_m
 
 ! Arguments:
@@ -758,5 +760,5 @@
     CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2(:j2), vtmp)
   ELSE
-    CALL grille_m   (ii, jj,   lon, lat, vari, i1-1, j1, lon1, lat1,     vtmp)
+    CALL grille_m   (lon, lat, vari, lon1, lat1, vtmp)
   END IF
   CALL gr_int_dyn(vtmp, varo, i1-1, j1)
Index: /trunk/LMDZ.COMMON/libf/misc/q_sat.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/misc/q_sat.F	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/misc/q_sat.F	(revision 1523)
@@ -0,0 +1,73 @@
+!
+! $Header$
+!
+!
+!
+
+      subroutine q_sat(np,temp,pres,qsat)
+!
+      IMPLICIT none
+!======================================================================
+! Autheur(s): Z.X. Li (LMD/CNRS)
+!  reecriture vectorisee par F. Hourdin.
+! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
+!======================================================================
+! Arguments:
+! kelvin---input-R: temperature en Kelvin
+! millibar--input-R: pression en mb
+!
+! q_sat----output-R: vapeur d'eau saturante en kg/kg
+!======================================================================
+!
+      integer np
+      REAL temp(np),pres(np),qsat(np)
+!
+      REAL r2es
+      PARAMETER (r2es=611.14 *18.0153/28.9644)
+!
+      REAL r3les, r3ies, r3es
+      PARAMETER (R3LES=17.269)
+      PARAMETER (R3IES=21.875)
+!
+      REAL r4les, r4ies, r4es
+      PARAMETER (R4LES=35.86)
+      PARAMETER (R4IES=7.66)
+!
+      REAL rtt
+      PARAMETER (rtt=273.16)
+!
+      REAL retv
+      PARAMETER (retv=28.9644/18.0153 - 1.0)
+
+      real zqsat
+      integer ip
+!
+!     ------------------------------------------------------------------
+!
+!
+
+      do ip=1,np
+
+!      write(*,*)'kelvin,millibar=',kelvin,millibar
+!       write(*,*)'temp,pres=',temp(ip),pres(ip)
+!
+         IF (temp(ip) .LE. rtt) THEN
+            r3es = r3ies
+            r4es = r4ies
+         ELSE
+            r3es = r3les
+            r4es = r4les
+         ENDIF
+!
+         zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
+         zqsat=MIN(0.5,ZQSAT)
+         zqsat=zqsat/(1.-retv *zqsat)
+!
+         qsat(ip)= zqsat
+!      write(*,*)'qsat=',qsat(ip)
+
+      enddo
+!
+      RETURN
+      END
+
Index: /trunk/LMDZ.COMMON/libf/misc/write_field.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/misc/write_field.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/misc/write_field.F90	(revision 1523)
@@ -0,0 +1,326 @@
+!
+! $Id: write_field.F90 1279 2009-12-10 09:02:56Z fairhead $
+!
+module write_field
+implicit none
+
+  integer, parameter :: MaxWriteField = 100
+  integer, dimension(MaxWriteField),save :: FieldId
+  integer, dimension(MaxWriteField),save :: FieldVarId
+  integer, dimension(MaxWriteField),save :: FieldIndex
+  character(len=255), dimension(MaxWriteField) ::  FieldName 
+   
+  integer,save :: NbField = 0
+  
+  interface WriteField
+    module procedure WriteField3d,WriteField2d,WriteField1d
+  end interface WriteField
+  contains
+  
+    function GetFieldIndex(name)
+    implicit none
+      integer          :: GetFieldindex
+      character(len=*) :: name
+    
+      character(len=255) :: TrueName
+      integer            :: i
+       
+      
+      TrueName=TRIM(ADJUSTL(name))
+    
+      GetFieldIndex=-1
+      do i=1,NbField
+        if (TrueName==FieldName(i)) then
+          GetFieldIndex=i
+          exit
+        endif
+      enddo
+    end function GetFieldIndex
+ 
+    subroutine WriteField3d(name,Field)
+    implicit none
+      character(len=*) :: name
+      real, dimension(:,:,:) :: Field 
+      integer, dimension(3) :: Dim
+      
+      Dim=shape(Field)
+      call WriteField_gen(name,Field,Dim(1),Dim(2),Dim(3))  
+  
+    end subroutine WriteField3d
+    
+    subroutine WriteField2d(name,Field)
+    implicit none
+      character(len=*) :: name
+      real, dimension(:,:) :: Field 
+      integer, dimension(2) :: Dim
+      
+      Dim=shape(Field)
+      call WriteField_gen(name,Field,Dim(1),Dim(2),1)  
+  
+    end subroutine WriteField2d
+    
+    subroutine WriteField1d(name,Field)
+    implicit none
+      character(len=*) :: name
+      real, dimension(:) :: Field 
+      integer, dimension(1) :: Dim
+      
+      Dim=shape(Field)
+      call WriteField_gen(name,Field,Dim(1),1,1)  
+  
+    end subroutine WriteField1d
+        
+    subroutine WriteField_gen(name,Field,dimx,dimy,dimz)
+    implicit none
+    include 'netcdf.inc'
+      character(len=*) :: name
+      integer :: dimx,dimy,dimz
+      real,dimension(dimx,dimy,dimz) :: Field
+      integer,dimension(dimx*dimy*dimz) :: ndex
+      integer :: status
+      integer :: index
+      integer :: start(4)
+      integer :: count(4)
+      
+           
+      Index=GetFieldIndex(name)
+      if (Index==-1) then
+        call CreateNewField(name,dimx,dimy,dimz)
+	Index=GetFieldIndex(name)
+      else
+        FieldIndex(Index)=FieldIndex(Index)+1.
+      endif
+      
+      start(1)=1
+      start(2)=1
+      start(3)=1
+      start(4)=FieldIndex(Index)
+
+      count(1)=dimx
+      count(2)=dimy
+      count(3)=dimz
+      count(4)=1
+
+      status = NF_PUT_VARA_DOUBLE(FieldId(Index),FieldVarId(Index),start,count,Field)
+      status = NF_SYNC(FieldId(Index))
+      
+    end subroutine WriteField_gen
+       
+    subroutine CreateNewField(name,dimx,dimy,dimz)
+    implicit none
+    include 'netcdf.inc'  
+      character(len=*) :: name
+      integer :: dimx,dimy,dimz
+      integer :: TabDim(4)
+      integer :: status
+      
+      
+      NbField=NbField+1
+      FieldName(NbField)=TRIM(ADJUSTL(name))
+      FieldIndex(NbField)=1
+      
+      
+      status = NF_CREATE(TRIM(ADJUSTL(name))//'.nc', NF_CLOBBER, FieldId(NbField))
+      status = NF_DEF_DIM(FieldId(NbField),'X',dimx,TabDim(1))
+      status = NF_DEF_DIM(FieldId(NbField),'Y',dimy,TabDim(2))
+      status = NF_DEF_DIM(FieldId(NbField),'Z',dimz,TabDim(3))
+      status = NF_DEF_DIM(FieldId(NbField),'iter',NF_UNLIMITED,TabDim(4))
+      status = NF_DEF_VAR(FieldId(NbField),FieldName(NbField),NF_DOUBLE,4,TabDim,FieldVarId(NbField))
+      status = NF_ENDDEF(FieldId(NbField))
+
+    end subroutine CreateNewField
+    
+    
+    
+  subroutine write_field1D(name,Field)
+    implicit none
+  
+    integer, parameter :: MaxDim=1
+    character(len=*)   :: name
+    real, dimension(:) :: Field
+    real, dimension(:),allocatable :: New_Field
+    character(len=20) :: str
+    integer, dimension(MaxDim) :: Dim
+    integer :: i,nb
+    integer, parameter :: id=10
+    integer, parameter :: NbCol=4
+    integer :: ColumnSize 
+    integer :: pos
+    character(len=255) :: form
+    character(len=255) :: MaxLen
+    
+    
+    open(unit=id,file=name//'.field',form='formatted',status='replace')
+    write (id,'("----- Field '//name//'",//)')
+    Dim=shape(Field)
+    MaxLen=int2str(len(trim(int2str(Dim(1)))))
+    ColumnSize=20+6+3+len(trim(int2str(Dim(1))))
+    Nb=0
+    Pos=2
+    do i=1,Dim(1)
+      nb=nb+1
+      
+      if (MOD(nb,NbCol)==0) then
+        form='(t'//trim(int2str(pos))// ',i'//trim(MaxLen) //'," ---> ",g22.16,/)'
+        Pos=2
+      else
+        form='(t'//trim(int2str(pos))// ',i'//trim(MaxLen) //'," ---> ",g22.16," | ",)'
+        Pos=Pos+ColumnSize
+      endif
+      write (id,form,advance='no') i,Field(i)
+    enddo
+     
+    close(id)
+
+  end subroutine write_field1D
+
+  subroutine write_field2D(name,Field)
+    implicit none
+  
+    integer, parameter :: MaxDim=2
+    character(len=*)   :: name
+    real, dimension(:,:) :: Field
+    real, dimension(:,:),allocatable :: New_Field
+    character(len=20) :: str
+    integer, dimension(MaxDim) :: Dim
+    integer :: i,j,nb
+    integer, parameter :: id=10
+    integer, parameter :: NbCol=4
+    integer :: ColumnSize 
+    integer :: pos,offset
+    character(len=255) :: form
+    character(len=255) :: spacing
+    
+    open(unit=id,file=name//'.field',form='formatted',status='replace')
+    write (id,'("----- Field '//name//'",//)')
+    
+    Dim=shape(Field)
+    offset=len(trim(int2str(Dim(1))))+len(trim(int2str(Dim(2))))+3
+    ColumnSize=20+6+3+offset
+
+    spacing='(t2,"'//repeat('-',ColumnSize*NbCol)//'")'
+    
+    do i=1,Dim(2)
+      nb=0
+      Pos=2
+      do j=1,Dim(1)
+        nb=nb+1
+      
+        if (MOD(nb,NbCol)==0) then
+          form='(t'//trim(int2str(pos))//            &
+               ',"('//trim(int2str(j))//','          &
+                    //trim(int2str(i))//')",t'       & 
+                    //trim(int2str(pos+offset))     &    
+                    //'," ---> ",g22.16,/)'
+          Pos=2
+        else
+          form='(t'//trim(int2str(pos))//            &
+               ',"('//trim(int2str(j))//','          &
+                    //trim(int2str(i))//')",t'       & 
+                    //trim(int2str(pos+offset))     &    
+                    //'," ---> ",g22.16," | ")'
+          Pos=Pos+ColumnSize
+        endif
+        write (id,form,advance='no') Field(j,i)
+      enddo
+      if (MOD(nb,NbCol)==0) then
+        write (id,spacing)
+      else
+        write (id,'("")')
+        write (id,spacing)
+      endif
+    enddo
+     
+  end subroutine write_field2D
+  
+  subroutine write_field3D(name,Field)
+    implicit none
+  
+    integer, parameter :: MaxDim=3
+    character(len=*)   :: name
+    real, dimension(:,:,:) :: Field
+    real, dimension(:,:,:),allocatable :: New_Field
+    integer, dimension(MaxDim) :: Dim
+    integer :: i,j,k,nb
+    integer, parameter :: id=10
+    integer, parameter :: NbCol=4
+    integer :: ColumnSize 
+    integer :: pos,offset
+    character(len=255) :: form
+    character(len=255) :: spacing
+
+    open(unit=id,file=name//'.field',form='formatted',status='replace')
+    write (id,'("----- Field '//name//'"//)')
+    
+    Dim=shape(Field)
+    offset=len(trim(int2str(Dim(1))))+len(trim(int2str(Dim(2))))+len(trim(int2str(Dim(3))))+4
+    ColumnSize=22+6+3+offset
+
+!    open(unit=id,file=name,form=formatted
+   
+    spacing='(t2,"'//repeat('-',ColumnSize*NbCol)//'")'
+    
+    do i=1,Dim(3)
+    
+      do j=1,Dim(2)
+        nb=0
+        Pos=2
+        
+        do k=1,Dim(1)
+        nb=nb+1
+      
+          if (MOD(nb,NbCol)==0) then
+            form='(t'//trim(int2str(pos))//            &
+                 ',"('//trim(int2str(k))//','          &
+                      //trim(int2str(j))//','          &
+                      //trim(int2str(i))//')",t'       & 
+                      //trim(int2str(pos+offset))      &    
+                      //'," ---> ",g22.16,/)'
+           Pos=2
+          else
+            form='(t'//trim(int2str(pos))//            &
+                 ',"('//trim(int2str(k))//','          &
+                      //trim(int2str(j))//','          &
+                      //trim(int2str(i))//')",t'       & 
+                      //trim(int2str(pos+offset))      &    
+                      //'," ---> ",g22.16," | ")'
+! dépent de l'implémention, sur compaq, c'est necessaire
+!            Pos=Pos+ColumnSize
+          endif
+          write (id,form,advance='no') Field(k,j,i)
+        enddo
+        if (MOD(nb,NbCol)==0) then
+          write (id,spacing)
+        else
+          write (id,'("")')
+          write (id,spacing)
+        endif
+      enddo
+      write (id,spacing)
+    enddo
+    
+    close(id)
+  
+  end subroutine write_field3D  
+  
+  function int2str(int)
+    implicit none
+    integer, parameter :: MaxLen=10
+    integer,intent(in) :: int
+    character(len=MaxLen) :: int2str
+    logical :: flag
+    integer :: i
+    flag=.true.
+    
+    i=int
+    
+    int2str=''
+    do while (flag)
+      int2str=CHAR(MOD(i,10)+48)//int2str
+      i=i/10
+      if (i==0) flag=.false.
+    enddo
+  end function int2str
+
+end module write_field
+  
Index: /trunk/LMDZ.COMMON/libf/phy_common/regular_lonlat_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/phy_common/regular_lonlat_mod.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/phy_common/regular_lonlat_mod.F90	(revision 1523)
@@ -0,0 +1,55 @@
+MODULE regular_lonlat_mod
+
+  ! Store information on the global physics grid
+  ! for a regular (longitude-latitude) grid
+  
+  INTEGER, PARAMETER :: north_east=1     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: north_west=2     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south_west=3     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south_east=4     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: east=1           ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: west=2           ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: north=1          ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south=2          ! boundaries of regular lontlat
+
+! Global definition, shared by all threads
+! Do not set threadprivate directives
+
+  REAL,SAVE,ALLOCATABLE :: lon_reg(:)      ! value of longitude cell (rad)
+
+  REAL,SAVE,ALLOCATABLE :: lat_reg(:)      ! value of longitude cell (rad)
+
+  REAL,SAVE,ALLOCATABLE :: boundslon_reg(:,:)      ! value of boundaries cell (1=>east, 2=>west)(rad)
+
+  REAL,SAVE,ALLOCATABLE :: boundslat_reg(:,:)      ! value of longitude cell (1=>north, 2=>south)(rad)
+
+
+CONTAINS
+  
+  SUBROUTINE init_regular_lonlat(nbp_lon, nbp_lat, lon_reg_, lat_reg_, boundslon_reg_, boundslat_reg_)
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: nbp_lon
+    INTEGER,INTENT(IN) :: nbp_lat
+    REAL,   INTENT(IN) :: lon_reg_(nbp_lon)  
+    REAL,   INTENT(IN) :: lat_reg_(nbp_lat) 
+    REAL,   INTENT(IN) :: boundslon_reg_(nbp_lon,2) 
+    REAL,   INTENT(IN) :: boundslat_reg_(nbp_lat,2)  
+
+    
+    ALLOCATE(lon_reg(nbp_lon))
+    lon_reg(:)=lon_reg_(:)
+
+    ALLOCATE(lat_reg(nbp_lat)) 
+    lat_reg(:)=lat_reg_(:)
+
+    ALLOCATE(boundslon_reg(nbp_lon,2)) 
+    boundslon_reg(:,:)=boundslon_reg_(:,:)
+
+    ALLOCATE(boundslat_reg(nbp_lat,2))  
+    boundslat_reg(:,:)=boundslat_reg_(:,:)
+
+  
+  END SUBROUTINE init_regular_lonlat    
+
+END MODULE regular_lonlat_mod
+
Index: /trunk/LMDZ.COMMON/libf/phy_common/write_field_phy.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/phy_common/write_field_phy.F90	(revision 1523)
+++ /trunk/LMDZ.COMMON/libf/phy_common/write_field_phy.F90	(revision 1523)
@@ -0,0 +1,36 @@
+!
+! $Header$
+!
+MODULE write_field_phy
+
+  CONTAINS 
+ 
+    SUBROUTINE WriteField_phy(name,Field,ll)
+    USE dimphy
+    USE mod_phys_lmdz_para
+    USE mod_grid_phy_lmdz
+    USE Write_Field
+    
+    IMPLICIT NONE
+    include 'dimensions.h'
+    include 'paramet.h'
+
+    character(len=*)   :: name
+    INTEGER :: ll
+    real, dimension(klon_omp,ll) :: Field
+    real,save,allocatable :: Field_tmp(:,:)
+    real, dimension(klon_glo,ll):: New_Field
+    real, dimension(iim,jjp1,ll):: Field_2d
+
+    CALL Gather(Field,New_Field)
+!$OMP MASTER
+    IF (is_mpi_root) THEN	
+      CALL Grid1Dto2D_glo(New_Field,Field_2D)
+      CALL WriteField(name,Field_2d)
+    ENDIF
+!$OMP END MASTER
+
+  
+   END SUBROUTINE WriteField_phy
+ 
+ END MODULE write_field_phy
Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 1522)
+++ /trunk/LMDZ.GENERIC/README	(revision 1523)
@@ -1135,2 +1135,8 @@
 - Adapted "create_make_gcm" and "makegcm*" scripts accordingly
   (deleted obsolete makegcm_g95)
+
+== 28/03/2016 == EM
+- Added module "regular_lonlat_mod.F90" (to store information on global
+  lon-lat grid) in phy_common.
+- Turned iniphysiq (in dynlonlat_phylonlat/phymars)into module
+  "iniphysiq_mod.F90".
Index: /trunk/LMDZ.GENERIC/libf/dyn3d/gcm.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/dyn3d/gcm.F	(revision 1522)
+++ /trunk/LMDZ.GENERIC/libf/dyn3d/gcm.F	(revision 1523)
@@ -12,4 +12,5 @@
      .		apdiss,purmats,physic,apphys
       USE temps_mod, ONLY: day_ini,day_end,itaufin,dt
+      USE iniphysiq_mod, ONLY: iniphysiq
       IMPLICIT NONE
 
@@ -260,6 +261,7 @@
 !         CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
          CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys,
-     &                rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,
-     &                1)
+     &                  rlatu,rlatv,rlonu,rlonv,
+     &                  aire,cu,cv,rad,g,r,cpp,
+     &                  1)
 !     &                iflag_phys)
 !#endif
Index: unk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq.F90	(revision 1522)
+++ 	(revision )
@@ -1,145 +1,0 @@
-subroutine iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep,           &
-                     rlatu,rlonv,aire,cu,cv,                             &
-                     prad,pg,pr,pcpp,iflag_phys)
-
-use dimphy, only : klev ! number of atmospheric levels
-use mod_grid_phy_lmdz, only : klon_glo ! number of atmospheric columns
-                                       ! (on full grid)
-use mod_phys_lmdz_para, only : klon_omp, & ! number of columns (on local omp grid)
-                               klon_omp_begin, & ! start index of local omp subgrid
-                               klon_omp_end, & ! end index of local omp subgrid
-                               klon_mpi_begin ! start indes of columns (on local mpi grid)
-
-use comgeomphy, only : initcomgeomphy, &
-                       airephy, & ! physics grid area (m2)
-                       cuphy, & ! cu coeff. (u_covariant = cu * u)
-                       cvphy, & ! cv coeff. (v_covariant = cv * v)
-                       rlond, & ! longitudes
-                       rlatd ! latitudes
-use infotrac, only : nqtot ! number of advected tracers
-use planete_mod, only: ini_planete_mod
-USE comvert_mod, ONLY: ap,bp,preff
-
-implicit none
-include "dimensions.h"
-include "iniprint.h"
-
-real,intent(in) :: prad ! radius of the planet (m)
-real,intent(in) :: pg ! gravitational acceleration (m/s2)
-real,intent(in) :: pr ! ! reduced gas constant R/mu
-real,intent(in) :: pcpp ! specific heat Cp
-real,intent(in) :: punjours ! length (in s) of a standard day
-!integer,intent(in) :: ngrid ! number of horizontal grid points in the physics (full grid)
-integer,intent(in) :: nlayer ! number of atmospheric layers
-integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes
-integer,intent(in) :: jj  ! number of atompsheric columns along latitudes
-real,intent(in) :: rlatu(jj+1) ! latitudes of the dynamics U grid
-real,intent(in) :: rlonv(ii+1) ! longitudes of the dynamics V grid
-real,intent(in) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2)
-real,intent(in) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
-real,intent(in) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
-integer,intent(in) :: pdayref ! reference day of for the simulation
-real,intent(in) :: ptimestep !physics time step (s)
-integer,intent(in) :: iflag_phys ! type of physics to be called
-
-integer :: ibegin,iend,offset
-integer :: i,j
-character(len=20) :: modname='iniphysiq'
-character(len=80) :: abort_message
-real :: total_area_phy, total_area_dyn
-
-
-! global array, on full physics grid:
-real,allocatable :: latfi(:)
-real,allocatable :: lonfi(:)
-real,allocatable :: cufi(:)
-real,allocatable :: cvfi(:)
-real,allocatable :: airefi(:)
-
-IF (nlayer.NE.klev) THEN
-  write(*,*) 'STOP in ',trim(modname)
-  write(*,*) 'Problem with dimensions :'
-  write(*,*) 'nlayer     = ',nlayer
-  write(*,*) 'klev   = ',klev
-  abort_message = ''
-  CALL abort_gcm (modname,abort_message,1)
-ENDIF
-
-!IF (ngrid.NE.klon_glo) THEN
-!  write(*,*) 'STOP in ',trim(modname)
-!  write(*,*) 'Problem with dimensions :'
-!  write(*,*) 'ngrid     = ',ngrid
-!  write(*,*) 'klon   = ',klon_glo
-!  abort_message = ''
-!  CALL abort_gcm (modname,abort_message,1)
-!ENDIF
-
-!call init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
-
-! Generate global arrays on full physics grid
-allocate(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
-latfi(1)=rlatu(1)
-lonfi(1)=0.
-cufi(1) = cu(1)
-cvfi(1) = cv(1)
-DO j=2,jj
-  DO i=1,ii
-    latfi((j-2)*ii+1+i)= rlatu(j)
-    lonfi((j-2)*ii+1+i)= rlonv(i)
-    cufi((j-2)*ii+1+i) = cu((j-1)*(ii+1)+i)
-    cvfi((j-2)*ii+1+i) = cv((j-1)*(ii+1)+i)
-  ENDDO
-ENDDO
-latfi(klon_glo)= rlatu(jj+1)
-lonfi(klon_glo)= 0.
-cufi(klon_glo) = cu((ii+1)*jj+1)
-cvfi(klon_glo) = cv((ii+1)*jj-ii)
-
-! build airefi(), mesh area on physics grid
-allocate(airefi(klon_glo))
-CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi)
-! Poles are single points on physics grid
-airefi(1)=sum(aire(1:ii,1))
-airefi(klon_glo)=sum(aire(1:ii,jj+1))
-
-! Sanity check: do total planet area match between physics and dynamics?
-total_area_dyn=sum(aire(1:ii,1:jj+1))
-total_area_phy=sum(airefi(1:klon_glo))
-IF (total_area_dyn/=total_area_phy) THEN
-  WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
-  WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
-  WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
-  IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
-    ! stop here if the relative difference is more than 0.001%
-    abort_message = 'planet total surface discrepancy'
-    CALL abort_gcm(modname, abort_message, 1)
-  ENDIF
-ENDIF
-
-
-!$OMP PARALLEL 
-! Now generate local lon/lat/cu/cv/area arrays
-call initcomgeomphy
-
-!!!!$OMP PARALLEL PRIVATE(ibegin,iend) & 
-!!!	!$OMP SHARED(airefi,cufi,cvfi,lonfi,latfi)
-
-offset=klon_mpi_begin-1
-airephy(1:klon_omp)=airefi(offset+klon_omp_begin:offset+klon_omp_end)
-cuphy(1:klon_omp)=cufi(offset+klon_omp_begin:offset+klon_omp_end)
-cvphy(1:klon_omp)=cvfi(offset+klon_omp_begin:offset+klon_omp_end)
-rlond(1:klon_omp)=lonfi(offset+klon_omp_begin:offset+klon_omp_end)
-rlatd(1:klon_omp)=latfi(offset+klon_omp_begin:offset+klon_omp_end)
-
-! copy over preff , ap() and bp() 
-call ini_planete_mod(nlayer,preff,ap,bp)
-
-! copy some fundamental parameters to physics 
-! and do some initializations 
-call inifis(klon_omp,nlayer,nqtot,pdayref,punjours,ptimestep, &
-            rlatd,rlond,airephy,prad,pg,pr,pcpp)
-
-!$OMP END PARALLEL
-
-
-end subroutine iniphysiq
Index: /trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90	(revision 1523)
+++ /trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90	(revision 1523)
@@ -0,0 +1,186 @@
+MODULE iniphysiq_mod
+
+CONTAINS
+
+subroutine iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep,           &
+                     rlatu,rlatv,rlonu,rlonv,aire,cu,cv,                 &
+                     prad,pg,pr,pcpp,iflag_phys)
+
+use dimphy, only : klev ! number of atmospheric levels
+use mod_grid_phy_lmdz, only : klon_glo ! number of atmospheric columns
+                                       ! (on full grid)
+use mod_phys_lmdz_para, only : klon_omp, & ! number of columns (on local omp grid)
+                               klon_omp_begin, & ! start index of local omp subgrid
+                               klon_omp_end, & ! end index of local omp subgrid
+                               klon_mpi_begin ! start indes of columns (on local mpi grid)
+
+use comgeomphy, only : initcomgeomphy, &
+                       airephy, & ! physics grid area (m2)
+                       cuphy, & ! cu coeff. (u_covariant = cu * u)
+                       cvphy, & ! cv coeff. (v_covariant = cv * v)
+                       rlond, & ! longitudes
+                       rlatd ! latitudes
+use infotrac, only : nqtot ! number of advected tracers
+use planete_mod, only: ini_planete_mod
+USE comvert_mod, ONLY: ap,bp,preff
+use regular_lonlat_mod, only: init_regular_lonlat, &
+                              east, west, north, south, &
+                              north_east, north_west, &
+                              south_west, south_east
+
+implicit none
+include "dimensions.h"
+include "iniprint.h"
+
+real,intent(in) :: prad ! radius of the planet (m)
+real,intent(in) :: pg ! gravitational acceleration (m/s2)
+real,intent(in) :: pr ! ! reduced gas constant R/mu
+real,intent(in) :: pcpp ! specific heat Cp
+real,intent(in) :: punjours ! length (in s) of a standard day
+!integer,intent(in) :: ngrid ! number of horizontal grid points in the physics (full grid)
+integer,intent(in) :: nlayer ! number of atmospheric layers
+integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes
+integer,intent(in) :: jj  ! number of atompsheric columns along latitudes
+real,intent(in) :: rlatu(jj+1) ! latitudes of the physics grid
+real,intent(in) :: rlatv(jj) ! latitude boundaries of the physics grid
+real,intent(in) :: rlonv(ii+1) ! longitudes of the physics grid
+real,intent(in) :: rlonu(ii+1) ! longitude boundaries of the physics grid
+real,intent(in) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2)
+real,intent(in) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
+real,intent(in) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
+integer,intent(in) :: pdayref ! reference day of for the simulation
+real,intent(in) :: ptimestep !physics time step (s)
+integer,intent(in) :: iflag_phys ! type of physics to be called
+
+integer :: ibegin,iend,offset
+integer :: i,j
+character(len=20) :: modname='iniphysiq'
+character(len=80) :: abort_message
+real :: total_area_phy, total_area_dyn
+real :: pi
+
+! boundaries, on global grid
+real,allocatable :: boundslon_reg(:,:)
+real,allocatable :: boundslat_reg(:,:)
+
+! global array, on full physics grid:
+real,allocatable :: latfi(:)
+real,allocatable :: lonfi(:)
+real,allocatable :: cufi(:)
+real,allocatable :: cvfi(:)
+real,allocatable :: airefi(:)
+
+pi=2.*asin(1.0)
+
+IF (nlayer.NE.klev) THEN
+  write(*,*) 'STOP in ',trim(modname)
+  write(*,*) 'Problem with dimensions :'
+  write(*,*) 'nlayer     = ',nlayer
+  write(*,*) 'klev   = ',klev
+  abort_message = ''
+  CALL abort_gcm (modname,abort_message,1)
+ENDIF
+
+!IF (ngrid.NE.klon_glo) THEN
+!  write(*,*) 'STOP in ',trim(modname)
+!  write(*,*) 'Problem with dimensions :'
+!  write(*,*) 'ngrid     = ',ngrid
+!  write(*,*) 'klon   = ',klon_glo
+!  abort_message = ''
+!  CALL abort_gcm (modname,abort_message,1)
+!ENDIF
+
+!call init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
+
+! init regular global longitude-latitude grid points and boundaries
+ALLOCATE(boundslon_reg(ii,2))
+ALLOCATE(boundslat_reg(jj+1,2))
+  
+DO i=1,ii
+   boundslon_reg(i,east)=rlonu(i) 
+   boundslon_reg(i,west)=rlonu(i+1) 
+ENDDO
+
+boundslat_reg(1,north)= PI/2 
+boundslat_reg(1,south)= rlatv(1)
+DO j=2,jj
+   boundslat_reg(j,north)=rlatv(j-1) 
+   boundslat_reg(j,south)=rlatv(j) 
+ENDDO
+boundslat_reg(jj+1,north)= rlatv(jj) 
+boundslat_reg(jj+1,south)= -PI/2
+
+! Write values in module regular_lonlat_mod
+CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, &
+                         boundslon_reg, boundslat_reg)
+
+! Generate global arrays on full physics grid
+allocate(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
+latfi(1)=rlatu(1)
+lonfi(1)=0.
+cufi(1) = cu(1)
+cvfi(1) = cv(1)
+DO j=2,jj
+  DO i=1,ii
+    latfi((j-2)*ii+1+i)= rlatu(j)
+    lonfi((j-2)*ii+1+i)= rlonv(i)
+    cufi((j-2)*ii+1+i) = cu((j-1)*(ii+1)+i)
+    cvfi((j-2)*ii+1+i) = cv((j-1)*(ii+1)+i)
+  ENDDO
+ENDDO
+latfi(klon_glo)= rlatu(jj+1)
+lonfi(klon_glo)= 0.
+cufi(klon_glo) = cu((ii+1)*jj+1)
+cvfi(klon_glo) = cv((ii+1)*jj-ii)
+
+! build airefi(), mesh area on physics grid
+allocate(airefi(klon_glo))
+CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi)
+! Poles are single points on physics grid
+airefi(1)=sum(aire(1:ii,1))
+airefi(klon_glo)=sum(aire(1:ii,jj+1))
+
+! Sanity check: do total planet area match between physics and dynamics?
+total_area_dyn=sum(aire(1:ii,1:jj+1))
+total_area_phy=sum(airefi(1:klon_glo))
+IF (total_area_dyn/=total_area_phy) THEN
+  WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
+  WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
+  WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
+  IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
+    ! stop here if the relative difference is more than 0.001%
+    abort_message = 'planet total surface discrepancy'
+    CALL abort_gcm(modname, abort_message, 1)
+  ENDIF
+ENDIF
+
+
+!$OMP PARALLEL 
+! Now generate local lon/lat/cu/cv/area arrays
+call initcomgeomphy
+
+!!!!$OMP PARALLEL PRIVATE(ibegin,iend) & 
+!!!	!$OMP SHARED(airefi,cufi,cvfi,lonfi,latfi)
+
+offset=klon_mpi_begin-1
+airephy(1:klon_omp)=airefi(offset+klon_omp_begin:offset+klon_omp_end)
+cuphy(1:klon_omp)=cufi(offset+klon_omp_begin:offset+klon_omp_end)
+cvphy(1:klon_omp)=cvfi(offset+klon_omp_begin:offset+klon_omp_end)
+rlond(1:klon_omp)=lonfi(offset+klon_omp_begin:offset+klon_omp_end)
+rlatd(1:klon_omp)=latfi(offset+klon_omp_begin:offset+klon_omp_end)
+
+! copy over preff , ap() and bp() 
+call ini_planete_mod(nlayer,preff,ap,bp)
+
+! copy some fundamental parameters to physics 
+! and do some initializations 
+call inifis(klon_omp,nlayer,nqtot,pdayref,punjours,ptimestep, &
+            rlatd,rlond,airephy,prad,pg,pr,pcpp)
+
+!$OMP END PARALLEL
+
+
+end subroutine iniphysiq
+
+
+END MODULE iniphysiq_mod
Index: /trunk/LMDZ.GENERIC/libf/phy_common/regular_lonlat_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phy_common/regular_lonlat_mod.F90	(revision 1523)
+++ /trunk/LMDZ.GENERIC/libf/phy_common/regular_lonlat_mod.F90	(revision 1523)
@@ -0,0 +1,55 @@
+MODULE regular_lonlat_mod
+
+  ! Store information on the global physics grid
+  ! for a regular (longitude-latitude) grid
+  
+  INTEGER, PARAMETER :: north_east=1     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: north_west=2     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south_west=3     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south_east=4     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: east=1           ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: west=2           ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: north=1          ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south=2          ! boundaries of regular lontlat
+
+! Global definition, shared by all threads
+! Do not set threadprivate directives
+
+  REAL,SAVE,ALLOCATABLE :: lon_reg(:)      ! value of longitude cell (rad)
+
+  REAL,SAVE,ALLOCATABLE :: lat_reg(:)      ! value of longitude cell (rad)
+
+  REAL,SAVE,ALLOCATABLE :: boundslon_reg(:,:)      ! value of boundaries cell (1=>east, 2=>west)(rad)
+
+  REAL,SAVE,ALLOCATABLE :: boundslat_reg(:,:)      ! value of longitude cell (1=>north, 2=>south)(rad)
+
+
+CONTAINS
+  
+  SUBROUTINE init_regular_lonlat(nbp_lon, nbp_lat, lon_reg_, lat_reg_, boundslon_reg_, boundslat_reg_)
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: nbp_lon
+    INTEGER,INTENT(IN) :: nbp_lat
+    REAL,   INTENT(IN) :: lon_reg_(nbp_lon)  
+    REAL,   INTENT(IN) :: lat_reg_(nbp_lat) 
+    REAL,   INTENT(IN) :: boundslon_reg_(nbp_lon,2) 
+    REAL,   INTENT(IN) :: boundslat_reg_(nbp_lat,2)  
+
+    
+    ALLOCATE(lon_reg(nbp_lon))
+    lon_reg(:)=lon_reg_(:)
+
+    ALLOCATE(lat_reg(nbp_lat)) 
+    lat_reg(:)=lat_reg_(:)
+
+    ALLOCATE(boundslon_reg(nbp_lon,2)) 
+    boundslon_reg(:,:)=boundslon_reg_(:,:)
+
+    ALLOCATE(boundslat_reg(nbp_lat,2))  
+    boundslat_reg(:,:)=boundslat_reg_(:,:)
+
+  
+  END SUBROUTINE init_regular_lonlat    
+
+END MODULE regular_lonlat_mod
+
Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 1522)
+++ /trunk/LMDZ.MARS/README	(revision 1523)
@@ -2219,2 +2219,8 @@
 - Adapted "create_make_gcm" and "makegcm*" scripts accordingly
   (deleted obsolete makegcm_g95)
+
+== 28/03/2016 == EM
+- Added module "regular_lonlat_mod.F90" (to store information on global
+  lon-lat grid) in phy_common.
+- Turned iniphysiq (in dynlonlat_phylonlat/phymars)into module
+  "iniphysiq_mod.F90".
Index: /trunk/LMDZ.MARS/libf/dyn3d/gcm.F
===================================================================
--- /trunk/LMDZ.MARS/libf/dyn3d/gcm.F	(revision 1522)
+++ /trunk/LMDZ.MARS/libf/dyn3d/gcm.F	(revision 1523)
@@ -12,4 +12,5 @@
      .		apdiss,purmats,physic,apphys
       USE temps_mod, ONLY: day_ini,day_end,dt,itaufin
+      USE iniphysiq_mod, ONLY: iniphysiq
       IMPLICIT NONE
 
@@ -254,6 +255,7 @@
 !         CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
          CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys,
-     &                rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,
-     &                1)
+     &                  rlatu,rlatv,rlonu,rlonv,
+     &                  aire,cu,cv,rad,g,r,cpp,
+     &                  1)
 !     &                iflag_phys)
 !#endif
Index: unk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/iniphysiq.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/iniphysiq.F90	(revision 1522)
+++ 	(revision )
@@ -1,142 +1,0 @@
-subroutine iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep,           &
-                     rlatu,rlonv,aire,cu,cv,                             &
-                     prad,pg,pr,pcpp,iflag_phys)
-
-use dimphy, only : klev ! number of atmospheric levels
-use mod_grid_phy_lmdz, only : klon_glo ! number of atmospheric columns
-                                       ! (on full grid)
-use mod_phys_lmdz_para, only : klon_omp, & ! number of columns (on local omp grid)
-                               klon_omp_begin, & ! start index of local omp subgrid
-                               klon_omp_end, & ! end index of local omp subgrid
-                               klon_mpi_begin ! start indes of columns (on local mpi grid)
-
-use comgeomphy, only : initcomgeomphy, &
-                       airephy, & ! physics grid area (m2)
-                       cuphy, & ! cu coeff. (u_covariant = cu * u)
-                       cvphy, & ! cv coeff. (v_covariant = cv * v)
-                       rlond, & ! longitudes
-                       rlatd ! latitudes
-use infotrac, only : nqtot ! number of advected tracers
-use comgeomfi_h, only: ini_fillgeom
-
-implicit none
-
-include "iniprint.h"
-
-real,intent(in) :: prad ! radius of the planet (m)
-real,intent(in) :: pg ! gravitational acceleration (m/s2)
-real,intent(in) :: pr ! ! reduced gas constant R/mu
-real,intent(in) :: pcpp ! specific heat Cp
-real,intent(in) :: punjours ! length (in s) of a standard day
-!integer,intent(in) :: ngrid ! number of horizontal grid points in the physics (full grid)
-integer,intent(in) :: nlayer ! number of atmospheric layers
-integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes
-integer,intent(in) :: jj  ! number of atompsheric columns along latitudes
-real,intent(in) :: rlatu(jj+1) ! latitudes of the dynamics U grid
-real,intent(in) :: rlonv(ii+1) ! longitudes of the dynamics V grid
-real,intent(in) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2)
-real,intent(in) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
-real,intent(in) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
-integer,intent(in) :: pdayref ! reference day of for the simulation
-real,intent(in) :: ptimestep !physics time step (s)
-integer,intent(in) :: iflag_phys ! type of physics to be called
-
-integer :: ibegin,iend,offset
-integer :: i,j
-character(len=20) :: modname='iniphysiq'
-character(len=80) :: abort_message
-real :: total_area_phy, total_area_dyn
-
-
-! global array, on full physics grid:
-real,allocatable :: latfi(:)
-real,allocatable :: lonfi(:)
-real,allocatable :: cufi(:)
-real,allocatable :: cvfi(:)
-real,allocatable :: airefi(:)
-
-IF (nlayer.NE.klev) THEN
-  write(*,*) 'STOP in ',trim(modname)
-  write(*,*) 'Problem with dimensions :'
-  write(*,*) 'nlayer     = ',nlayer
-  write(*,*) 'klev   = ',klev
-  abort_message = ''
-  CALL abort_gcm (modname,abort_message,1)
-ENDIF
-
-!IF (ngrid.NE.klon_glo) THEN
-!  write(*,*) 'STOP in ',trim(modname)
-!  write(*,*) 'Problem with dimensions :'
-!  write(*,*) 'ngrid     = ',ngrid
-!  write(*,*) 'klon   = ',klon_glo
-!  abort_message = ''
-!  CALL abort_gcm (modname,abort_message,1)
-!ENDIF
-
-! Generate global arrays on full physics grid
-allocate(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
-latfi(1)=rlatu(1)
-lonfi(1)=0.
-cufi(1) = cu(1)
-cvfi(1) = cv(1)
-DO j=2,jj
-  DO i=1,ii
-    latfi((j-2)*ii+1+i)= rlatu(j)
-    lonfi((j-2)*ii+1+i)= rlonv(i)
-    cufi((j-2)*ii+1+i) = cu((j-1)*(ii+1)+i)
-    cvfi((j-2)*ii+1+i) = cv((j-1)*(ii+1)+i)
-  ENDDO
-ENDDO
-latfi(klon_glo)= rlatu(jj+1)
-lonfi(klon_glo)= 0.
-cufi(klon_glo) = cu((ii+1)*jj+1)
-cvfi(klon_glo) = cv((ii+1)*jj-ii)
-
-! build airefi(), mesh area on physics grid
-allocate(airefi(klon_glo))
-CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi)
-! Poles are single points on physics grid
-airefi(1)=sum(aire(1:ii,1))
-airefi(klon_glo)=sum(aire(1:ii,jj+1))
-
-! Sanity check: do total planet area match between physics and dynamics?
-total_area_dyn=sum(aire(1:ii,1:jj+1))
-total_area_phy=sum(airefi(1:klon_glo))
-IF (total_area_dyn/=total_area_phy) THEN
-  WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
-  WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
-  WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
-  IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
-    ! stop here if the relative difference is more than 0.001%
-    abort_message = 'planet total surface discrepancy'
-    CALL abort_gcm(modname, abort_message, 1)
-  ENDIF
-ENDIF
-
-
-
-!$OMP PARALLEL 
-! Now generate local lon/lat/cu/cv/area arrays
-call initcomgeomphy
-
-!!!!$OMP PARALLEL PRIVATE(ibegin,iend) 
-!!!$OMP+         SHARED(parea,pcu,pcv,plon,plat)
-      
-offset=klon_mpi_begin-1
-airephy(1:klon_omp)=airefi(offset+klon_omp_begin:offset+klon_omp_end)
-cuphy(1:klon_omp)=cufi(offset+klon_omp_begin:offset+klon_omp_end)
-cvphy(1:klon_omp)=cvfi(offset+klon_omp_begin:offset+klon_omp_end)
-rlond(1:klon_omp)=lonfi(offset+klon_omp_begin:offset+klon_omp_end)
-rlatd(1:klon_omp)=latfi(offset+klon_omp_begin:offset+klon_omp_end)
-
-! copy some fundamental parameters to physics 
-! and do some initializations 
-call phys_state_var_init(klon_omp,nlayer,nqtot, &
-                         punjours,ptimestep,prad,pg,pr,pcpp)
-call ini_fillgeom(klon_omp,rlatd,rlond,airephy)
-call conf_phys(klon_omp,nlayer,nqtot)
-
-!$OMP END PARALLEL
-
-
-end subroutine iniphysiq
Index: /trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/iniphysiq_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/iniphysiq_mod.F90	(revision 1523)
+++ /trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/iniphysiq_mod.F90	(revision 1523)
@@ -0,0 +1,183 @@
+MODULE iniphysiq_mod
+
+CONTAINS
+
+subroutine iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep,           &
+                     rlatu,rlatv,rlonu,rlonv,aire,cu,cv,                 &
+                     prad,pg,pr,pcpp,iflag_phys)
+
+use dimphy, only : klev ! number of atmospheric levels
+use mod_grid_phy_lmdz, only : klon_glo ! number of atmospheric columns
+                                       ! (on full grid)
+use mod_phys_lmdz_para, only : klon_omp, & ! number of columns (on local omp grid)
+                               klon_omp_begin, & ! start index of local omp subgrid
+                               klon_omp_end, & ! end index of local omp subgrid
+                               klon_mpi_begin ! start indes of columns (on local mpi grid)
+
+use comgeomphy, only : initcomgeomphy, &
+                       airephy, & ! physics grid area (m2)
+                       cuphy, & ! cu coeff. (u_covariant = cu * u)
+                       cvphy, & ! cv coeff. (v_covariant = cv * v)
+                       rlond, & ! longitudes
+                       rlatd ! latitudes
+use infotrac, only : nqtot ! number of advected tracers
+use comgeomfi_h, only: ini_fillgeom
+use regular_lonlat_mod, only: init_regular_lonlat, &
+                              east, west, north, south, &
+                              north_east, north_west, &
+                              south_west, south_east
+
+implicit none
+
+include "iniprint.h"
+
+real,intent(in) :: prad ! radius of the planet (m)
+real,intent(in) :: pg ! gravitational acceleration (m/s2)
+real,intent(in) :: pr ! ! reduced gas constant R/mu
+real,intent(in) :: pcpp ! specific heat Cp
+real,intent(in) :: punjours ! length (in s) of a standard day
+!integer,intent(in) :: ngrid ! number of horizontal grid points in the physics (full grid)
+integer,intent(in) :: nlayer ! number of atmospheric layers
+integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes
+integer,intent(in) :: jj  ! number of atompsheric columns along latitudes
+real,intent(in) :: rlatu(jj+1) ! latitudes of the physics grid
+real,intent(in) :: rlatv(jj) ! latitude boundaries of the physics grid
+real,intent(in) :: rlonv(ii+1) ! longitudes of the physics grid
+real,intent(in) :: rlonu(ii+1) ! longitude boundaries of the physics grid
+real,intent(in) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2)
+real,intent(in) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
+real,intent(in) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
+integer,intent(in) :: pdayref ! reference day of for the simulation
+real,intent(in) :: ptimestep !physics time step (s)
+integer,intent(in) :: iflag_phys ! type of physics to be called
+
+integer :: ibegin,iend,offset
+integer :: i,j
+character(len=20) :: modname='iniphysiq'
+character(len=80) :: abort_message
+real :: total_area_phy, total_area_dyn
+real :: pi
+
+! boundaries, on global grid
+real,allocatable :: boundslon_reg(:,:)
+real,allocatable :: boundslat_reg(:,:)
+
+! global array, on full physics grid:
+real,allocatable :: latfi(:)
+real,allocatable :: lonfi(:)
+real,allocatable :: cufi(:)
+real,allocatable :: cvfi(:)
+real,allocatable :: airefi(:)
+
+pi=2.*asin(1.0)
+
+IF (nlayer.NE.klev) THEN
+  write(*,*) 'STOP in ',trim(modname)
+  write(*,*) 'Problem with dimensions :'
+  write(*,*) 'nlayer     = ',nlayer
+  write(*,*) 'klev   = ',klev
+  abort_message = ''
+  CALL abort_gcm (modname,abort_message,1)
+ENDIF
+
+!IF (ngrid.NE.klon_glo) THEN
+!  write(*,*) 'STOP in ',trim(modname)
+!  write(*,*) 'Problem with dimensions :'
+!  write(*,*) 'ngrid     = ',ngrid
+!  write(*,*) 'klon   = ',klon_glo
+!  abort_message = ''
+!  CALL abort_gcm (modname,abort_message,1)
+!ENDIF
+
+! init regular global longitude-latitude grid points and boundaries
+ALLOCATE(boundslon_reg(ii,2))
+ALLOCATE(boundslat_reg(jj+1,2))
+  
+DO i=1,ii
+   boundslon_reg(i,east)=rlonu(i) 
+   boundslon_reg(i,west)=rlonu(i+1) 
+ENDDO
+
+boundslat_reg(1,north)= PI/2 
+boundslat_reg(1,south)= rlatv(1)
+DO j=2,jj
+   boundslat_reg(j,north)=rlatv(j-1) 
+   boundslat_reg(j,south)=rlatv(j) 
+ENDDO
+boundslat_reg(jj+1,north)= rlatv(jj) 
+boundslat_reg(jj+1,south)= -PI/2
+
+! Write values in module regular_lonlat_mod
+CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, &
+                         boundslon_reg, boundslat_reg)
+
+! Generate global arrays on full physics grid
+allocate(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
+latfi(1)=rlatu(1)
+lonfi(1)=0.
+cufi(1) = cu(1)
+cvfi(1) = cv(1)
+DO j=2,jj
+  DO i=1,ii
+    latfi((j-2)*ii+1+i)= rlatu(j)
+    lonfi((j-2)*ii+1+i)= rlonv(i)
+    cufi((j-2)*ii+1+i) = cu((j-1)*(ii+1)+i)
+    cvfi((j-2)*ii+1+i) = cv((j-1)*(ii+1)+i)
+  ENDDO
+ENDDO
+latfi(klon_glo)= rlatu(jj+1)
+lonfi(klon_glo)= 0.
+cufi(klon_glo) = cu((ii+1)*jj+1)
+cvfi(klon_glo) = cv((ii+1)*jj-ii)
+
+! build airefi(), mesh area on physics grid
+allocate(airefi(klon_glo))
+CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi)
+! Poles are single points on physics grid
+airefi(1)=sum(aire(1:ii,1))
+airefi(klon_glo)=sum(aire(1:ii,jj+1))
+
+! Sanity check: do total planet area match between physics and dynamics?
+total_area_dyn=sum(aire(1:ii,1:jj+1))
+total_area_phy=sum(airefi(1:klon_glo))
+IF (total_area_dyn/=total_area_phy) THEN
+  WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
+  WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
+  WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
+  IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
+    ! stop here if the relative difference is more than 0.001%
+    abort_message = 'planet total surface discrepancy'
+    CALL abort_gcm(modname, abort_message, 1)
+  ENDIF
+ENDIF
+
+
+
+!$OMP PARALLEL 
+! Now generate local lon/lat/cu/cv/area arrays
+call initcomgeomphy
+
+!!!!$OMP PARALLEL PRIVATE(ibegin,iend) 
+!!!$OMP+         SHARED(parea,pcu,pcv,plon,plat)
+      
+offset=klon_mpi_begin-1
+airephy(1:klon_omp)=airefi(offset+klon_omp_begin:offset+klon_omp_end)
+cuphy(1:klon_omp)=cufi(offset+klon_omp_begin:offset+klon_omp_end)
+cvphy(1:klon_omp)=cvfi(offset+klon_omp_begin:offset+klon_omp_end)
+rlond(1:klon_omp)=lonfi(offset+klon_omp_begin:offset+klon_omp_end)
+rlatd(1:klon_omp)=latfi(offset+klon_omp_begin:offset+klon_omp_end)
+
+! copy some fundamental parameters to physics 
+! and do some initializations 
+call phys_state_var_init(klon_omp,nlayer,nqtot, &
+                         punjours,ptimestep,prad,pg,pr,pcpp)
+call ini_fillgeom(klon_omp,rlatd,rlond,airephy)
+call conf_phys(klon_omp,nlayer,nqtot)
+
+!$OMP END PARALLEL
+
+
+end subroutine iniphysiq
+
+
+END MODULE iniphysiq_mod
Index: /trunk/LMDZ.MARS/libf/phy_common/regular_lonlat_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phy_common/regular_lonlat_mod.F90	(revision 1523)
+++ /trunk/LMDZ.MARS/libf/phy_common/regular_lonlat_mod.F90	(revision 1523)
@@ -0,0 +1,55 @@
+MODULE regular_lonlat_mod
+
+  ! Store information on the global physics grid
+  ! for a regular (longitude-latitude) grid
+  
+  INTEGER, PARAMETER :: north_east=1     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: north_west=2     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south_west=3     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south_east=4     ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: east=1           ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: west=2           ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: north=1          ! boundaries of regular lontlat
+  INTEGER, PARAMETER :: south=2          ! boundaries of regular lontlat
+
+! Global definition, shared by all threads
+! Do not set threadprivate directives
+
+  REAL,SAVE,ALLOCATABLE :: lon_reg(:)      ! value of longitude cell (rad)
+
+  REAL,SAVE,ALLOCATABLE :: lat_reg(:)      ! value of longitude cell (rad)
+
+  REAL,SAVE,ALLOCATABLE :: boundslon_reg(:,:)      ! value of boundaries cell (1=>east, 2=>west)(rad)
+
+  REAL,SAVE,ALLOCATABLE :: boundslat_reg(:,:)      ! value of longitude cell (1=>north, 2=>south)(rad)
+
+
+CONTAINS
+  
+  SUBROUTINE init_regular_lonlat(nbp_lon, nbp_lat, lon_reg_, lat_reg_, boundslon_reg_, boundslat_reg_)
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: nbp_lon
+    INTEGER,INTENT(IN) :: nbp_lat
+    REAL,   INTENT(IN) :: lon_reg_(nbp_lon)  
+    REAL,   INTENT(IN) :: lat_reg_(nbp_lat) 
+    REAL,   INTENT(IN) :: boundslon_reg_(nbp_lon,2) 
+    REAL,   INTENT(IN) :: boundslat_reg_(nbp_lat,2)  
+
+    
+    ALLOCATE(lon_reg(nbp_lon))
+    lon_reg(:)=lon_reg_(:)
+
+    ALLOCATE(lat_reg(nbp_lat)) 
+    lat_reg(:)=lat_reg_(:)
+
+    ALLOCATE(boundslon_reg(nbp_lon,2)) 
+    boundslon_reg(:,:)=boundslon_reg_(:,:)
+
+    ALLOCATE(boundslat_reg(nbp_lat,2))  
+    boundslat_reg(:,:)=boundslat_reg_(:,:)
+
+  
+  END SUBROUTINE init_regular_lonlat    
+
+END MODULE regular_lonlat_mod
+
Index: unk/LMDZ.TITAN/libf/phytitan/dyn1d/iniphysiq.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/iniphysiq.F90	(revision 1522)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dynlonlat_phylonlat/phytitan/iniphysiq.F90
Index: /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/iniphysiq_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/iniphysiq_mod.F90	(revision 1523)
+++ /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/iniphysiq_mod.F90	(revision 1523)
@@ -0,0 +1,1 @@
+link ../../dynlonlat_phylonlat/phytitan/iniphysiq_mod.F90
Index: /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F	(revision 1522)
+++ /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F	(revision 1523)
@@ -11,4 +11,5 @@
      .                       aps,bps,scaleheight,pseudoalt,
      .                       disvert_type,pressure_exner
+      USE iniphysiq_mod, ONLY: iniphysiq
       IMPLICIT NONE
 
@@ -191,8 +192,4 @@
 
       CALL init_phys_lmdz(iim,jjm,llm,1,(/1/))
-      CALL iniphysiq(1,1,llm,daysec,day0,dtphys,
-     .            lati,long,area,cufi,cvfi,rad,g,r,cpp,1)
-
-      call ini_cpdet
 
 c   la surface de chaque maille est inutile en 1D --->
@@ -201,4 +198,17 @@
       cufi(1)=1.E+0
       cvfi(1)=1.E+0
+
+c Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid,
+c e.g. for cell boundaries, which are meaningless in 1D; so pad these 
+c with '0.' when necessary
+      CALL iniphysiq(1,1,llm,daysec,day0,dtphys,
+     &            (/lati(1),0./),(/0./),
+     &            (/0.,0./),(/long(1),0./),
+     &            (/ (/area,0./),(/0.,0./) /),
+     &            (/cufi,0.,0.,0./),
+     &            (/cvfi,0./),
+     &            rad,g,r,cpp,1)
+
+      call ini_cpdet
 
 c   le geopotentiel au sol est inutile en 1D car tout est controle
Index: unk/LMDZ.TITAN/libf/phytitan/dyn1d/write_field.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/write_field.F90	(revision 1522)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d_common/write_field.F90
Index: unk/LMDZ.TITAN/libf/phytitan/write_field_phy.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/write_field_phy.F90	(revision 1522)
+++ 	(revision )
@@ -1,36 +1,0 @@
-!
-! $Header$
-!
-MODULE write_field_phy
-
-  CONTAINS 
- 
-    SUBROUTINE WriteField_phy(name,Field,ll)
-    USE dimphy
-    USE mod_phys_lmdz_para
-    USE mod_grid_phy_lmdz
-    USE Write_Field
-    
-    IMPLICIT NONE
-    include 'dimensions.h'
-    include 'paramet.h'
-
-    character(len=*)   :: name
-    INTEGER :: ll
-    real, dimension(klon_omp,ll) :: Field
-    real,save,allocatable :: Field_tmp(:,:)
-    real, dimension(klon_glo,ll):: New_Field
-    real, dimension(iim,jjp1,ll):: Field_2d
-
-    CALL Gather(Field,New_Field)
-!$OMP MASTER
-    IF (is_mpi_root) THEN	
-      CALL Grid1Dto2D_glo(New_Field,Field_2D)
-      CALL WriteField(name,Field_2d)
-    ENDIF
-!$OMP END MASTER
-
-  
-   END SUBROUTINE WriteField_phy
- 
- END MODULE write_field_phy
Index: unk/LMDZ.VENUS/libf/phyvenus/dyn1d/iniphysiq.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/iniphysiq.F90	(revision 1522)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dynlonlat_phylonlat/phyvenus/iniphysiq.F90
Index: /trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/iniphysiq_mod.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/iniphysiq_mod.F90	(revision 1523)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/iniphysiq_mod.F90	(revision 1523)
@@ -0,0 +1,1 @@
+link ../../dynlonlat_phylonlat/phyvenus/iniphysiq_mod.F90
Index: /trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/rcm1d.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/rcm1d.F	(revision 1522)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/rcm1d.F	(revision 1523)
@@ -13,4 +13,5 @@
      .                       disvert_type,pressure_exner
       use conc, only: rho
+      USE iniphysiq_mod, ONLY: iniphysiq
 
       IMPLICIT NONE
@@ -199,8 +200,4 @@
 
       CALL init_phys_lmdz(iim,jjm,llm,1,(/1/))
-      CALL iniphysiq(1,1,llm,daysec,day0,dtphys,
-     .            lati,long,area,cufi,cvfi,rad,g,r,cpp,1)
-
-      call ini_cpdet
 
 c   la surface de chaque maille est inutile en 1D --->
@@ -209,4 +206,17 @@
       cufi(1)=1.E+0
       cvfi(1)=1.E+0
+
+c Ehouarn: iniphysiq requires arrays related to (3D) dynamics grid,
+c e.g. for cell boundaries, which are meaningless in 1D; so pad these 
+c with '0.' when necessary
+      CALL iniphysiq(1,1,llm,daysec,day0,dtphys,
+     &            (/lati(1),0./),(/0./),
+     &            (/0.,0./),(/long(1),0./),
+     &            (/ (/area,0./),(/0.,0./) /),
+     &            (/cufi,0.,0.,0./),
+     &            (/cvfi,0./),
+     &            rad,g,r,cpp,1)
+
+      call ini_cpdet
 
 c   le geopotentiel au sol est inutile en 1D car tout est controle
Index: unk/LMDZ.VENUS/libf/phyvenus/dyn1d/write_field.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/write_field.F90	(revision 1522)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d_common/write_field.F90
Index: unk/LMDZ.VENUS/libf/phyvenus/write_field_phy.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/write_field_phy.F90	(revision 1522)
+++ 	(revision )
@@ -1,36 +1,0 @@
-!
-! $Header$
-!
-MODULE write_field_phy
-
-  CONTAINS 
- 
-    SUBROUTINE WriteField_phy(name,Field,ll)
-    USE dimphy
-    USE mod_phys_lmdz_para
-    USE mod_grid_phy_lmdz
-    USE Write_Field
-    
-    IMPLICIT NONE
-    include 'dimensions.h'
-    include 'paramet.h'
-
-    character(len=*)   :: name
-    INTEGER :: ll
-    real, dimension(klon_omp,ll) :: Field
-    real,save,allocatable :: Field_tmp(:,:)
-    real, dimension(klon_glo,ll):: New_Field
-    real, dimension(iim,jjp1,ll):: Field_2d
-
-    CALL Gather(Field,New_Field)
-!$OMP MASTER
-    IF (is_mpi_root) THEN	
-      CALL Grid1Dto2D_glo(New_Field,Field_2D)
-      CALL WriteField(name,Field_2d)
-    ENDIF
-!$OMP END MASTER
-
-  
-   END SUBROUTINE WriteField_phy
- 
- END MODULE write_field_phy
