Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/etat0_netcdf.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/etat0_netcdf.F	(revision 1292)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/etat0_netcdf.F	(revision 1293)
@@ -263,6 +263,6 @@
         varname = 'masque'
         masque(:,:) = 0.0
-        CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
-     ,  jjm ,rlonu,rlatv , interbar )
+        CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, masque,
+     $       0.0, jjm ,rlonu,rlatv , interbar )
         WRITE(*,*) 'MASQUE construit : Masque'
         WRITE(*,'(97I1)') nINT(masque(:,:))
@@ -324,6 +324,6 @@
       ! This line needs to be replaced by a call to restget to get the values in the restart file
       orog(:,:) = 0.0
-       CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
-     , jjm ,rlonu,rlatv , interbar, masque )
+       CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, orog,
+     $     0.0 , jjm ,rlonu,rlatv , interbar, masque )
       !
       WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
@@ -333,6 +333,6 @@
       ! This line needs to be replaced by a call to restget to get the values in the restart file
       rugo(:,:) = 0.0
-       CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
-     , jjm, rlonu,rlatv , interbar )
+       CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, rugo,
+     $     0.0 , jjm, rlonu,rlatv , interbar )
       !
       WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite' 
@@ -346,6 +346,6 @@
       varname = 'psol'
       psol(:,:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
-     , jjm ,rlonu,rlatv , interbar )
+      CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, psol,
+     $     0.0 , jjm ,rlonu,rlatv , interbar )
       !
       !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM 
@@ -375,21 +375,21 @@
       varname = 'surfgeo'
       phis(:,:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
-     , jjm ,rlonu,rlatv, interbar )
+      CALL startget_phys2d(varname, iip1, jjp1, rlonv, rlatu, phis,
+     $     0.0 , jjm ,rlonu,rlatv, interbar )
       !
       varname = 'u'
       uvent(:,:,:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
-     . workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
+      CALL startget_dyn(varname, rlonu, rlatu, pls, workvar, uvent, 0.,
+     $     rlonv, rlatv, interbar )
       !  
       varname = 'v'
       vvent(:,:,:) = 0.0
-      CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
-     . workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
+      CALL startget_dyn(varname, rlonv, rlatv, pls(:, :jjm, :),
+     . workvar(:, :jjm, :), vvent, 0., rlonu, rlatu(:jjm), interbar )
       !
       varname = 't'
       t3d(:,:,:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
-     . workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
+      CALL startget_dyn(varname, rlonv, rlatu, pls, workvar, t3d, 0.,
+     $     rlonu, rlatv , interbar )
       !
       WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
@@ -397,6 +397,6 @@
       varname = 'tpot'
       tpot(:,:,:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
-     . pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
+      CALL startget_dyn(varname, rlonv, rlatu, pls, pk, tpot, 0., rlonu,
+     $     rlatv, interbar)
       !
       WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
@@ -420,6 +420,6 @@
       WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
      .                           maxval(qsat(:,:,:))
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
-     . qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
+      CALL startget_dyn(varname, rlonv, rlatu, pls, qsat, qd, 0., rlonu,
+     $     rlatv , interbar )
       q3d(:,:,:,1) = qd(:,:,:)
       !
@@ -431,6 +431,6 @@
       ! This line needs to be replaced by a call to restget to get the values in the restart file
       tsol(:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
-     .    jjm, rlonu, rlatv , interbar )
+      CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon,
+     $     tsol, 0.0, jjm, rlonu, rlatv , interbar )
       !
       WRITE(*,*) 'TSOL construit :'
@@ -439,55 +439,55 @@
       varname = 'qsol'
       qsol(:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0,
-     .   jjm, rlonu, rlatv , interbar )
+      CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon,
+     $     qsol, 0.0, jjm, rlonu, rlatv , interbar )
       !
       varname = 'snow'
       sn(:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0,
-     .    jjm, rlonu, rlatv , interbar )
+      CALL startget_phys1d(varname, iip1, jjp1, rlonv, rlatu, klon, sn,
+     $     0.0, jjm, rlonu, rlatv , interbar )
       !
       varname = 'rads'
       radsol(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
-     .    jjm, rlonu, rlatv , interbar )
+      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,
+     $     0.0, jjm, rlonu, rlatv , interbar )
       !
       varname = 'rugmer'
       rugmer(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
-     .     jjm, rlonu, rlatv , interbar )
+      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,
+     $     0.0, jjm, rlonu, rlatv , interbar )
       !
 !      varname = 'agesno'
 !      agesno(:) = 0.0
-!      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
+!      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
 !     .     jjm, rlonu, rlatv , interbar )
 
       varname = 'zmea'
       zmea(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
+      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
      .     jjm, rlonu, rlatv , interbar )
 
       varname = 'zstd'
       zstd(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
+      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
      .     jjm, rlonu, rlatv , interbar )
       varname = 'zsig'
       zsig(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
+      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
      .     jjm, rlonu, rlatv , interbar )
       varname = 'zgam'
       zgam(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
+      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
      .     jjm, rlonu, rlatv , interbar )
       varname = 'zthe'
       zthe(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
+      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
      .     jjm, rlonu, rlatv , interbar )
       varname = 'zpic'
       zpic(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
+      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
      .     jjm, rlonu, rlatv , interbar )
       varname = 'zval'
       zval(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
+      CALL startget_phys1d(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
      .     jjm, rlonu, rlatv , interbar )
 c
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/inter_barx.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/inter_barx.F	(revision 1292)
+++ 	(revision )
@@ -1,218 +1,0 @@
-!
-! $Header$
-!
-       SUBROUTINE inter_barx ( idatmax,xidat,fdat,imodmax,ximod,fmod ) 
-
-c      .... Auteurs :  Robert Sadourny ,  P. Le Van  .....
-c
-       IMPLICIT NONE
-c    ----------------------------------------------------------
-c        INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES
-c            VERSION UNIDIMENSIONNELLE  ,   EN  LONGITUDE .
-c    ----------------------------------------------------------
-c
-c     idat : indice du champ de donnees, de 1 a idatmax
-c     imod : indice du champ du modele,  de 1 a  imodmax
-c     fdat(idat) : champ de donnees (entrees)
-c     fmod(imod) : champ du modele (sorties)
-c     xidat(idat): abscisses des interfaces des mailles donnees
-c     ximod(imod): abscisses des interfaces des mailles modele
-c      ( L'indice 1 correspond a l'interface mailLE 1 / maille 2)
-c      ( Les abscisses sont exprimes en degres)
-
-
-      INTEGER idatmax, imodmax
-      REAL xidat(idatmax),fdat(idatmax),ximod(imodmax),fmod(imodmax)
-
-c    ...  Variables locales ... 
-    
-      REAL xxid(idatmax+1), xxd(idatmax+1), fdd(idatmax+1)
-      REAL  fxd(idatmax+1), xchan(idatmax+1), fdchan(idatmax+1) 
-      REAL  xxim(imodmax)
-
-      REAL x0,xim0,dx,dxm
-      REAL chmin,chmax,pi
-
-      INTEGER imod,idat,i,ichang,id0,id1,nid,idatmax1
-      
-      pi = 2. * ASIN(1.)
-
-c  -----------------------------------------------------
-c   REDEFINITION DE L'ORIGINE DES ABSCISSES
-c    A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE  
-c  -----------------------------------------------------
-      DO imod = 1,imodmax
-       xxim(imod) = ximod(imod)
-      ENDDO
-
-      CALL minmax( imodmax,xxim,chmin,chmax)
-       IF( chmax.LT.6.50 )   THEN
-c        PRINT  3
-c        PRINT *,'   conversion des longit. ximod (donnees en radians)'
-c     , ,' en degres  .' 
-c        PRINT  3
-        DO imod = 1, imodmax
-        xxim(imod) = xxim(imod) * 180./pi
-        ENDDO  
-       ENDIF
-
-      xim0 = xxim(imodmax) - 360.
-
-      DO imod = 1, imodmax
-       xxim(imod) = xxim(imod) - xim0
-      ENDDO
-
-      idatmax1 = idatmax +1
-
-      DO idat = 1, idatmax
-       xxd(idat) = xidat(idat)
-      ENDDO
-
-      CALL minmax( idatmax,xxd,chmin,chmax)
-       IF( chmax.LT.6.50 )  THEN
-c        PRINT  3
-c        PRINT *,'   conversion des longit. ximod (donnees en radians)'
-c     , ,' en degres  .' 
-c        PRINT  3
-        DO idat = 1, idatmax
-        xxd(idat) = xxd(idat) * 180./pi
-        ENDDO  
-       ENDIF
-
-      DO idat = 1, idatmax
-       xxd(idat) = MOD( xxd(idat) - xim0, 360. )
-       fdd(idat) = fdat (idat)
-      ENDDO
-c       PRINT *,' xxd redef. origine abscisses '
-c       PRINT 2,(xxd(i),i=1,idatmax)
-
-      DO i = 2, idatmax
-        IF( ( xxd(i) - xxd(i-1)).LT.0. )  THEN
-         ichang = i
-         GO TO 5
-        ENDIF
-      ENDDO
-      GO TO 6
-c
-c  ***  reorganisation  des longitudes entre 0. et 360. degres ****
-c
- 5    nid = idatmax - ichang +1
-      DO i = 1, nid
-        xchan (i) = xxd(i+ichang -1 )
-        fdchan(i) = fdd(i+ichang -1 )
-      ENDDO
-       DO i=1,ichang -1
-        xchan (i+ nid) = xxd(i)
-        fdchan(i+nid) = fdd(i) 
-       ENDDO
-      DO i =1,idatmax
-       xxd(i) = xchan(i)
-       fdd(i) = fdchan(i)
-      ENDDO
-
- 6    continue
-
-
-c  ------------------------------------------------
-c    translation des champs de donnees par rapport
-c    a la nouvelle origine, avec redondance de la
-c       maille a cheval sur les bords
-c -----------------------------------------------
-
-      id0 = 0
-      id1 = 0
-
-      DO idat = 1, idatmax
-       IF ( xxd( idatmax1- idat ).LT.360.)   GO TO 10
-       id1 = id1 + 1
-      ENDDO
-
- 10   DO idat = 1, idatmax
-       IF (xxd(idat).GT.0.) GO TO 20
-       id0 = id0 + 1
-      END DO
-
- 20   IF( id1.EQ.0 ) GO TO 30
-      DO idat = 1, id1
-       xxid(idat) = xxd(idatmax - id1 + idat) - 360.
-       fxd (idat) = fdd(idatmax - id1 + idat)     
-      END DO
-      DO idat = 1, idatmax - id1
-       xxid(idat + id1) = xxd(idat)
-       fxd (idat + id1) = fdd(idat)
-      END DO
-
-  30  IF(id0.EQ.0) GO TO 40
-      DO idat = 1, idatmax - id0
-       xxid(idat) = xxd(idat + id0)
-       fxd (idat) = fdd(idat + id0)
-      END DO
-
-      DO idat = 1, id0
-       xxid (idatmax - id0 + idat) =  xxd(idat) + 360.
-       fxd  (idatmax - id0 + idat) =  fdd(idat)   
-      END DO
-      GO TO 50
- 
- 40   DO idat = 1, idatmax
-       xxid(idat)  = xxd(idat)
-       fxd (idat)  = fdd(idat)
-      ENDDO
-
- 50   xxid(idatmax1) = xxid(1) + 360.
-      fxd (idatmax1) = fxd(1)
-
-c  ------------------------------------
-c   initialisation du champ du modele
-
-      DO imod = 1, imodmax
-       fmod(imod) = 0.
-      END DO
- 
-c      PRINT *,' id0 id1 ',id0,id1
-c      PRINT *,' xxim apres translation  '
-c      PRINT 2,(xxim(i),i=1,imodmax)
-c      PRINT *,' xxid apres translation '
-c      PRINT 2,(xxid(i),i=1,idatmax)
-c ---------------------------------------
-c iteration
-
-      x0   = xim0
-      dxm  = 0.
-      imod = 1
-      idat = 1
- 
- 100  IF (xxim(imod).LT.xxid(idat)) THEN
-       dx   = xxim(imod) - x0
-       dxm  = dxm + dx
-       fmod(imod) = (fmod(imod) + dx * fxd(idat)) / dxm
-       x0   = xxim(imod)
-       dxm  = 0.
-       imod = imod + 1
-       IF (imod.LE.imodmax) GO TO 100
-  
-      ELSE IF (xxim(imod).GT.xxid(idat)) THEN
-       dx   = xxid(idat) - x0
-       dxm  = dxm + dx
-       fmod(imod) = fmod(imod) + dx * fxd(idat)
-       x0   = xxid(idat)
-       idat = idat + 1
-       GO TO 100
- 
-      ELSE
-       dx   = xxim(imod) - x0
-       dxm  = dxm + dx
-       fmod(imod) = (fmod(imod) + dx * fxd(idat)) / dxm
-       x0   = xxim(imod)
-       dxm  = 0.
-       imod = imod + 1
-       idat = idat + 1
-       IF (imod.LE.imodmax) GO TO 100
-      END IF
-     
-
-3      FORMAT(1x,70("-"))
-2      FORMAT(1x,8f8.2)
-
-       RETURN
-       END
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/inter_barxy.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/inter_barxy.F	(revision 1292)
+++ 	(revision )
@@ -1,59 +1,0 @@
-!
-! $Header$
-!
-       SUBROUTINE inter_barxy ( interfd,jnterfd,dlonid,dlatid ,
-     ,        champ,imod,jmod,rlonimod,rlatimod, jsort,champint )
-
-c    Auteur :   P. Le Van
-c
-       INTEGER interfd,jnterfd,imod,jmod
-       REAL champ(interfd,jnterfd +1 ),dlonid(interfd),dlatid(jnterfd),
-     ,      champint(imod,jsort)
-       REAL rlonimod(imod),rlatimod(jmod)
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom2.h"
-
-       REAL champx(imod),champy(jnterfd +1,imod),chpn(imod),chps(imod)
-       REAL chhpn,chhps
-       REAL fmody(jjp1)
-c
-
-       DO j = 1, jnterfd + 1
-        CALL inter_barx( interfd, dlonid, champ( 1,j ),
-     ,                       imod, rlonimod , champx )
-         DO i = 1,imod
-           champy(j,i) = champx(i)
-         ENDDO
-       ENDDO
-
-       DO i = 1, imod
-        CALL inter_bary( jjm,jnterfd,dlatid,champy(1,i),
-     ,                     jmod ,rlatimod,  fmody     )
-          DO j = 1, jsort
-           champint(i,j) = fmody(j)
-          ENDDO
-       ENDDO
-
-       IF( jsort.EQ.jjp1)  THEN
-
-c   ....  Valeurs uniques  aux  poles ....
-c
-         DO i =  1,imod
-          chpn(i)  = aire( i,  1   ) * champint( i, 1   )
-          chps(i)  = aire( i, jjp1 ) * champint( i,jjp1 )
-         ENDDO
-          chhpn  = SSUM(imod,chpn,1)/apoln
-          chhps  = SSUM(imod,chps,1)/apols
-
-         DO i = 1, imod
-          champint( i,  1  ) = chhpn
-          champint( i, jjp1) = chhps
-         ENDDO
-c
-       ENDIF
-
-       RETURN
-       END
-
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/inter_barxy_m.F90
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/inter_barxy_m.F90	(revision 1293)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/inter_barxy_m.F90	(revision 1293)
@@ -0,0 +1,450 @@
+module inter_barxy_m
+
+  ! Authors: Robert SADOURNY, Phu LE VAN, Lionel GUEZ
+
+  implicit none
+
+  private
+  public inter_barxy
+
+contains
+
+  SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint)
+
+    use assert_eq_m, only: assert_eq
+    use assert_m, only: assert
+
+    include "dimensions.h"
+    ! (for "iim", "jjm")
+
+    include "paramet.h"
+    ! (for other included files)
+
+    include "comgeom2.h"
+    ! (for "aire", "apoln", "apols")
+
+    REAL, intent(in):: dlonid(:)
+    ! (longitude from input file, in rad, from -pi to pi)
+
+    REAL, intent(in):: dlatid(:), champ(:, :), rlonimod(:)
+
+    REAL, intent(in):: rlatimod(:)
+    ! (latitude angle, in degrees or rad, in strictly decreasing order)
+
+    real, intent(out):: champint(:, :)
+    ! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les
+    ! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U)
+    ! Si taille de la seconde dim = jjm, on veut interpoler sur les
+    ! jjm latitudes rlatv du modèle (latitudes de V) 
+
+    ! Variables local to the procedure:
+
+    REAL champy(iim, size(champ, 2))
+    integer j, i, jnterfd, jmods
+
+    REAL yjmod(size(champint, 2))
+    ! (angle, in degrees, in strictly increasing order)
+
+    REAL   yjdat(size(dlatid) + 1) ! angle, in degrees, in increasing order
+    LOGICAL decrois ! "dlatid" is in decreasing order
+
+    !-----------------------------------
+
+    jnterfd = assert_eq(size(champ, 2) - 1, size(dlatid), &
+         "inter_barxy jnterfd")
+    jmods = size(champint, 2)
+    call assert(size(champ, 1) == size(dlonid), "inter_barxy size(champ, 1)")
+    call assert((/size(rlonimod), size(champint, 1)/) == iim, &
+         "inter_barxy iim")
+    call assert(any(jmods == (/jjm, jjm + 1/)), 'inter_barxy jmods')
+    call assert(size(rlatimod) == jjm, "inter_barxy size(rlatimod)")
+
+    ! Check decreasing order for "rlatimod":
+    DO i = 2, jjm
+       IF (rlatimod(i) >= rlatimod(i-1)) stop &
+            '"inter_barxy": "rlatimod" should be strictly decreasing'
+    ENDDO
+
+    yjmod(:jjm) = ord_coordm(rlatimod)
+    IF (jmods == jjm + 1) THEN
+       IF (90. - yjmod(jjm) < 0.01) stop &
+            '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.'
+    ELSE
+       ! jmods = jjm
+       IF (ABS(yjmod(jjm) - 90.) > 0.01) stop &
+            '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.'
+    ENDIF
+
+    if (jmods == jjm + 1) yjmod(jjm + 1) = 90.
+
+    DO j = 1, jnterfd + 1
+       champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod)
+    ENDDO
+
+    CALL ord_coord(dlatid, yjdat, decrois) 
+    IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1)
+    DO i = 1, iim
+       champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod)
+    ENDDO
+    champint(:, :) = champint(:, jmods:1:-1)
+
+    IF (jmods == jjm + 1) THEN
+       ! Valeurs uniques aux poles
+       champint(:, 1) = SUM(aire(:iim,  1) * champint(:, 1)) / apoln
+       champint(:, jjm + 1) = SUM(aire(:iim, jjm + 1) &
+            * champint(:, jjm + 1)) / apols
+    ENDIF
+
+  END SUBROUTINE inter_barxy
+
+  !******************************
+
+  function inter_barx(dlonid, fdat, rlonimod) 
+
+    !        INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES
+    !            VERSION UNIDIMENSIONNELLE  ,   EN  LONGITUDE .
+
+    !     idat : indice du champ de donnees, de 1 a idatmax
+    !     imod : indice du champ du modele,  de 1 a  imodmax
+    !     fdat(idat) : champ de donnees (entrees)
+    !     inter_barx(imod) : champ du modele (sorties)
+    !     dlonid(idat): abscisses des interfaces des mailles donnees
+    !     rlonimod(imod): abscisses des interfaces des mailles modele
+    !      ( L'indice 1 correspond a l'interface mailLE 1 / maille 2)
+    !      ( Les abscisses sont exprimées en degres)
+
+    use assert_eq_m, only: assert_eq
+
+    IMPLICIT NONE
+
+    REAL, intent(in):: dlonid(:)
+    real, intent(in):: fdat(:)
+    real, intent(in):: rlonimod(:)
+
+    real inter_barx(size(rlonimod))
+
+    !    ...  Variables locales ... 
+
+    INTEGER idatmax, imodmax
+    REAL xxid(size(dlonid)+1), xxd(size(dlonid)+1), fdd(size(dlonid)+1)
+    REAL  fxd(size(dlonid)+1), xchan(size(dlonid)+1), fdchan(size(dlonid)+1) 
+    REAL  xxim(size(rlonimod))
+
+    REAL x0, xim0, dx, dxm
+    REAL chmin, chmax, pi
+
+    INTEGER imod, idat, i, ichang, id0, id1, nid, idatmax1
+
+    !-----------------------------------------------------
+
+    idatmax = assert_eq(size(dlonid), size(fdat), "inter_barx idatmax")
+    imodmax = size(rlonimod)
+
+    pi = 2. * ASIN(1.)
+
+    !   REDEFINITION DE L'ORIGINE DES ABSCISSES
+    !    A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE  
+    DO imod = 1, imodmax
+       xxim(imod) = rlonimod(imod)
+    ENDDO
+
+    CALL minmax( imodmax, xxim, chmin, chmax)
+    IF( chmax.LT.6.50 )   THEN
+       DO imod = 1, imodmax
+          xxim(imod) = xxim(imod) * 180./pi
+       ENDDO
+    ENDIF
+
+    xim0 = xxim(imodmax) - 360.
+
+    DO imod = 1, imodmax
+       xxim(imod) = xxim(imod) - xim0
+    ENDDO
+
+    idatmax1 = idatmax +1
+
+    DO idat = 1, idatmax
+       xxd(idat) = dlonid(idat)
+    ENDDO
+
+    CALL minmax( idatmax, xxd, chmin, chmax)
+    IF( chmax.LT.6.50 )  THEN
+       DO idat = 1, idatmax
+          xxd(idat) = xxd(idat) * 180./pi
+       ENDDO
+    ENDIF
+
+    DO idat = 1, idatmax
+       xxd(idat) = AMOD( xxd(idat) - xim0, 360. )
+       fdd(idat) = fdat (idat)
+    ENDDO
+
+    i = 2
+    DO while (xxd(i) >= xxd(i-1) .and. i < idatmax)
+       i = i + 1
+    ENDDO
+    IF (xxd(i) < xxd(i-1)) THEN
+       ichang = i
+       !  ***  reorganisation  des longitudes entre 0. et 360. degres ****
+       nid = idatmax - ichang +1
+       DO i = 1, nid
+          xchan (i) = xxd(i+ichang -1 )
+          fdchan(i) = fdd(i+ichang -1 )
+       ENDDO
+       DO i=1, ichang -1
+          xchan (i+ nid) = xxd(i)
+          fdchan(i+nid) = fdd(i) 
+       ENDDO
+       DO i =1, idatmax
+          xxd(i) = xchan(i)
+          fdd(i) = fdchan(i)
+       ENDDO
+    end IF
+
+    !    translation des champs de donnees par rapport
+    !    a la nouvelle origine, avec redondance de la
+    !       maille a cheval sur les bords
+
+    id0 = 0
+    id1 = 0
+
+    DO idat = 1, idatmax
+       IF ( xxd( idatmax1- idat ).LT.360.) exit
+       id1 = id1 + 1
+    ENDDO
+
+    DO idat = 1, idatmax
+       IF (xxd(idat).GT.0.) exit
+       id0 = id0 + 1
+    END DO
+
+    IF( id1 /= 0 ) then
+       DO idat = 1, id1
+          xxid(idat) = xxd(idatmax - id1 + idat) - 360.
+          fxd (idat) = fdd(idatmax - id1 + idat)     
+       END DO
+       DO idat = 1, idatmax - id1
+          xxid(idat + id1) = xxd(idat)
+          fxd (idat + id1) = fdd(idat)
+       END DO
+    end IF
+
+    IF(id0 /= 0) then
+       DO idat = 1, idatmax - id0
+          xxid(idat) = xxd(idat + id0)
+          fxd (idat) = fdd(idat + id0)
+       END DO
+
+       DO idat = 1, id0
+          xxid (idatmax - id0 + idat) =  xxd(idat) + 360.
+          fxd  (idatmax - id0 + idat) =  fdd(idat)   
+       END DO
+    else 
+       DO idat = 1, idatmax
+          xxid(idat)  = xxd(idat)
+          fxd (idat)  = fdd(idat)
+       ENDDO
+    end IF
+    xxid(idatmax1) = xxid(1) + 360.
+    fxd (idatmax1) = fxd(1)
+
+    !   initialisation du champ du modele
+
+    inter_barx(:) = 0.
+
+    ! iteration
+
+    x0   = xim0
+    dxm  = 0.
+    imod = 1
+    idat = 1
+
+    do while (imod <= imodmax)
+       do while (xxim(imod).GT.xxid(idat))
+          dx   = xxid(idat) - x0
+          dxm  = dxm + dx
+          inter_barx(imod) = inter_barx(imod) + dx * fxd(idat)
+          x0   = xxid(idat)
+          idat = idat + 1
+       end do
+       IF (xxim(imod).LT.xxid(idat)) THEN
+          dx   = xxim(imod) - x0
+          dxm  = dxm + dx
+          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
+          x0   = xxim(imod)
+          dxm  = 0.
+          imod = imod + 1
+       ELSE
+          dx   = xxim(imod) - x0
+          dxm  = dxm + dx
+          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
+          x0   = xxim(imod)
+          dxm  = 0.
+          imod = imod + 1
+          idat = idat + 1
+       END IF
+    end do
+
+  END function inter_barx
+
+  !******************************
+
+  function inter_bary(yjdat, fdat, yjmod)
+
+    ! Interpolation barycentrique basée sur les aires.
+    ! Version unidimensionnelle, en latitude.
+    ! L'indice 1 correspond à l'interface maille 1 -- maille 2.
+
+    use assert_m, only: assert
+
+    IMPLICIT NONE
+
+    REAL, intent(in):: yjdat(:)
+    ! (angles, ordonnées des interfaces des mailles des données, in
+    ! degrees, in increasing order)
+
+    REAL, intent(in):: fdat(:) ! champ de données
+
+    REAL, intent(in):: yjmod(:)
+    ! (ordonnées des interfaces des mailles du modèle)
+    ! (in degrees, in strictly increasing order)
+
+    REAL inter_bary(size(yjmod)) ! champ du modèle
+
+    ! Variables local to the procedure:
+
+    REAL y0, dy, dym 
+    INTEGER jdat ! indice du champ de données
+    integer jmod ! indice du champ du modèle
+
+    !------------------------------------
+
+    call assert(size(yjdat) == size(fdat), "inter_bary")
+
+    ! Initialisation des variables
+    inter_bary(:) = 0.
+    y0    = -90.
+    dym   = 0.
+    jmod  = 1
+    jdat  = 1
+
+    do while (jmod <= size(yjmod))
+       do while (yjmod(jmod) > yjdat(jdat))
+          dy         = yjdat(jdat) - y0
+          dym        = dym + dy
+          inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat)
+          y0         = yjdat(jdat)
+          jdat       = jdat + 1
+       end do
+       IF (yjmod(jmod) < yjdat(jdat)) THEN
+          dy         = yjmod(jmod) - y0
+          dym        = dym + dy
+          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
+          y0         = yjmod(jmod)
+          dym        = 0.
+          jmod       = jmod + 1
+       ELSE
+          ! {yjmod(jmod) == yjdat(jdat)}
+          dy         = yjmod(jmod) - y0
+          dym        = dym + dy
+          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
+          y0         = yjmod(jmod)
+          dym        = 0.
+          jmod       = jmod + 1
+          jdat       = jdat + 1
+       END IF
+    end do
+    ! Le test de fin suppose que l'interface 0 est commune aux deux
+    ! grilles "yjdat" et "yjmod".
+
+  END function inter_bary
+
+  !******************************
+
+  SUBROUTINE ord_coord(xi, xo, decrois)
+
+    ! This procedure receives an array of latitudes.
+    ! It converts them to degrees if they are in radians.
+    ! If the input latitudes are in decreasing order, the procedure
+    ! reverses their order.
+    ! Finally, the procedure adds 90° as the last value of the array.
+
+    use assert_eq_m, only: assert_eq
+
+    IMPLICIT NONE
+
+    include "comconst.h"
+    ! (for "pi")
+
+    REAL, intent(in):: xi(:)
+    ! (latitude, in degrees or radians, in increasing or decreasing order)
+    ! ("xi" should contain latitudes from pole to pole.
+    ! "xi" should contain the latitudes of the boundaries of grid
+    ! cells, not the centers of grid cells.
+    ! So the extreme values should not be 90° and -90°.)
+
+    REAL, intent(out):: xo(:) ! angles in degrees
+    LOGICAL, intent(out):: decrois
+
+    ! Variables  local to the procedure:
+    INTEGER nmax, i
+
+    !--------------------
+
+    nmax = assert_eq(size(xi), size(xo) - 1, "ord_coord")
+
+    ! Check monotonicity:
+    decrois = xi(2) < xi(1)
+    DO i = 3, nmax
+       IF (decrois .neqv. xi(i) < xi(i-1)) stop &
+            '"ord_coord":  latitudes are not monotonic'
+    ENDDO
+
+    IF (abs(xi(1)) < pi) then
+       ! "xi" contains latitudes in radians
+       xo(:nmax) = xi(:) * 180. / pi
+    else
+       ! "xi" contains latitudes in degrees
+       xo(:nmax) = xi(:)
+    end IF
+
+    IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN
+       print *, "ord_coord"
+       PRINT *, '"xi" should contain the latitudes of the boundaries of ' &
+            // 'grid cells, not the centers of grid cells.'
+       STOP
+    ENDIF
+
+    IF (decrois) xo(:nmax) = xo(nmax:1:- 1)
+    xo(nmax + 1) = 90.
+
+  END SUBROUTINE ord_coord
+
+  !***********************************
+
+  function ord_coordm(xi)
+
+    ! This procedure converts to degrees, if necessary, and inverts the
+    ! order.
+
+    IMPLICIT NONE
+
+    include "comconst.h"
+    ! (for "pi")
+
+    REAL, intent(in):: xi(:) ! angle, in rad or degrees
+    REAL ord_coordm(size(xi)) ! angle, in degrees
+
+    !-----------------------------
+
+    IF (xi(1) < 6.5) THEN
+       ! "xi" is in rad
+       ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi
+    else
+       ! "xi" is in degrees
+       ord_coordm(:) = xi(size(xi):1:-1)
+    ENDIF
+
+  END function ord_coordm
+
+end module inter_barxy_m
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/inter_bary.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/inter_bary.F	(revision 1292)
+++ 	(revision )
@@ -1,135 +1,0 @@
-!
-! $Header$
-!
-       SUBROUTINE inter_bary( jjm, jdatmax, yjdatt, fdatt  ,
-     ,                       jmodmax, yjmodd,  fmod      )
-c
-c    ...  Auteurs :  Robert Sadourny  , P. Le Van ...
-c
-       IMPLICIT NONE
-
-c  ----------------------------------------------------------
-c       INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES  .
-c         VERSION UNIDIMENSIONNELLE  ,    EN LATITUDE  .
-c  ----------------------------------------------------------
-c
-c     jdat : indice du champ de donnees, de 1 a jdatmax
-c     jmod : indice du champ du modele,  de 1 a jmodmax
-c     fdatt(jdatmax) : champ de donnees (entrees)
-c     yjdatt(jdatmax): ordonnees des interfaces des mailles donnees
-c     yjmodd(jmodmax): ordonnees des interfaces des mailles modele
-c     fmod(jmodmax)  : champ du modele  (sorties)
-c
-c      ( L'indice 1 correspond a l'interface maille 1 / maille 2)
-c      ( Les ordonnees sont exprimees en degres)
-c
-c     jdatmax = nb. d'interfaces  donnees =  nombre de donnees - 1 
-c     jmodmax = nb. d'interfaces  modele
-
-c     Si jmodmax = jjm , on veut interpoler sur les jjm+1 latitudes 
-c       rlatu   du modele ( lat.  des scalaires et de U ) 
-c
-c     Si jmodmax = jjp1 , on veut interpoler sur les jjm latitudes 
-c       rlatv du modele  ( lat.  de  V ) 
-
-c  ....  Arguments  en entree  .......
-
-       INTEGER jjm , jdatmax, jmodmax
-       REAL    yjdatt( jdatmax ) , fdatt( jdatmax +1 )
-       REAL    yjmodd( jmodmax )     
-
-c  ....  Arguments  en sortie  .......
-c
-       REAL    fmod( jmodmax + 1 )
-c
-c   ...... Variables locales  ......
-
-       INTEGER      jmods
-
-       REAL       yjdat ( jdatmax +1 ), fdat( jdatmax +1) 
-       REAL       fscrat( jdatmax +1 )
-       REAL       yjmod ( jmodmax +1 )
-       LOGICAL    decrois
-       SAVE       decrois
-c
-       REAL y0,dy,dym 
-       INTEGER jdat, jmod,i
-c
-
-        DO i = 1, jdatmax +1
-         fdat (i) = fdatt (i)
-        ENDDO
-
-       CALL ord_coord (  jdatmax , yjdatt(1), yjdat(1), decrois ) 
-
-       IF( decrois )   THEN
-         DO i = 1,jdatmax + 1
-          fscrat(i) = fdat(i)
-         ENDDO
-         DO i = 1, jdatmax + 1
-          fdat(i) = fscrat( jdatmax + 2 -i )
-         ENDDO
-
-       ENDIF
-
-       CALL ord_coordm (jmodmax,yjmodd(1),yjmod(1),jjm,jmods,decrois ) 
-c
-c      Initialisation des variables
-c    --------------------------------
-
-       DO jmod = 1, jmods
-        fmod(jmod) = 0.
-       END DO
-
-       y0    = 0.
-       dym   = 0.
-       jmod  = 1
-       jdat  = 1
-c  --------------------
-c      Iteration
-c  --------------------
-
-100    IF ( yjmod(jmod).LT.yjdat(jdat) ) THEN
-        dy         = yjmod(jmod) - y0
-        dym        = dym + dy
-        fmod(jmod) = (fmod(jmod) + dy * fdat(jdat)) / dym
-        y0         = yjmod(jmod)
-        dym        = 0.
-        jmod       = jmod + 1
-        GO TO 100
-
-       ELSE IF ( yjmod(jmod).GT.yjdat(jdat) ) THEN
-        dy         = yjdat(jdat) - y0
-        dym        = dym + dy
-        fmod(jmod) = fmod(jmod) + dy * fdat(jdat)
-        y0         = yjdat(jdat)
-        jdat       = jdat + 1
-
-       GO TO 100
-
-       ELSE
-        dy         = yjmod(jmod) - y0
-        dym        = dym + dy
-        fmod(jmod) = (fmod(jmod) + dy * fdat(jdat)) / dym
-        y0         = yjmod(jmod)
-        dym        = 0.
-        jmod       = jmod + 1
-        jdat       = jdat + 1
-
-        IF ( jmod.LE.jmods ) GO TO 100
-       END IF
-c   ---------------------------------------------
-c    Le test de fin suppose que l'interface 0
-c    est commune aux deux grilles yjdat et yjmod.
-c   ----------------------------------------------
-       IF( decrois )  THEN
-         DO i = 1,jmods
-          fscrat(i) = fmod(i)
-         ENDDO
-         DO i = 1, jmods
-          fmod(i) = fscrat( jmods + 1 -i )
-         ENDDO
-       ENDIF
-
-       RETURN
-       END
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/limit_netcdf.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/limit_netcdf.F	(revision 1292)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/limit_netcdf.F	(revision 1293)
@@ -9,4 +9,6 @@
       USE dimphy
       use phys_state_var_mod , ONLY : pctsrf
+      use inter_barxy_m, only: inter_barxy
+
       IMPLICIT none
 c
@@ -70,5 +72,5 @@
       REAL , ALLOCATABLE :: dlon(:), dlat(:)
       REAL , ALLOCATABLE :: dlon_ini(:), dlat_ini(:)
-      REAL , ALLOCATABLE :: champ_msk(:), champ(:)
+      REAL , ALLOCATABLE :: champ_msk(:), champ(:, :)
       REAL , ALLOCATABLE :: work(:,:)
 
@@ -335,5 +337,5 @@
       ENDIF
 c
-      ALLOCATE( champ(imdep*jmdep) )
+      ALLOCATE( champ(imdep, jmdep) )
 
       DO  200 l = 1, lmdep
@@ -364,8 +366,4 @@
 
        IF ( interbar )   THEN
-         DO j = 1, imdep * jmdep
-           champ(j) = LOG(champ(j))
-         ENDDO
-
         IF( l.EQ.1 )  THEN
          WRITE(6,*) '-------------------------------------------------',
@@ -376,6 +374,6 @@
      ,'------------------------'
         ENDIF
-        CALL inter_barxy ( imdep,jmdep -1,dlon,dlat,champ ,
-     ,                  iim,jjm,rlonu,rlatv, jjp1,champint )
+        CALL inter_barxy(dlon, dlat(:jmdep -1), log(champ), rlonu(:iim),
+     $       rlatv, champint)
          DO j=1,jjp1
           DO i=1,iim
@@ -556,5 +554,5 @@
       ENDIF 
 c
-      ALLOCATE ( champ(imdep*jmdep) )
+      ALLOCATE ( champ(imdep, jmdep) )
 
       DO l = 1, lmdep
@@ -604,6 +602,6 @@
 cIM   ENDDO
 cIM22/02/2002
-         CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
-     ,     champ, iim, jjm, rlonu, rlatv, jjp1, champint )
+       CALL inter_barxy (dlon, dlat(:jmdep -1), champ, rlonu(:iim),
+     $      rlatv, champint)
       ELSE
          CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
@@ -881,5 +879,5 @@
       ENDIF 
 
-       ALLOCATE( champ(imdep*jmdep) )
+       ALLOCATE( champ(imdep, jmdep) )
        IF( extrap )   THEN
          ALLOCATE ( work(imdep,jmdep) )
@@ -924,6 +922,6 @@
      ,'------------------------'
         ENDIF
-       CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
-     , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
+        CALL inter_barxy (dlon, dlat(:jmdep -1), champ, rlonu(:iim),
+     $       rlatv, champint)
       ELSE
        CALL grille_m(imdep, jmdep, dlon, dlat, champ,
@@ -1082,5 +1080,5 @@
       ENDIF 
 c
-      ALLOCATE ( champ(imdep*jmdep) )
+      ALLOCATE ( champ(imdep, jmdep) )
 
       DO l = 1, lmdep
@@ -1120,6 +1118,6 @@
         ENDIF
 
-       CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
-     , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
+       CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim),
+     $       rlatv, champint)
       ELSE
        CALL grille_m(imdep, jmdep, dlon, dlat, champ,
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/ord_coord.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/ord_coord.F	(revision 1292)
+++ 	(revision )
@@ -1,95 +1,0 @@
-!
-! $Header$
-!
-       SUBROUTINE ord_coord ( nmax, xi, xo, decrois )
-
-c   .... Auteur :  P. Le Van  .... 
-c
-c   ... Reordonne eventuellement les coordonnees de la grille donnees ...
-c
-       IMPLICIT NONE
-
-c    .....  Arguments  en entree  .....
-
-       INTEGER nmax
-       REAL xi(nmax)
-
-c    .....  Arguments  en sortie  .....
-c
-       REAL xo(nmax+1) 
-       LOGICAL decrois
-
-c    .... Variables  locales  ....
-
-       REAL xscr(nmax)
-       INTEGER i,ii
-       REAL pi, degres, chmin, chmax, mult
-c
-
-       pi     = 2.*ASIN(1.)
-       degres = 180./pi
-       decrois = .FALSE.
-
-       DO i = 1, nmax 
-        xo(i) = xi(i)
-       ENDDO
-     
-       mult = 1.
-       IF( xo(1).GT.xo(nmax) ) mult = -1.
-
-       CALL minmax(nmax,xo(1),chmin,chmax)
-
-       IF(chmax.LT.6.5 )  THEN
-          DO i = 1,nmax
-           xo(i) = xo(i) * degres 
-          ENDDO
-       ENDIF
-
-       IF( ABS( xo( 1  ) + mult* 90. ). LT .0.001. OR .
-     ,     ABS( xo(nmax) - mult* 90. ). LT .0.001      )  THEN
-        PRINT *,' Reverifier les valeurs de  xidat  pour les donnees .'
-        PRINT *,' Elles doivent correspondre aux interfaces et non aux',
-     , 'ordonnees des donnees,egales a -90. et 90.deg aux 2 extremites '
-         CALL ABORT
-       ENDIF
-
-       IF( xo(1).GT.xo(nmax) )   THEN
-          DO i = 1, nmax 
-            xscr(i) = xo(i)
-          ENDDO
-          DO i = 1, nmax 
-            xo(i+1) = xscr(i)
-          ENDDO
-            xo (   1    ) =   90.
-       ELSE
-            xo ( nmax +1) =   90.
-       ENDIF
-
-       IF ( xo(2).LT.xo(1) ) decrois =.TRUE.
-
-       DO i = 3, nmax 
-
-        IF(decrois.AND.xo(i).GT.xo(i-1) ) THEN
-         PRINT 1
-         PRINT 2,(xo(ii),ii=1,nmax)
-         CALL ABORT 
-        ENDIF
-        IF(.NOT.decrois.AND.xo(i).LT.xo(i-1) ) THEN
-         PRINT 1
-         PRINT 2,(xo(ii),ii=1,nmax)
-         CALL ABORT
-        ENDIF
-
-       ENDDO
-        
-       IF( decrois )  THEN
-c         CALL sort(nmax+1,xo(1))
-        CALL sort(nmax+1,xo)
-       ENDIF
-
-1      FORMAT(5x,' Incoherence dans les valeurs des latitudes de la ',
-     ,  'grille du modele ')
-2      FORMAT(1x,8f8.2)
-
-       RETURN
-       END
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/ord_coordm.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/ord_coordm.F	(revision 1292)
+++ 	(revision )
@@ -1,110 +1,0 @@
-!
-! $Header$
-!
-       SUBROUTINE ord_coordm ( nmax, xi, xo, jjm, jmods, decrois )
-
-c   ....  Auteur :  P. Le Van  .... 
-
-c   ... Reordonne eventuellement les coordonnees de la grille modele ...
-c
-       IMPLICIT NONE
-
-c    .....  Arguments  en entree  .....
-
-       INTEGER nmax,jjm
-       REAL xi(nmax)
-
-c    .....  Arguments  en sortie  .....
-c
-       REAL xo(nmax+1) 
-       LOGICAL decrois
-       INTEGER jmods
-
-c    .... Variables  locales  ....
-
-       REAL xscr(nmax)
-       INTEGER i
-       REAL pi, degres, chmin, chmax,mult
-c
-       DO i = 1, nmax 
-        xo(i) = xi(i)
-       ENDDO
-
-       mult = 1.
-       IF( xo(1).GT.xo(nmax) )  mult = - 1.
-       IF( nmax.EQ.jjm    ) jmods = nmax +1
-       IF( nmax.EQ.jjm +1 ) jmods = nmax -1 
-     
-       pi     = 2.*ASIN(1.)
-       degres = 180./pi
-       decrois = .FALSE.
-
-       CALL minmax(nmax,xo(1),chmin,chmax)
-
-       IF(chmax.LT.6.5 )  THEN
-          DO i = 1,nmax
-           xo(i) = xo(i) * degres 
-          ENDDO
-       ENDIF
-
-       IF( nmax.EQ.jjm )   THEN
-         IF( xo(1).GT.xo(nmax) )   THEN
-           DO i = 1, nmax 
-            xscr(i) = xo(i)
-           ENDDO
-           DO i = 1, nmax 
-            xo(i+1) = xscr(i)
-           ENDDO
-            xo (   1    ) =   90.
-         ELSE
-            xo ( nmax+1 ) =   90.
-         ENDIF
-       ELSE
-          IF( nmax.NE.jjm +1 )   THEN
-             PRINT *,'  Dans la routine ord_coordm , l argument nmax '
-             PRINT *,'  n est pas egal a jjm ni a jjm +1 . Corriger !'
-             CALL ABORT
-          ELSE
-            IF( ABS( xo(1)+ mult * 90.).GT.0.01 )  THEN
-              PRINT *,' Avec nmax =',nmax,'on devrait avoir des',
-     ,    ' ordonnees = 90. deg pour j=1 ou jjm+1 ! '
-             CALL ABORT
-            ELSE
-               IF( xo(1).LT.xo(nmax) )  THEN
-                 DO i = 1, nmax
-                  xscr(i) = xo(i)
-                 ENDDO
-                 DO i = 1, nmax -1
-                  xo(i) = xscr(i+1)
-                 ENDDO
-               ENDIF
-            ENDIF
-          ENDIF
-       ENDIF
-
-       IF ( xo(2).LT.xo(1) ) decrois =.TRUE.
-
-       DO i = 3, nmax 
-
-        IF(decrois.AND.xo(i).GT.xo(i-1) ) THEN
-         PRINT 1
-         CALL ABORT 
-        ENDIF
-        IF(.NOT.decrois.AND.xo(i).LT.xo(i-1) ) THEN
-         PRINT 1
-         CALL ABORT
-        ENDIF
-
-       ENDDO
-        
-       IF( decrois )  THEN
-         CALL sort(jmods,xo(1))
-       ENDIF
-
-
-1      FORMAT(5x,' Incoherence dans les valeurs des latitudes de la ',
-     ,  'grille du modele ')
-2      FORMAT(1x,8f8.2)
-
-       RETURN
-       END
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/startvar.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/startvar.F	(revision 1292)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/startvar.F	(revision 1293)
@@ -9,18 +9,16 @@
     !      There are three ways to access data from the database of atmospheric data which 
     !       can be used to initialize the model. This depends on the type of field which needs 
-    !       to be extracted. In any case the call should come after a restget and should be of the type :
-    !                CALL startget(...)
-    !
+    !       to be extracted. 
     !       We will details the possible arguments to startget here :
     !
     !        - A 2D variable on the dynamical grid :
-    !           CALL startget(varname, iml, jml, lon_in, lat_in, champ, val_ex, jml2, lon_in2, lat_in2, interbar )             
+    !           CALL startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_ex, jml2, lon_in2, lat_in2, interbar )             
     !
     !        - A 1D variable on the physical grid :
-    !            CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
+    !            CALL startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
     !
     !
     !         - A 3D variable on the dynamical grid :
-    !            CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
+    !            CALL startget_dyn(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
     !
     !
@@ -46,10 +44,5 @@
     !
       PRIVATE
-      PUBLIC startget
-    !
-    !
-      INTERFACE startget
-        MODULE PROCEDURE startget_phys2d, startget_phys1d, startget_dyn
-      END INTERFACE
+      public startget_phys2d, startget_phys1d, startget_dyn
     !
       INTEGER, SAVE :: fid_phys, fid_dyn
@@ -553,4 +546,6 @@
       SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in, jml2,
      .                 lon_in2, lat_in2 , interbar )
+
+      use inter_barxy_m, only: inter_barxy
     !
       INTEGER, INTENT(in) :: iml, jml ,jml2
@@ -648,6 +643,6 @@
         WRITE(6,*) '-------------------------------------------------',
      ,'--------------'
-        CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad ,
-     ,   var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var   ) 
+        CALL inter_barxy(lon_rad, lat_rad(:jml_phys -1), var_ana,
+     $       lon_in2(:iml-1), lat_in2(:jml-1), tmp_var) 
       ELSE
         CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
@@ -674,6 +669,6 @@
         WRITE(6,*) '-------------------------------------------------',
      ,'--------------'
-        CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad ,
-     ,   var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var  )
+        CALL inter_barxy(lon_rad, lat_rad(:jml_phys -1), var_ana,
+     $       lon_in2(:iml-1), lat_in2(:jml-1), tmp_var) 
       ELSE
         CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
@@ -691,17 +686,20 @@
     !
     !
-      SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in,
-     . lml, pls, workvar, champ, val_exp,jml2, lon_in2, lat_in2 ,
+      SUBROUTINE startget_dyn(varname, lon_in, lat_in,
+     . pls, workvar, champ, val_exp, lon_in2, lat_in2 ,
      ,  interbar )
+
+      use assert_eq_m, only: assert_eq
     !
     !   ARGUMENTS
     !
-      CHARACTER*(*), INTENT(in) :: varname
-      INTEGER, INTENT(in) :: iml, jml, lml, jml2
-      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
-      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
-      REAL, INTENT(in) :: pls(iml, jml, lml)
-      REAL, INTENT(in) :: workvar(iml, jml, lml)
-      REAL, INTENT(inout) :: champ(iml, jml, lml)
+      CHARACTER(len=*), INTENT(in) :: varname
+      REAL, INTENT(in) :: lon_in(:) ! dim(iml)
+      REAL, INTENT(in) :: lat_in(:) ! dim(jml)
+      REAL, INTENT(in) :: lon_in2(:) ! dim(iml)
+      REAL, INTENT(in) :: lat_in2(:) ! dim(jml2)
+      REAL, INTENT(in) :: pls(:, :, :) ! dim(iml, jml, lml)
+      REAL, INTENT(in) :: workvar(:, :, :) ! dim(iml, jml, lml)
+      REAL, INTENT(inout) :: champ(:, :, :) ! dim(iml, jml, lml)
       REAL, INTENT(in) :: val_exp
       LOGICAL interbar
@@ -709,5 +707,5 @@
     !    LOCAL
     !
-      INTEGER :: il, ij, ii
+      INTEGER :: il, ij, ii, iml, jml, lml, jml2
       REAL :: xppn, xpps
     !
@@ -719,4 +717,14 @@
     !   This routine only works if the variable does not exist or is constant
     !
+C     -----------------------------
+
+      iml = assert_eq((/size(lon_in), size(pls, 1), size(workvar, 1),
+     $     size(champ, 1), size(lon_in2)/), "startget_dyn iml")
+      jml = assert_eq(size(lat_in), size(pls, 2), size(workvar, 2),
+     $     size(champ, 2), "startget_dyn jml")
+      lml = assert_eq(size(pls, 3), size(workvar, 3), size(champ, 3),
+     $     "startget_dyn lml")
+      jml2 = size(lat_in2)
+
       IF ( MINVAL(champ(:,:,:)).EQ.MAXVAL(champ(:,:,:)) .AND.
      . MINVAL(champ(:,:,:)).EQ.val_exp ) THEN
@@ -834,4 +842,6 @@
      ,             lat_in2 , interbar )
     !
+      use inter_barxy_m, only: inter_barxy
+
       INTEGER, INTENT(in) :: iml, jml, jml2
       REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
@@ -929,6 +939,6 @@
         WRITE(6,*) '-------------------------------------------------',
      ,'--------------'
-        CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad ,
-     ,    var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var) 
+        CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana,
+     $       lon_in2(:iml-1), lat_in2(:jml-1), tmp_var) 
       ELSE
         CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
@@ -954,6 +964,6 @@
         WRITE(6,*) '-------------------------------------------------',
      ,'--------------'
-        CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad ,
-     ,    var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var) 
+        CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1), var_ana,
+     $       lon_in2(:iml-1), lat_in2(:jml-1), tmp_var) 
       ELSE
         CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
@@ -1035,4 +1045,6 @@
     !
     !
+      use inter_barxy_m, only: inter_barxy
+
     !    ARGUMENTS
     !
@@ -1120,6 +1132,6 @@
      ,'--------------'
        ENDIF
-       CALL inter_barxy ( iml_dyn, jml_dyn -1,lon_rad, lat_rad, 
-     , var_ana3d(:,:,il),iml-1, jml2, lon_in2, lat_in2,jml,var_tmp2d ) 
+       CALL inter_barxy(lon_rad, lat_rad(:jml_dyn -1),
+     $      var_ana3d(:,:,il), lon_in2(:iml-1), lat_in2, var_tmp2d) 
       ELSE
        CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad, 
