Index: LMDZ.3.3/trunk/libf/dyn3d/conf_dat2d.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/conf_dat2d.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/conf_dat2d.F	(revision 259)
@@ -0,0 +1,217 @@
+C
+C $Header$
+C
+      SUBROUTINE conf_dat2d( title,lons,lats,xd,yd,xf,yf,champd ,
+     ,                           interbar                        )
+c
+c     Auteur :  P. Le Van
+
+c    Ce s-pr. configure le champ de donnees 2D 'champd' de telle facon que
+c       qu'on ait     - pi    a    pi    en longitude
+c       et qu'on ait   pi/2.  a - pi/2.  en latitude
+c
+c      xd et yd  sont les longitudes et latitudes initiales
+c      xf et yf  sont les longitudes et latitudes en sortie , eventuellement
+c      modifiees pour etre configurees comme ci-dessus .
+
+      IMPLICIT NONE
+ 
+c    ***       Arguments en  entree      ***
+      INTEGER lons,lats
+      CHARACTER*25 title
+      REAL xd(lons),yd(lats)
+      LOGICAL interbar
+c
+c    ***       Arguments en  sortie      ***
+      REAL xf(lons),yf(lats)
+c
+c    ***  Arguments en entree et  sortie ***
+      REAL champd(lons,lats)
+
+c   ***     Variables  locales  ***
+c
+      REAL pi,pis2,depi
+      LOGICAL radianlon, invlon ,radianlat, invlat, alloc
+      REAL rlatmin,rlatmax,oldxd1
+      INTEGER i,j,ip180,ind
+
+      REAL, ALLOCATABLE :: xtemp(:) 
+      REAL, ALLOCATABLE :: ytemp(:) 
+      REAL, ALLOCATABLE :: champf(:,:)
+     
+c
+c      WRITE(6,*) ' conf_dat2d  pour la variable ', title
+
+      ALLOCATE( xtemp(lons) )
+      ALLOCATE( ytemp(lats) )
+
+      DO i = 1, lons
+       xtemp(i) = xd(i)
+      ENDDO
+      DO j = 1, lats
+       ytemp(j) = yd(j)
+      ENDDO
+
+      pi   = 2. * ASIN(1.) 
+      pis2 = pi/2.
+      depi = 2. * pi
+
+            radianlon = .FALSE.
+      IF( xtemp(1).GE.-pi-0.5.AND. xtemp(lons).LE.pi+0.5 )  THEN
+            radianlon = .TRUE.
+            invlon    = .FALSE.
+      ELSE IF (xtemp(1).GE.-0.5.AND.xtemp(lons).LE.depi+0.5 ) THEN
+            radianlon = .TRUE.
+            invlon    = .TRUE.
+      ELSE IF ( xtemp(1).GE.-180.5.AND. xtemp(lons).LE.180.5 )   THEN
+            radianlon = .FALSE.
+            invlon    = .FALSE.
+      ELSE IF ( xtemp(1).GE.-0.5.AND.xtemp(lons).LE.360.5 )   THEN
+            radianlon = .FALSE.
+            invlon    = .TRUE.
+      ELSE
+        WRITE(6,*) 'Pbs. sur les longitudes des donnees pour le fichier'
+     ,  , title
+      ENDIF
+
+      invlat = .FALSE.
+      
+      IF( ytemp(1).LT.ytemp(lats) ) THEN
+        invlat = .TRUE.
+      ENDIF
+
+      rlatmin = MIN( ytemp(1), ytemp(lats) )
+      rlatmax = MAX( ytemp(1), ytemp(lats) )
+      
+      IF( rlatmin.GE.-pis2-0.5.AND.rlatmax.LE.pis2+0.5)THEN
+             radianlat = .TRUE.
+      ELSE IF ( rlatmin.GE.-90.-0.5.AND.rlatmax.LE.90.+0.5 ) THEN
+             radianlat = .FALSE.
+      ELSE
+        WRITE(6,*) ' Pbs. sur les latitudes des donnees pour le fichier'
+     ,  , title
+      ENDIF
+
+       IF( .NOT. radianlon )  THEN
+         DO i = 1, lons
+          xtemp(i) = xtemp(i) * pi/180.
+         ENDDO
+       ENDIF
+
+       IF( .NOT. radianlat )  THEN
+         DO j = 1, lats
+          ytemp(j) = ytemp(j) * pi/180.
+         ENDDO   
+       ENDIF
+
+
+        IF ( invlon )   THEN
+
+           ALLOCATE(champf(lons,lats))
+           alloc = .TRUE.
+
+           DO i = 1 ,lons
+            xf(i) = xtemp(i)
+           ENDDO
+c
+c    ***  On tourne les longit.  pour avoir  - pi  a  +  pi  ****
+c
+           DO i=1,lons
+            IF( xf(i).GT. pi )  THEN
+            GO TO 88
+            ENDIF
+           ENDDO
+
+88         CONTINUE
+c
+           ip180 = i
+
+           DO i = 1,lons
+            IF (xf(i).GT. pi)  THEN
+             xf(i) = xf(i) - depi
+            ENDIF
+           ENDDO
+
+           DO i= ip180,lons
+            ind = i-ip180 +1
+            xtemp(ind) = xf(i)
+           ENDDO
+
+           DO i= ind +1,lons
+            xtemp(i) = xf(i-ind)
+           ENDDO
+
+c   .....    on tourne les longitudes  pour  champf ....
+c
+           DO j = 1,lats
+
+             DO i = ip180,lons
+              ind  = i-ip180 +1
+              champd (ind,j) = champf (i,j)
+             ENDDO
+   
+             DO i= ind +1,lons
+              champd (i,j)  = champf (i-ind,j)
+             ENDDO
+
+           ENDDO
+
+
+        ENDIF
+c
+c    *****   fin  de   IF(invlon)   ****
+
+         IF ( invlat )    THEN
+
+           IF(.NOT.alloc)  THEN
+            ALLOCATE(champf(lons,lats))
+            alloc = .TRUE.
+           ENDIF
+
+           DO j = 1,lats
+            yf(j) = ytemp(j)
+           ENDDO
+
+           DO j = 1, lats
+             DO i = 1,lons
+              champf(i,j) = champd(i,j)
+             ENDDO
+           ENDDO
+
+           DO j = 1, lats
+              ytemp( lats-j+1 ) = yf(j)
+              DO i = 1, lons
+               champd (i,lats-j+1) = champf (i,j)
+              ENDDO
+           ENDDO
+
+
+         ENDIF
+
+c    *****  fin  de  IF(invlat)   ****
+
+c        
+      IF( interbar )  THEN
+        oldxd1 = xtemp(1)
+        DO i = 1, lons -1
+          xtemp(i) = 0.5 * ( xtemp(i) + xtemp(i+1) )
+        ENDDO
+          xtemp(lons) = 0.5 * ( xtemp(lons) + oldxd1 + depi )
+
+        DO j = 1, lats -1
+          ytemp(j) = 0.5 * ( ytemp(j) + ytemp(j+1) )
+        ENDDO
+
+      ENDIF
+c
+       IF(alloc) DEALLOCATE(champf)
+
+       DO i = 1, lons
+        xf(i) = xtemp(i)
+       ENDDO
+       DO j = 1, lats
+        yf(j) = ytemp(j)
+       ENDDO
+
+      RETURN
+      END
Index: LMDZ.3.3/trunk/libf/dyn3d/conf_dat3d.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/conf_dat3d.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/conf_dat3d.F	(revision 259)
@@ -0,0 +1,293 @@
+C
+C $Header$
+C
+      SUBROUTINE conf_dat3d( title, lons,lats,levs,xd,yd,zd,xf,yf,zf,
+     ,                                 champd , interbar             )
+c
+c     Auteur : P. Le Van
+c
+c    Ce s-pr. configure le champ de donnees 3D 'champd' de telle facon 
+c       qu'on ait     - pi    a    pi    en longitude
+c       qu'on ait      pi/2.  a - pi/2.  en latitude
+c      et qu'on ait les niveaux verticaux variant du sol vers le ht de l'atmos.
+c           (     en Pascals   ) .
+c
+c      xd et yd  sont les longitudes et latitudes initiales
+c      zd  les pressions initiales
+c
+c      xf et yf  sont les longitudes et latitudes en sortie , eventuellement
+c       modifiees pour etre configurees comme ci-dessus .
+c      zf  les pressions en sortie
+c
+c      champd   en meme temps le champ initial et  final
+c
+c      interbar = .TRUE.  si on appelle l'interpo. barycentrique inter_barxy
+c          sinon , l'interpolation   grille_m  ( grid_atob ) .
+c
+
+      IMPLICIT NONE
+ 
+c    ***       Arguments en  entree      ***
+      CHARACTER*(*) :: title
+      INTEGER lons, lats, levs
+      REAL xd(lons), yd(lats), zd(levs)
+      LOGICAL interbar
+c
+c    ***       Arguments en  sortie      ***
+      REAL xf(lons), yf(lats), zf(levs)
+
+c    ***  Arguments en entree et  sortie ***
+      REAL  champd(lons,lats,levs)
+
+c    ***  Variables locales  ***
+c
+      REAL pi,pis2,depi,presmax
+      LOGICAL radianlon, invlon ,radianlat, invlat, invlev, alloc
+      REAL rlatmin,rlatmax,oldxd1
+      INTEGER i,j,ip180,ind,l
+
+      REAL, ALLOCATABLE :: xtemp(:)
+      REAL, ALLOCATABLE :: ytemp(:)
+      REAL, ALLOCATABLE :: ztemp(:)
+      REAL, ALLOCATABLE :: champf(:,:,:)
+     
+
+c      WRITE(6,*) '  Conf_dat3d  pour  ',title
+
+      ALLOCATE(xtemp(lons))
+      ALLOCATE(ytemp(lats))
+      ALLOCATE(ztemp(levs))
+
+      DO i = 1, lons
+       xtemp(i) = xd(i)
+      ENDDO
+      DO j = 1, lats
+       ytemp(j) = yd(j)
+      ENDDO
+      DO l = 1, levs
+       ztemp(l) = zd(l)
+      ENDDO
+
+      pi   = 2. * ASIN(1.) 
+      pis2 = pi/2.
+      depi = 2. * pi
+
+      IF( xtemp(1).GE.-pi-0.5.AND. xtemp(lons).LE.pi+0.5 )  THEN
+            radianlon = .TRUE.
+            invlon    = .FALSE.
+      ELSE IF (xtemp(1).GE.-0.5.AND.xtemp(lons).LE.depi+0.5 ) THEN
+            radianlon = .TRUE.
+            invlon    = .TRUE.
+      ELSE IF ( xtemp(1).GE.-180.5.AND. xtemp(lons).LE.180.5 )   THEN
+            radianlon = .FALSE.
+            invlon    = .FALSE.
+      ELSE IF ( xtemp(1).GE.-0.5.AND.xtemp(lons).LE.360.5 )   THEN
+            radianlon = .FALSE.
+            invlon    = .TRUE.
+      ELSE
+        WRITE(6,*) 'Pbs. sur les longitudes des donnees pour le fichier'
+     ,  , title
+      ENDIF
+
+      invlat = .FALSE.
+      
+      IF( ytemp(1).LT.ytemp(lats) ) THEN
+        invlat = .TRUE.
+      ENDIF
+
+      rlatmin = MIN( ytemp(1), ytemp(lats) )
+      rlatmax = MAX( ytemp(1), ytemp(lats) )
+      
+      IF( rlatmin.GE.-pis2-0.5.AND.rlatmax.LE.pis2+0.5)THEN
+             radianlat = .TRUE.
+      ELSE IF ( rlatmin.GE.-90.-0.5.AND.rlatmax.LE.90.+0.5 ) THEN
+             radianlat = .FALSE.
+      ELSE
+        WRITE(6,*) ' Pbs. sur les latitudes des donnees pour le fichier'
+     ,  , title
+      ENDIF
+
+       IF( .NOT. radianlon )  THEN
+         DO i = 1, lons
+          xtemp(i) = xtemp(i) * pi/180.
+         ENDDO
+       ENDIF
+
+       IF( .NOT. radianlat )  THEN
+         DO j = 1, lats
+          ytemp(j) = ytemp(j) * pi/180.
+         ENDDO   
+       ENDIF
+
+
+        alloc =.FALSE.
+
+        IF ( invlon )   THEN
+
+            ALLOCATE(champf(lons,lats,levs))
+            alloc = .TRUE.
+
+            DO i = 1 ,lons
+             xf(i) = xtemp(i)
+            ENDDO
+
+            DO l = 1, levs
+             DO j = 1, lats
+              DO i= 1, lons
+               champf (i,j,l)  = champd (i,j,l)
+              ENDDO
+             ENDDO
+            ENDDO
+c
+c    ***  On tourne les longit.  pour avoir  - pi  a  +  pi  ****
+c
+            DO i=1,lons
+             IF( xf(i).GT. pi )  THEN
+              GO TO 88
+             ENDIF
+            ENDDO
+
+88          CONTINUE
+c
+            ip180 = i
+
+            DO i = 1,lons
+             IF (xf(i).GT. pi)  THEN
+              xf(i) = xf(i) - depi
+             ENDIF
+            ENDDO
+
+            DO i= ip180,lons
+             ind = i-ip180 +1
+             xtemp(ind) = xf(i)
+            ENDDO
+
+            DO i= ind +1,lons
+             xtemp(i) = xf(i-ind)
+            ENDDO
+
+c   .....    on tourne les longitudes  pour champf  ....
+c
+            DO l = 1,levs
+              DO j = 1,lats
+               DO i = ip180,lons
+                ind  = i-ip180 +1
+                champd (ind,j,l) = champf (i,j,l)
+               ENDDO
+   
+               DO i= ind +1,lons
+                champd (i,j,l)  = champf (i-ind,j,l)
+               ENDDO
+              ENDDO
+            ENDDO
+
+            DEALLOCATE(xtemp)
+        ENDIF
+c
+c    *****   fin  de   IF(invlon)   ****
+         
+         IF ( invlat )    THEN
+
+           IF(.NOT.alloc)  THEN 
+            ALLOCATE(champf(lons,lats,levs))
+            alloc = .TRUE.
+           ENDIF
+
+           DO j = 1, lats
+            yf(j) = ytemp(j)
+           ENDDO
+         
+           DO l = 1,levs
+            DO j = 1, lats
+             DO i = 1,lons
+              champf(i,j,l) = champd(i,j,l)
+             ENDDO
+            ENDDO
+
+            DO j = 1, lats
+              ytemp( lats-j+1 ) = yf(j)
+              DO i = 1, lons
+               champd (i,lats-j+1,l) = champf (i,j,l)
+              ENDDO
+            ENDDO
+          ENDDO
+
+          DEALLOCATE(ytemp)
+
+         ENDIF
+
+c    *****  fin  de  IF(invlat)   ****
+c
+c
+      IF( interbar )  THEN
+        oldxd1 = xtemp(1)
+        DO i = 1, lons -1
+          xtemp(i) = 0.5 * ( xtemp(i) + xtemp(i+1) )
+        ENDDO
+          xtemp(lons) = 0.5 * ( xtemp(lons) + oldxd1 + depi )
+
+        DO j = 1, lats -1
+          ytemp(j) = 0.5 * ( ytemp(j) + ytemp(j+1) )
+        ENDDO
+      ENDIF
+c
+
+      invlev = .FALSE.
+      IF( ztemp(1).LT.ztemp(levs) )  invlev = .TRUE.
+
+      presmax = MAX( ztemp(1), ztemp(levs) )
+      IF( presmax.LT.1200. ) THEN
+         DO l = 1,levs
+           ztemp(l) = ztemp(l) * 100.
+         ENDDO
+      ENDIF
+
+      IF( invlev )  THEN
+
+          IF(.NOT.alloc)  THEN
+            ALLOCATE(champf(lons,lats,levs))
+            alloc = .TRUE.
+          ENDIF
+
+          DO l = 1,levs
+            zf(l) = ztemp(l)
+          ENDDO
+
+          DO l = 1,levs
+            DO j = 1, lats
+             DO i = 1,lons
+              champf(i,j,l) = champd(i,j,l)
+             ENDDO
+            ENDDO
+          ENDDO
+
+          DO l = 1,levs
+            ztemp(levs+1-l) = zf(l)
+          ENDDO
+
+          DO l = 1,levs
+            DO j = 1, lats
+             DO i = 1,lons
+              champd(i,j,levs+1-l) = champf(i,j,l)
+             ENDDO
+            ENDDO
+          ENDDO
+
+          DEALLOCATE(ztemp)
+
+      ENDIF
+
+         IF(alloc)  DEALLOCATE(champf)
+
+         DO i = 1, lons
+           xf(i) = xtemp(i)
+         ENDDO
+         DO j = 1, lats
+           yf(j) = ytemp(j)
+         ENDDO
+         DO l = 1, levs
+           zf(l) = ztemp(l)
+         ENDDO
+
+      RETURN
+      END
Index: LMDZ.3.3/trunk/libf/dyn3d/create_etat0_limit.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/create_etat0_limit.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/create_etat0_limit.F	(revision 259)
@@ -0,0 +1,28 @@
+C
+C $Header$
+C
+       PROGRAM create_etat0_limit
+c
+c
+c     Programme d'appel a etat0, creation des etats initiaux et limit_netcdf
+c   
+c
+c     interbar = .T . si appel a  interpol. barycentrique inter_barxy
+c
+c     extrap   = .T . si on fait une extrapolation de donnees , comme pour
+c       les  SST  lorsque  le fichier ne contient pas uniquement  des points 
+c     oceaniques .
+c
+c     oldice   = .T. si l'on veut garder les anciennes glaces , obtenues
+c     par  grille_m  ( grid_atob ) .
+
+
+      LOGICAL interbar, extrap , oldice
+      PARAMETER ( interbar = .TRUE. , extrap = .FALSE. , oldice=.TRUE.)
+
+      CALL etat0_netcdf ( interbar )
+c
+      CALL  limit_netcdf ( interbar, extrap , oldice )
+
+      STOP
+      END
Index: LMDZ.3.3/trunk/libf/dyn3d/etat0_netcdf.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/etat0_netcdf.F	(revision 256)
+++ LMDZ.3.3/trunk/libf/dyn3d/etat0_netcdf.F	(revision 259)
@@ -1,3 +1,6 @@
-      SUBROUTINE etat0_netcdf
+C
+C $Header$
+C
+      SUBROUTINE etat0_netcdf ( interbar )
     
       USE startvar
@@ -17,4 +20,5 @@
 #include "comconst.h"
       !
+      LOGICAL interbar
       REAL :: latfi(klon), lonfi(klon)
       REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
@@ -34,5 +38,5 @@
       CHARACTER*80 :: varname
       !
-      INTEGER :: i,j, ig, l
+      INTEGER :: i,j, ig, l,ii1,ii2
       REAL :: xpi
       !
@@ -65,4 +69,5 @@
       REAL ::phystep,co2_ppm,solaire
       INTEGER :: radpas
+       real zrel(iip1*jjp1),chmin,chmax
 
       CHARACTER*80 :: visu_file
@@ -122,5 +127,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)
+       CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
+     , jjm ,rlonu,rlatv , interbar )
       !
       WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
@@ -130,5 +136,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)
+       CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
+     , jjm, rlonu,rlatv , interbar )
       !
       WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite' 
@@ -138,13 +145,15 @@
       ! This line needs to be replaced by a call to restget to get the values in the restart file
       masque(:,:) = 0.0
-       CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0)
-      !
-      WRITE(*,*) 'MASQUE construit : Masque'
-      WRITE(*,'(49I1)') INT(masque(:,:))
+       CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
+     , jjm ,rlonu,rlatv , interbar )
+      !
+c      WRITE(*,*) 'MASQUE construit : Masque'
+c      WRITE(*,'(49I1)') INT(masque(:,:))
       !
       !
       varname = 'psol'
       psol(:,:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0)
+      CALL startget(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 
@@ -174,20 +183,21 @@
       varname = 'surfgeo'
       phis(:,:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0)
+      CALL startget(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)
+     . workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
       !  
       varname = 'v'
       vvent(:,:,:) = 0.0
       CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
-     . workvar, vvent, 0.0)
+     . workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
       !
       varname = 't'
       t3d(:,:,:) = 0.0
       CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
-     . workvar, t3d, 0.0)
+     . workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
       !
       WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
@@ -196,5 +206,5 @@
       tpot(:,:,:) = 0.0
       CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
-     . pk, tpot, 0.0)
+     . pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
       !
       WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
@@ -221,5 +231,5 @@
      .                           maxval(qsat(:,:,:))
       CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
-     . qsat, qd, 0.0)
+     . qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
       q3d(:,:,:,1) = qd(:,:,:)
       !
@@ -227,5 +237,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)
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
+     .    jjm, rlonu, rlatv , interbar )
       !
       WRITE(*,*) 'TSOL construit :'
@@ -234,54 +245,77 @@
       varname = 'qsol'
       qsol(:) = 0.0
-      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0)
+      CALL startget(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)
+      CALL startget(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)
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
+     .    jjm, rlonu, rlatv , interbar )
       !
       varname = 'deltat'
       deltat(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0)
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0,
+     .     jjm, rlonu, rlatv , interbar )
       !
       varname = 'rugmer'
       rugmer(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0)
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
+     .     jjm, rlonu, rlatv , interbar )
       !
       varname = 'agsno'
       agesno(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0)
+      CALL startget(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(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(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(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(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(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(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(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+c
       rugsrel(:) = 0.0
-C Calcul intermediaire
+      IF(ok_orodr)  THEN
+        DO i = 1, iip1* jjp1
+         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
+        ENDDO
+
+      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite  rugsrel'
+      WRITE(*,'(49I1)') INT(rugsrel(:)*10)
+
+      ENDIF
+c
+c  Calcul intermediaire
 c 
       CALL massdair( p3d, masse  )
 c
-
-      print *,' ALPHAX ',alphax
 
       DO  l = 1, llm
@@ -305,5 +339,4 @@
       CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
      *                tetagdiv, tetagrot , tetatemp              )
-      print*,'sortie inidissip'
       itau = 0
       iday = dayref +itau/day_step
@@ -315,9 +348,8 @@
       ENDIF
       CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
-      print*,'sortie geopot'
       
       CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
      *                phi,w, pbaru,pbarv,time+iday-dayref   )
-       print*,'sortie caldyn0'     
+      print*,'sortie caldyn0'     
       CALL dynredem0("start.nc",dayref,anneeref,phis,nqmx)
       print*,'sortie dynredem0'
@@ -339,4 +371,5 @@
      .                   zgam, zthe, zpic, zval,
      .                   rugsrel)
+cc
 
 C     Sortie Visu pour les champs dynamiques
Index: LMDZ.3.3/trunk/libf/dyn3d/extrapol.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/extrapol.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/extrapol.F	(revision 259)
@@ -0,0 +1,198 @@
+C
+C $Header$
+C
+      SUBROUTINE extrapol (pfild, kxlon, kylat, pmask,
+     .                   norsud, ldper, knbor, pwork)
+      IMPLICIT none
+c
+c OASIS routine (Adaptation: Laurent Li, le 14 mars 1997)
+c Fill up missed values by using the neighbor points
+c
+      INTEGER kxlon, kylat ! longitude and latitude dimensions (Input)
+      INTEGER knbor ! minimum neighbor number (Input)
+      LOGICAL norsud ! True if field is from North to South (Input)
+      LOGICAL ldper ! True if take into account the periodicity (Input)
+      REAL pmask ! mask value (Input)
+      REAL pfild(kxlon,kylat) ! field to be extrapolated (Input/Output)
+      REAL pwork(kxlon,kylat) ! working space
+c
+      REAL zwmsk
+      INTEGER incre, idoit, i, j, k, inbor, ideb, ifin, ilon, jlat
+      INTEGER ix(9), jy(9) ! index arrays for the neighbors coordinates
+      REAL zmask(9)
+C
+C  We search over the eight closest neighbors
+C
+C            j+1  7  8  9
+C              j  4  5  6    Current point 5 --> (i,j)
+C            j-1  1  2  3
+C                i-1 i i+1
+c
+c
+      IF (norsud) THEN
+         DO j = 1, kylat
+         DO i = 1, kxlon
+            pwork(i,j) = pfild(i,kylat-j+1)
+         ENDDO
+         ENDDO
+         DO j = 1, kylat
+         DO i = 1, kxlon
+            pfild(i,j) = pwork(i,j)
+         ENDDO
+         ENDDO
+      ENDIF
+c
+      incre = 0
+c
+      DO j = 1, kylat
+      DO i = 1, kxlon
+         pwork(i,j) = pfild(i,j)
+      ENDDO
+      ENDDO
+c
+C* To avoid problems in floating point tests
+      zwmsk = pmask - 1.0
+c
+200   CONTINUE
+      incre = incre + 1
+      DO 99999 j = 1, kylat
+      DO 99999 i = 1, kxlon
+      IF (pfild(i,j).GT. zwmsk) THEN
+         pwork(i,j) = pfild(i,j)
+         inbor = 0
+         ideb = 1
+         ifin = 9
+C
+C* Fill up ix array
+         ix(1) = MAX (1,i-1)
+         ix(2) = i
+         ix(3) = MIN (kxlon,i+1)
+         ix(4) = MAX (1,i-1)
+         ix(5) = i
+         ix(6) = MIN (kxlon,i+1)
+         ix(7) = MAX (1,i-1)
+         ix(8) = i
+         ix(9) = MIN (kxlon,i+1)
+C
+C* Fill up iy array
+         jy(1) = MAX (1,j-1)
+         jy(2) = MAX (1,j-1)
+         jy(3) = MAX (1,j-1)
+         jy(4) = j
+         jy(5) = j
+         jy(6) = j
+         jy(7) = MIN (kylat,j+1)
+         jy(8) = MIN (kylat,j+1)
+         jy(9) = MIN (kylat,j+1)
+C
+C* Correct latitude bounds if southernmost or northernmost points
+         IF (j .EQ. 1) ideb = 4
+         IF (j .EQ. kylat) ifin = 6
+C
+C* Account for periodicity in longitude
+C
+         IF (ldper) THEN 
+            IF (i .EQ. kxlon) THEN
+               ix(3) = 1
+               ix(6) = 1
+               ix(9) = 1
+            ELSE IF (i .EQ. 1) THEN
+               ix(1) = kxlon
+               ix(4) = kxlon
+               ix(7) = kxlon
+            ENDIF
+         ELSE
+            IF (i .EQ. 1) THEN
+               ix(1) = i
+               ix(2) = i + 1
+               ix(3) = i
+               ix(4) = i + 1
+               ix(5) = i
+               ix(6) = i + 1
+            ENDIF 
+            IF (i .EQ. kxlon) THEN
+               ix(1) = i -1
+               ix(2) = i
+               ix(3) = i - 1
+               ix(4) = i
+               ix(5) = i - 1
+               ix(6) = i
+            ENDIF
+C
+            IF (i .EQ. 1 .OR. i .EQ. kxlon) THEN 
+               jy(1) = MAX (1,j-1)
+               jy(2) = MAX (1,j-1)
+               jy(3) = j
+               jy(4) = j
+               jy(5) = MIN (kylat,j+1)
+               jy(6) = MIN (kylat,j+1)
+C
+               ideb = 1
+               ifin = 6
+               IF (j .EQ. 1) ideb = 3
+               IF (j .EQ. kylat) ifin = 4
+            ENDIF
+         ENDIF ! end for ldper test
+C
+C* Find unmasked neighbors
+C
+         DO 230 k = ideb, ifin
+            zmask(k) = 0.
+            ilon = ix(k)
+            jlat = jy(k)
+            IF (pfild(ilon,jlat) .LT. zwmsk) THEN
+               zmask(k) = 1.
+               inbor = inbor + 1
+            ENDIF
+ 230     CONTINUE
+C
+C* Not enough points around point P are unmasked; interpolation on P 
+C  will be done in a future call to extrap.
+C
+         IF (inbor .GE. knbor) THEN
+            pwork(i,j) = 0.
+            DO k = ideb, ifin
+               ilon = ix(k)
+               jlat = jy(k)
+               pwork(i,j) = pwork(i,j)
+     $                      + pfild(ilon,jlat) * zmask(k)/FLOAT(inbor)
+            ENDDO
+         ENDIF
+C
+      ENDIF
+99999 CONTINUE
+C
+C*    3. Writing back unmasked field in pfild
+C        ------------------------------------
+C
+C* pfild then contains:
+C     - Values which were not masked
+C     - Interpolated values from the inbor neighbors
+C     - Values which are not yet interpolated
+C
+      idoit = 0
+      DO j = 1, kylat
+      DO i = 1, kxlon
+         IF (pwork(i,j) .GT. zwmsk) idoit = idoit + 1
+         pfild(i,j) = pwork(i,j)
+      ENDDO
+      ENDDO
+c
+      IF (idoit .ne. 0) GOTO 200
+ccc      PRINT*, "Number of extrapolation steps incre =", incre
+c
+      IF (norsud) THEN
+         DO j = 1, kylat
+         DO i = 1, kxlon
+            pwork(i,j) = pfild(i,kylat-j+1)
+         ENDDO
+         ENDDO
+         DO j = 1, kylat
+         DO i = 1, kxlon
+            pfild(i,j) = pwork(i,j)
+         ENDDO
+         ENDDO
+      ENDIF
+c
+      RETURN
+      END
Index: LMDZ.3.3/trunk/libf/dyn3d/inter_barx.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/inter_barx.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/inter_barx.F	(revision 259)
@@ -0,0 +1,228 @@
+C
+C $Header$
+C
+       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)
+      INTEGER imod,idat,i,ichang,id0,id1,nid,idatmax1
+      REAL x0,xim0,dx,dxm
+      REAL chmin,chmax,pi
+
+    
+      REAL,ALLOCATABLE :: xxid(:)
+      REAL,ALLOCATABLE :: xxd(:)
+      REAL,ALLOCATABLE :: fdd(:)
+      REAL,ALLOCATABLE :: fxd(:)
+      REAL,ALLOCATABLE :: xchan(:)
+      REAL,ALLOCATABLE :: fdchan(:)
+      REAL,ALLOCATABLE :: xxim(:)
+
+      ALLOCATE( xxid(idatmax+1))
+      ALLOCATE( xxd(idatmax+1))
+      ALLOCATE( fdd(idatmax+1))
+      ALLOCATE( fxd(idatmax+1))
+      ALLOCATE( xchan(idatmax+1))
+      ALLOCATE( fdchan(idatmax+1))
+      ALLOCATE( xxim(imodmax))
+
+      
+      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) = AMOD( 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(1h-))
+2      FORMAT(1x,8f8.2)
+
+       RETURN
+       END
Index: LMDZ.3.3/trunk/libf/dyn3d/inter_barxy.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/inter_barxy.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/inter_barxy.F	(revision 259)
@@ -0,0 +1,66 @@
+C
+C $Header$
+C
+       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(iip1),chpn(iip1),chps(iip1)
+       REAL chhpn,chhps
+       REAL fmody(jjp1)
+       INTEGER jmp1
+    
+       REAL,ALLOCATABLE :: champy(:,:)
+       ALLOCATE(champy(jnterfd+1,iip1))
+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
+
+         DEALLOCATE(champy)
+
+       RETURN
+       END
+
Index: LMDZ.3.3/trunk/libf/dyn3d/inter_bary.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/inter_bary.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/inter_bary.F	(revision 259)
@@ -0,0 +1,144 @@
+C
+C $Header$
+C
+       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  .......
+c
+       INTEGER jjm , jdatmax, jmodmax
+       REAL    yjdatt( 1 ) , fdatt( 1 )
+       REAL    yjmodd( 1 )     
+
+c  ....  Arguments  en sortie  .......
+c
+       REAL    fmod( 1 )
+c
+c   ...... Variables locales  ......
+       LOGICAL    decrois
+       SAVE       decrois
+
+       REAL y0,dy,dym 
+       INTEGER jdat, jmod,i
+
+       INTEGER     jmods
+       REAL , ALLOCATABLE :: fdat(:)
+       REAL , ALLOCATABLE :: yjdat(:)
+       REAL , ALLOCATABLE :: yjmod(:)
+       REAL , ALLOCATABLE :: fscrat(:)
+       ALLOCATE (fdat(jdatmax+1))
+       ALLOCATE (yjdat(jdatmax+1))
+       ALLOCATE (yjmod(jmodmax))
+       ALLOCATE (fscrat(jdatmax+1))
+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
+
+       DEALLOCATE(fdat)
+       DEALLOCATE(yjdat)
+       DEALLOCATE(yjmod)
+       DEALLOCATE(fscrat)
+
+       RETURN
+       END
Index: LMDZ.3.3/trunk/libf/dyn3d/limit_netcdf.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/limit_netcdf.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/limit_netcdf.F	(revision 259)
@@ -0,0 +1,1062 @@
+C
+C $Header$
+C
+      SUBROUTINE limit_netcdf ( interbar, extrap, oldice )
+c
+      IMPLICIT none
+c
+c-------------------------------------------------------------
+C Author : L. Fairhead
+C Date   : 27/01/94
+C Objet  : Construction des fichiers de conditions aux limites
+C          pour le nouveau
+C          modele a partir de fichiers de climatologie. Les deux
+C          grilles doivent etre regulieres
+c
+c Modifie par z.x.li (le23mars1994)
+c Modifie par L. Fairhead (fairhead@lmd.jussieu.fr) septembre 1999
+c                         pour lecture netcdf dans LMDZ.3.3
+c Modifie par P;Le Van  ,  juillet 2001
+c-------------------------------------------------------------
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "control.h"
+#include "logic.h"
+#include "netcdf.inc"
+#include "comvert.h"
+#include "comgeom2.h"
+#include "comconst.h"
+c
+c-----------------------------------------------------------------------
+      LOGICAL interbar, extrap, oldice
+
+      INTEGER KIDIA, KFDIA, KLON, KLEV
+      PARAMETER (KIDIA=1,KFDIA=iim*(jjm-1)+2,
+     .           KLON=KFDIA-KIDIA+1,KLEV=llm)
+c-----------------------------------------------------------------------
+      REAL phy_nat(klon,360), phy_nat0(klon)
+      REAL phy_alb(klon,360)
+      REAL phy_sst(klon,360)
+      REAL phy_bil(klon,360)
+      REAL phy_rug(klon,360)
+      REAL phy_ice(klon,360)
+c
+      REAL masque(iip1,jjp1)
+      REAL mask(iim,jjp1)
+
+C Declarations pour le champ de depart
+      INTEGER imdep, jmdep,lmdep
+      INTEGER  tbid
+      PARAMETER ( tbid = 60 )        ! >52 semaines
+      REAL  timecoord(tbid)
+c
+      REAL , ALLOCATABLE :: dlon_msk(:), dlat_msk(:)
+      REAL , ALLOCATABLE :: lonmsk_ini(:), latmsk_ini(:)
+      REAL , ALLOCATABLE :: dlon(:), dlat(:)
+      REAL , ALLOCATABLE :: dlon_ini(:), dlat_ini(:)
+      REAL , ALLOCATABLE :: champ_msk(:), champ(:)
+      REAL , ALLOCATABLE :: work(:,:)
+
+      CHARACTER*25 title
+
+C Declarations pour le champ interpole 2D
+      REAL champint(iim,jjp1)
+      real chmin,chmax
+
+C Declarations pour le champ interpole 3D
+      REAL champtime(iim,jjp1,tbid)
+      REAL timeyear(tbid)
+      REAL champan(iip1,jjp1,366)
+
+C Declarations pour l'inteprolation verticale
+      REAL ax(tbid), ay(tbid)
+      REAL by
+      REAL yder(tbid)
+
+
+      INTEGER ierr
+      INTEGER dimfirst(3)
+      INTEGER dimlast(3)
+c
+      INTEGER nid, ndim, ntim
+      INTEGER dims(2), debut(2), epais(2)
+      INTEGER id_tim
+      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
+
+      INTEGER i, j, k, l
+c Diverses variables locales
+      REAL time
+
+      INTEGER          longcles
+      PARAMETER      ( longcles = 20 )
+      REAL  clesphy0 ( longcles      )
+#include "serre.h"
+      INTEGER ncid,varid,ndimid(4),dimid
+      character*30 namedim
+
+      pi     = 4. * ATAN(1.)
+      rad    = 6 371 229.
+      omeg   = 4.* ASIN(1.)/(24.*3600.)
+      g      = 9.8
+      daysec = 86400.
+      kappa  = 0.2857143
+      cpp    = 1004.70885
+      dtvr    = daysec/FLOAT(day_step)
+c
+C Traitement du relief au sol
+c
+      write(*,*) 'Traitement du relief au sol pour fabriquer masque'
+      ierr = NF_OPEN('Relief.nc', NF_NOWRITE, ncid)
+
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ierr = NF_INQ_VARID(ncid,'RELIEF',varid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      print*,'variable ', namedim, 'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ALLOCATE( lonmsk_ini(imdep) )
+      ALLOCATE(   dlon_msk(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,lonmsk_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,lonmsk_ini)
+#endif
+
+c
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      print*,'variable ', namedim, 'dimension ', jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ALLOCATE( latmsk_ini(jmdep) )
+      ALLOCATE(   dlat_msk(jmdep) )
+      ALLOCATE(  champ_msk(imdep*jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,latmsk_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,latmsk_ini)
+#endif
+c
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,varid,champ_msk)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk)
+#endif
+c
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+c
+      title='RELIEF'
+
+      CALL conf_dat2d(title,imdep, jmdep, lonmsk_ini, latmsk_ini,
+     . dlon_msk, dlat_msk, champ_msk, interbar  )
+
+      CALL mask_c_o(imdep, jmdep, dlon_msk, dlat_msk,champ_msk,
+     .             iim, jjp1, rlonv, rlatu, champint)
+      CALL gr_int_dyn(champint, masque, iim, jjp1)
+      DO i = 1, iim
+         masque(i,1) = FLOAT(NINT(masque(i,1)))
+         masque(i,jjp1) = FLOAT(NINT(masque(i,jjp1)))
+      ENDDO
+      DO i = 1, iim
+      DO j = 1, jjp1
+         mask(i,j) = champint(i,j)
+      ENDDO
+      ENDDO
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0)
+      ierr = NF_CLOSE(ncid)
+c
+c
+C Traitement de la rugosite
+c
+      PRINT*, 'Traitement de la rugosite'
+      ierr = NF_OPEN('Rugos.nc', NF_NOWRITE, ncid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ierr = NF_INQ_VARID(ncid,'RUGOS',varid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE( dlon_ini(imdep) )
+      ALLOCATE(     dlon(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE( dlat_ini(jmdep) )
+      ALLOCATE(     dlat(jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      print*,'variable ', namedim, 'dimension ', lmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+c
+      ALLOCATE( champ(imdep*jmdep) )
+
+      DO  200 l = 1, lmdep
+         dimfirst(1) = 1
+         dimfirst(2) = 1
+         dimfirst(3) = l
+c
+         dimlast(1) = imdep
+         dimlast(2) = jmdep
+         dimlast(3) = 1
+c
+         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+         print*,dimfirst,dimlast
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
+         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
+         if (ierr.ne.0) then
+           print *, NF_STRERROR(ierr)
+           STOP
+         ENDIF 
+   
+        title = 'Rugosite Amip '
+c
+        CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
+     .                      dlon, dlat, champ, interbar          )
+
+       IF ( interbar )   THEN
+         DO j = 1, imdep * jmdep
+           champ(j) = LOG(champ(j))
+         ENDDO
+
+        IF( l.EQ.1 )  THEN
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour la rugosite $$$ '
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+        ENDIF
+        CALL inter_barxy ( imdep,jmdep -1,dlon,dlat,champ ,
+     ,                  iim,jjm,rlonu,rlatv, jjp1,champint )
+         DO j=1,jjp1
+          DO i=1,iim
+           champint(i,j)=EXP(champint(i,j))
+          ENDDO
+         ENDDO
+
+         DO j = 1, jjp1
+           DO i = 1, iim
+             IF(NINT(mask(i,j)).NE.1)  THEN
+               champint( i,j ) = 0.001
+             ENDIF
+           ENDDO
+         ENDDO
+      ELSE
+         CALL rugosite(imdep, jmdep, dlon, dlat, champ,
+     .             iim, jjp1, rlonv, rlatu, champint, mask)
+      ENDIF
+         DO j = 1,jjp1
+         DO i = 1, iim
+            champtime (i,j,l) = champint(i,j)
+         ENDDO
+         ENDDO
+200      CONTINUE
+c
+      DO l = 1, lmdep
+         timeyear(l) = timecoord(l)
+      ENDDO
+
+      PRINT 222, timeyear
+222   FORMAT(2x,' Time year ',10f6.1)
+c
+        
+      PRINT*, 'Interpolation temporelle dans l annee'
+
+      DO j = 1, jjp1
+      DO i = 1, iim
+          DO l = 1, lmdep
+            ax(l) = timeyear(l)
+            ay(l) = champtime (i,j,l)
+          ENDDO
+          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
+          DO k = 1, 360
+            time = FLOAT(k-1)
+            CALL SPLINT(ax,ay,yder,lmdep,time,by)
+            champan(i,j,k) = by
+          ENDDO
+      ENDDO
+      ENDDO
+      DO k = 1, 360
+      DO j = 1, jjp1
+         champan(iip1,j,k) = champan(1,j,k)
+      ENDDO
+        IF ( k.EQ.10 )  THEN
+          DO j = 1, jjp1
+            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
+            PRINT *,' Rugosite au temps 10 ', chmin,chmax,j
+          ENDDO
+        ENDIF
+      ENDDO
+c
+      DO k = 1, 360
+         CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k), phy_rug(1,k))
+      ENDDO
+c
+      ierr = NF_CLOSE(ncid)
+
+       DEALLOCATE( dlon      )
+       DEALLOCATE( dlon_ini  )
+       DEALLOCATE( dlat      )
+       DEALLOCATE( dlat_ini  )
+       DEALLOCATE( champ     )
+c
+c
+C Traitement de la glace oceanique
+c
+      PRINT*, 'Traitement de la glace oceanique'
+      ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE ( dlon_ini(imdep) )
+      ALLOCATE (     dlon(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE ( dlat_ini(jmdep) )
+      ALLOCATE (     dlat(jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, lmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+c
+      ALLOCATE ( champ(imdep*jmdep) )
+
+      DO l = 1, lmdep
+         dimfirst(1) = 1
+         dimfirst(2) = 1
+         dimfirst(3) = l
+c
+         dimlast(1) = imdep
+         dimlast(2) = jmdep
+         dimlast(3) = 1
+c
+         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
+         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
+         if (ierr.ne.0) then
+           print *, NF_STRERROR(ierr)
+           STOP
+         ENDIF
+ 
+         title = 'Sea-ice Amip '
+c
+         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
+     .                        dlon, dlat, champ, interbar          )
+c
+      IF( oldice )  THEN
+                 CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
+     .             iim, jjp1, rlonv, rlatu, champint )
+      ELSEIF ( interbar )  THEN
+       IF( l.EQ.1 )  THEN
+        WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour Sea-ice Amip  $$$ '
+        WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+       ENDIF
+
+         CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
+     ,     champ, iim, jjm, rlonu, rlatv, jjp1, champint )
+      ELSE
+         CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
+     .             iim, jjp1, rlonv, rlatu, champint )
+      ENDIF
+         DO j = 1,jjp1
+         DO i = 1, iim
+            champtime (i,j,l) = champint(i,j)
+         ENDDO
+         ENDDO
+      ENDDO
+c
+      DO l = 1, lmdep
+         timeyear(l) = timecoord(l)
+      ENDDO
+      PRINT 222,  timeyear
+c
+      PRINT*, 'Interpolation temporelle'
+      DO j = 1, jjp1
+      DO i = 1, iim
+          DO l = 1, lmdep
+            ax(l) = timeyear(l)
+            ay(l) = champtime (i,j,l)
+          ENDDO
+          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
+          DO k = 1, 360
+            time = FLOAT(k-1)
+            CALL SPLINT(ax,ay,yder,lmdep,time,by)
+            champan(i,j,k) = by
+          ENDDO
+      ENDDO
+      ENDDO
+      DO k = 1, 360
+      DO j = 1, jjp1
+         champan(iip1, j, k) = champan(1, j, k)
+      ENDDO
+        IF ( k.EQ.10 )  THEN
+          DO j = 1, jjp1
+            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
+            PRINT *,' Sea ice au temps 10 ', chmin,chmax,j
+          ENDDO
+        ENDIF
+      ENDDO
+c
+
+      DO k = 1, 360
+         CALL gr_dyn_fi(1, iip1, jjp1, klon,
+     .                  champan(1,1,k), phy_ice(1,k))
+         DO i = 1, klon
+            phy_nat(i,k) = phy_nat0(i)
+            IF ( (phy_ice(i,k) - 0.5).GE.1.e-5 ) THEN
+               IF (NINT(phy_nat0(i)).EQ.0) THEN
+                  phy_nat(i,k) = 3.0
+               ELSE
+                  phy_nat(i,k) = 2.0
+               ENDIF
+            ENDIF
+            IF( NINT(phy_nat(i,k)).EQ.0 ) THEN
+              IF ( phy_rug(i,k).NE.0.001 ) phy_rug(i,k) = 0.001
+            ENDIF
+         ENDDO
+
+      ENDDO
+c
+
+      ierr = NF_CLOSE(ncid)
+c
+       DEALLOCATE( dlon      )
+       DEALLOCATE( dlon_ini  )
+       DEALLOCATE( dlat      )
+       DEALLOCATE( dlat_ini  )
+       DEALLOCATE( champ     )
+
+477    continue
+c
+C Traitement de la sst
+c
+      PRINT*, 'Traitement de la sst'
+      ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ierr = NF_INQ_VARID(ncid,'SST',varid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable SST ', namedim,'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+ 
+      ALLOCATE( dlon_ini(imdep) )
+      ALLOCATE(     dlon(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
+#endif
+
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable SST ', namedim, 'dimension ', jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE( dlat_ini(jmdep) )
+      ALLOCATE(     dlat(jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', lmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+       ALLOCATE( champ(imdep*jmdep) )
+       IF( extrap )   THEN
+         ALLOCATE ( work(imdep,jmdep) )
+       ENDIF
+c
+      DO l = 1, lmdep
+         dimfirst(1) = 1
+         dimfirst(2) = 1
+         dimfirst(3) = l
+c
+         dimlast(1) = imdep
+         dimlast(2) = jmdep
+         dimlast(3) = 1
+c
+         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
+         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
+         if (ierr.ne.0) then
+           print *, NF_STRERROR(ierr)
+           STOP
+         ENDIF
+
+         title='Sst Amip'
+c
+         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
+     .                            dlon, dlat, champ, interbar     )
+       IF ( extrap )  THEN
+        CALL extrapol(champ, imdep, jmdep, 999999.,.TRUE.,.TRUE.,2,work)
+       ENDIF
+c
+
+      IF ( interbar )  THEN
+        IF( l.EQ.1 )  THEN
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour la Sst Amip $$$ '
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+        ENDIF
+       CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
+     , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
+      ELSE
+       CALL grille_m(imdep, jmdep, dlon, dlat, champ,
+     .          iim, jjp1, rlonv, rlatu, champint   )
+      ENDIF
+
+         DO j = 1,jjp1
+         DO i = 1, iim
+            champtime (i,j,l) = champint(i,j)
+         ENDDO
+         ENDDO
+      ENDDO
+c
+      DO l = 1, lmdep
+         timeyear(l) = timecoord(l)
+      ENDDO
+      print 222,  timeyear
+c
+C interpolation temporelle
+      DO j = 1, jjp1
+      DO i = 1, iim
+          DO l = 1, lmdep
+            ax(l) = timeyear(l)
+            ay(l) = champtime (i,j,l)
+          ENDDO
+          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
+          DO k = 1, 360
+            time = FLOAT(k-1)
+            CALL SPLINT(ax,ay,yder,lmdep,time,by)
+            champan(i,j,k) = by
+          ENDDO
+      ENDDO
+      ENDDO
+      DO k = 1, 360
+      DO j = 1, jjp1
+         champan(iip1,j,k) = champan(1,j,k)
+      ENDDO
+        IF ( k.EQ.10 )  THEN
+          DO j = 1, jjp1
+            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
+            PRINT *,' SST au temps 10 ', chmin,chmax,j
+          ENDDO
+        ENDIF
+      ENDDO
+c
+      DO k = 1, 360
+         CALL gr_dyn_fi(1, iip1, jjp1, klon, 
+     .                  champan(1,1,k), phy_sst(1,k))
+      ENDDO
+c
+      ierr = NF_CLOSE(ncid)
+c
+       DEALLOCATE( dlon      )
+       DEALLOCATE( dlon_ini  )
+       DEALLOCATE( dlat      )
+       DEALLOCATE( dlat_ini  )
+       DEALLOCATE( champ     )
+c
+C Traitement de l'albedo
+c
+      PRINT*, 'Traitement de l albedo'
+      ierr = NF_OPEN('Albedo.nc', NF_NOWRITE, ncid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARID(ncid,'ALBEDO',varid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE ( dlon_ini(imdep) )
+      ALLOCATE (     dlon(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE ( dlat_ini(jmdep) )
+      ALLOCATE (     dlat(jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', lmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+c
+      ALLOCATE ( champ(imdep*jmdep) )
+
+      DO l = 1, lmdep
+         dimfirst(1) = 1
+         dimfirst(2) = 1
+         dimfirst(3) = l
+c
+         dimlast(1) = imdep
+         dimlast(2) = jmdep
+         dimlast(3) = 1
+c
+         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
+         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
+         if (ierr.ne.0) then
+           print *, NF_STRERROR(ierr)
+           STOP
+         ENDIF
+
+         title='Albedo Amip'
+c
+         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
+     .                            dlon, dlat, champ, interbar      )
+c
+c
+      IF ( interbar )  THEN
+        IF( l.EQ.1 )  THEN
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour l Albedo Amip $$$ '
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+        ENDIF
+
+       CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
+     , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
+      ELSE
+       CALL grille_m(imdep, jmdep, dlon, dlat, champ,
+     .          iim, jjp1, rlonv, rlatu, champint   )
+      ENDIF
+c
+         DO j = 1,jjp1
+         DO i = 1, iim
+            champtime (i, j, l) = champint(i, j)
+         ENDDO
+         ENDDO
+      ENDDO
+c
+      DO l = 1, lmdep
+         timeyear(l) = timecoord(l)
+      ENDDO
+      print 222,  timeyear
+c
+C interpolation temporelle
+      DO j = 1, jjp1
+      DO i = 1, iim
+          DO l = 1, lmdep
+            ax(l) = timeyear(l)
+            ay(l) = champtime (i, j, l)
+          ENDDO
+          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
+          DO k = 1, 360
+            time = FLOAT(k-1)
+            CALL SPLINT(ax,ay,yder,lmdep,time,by)
+            champan(i,j,k) = by
+          ENDDO
+      ENDDO
+      ENDDO
+      DO k = 1, 360
+      DO j = 1, jjp1
+         champan(iip1, j, k) = champan(1, j, k)
+      ENDDO
+        IF ( k.EQ.10 )  THEN
+          DO j = 1, jjp1
+            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
+            PRINT *,' Albedo au temps 10 ', chmin,chmax,j
+          ENDDO
+        ENDIF
+      ENDDO
+c
+      DO k = 1, 360
+         CALL gr_dyn_fi(1, iip1, jjp1, klon,
+     .                  champan(1,1,k), phy_alb(1,k))
+      ENDDO
+c
+      ierr = NF_CLOSE(ncid)
+c
+c
+      DO k = 1, 360
+      DO i = 1, klon
+         phy_bil(i,k) = 0.0
+      ENDDO
+      ENDDO
+c
+      PRINT*, 'Ecriture du fichier limit'
+c
+      ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
+c
+      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
+     .                       "Fichier conditions aux limites")
+      ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
+      ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
+c
+      dims(1) = ndim
+      dims(2) = ntim
+c
+ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
+      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
+      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
+     .                        "Jour dans l annee")
+ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
+      ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
+      ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
+     .                        "Nature du sol (0,1,2,3)")
+ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
+      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
+      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
+     .                        "Temperature superficielle de la mer")
+ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
+      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
+      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
+     .                        "Reference flux de chaleur au sol")
+ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
+      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
+      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
+     .                        "Albedo a la surface")
+ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
+      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
+      ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
+     .                        "Rugosite")
+c
+      ierr = NF_ENDDEF(nid)
+c
+      DO k = 1, 360
+c
+      debut(1) = 1
+      debut(2) = k
+      epais(1) = klon
+      epais(2) = 1
+c
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
+#else
+      ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
+      ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
+#endif
+c
+      ENDDO
+c
+      ierr = NF_CLOSE(nid)
+c
+      STOP
+      END
Index: LMDZ.3.3/trunk/libf/dyn3d/ord_coord.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/ord_coord.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/ord_coord.F	(revision 259)
@@ -0,0 +1,102 @@
+C
+C $Header$
+C
+       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(1)
+       INTEGER nscr
+       PARAMETER ( nscr = 4000 )
+       REAL xscr(nscr)
+
+c    .....  Arguments  en sortie  .....
+c
+       REAL xo(1) 
+       LOGICAL decrois
+
+c    .... Variables  locales  ....
+
+       INTEGER i,ii
+       REAL pi, degres, chmin, chmax, mult
+c
+
+       IF( nmax.GT.nscr )   THEN
+         PRINT *,' Dans ord_coord  , nmax = ',nmax, 'est plus grand que'
+     ,, ' nscr !  Mettre nscr a la valeur de nmax et repasser ! '  
+         CALL ABORT
+       ENDIF
+
+       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
+         CALL sort(nmax+1,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: LMDZ.3.3/trunk/libf/dyn3d/ord_coordm.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/ord_coordm.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/ord_coordm.F	(revision 259)
@@ -0,0 +1,118 @@
+C
+C $Header$
+C
+       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
+       REAL xi(1)
+       INTEGER nscr
+       PARAMETER ( nscr = 4000 )
+       REAL xscr(nscr)
+
+c    .....  Arguments  en sortie  .....
+c
+       REAL xo(1) 
+       LOGICAL decrois
+       INTEGER jmods
+
+c    .... Variables  locales  ....
+
+       INTEGER i,jjm
+       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 
+     
+       IF( nmax.GT.nscr )   THEN
+         PRINT *,' Dans ord_coord  , nmax = ',nmax, 'est plus grand que'
+     ,, ' nscr !  Mettre nscr a la valeur de nmax et repasser ! '  
+         CALL ABORT
+       ENDIF
+
+       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: LMDZ.3.3/trunk/libf/dyn3d/sort.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/sort.F	(revision 259)
+++ LMDZ.3.3/trunk/libf/dyn3d/sort.F	(revision 259)
@@ -0,0 +1,35 @@
+C
+C $Header$
+C
+      SUBROUTINE sort(n,d)
+c
+c     P.Le Van
+c      
+c...  cette routine met le tableau d  dans l'ordre croissant  ....
+cc   ( pour avoir l'ordre decroissant,il suffit de remplacer l'instruc
+c      tion  situee + bas  IF(d(j).LE.p)  THEN     par
+c                           IF(d(j).GE.p)  THEN
+c
+
+      INTEGER n
+      REAL d(1) , p
+      INTEGER i,j,k
+
+      DO i=1,n-1
+        k=i
+        p=d(i)
+        DO j=i+1,n
+         IF(d(j).LE.p) THEN
+           k=j
+           p=d(j)
+         ENDIF
+        ENDDO
+
+       IF(k.ne.i) THEN
+         d(k)=d(i)
+         d(i)=p
+       ENDIF
+      ENDDO
+
+       RETURN
+       END
Index: LMDZ.3.3/trunk/libf/dyn3d/startvar.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/startvar.F	(revision 256)
+++ LMDZ.3.3/trunk/libf/dyn3d/startvar.F	(revision 259)
@@ -1,2 +1,5 @@
+C
+C $Header$
+C
       MODULE startvar
     !
@@ -10,13 +13,12 @@
     !
     !        - A 2D variable on the dynamical grid :
-    !           CALL startget(varname, iml, jml, lon_in, lat_in, champ, val_ex)
-    !            
+    !           CALL startget(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)
+    !            CALL startget(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)
+    !            CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
     !
     !
@@ -57,5 +59,5 @@
       REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lat_phys, lat_rug,
      . lat_alb, lat_rel, lat_dyn
-      REAL, ALLOCATABLE, SAVE, DIMENSION (:)  :: lev_dyn
+      REAL, ALLOCATABLE, SAVE, DIMENSION (:)  :: levdyn_ini
       REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: relief, zstd, zsig,
      . zgam, zthe, zpic, zval
@@ -70,5 +72,5 @@
     !
       SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in,
-     . champ, val_exp)
+     . champ, val_exp, jml2, lon_in2, lat_in2 , interbar )
     !
     !    There is a big mess with the size in logitude, should it be iml or iml+1.
@@ -79,8 +81,10 @@
     !
       CHARACTER*(*), INTENT(in) :: varname
-      INTEGER, INTENT(in) :: iml, jml
+      INTEGER, INTENT(in) :: iml, jml ,jml2
       REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
+      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
       REAL, INTENT(inout) :: champ(iml,jml)
       REAL, INTENT(in) :: val_exp
+      LOGICAL interbar
     !
     !   This routine only works if the variable does not exist or is constant
@@ -97,5 +101,6 @@
                   IF ( .NOT.ALLOCATED(relief)) THEN
                       !
-                      CALL start_init_orog( iml, jml, lon_in, lat_in)
+                      CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .                    jml2,lon_in2,lat_in2, interbar )
                       !
                   ENDIF
@@ -117,5 +122,6 @@
                   IF ( .NOT.ALLOCATED(rugo)) THEN
                       !
-                      CALL start_init_orog( iml, jml, lon_in, lat_in)
+                      CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .                    jml2,lon_in2,lat_in2 , interbar )
                       !
                   ENDIF
@@ -137,5 +143,6 @@
                   IF ( .NOT.ALLOCATED(masque)) THEN
                       !
-                      CALL start_init_orog( iml, jml, lon_in, lat_in)
+                      CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .                     jml2,lon_in2,lat_in2 , interbar )
                       !
                   ENDIF
@@ -157,5 +164,6 @@
                   IF ( .NOT.ALLOCATED(phis)) THEN
                       !
-                      CALL start_init_orog( iml, jml, lon_in, lat_in)
+                      CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .                   jml2,lon_in2, lat_in2 , interbar )
                       !
                   ENDIF
@@ -177,5 +185,6 @@
                   IF ( .NOT.ALLOCATED(psol_dyn)) THEN
                       !
-                      CALL start_init_dyn( iml, jml, lon_in, lat_in)
+                      CALL start_init_dyn( iml, jml, lon_in, lat_in,
+     .                   jml2,lon_in2, lat_in2 , interbar )
                       !
                   ENDIF
@@ -231,18 +240,24 @@
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !
-      SUBROUTINE start_init_orog( iml, jml, lon_in, lat_in)
-    !
-      INTEGER, INTENT(in) :: iml, jml
+      SUBROUTINE start_init_orog ( iml,jml,lon_in, lat_in,jml2,lon_in2 ,
+     ,   lat_in2 , interbar )
+    !
+      INTEGER, INTENT(in) :: iml, jml, jml2
       REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
+      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
+      LOGICAL interbar
     !
     !  LOCAL
     !
-      REAL :: lev(1), date, dt
+      LOGICAL interbar2
+      REAL :: lev(1), date, dt,chmin,chmax
       INTEGER :: itau(1), fid
       INTEGER ::  llm_tmp, ttm_tmp
       INTEGER :: i, j
       INTEGER :: iret
+      CHARACTER*25 title
       REAL, ALLOCATABLE :: relief_hi(:,:)
       REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
+      REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:)
       REAL, ALLOCATABLE :: tmp_var(:,:)
       INTEGER, ALLOCATABLE :: tmp_int(:,:)
@@ -274,17 +289,29 @@
     !
       ALLOCATE(lon_rad(iml_rel))
+      ALLOCATE(lon_ini(iml_rel))
+
       IF ( MAXVAL(lon_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
-          lon_rad(:) = lon_rel(:,1) * 2.0 * ASIN(1.0) / 180.0
-      ELSE
-          lon_rad(:) = lon_rel(:,1) 
-      ENDIF
+          lon_ini(:) = lon_rel(:,1) * 2.0 * ASIN(1.0) / 180.0
+      ELSE
+          lon_ini(:) = lon_rel(:,1) 
+      ENDIF
+
       ALLOCATE(lat_rad(jml_rel))
+      ALLOCATE(lat_ini(jml_rel))
+
       IF ( MAXVAL(lat_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
-          lat_rad(:) = lat_rel(1,:) * 2.0 * ASIN(1.0) / 180.0
-      ELSE
-          lat_rad(:) = lat_rel(1,:) 
-      ENDIF
-    !
-    !
+          lat_ini(:) = lat_rel(1,:) * 2.0 * ASIN(1.0) / 180.0
+      ELSE
+          lat_ini(:) = lat_rel(1,:) 
+      ENDIF
+    !
+    !
+
+      title='RELIEF'
+
+      interbar2 = .FALSE.
+      CALL conf_dat2d(title,iml_rel, jml_rel, lon_ini, lat_ini,
+     . lon_rad, lat_rad, relief_hi , interbar2  )
+
       IF ( check ) WRITE(*,*) 'Computes all the parameters needed',
      .' for the gravity wave drag code'
@@ -337,4 +364,6 @@
         rugo(iml,j) = tmp_var(1,j)
       ENDDO
+c
+cc   ***   rugo  n'est pas utilise pour l'instant  ******
     !
     !   Build land-sea mask
@@ -348,11 +377,13 @@
     !
       SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, 
-     .lat_in, nbindex, champ, val_exp)
+     .lat_in, nbindex, champ, val_exp ,jml2, lon_in2, lat_in2,interbar)
     !
       CHARACTER*(*), INTENT(in) :: varname
-      INTEGER, INTENT(in) :: iml, jml, nbindex
+      INTEGER, INTENT(in) :: iml, jml, nbindex, jml2
       REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
+      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
       REAL, INTENT(inout) :: champ(nbindex)
       REAL, INTENT(in) :: val_exp
+      LOGICAL interbar
     !
     !
@@ -364,5 +395,6 @@
             CASE ('tsol')
               IF ( .NOT.ALLOCATED(tsol)) THEN
-                CALL start_init_phys( iml, jml, lon_in, lat_in)
+                CALL start_init_phys( iml, jml, lon_in, lat_in,
+     .              jml2, lon_in2, lat_in2, interbar )
               ENDIF
               IF ( SIZE(tsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
@@ -374,5 +406,6 @@
             CASE ('qsol')
               IF ( .NOT.ALLOCATED(qsol)) THEN
-                CALL start_init_phys( iml, jml, lon_in, lat_in)
+                CALL start_init_phys( iml, jml, lon_in, lat_in,
+     .              jml2, lon_in2,lat_in2 , interbar )
               ENDIF
               IF ( SIZE(qsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
@@ -384,5 +417,6 @@
             CASE ('psol')
               IF ( .NOT.ALLOCATED(psol_dyn)) THEN
-                CALL start_init_dyn( iml, jml, lon_in, lat_in)
+                CALL start_init_dyn( iml, jml, lon_in, lat_in,
+     .              jml2, lon_in2,lat_in2 , interbar )
               ENDIF
               IF (SIZE(psol_dyn) .NE. SIZE(lon_in)*SIZE(lat_in)) THEN
@@ -394,5 +428,6 @@
             CASE ('zmea')
               IF ( .NOT.ALLOCATED(relief)) THEN
-                CALL start_init_orog( iml, jml, lon_in, lat_in)
+                CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .            jml2, lon_in2,lat_in2 , interbar )
               ENDIF
               IF ( SIZE(relief) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
@@ -404,5 +439,6 @@
             CASE ('zstd')
               IF ( .NOT.ALLOCATED(zstd)) THEN
-                CALL start_init_orog( iml, jml, lon_in, lat_in)
+                CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .              jml2, lon_in2,lat_in2 , interbar )
               ENDIF
               IF ( SIZE(zstd) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
@@ -414,5 +450,6 @@
             CASE ('zsig')
               IF ( .NOT.ALLOCATED(zsig)) THEN
-                CALL start_init_orog( iml, jml, lon_in, lat_in)
+                CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .               jml2, lon_in2,lat_in2 , interbar )
               ENDIF
               IF ( SIZE(zsig) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
@@ -424,5 +461,6 @@
             CASE ('zgam')
               IF ( .NOT.ALLOCATED(zgam)) THEN
-                CALL start_init_orog( iml, jml, lon_in, lat_in)
+                CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .            jml2, lon_in2,lat_in2 , interbar )
               ENDIF
               IF ( SIZE(zgam) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
@@ -434,5 +472,6 @@
             CASE ('zthe')
               IF ( .NOT.ALLOCATED(zthe)) THEN
-                CALL start_init_orog( iml, jml, lon_in, lat_in)
+                CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .            jml2, lon_in2,lat_in2 , interbar )
               ENDIF
               IF ( SIZE(zthe) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
@@ -444,5 +483,6 @@
             CASE ('zpic')
               IF ( .NOT.ALLOCATED(zpic)) THEN
-                CALL start_init_orog( iml, jml, lon_in, lat_in)
+                CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .            jml2, lon_in2,lat_in2 , interbar )
               ENDIF
               IF ( SIZE(zpic) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
@@ -454,5 +494,6 @@
             CASE ('zval')
               IF ( .NOT.ALLOCATED(zval)) THEN
-                CALL start_init_orog( iml, jml, lon_in, lat_in)
+                CALL start_init_orog( iml, jml, lon_in, lat_in,
+     .            jml2, lon_in2,lat_in2 , interbar )
               ENDIF
               IF ( SIZE(zval) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
@@ -462,5 +503,5 @@
               ENDIF
               CALL gr_dyn_fi(1, iml, jml, nbindex,zval, champ)
-             CASE ('rads')
+            CASE ('rads')
                   champ(:) = 0.0
             CASE ('snow')
@@ -493,8 +534,11 @@
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !
-      SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in)
-    !
-      INTEGER, INTENT(in) :: iml, jml
+      SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in, jml2,
+     .                 lon_in2, lat_in2 , interbar )
+    !
+      INTEGER, INTENT(in) :: iml, jml ,jml2
       REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
+      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
+      LOGICAL interbar
     !
     !  LOCAL
@@ -505,8 +549,10 @@
       INTEGER :: i, j
     !
+      CHARACTER*25 title
       CHARACTER*120 :: physfname
       LOGICAL :: check=.TRUE.
     !
       REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
+      REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:)
       REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:)
     !
@@ -533,15 +579,22 @@
     !
       ALLOCATE(lon_rad(iml_phys))
+      ALLOCATE(lon_ini(iml_phys))
+
       IF ( MAXVAL(lon_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
-          lon_rad(:) = lon_phys(:,1) * 2.0 * ASIN(1.0) / 180.0
-      ELSE
-          lon_rad(:) = lon_phys(:,1) 
-      ENDIF
+          lon_ini(:) = lon_phys(:,1) * 2.0 * ASIN(1.0) / 180.0
+      ELSE
+          lon_ini(:) = lon_phys(:,1) 
+      ENDIF
+
       ALLOCATE(lat_rad(jml_phys))
+      ALLOCATE(lat_ini(jml_phys))
+
       IF ( MAXVAL(lat_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
-          lat_rad(:) = lat_phys(1,:) * 2.0 * ASIN(1.0) / 180.0
-      ELSE
-          lat_rad(:) = lat_phys(1,:) 
-      ENDIF
+          lat_ini(:) = lat_phys(1,:) * 2.0 * ASIN(1.0) / 180.0
+      ELSE
+          lat_ini(:) = lat_phys(1,:) 
+      ENDIF
+
+
     !
     !   We get the two standard varibales
@@ -552,8 +605,26 @@
     !
     !
+
       CALL flinget(fid_phys, 'ST', iml_phys, jml_phys, 
      .llm_tmp, ttm_tmp, 1, 1, var_ana)
-      CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
-     . var_ana, iml-1, jml, lon_in, lat_in, tmp_var)
+
+      title='ST'
+      CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini,
+     . lon_rad, lat_rad, var_ana , interbar  )
+
+      IF ( interbar )   THEN
+        WRITE(6,*) '-------------------------------------------------',
+     ,'--------------'
+        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour  ST $$$ '
+        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   ) 
+      ELSE
+        CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
+     .    var_ana, iml-1, jml, lon_in, lat_in, tmp_var     )
+      ENDIF
+
       CALL gr_int_dyn(tmp_var, tsol, iml-1, jml)
     !
@@ -563,9 +634,26 @@
       CALL flinget(fid_phys, 'CDSW', iml_phys, jml_phys,
      . llm_tmp, ttm_tmp, 1, 1, var_ana)
-      CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
-     . var_ana, iml-1, jml, lon_in, lat_in, tmp_var)
-      CALL gr_int_dyn(tmp_var, qsol, iml-1, jml)
-    !
-      CALL flinclo(fid_phys)
+
+      title='CDSW'
+      CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini,
+     . lon_rad, lat_rad, var_ana, interbar  )
+
+      IF ( interbar )   THEN
+        WRITE(6,*) '-------------------------------------------------',
+     ,'--------------'
+        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour  CDSW $$$ '
+        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  )
+      ELSE
+        CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
+     .    var_ana, iml-1, jml, lon_in, lat_in, tmp_var     )
+      ENDIF
+c
+        CALL gr_int_dyn(tmp_var, qsol, iml-1, jml)
+    !
+       CALL flinclo(fid_phys)
     !
       END SUBROUTINE start_init_phys
@@ -576,15 +664,18 @@
     !
       SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in,
-     . lml, pls, workvar, champ, val_exp)
+     . lml, pls, workvar, champ, val_exp,jml2, lon_in2, lat_in2 ,
+     ,  interbar )
     !
     !   ARGUMENTS
     !
       CHARACTER*(*), INTENT(in) :: varname
-      INTEGER, INTENT(in) :: iml, jml, lml
+      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)
       REAL, INTENT(in) :: val_exp
+      LOGICAL interbar
     !
     !    LOCAL
@@ -592,7 +683,4 @@
       INTEGER :: il, ij, ii
       REAL :: xppn, xpps
-    !
-    !   C'est vraiment une galere de devoir rajouter tant de commons just pour avoir les aires.
-    !   Il faudrait mettre une structure plus flexible et moins dangereuse.
     !
 #include "dimensions.h"
@@ -609,8 +697,9 @@
           CASE ('u')
             IF ( .NOT.ALLOCATED(psol_dyn)) THEN
-              CALL start_init_dyn( iml, jml, lon_in, lat_in)
+              CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 ,
+     .          lon_in2,lat_in2 , interbar )
             ENDIF
             CALL start_inter_3d('U', iml, jml, lml, lon_in,
-     .                           lat_in, pls, champ)
+     .       lat_in, jml2, lon_in2, lat_in2,  pls, champ,interbar )
             DO il=1,lml
               DO ij=1,jml
@@ -623,8 +712,9 @@
           CASE ('v')
             IF ( .NOT.ALLOCATED(psol_dyn)) THEN
-              CALL start_init_dyn( iml, jml, lon_in, lat_in)
+              CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2, 
+     .           lon_in2, lat_in2 , interbar )
             ENDIF
-            CALL start_inter_3d('V', iml, jml, lml, lon_in,
-     .                          lat_in, pls, champ)
+            CALL start_inter_3d('V', iml, jml, lml, lon_in, 
+     .       lat_in, jml2, lon_in2, lat_in2,  pls, champ, interbar )
             DO il=1,lml
               DO ij=1,jml
@@ -637,15 +727,17 @@
           CASE ('t')
             IF ( .NOT.ALLOCATED(psol_dyn)) THEN
-              CALL start_init_dyn( iml, jml, lon_in, lat_in)
+              CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 ,
+     .           lon_in2, lat_in2 ,interbar )
             ENDIF
             CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
-     .                           lat_in, pls, champ)
+     .       lat_in, jml2, lon_in2, lat_in2,  pls, champ, interbar )
  
           CASE ('tpot')
             IF ( .NOT.ALLOCATED(psol_dyn)) THEN
-              CALL start_init_dyn( iml, jml, lon_in, lat_in)
+              CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2 ,
+     .            lon_in2, lat_in2 , interbar )
             ENDIF
             CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
-     .                           lat_in, pls, champ)
+     .       lat_in, jml2, lon_in2, lat_in2,  pls, champ, interbar )
             IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
      .                                    THEN
@@ -672,8 +764,9 @@
           CASE ('q')
             IF ( .NOT.ALLOCATED(psol_dyn)) THEN
-              CALL start_init_dyn( iml, jml, lon_in, lat_in)
+              CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 ,
+     .           lon_in2, lat_in2 , interbar )
             ENDIF
             CALL start_inter_3d('R', iml, jml, lml, lon_in, lat_in,
-     .                           pls, champ)
+     .        jml2, lon_in2, lat_in2,  pls, champ, interbar )
             IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) ) 
      .                                     THEN
@@ -710,8 +803,11 @@
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !
-      SUBROUTINE start_init_dyn( iml, jml, lon_in, lat_in)
-    !
-      INTEGER, INTENT(in) :: iml, jml
+      SUBROUTINE start_init_dyn( iml, jml, lon_in, lat_in,jml2,lon_in2 ,
+     ,             lat_in2 , interbar )
+    !
+      INTEGER, INTENT(in) :: iml, jml, jml2
       REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
+      REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
+      LOGICAL interbar
     !
     !  LOCAL
@@ -726,13 +822,16 @@
     !
       REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
+      REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:)
       REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:), z(:,:)
       REAL, ALLOCATABLE :: xppn(:), xpps(:)
       LOGICAL :: allo
     !
-    !   Ce n'est pas tres pratique d'avoir a charger 3 include pour avoir la grille du modele
     !
 #include "dimensions.h"
 #include "paramet.h"
 #include "comgeom2.h"
+
+      CHARACTER*25 title
+
     !
       physfname = 'ECDYN.nc'
@@ -747,8 +846,8 @@
       ALLOCATE (lat_dyn(iml_dyn,jml_dyn), stat=iret)
       ALLOCATE (lon_dyn(iml_dyn,jml_dyn), stat=iret)
-      ALLOCATE (lev_dyn(llm_dyn), stat=iret)
+      ALLOCATE (levdyn_ini(llm_dyn), stat=iret)
     !
       CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,
-     . lon_dyn, lat_dyn, lev_dyn, ttm_dyn, 
+     . lon_dyn, lat_dyn, levdyn_ini, ttm_dyn, 
      . itau, date, dt, fid_dyn)
     !
@@ -764,18 +863,25 @@
         DEALLOCATE(lon_rad, stat=iret)
       endif
-        ALLOCATE(lon_rad(iml_dyn), stat=iret)
+
+      ALLOCATE(lon_rad(iml_dyn), stat=iret)
+      ALLOCATE(lon_ini(iml_dyn))
        
       IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
-          lon_rad(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
-      ELSE
-          lon_rad(:) = lon_dyn(:,1) 
-      ENDIF
+          lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
+      ELSE
+          lon_ini(:) = lon_dyn(:,1) 
+      ENDIF
+
       ALLOCATE(lat_rad(jml_dyn))
+      ALLOCATE(lat_ini(jml_dyn))
+
       IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
-          lat_rad(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
-      ELSE
-          lat_rad(:) = lat_dyn(1,:) 
-      ENDIF
-    !
+          lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
+      ELSE
+          lat_ini(:) = lat_dyn(1,:) 
+      ENDIF
+    !
+
+
       ALLOCATE(z(iml, jml))
       ALLOCATE(tmp_var(iml-1,jml))
@@ -783,6 +889,23 @@
       CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn,
      .              1, 1, var_ana)
-      CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
+c
+      title='Z'
+      CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini,
+     . lon_rad, lat_rad, var_ana, interbar  )
+c
+      IF ( interbar )   THEN
+        WRITE(6,*) '-------------------------------------------------',
+     ,'--------------'
+        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour  Z  $$$ '
+        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) 
+      ELSE
+        CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
      .               iml-1, jml, lon_in, lat_in, tmp_var)
+      ENDIF
+
       CALL gr_int_dyn(tmp_var, z, iml-1, jml)
     !
@@ -791,6 +914,23 @@
       CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn,
      .              1, 1, var_ana)
-      CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
-     .               iml-1, jml, lon_in, lat_in, tmp_var)
+
+       title='SP'
+      CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini,
+     . lon_rad, lat_rad, var_ana, interbar  )
+
+      IF ( interbar )   THEN
+        WRITE(6,*) '-------------------------------------------------',
+     ,'--------------'
+        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour  SP  $$$ '
+        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) 
+      ELSE
+        CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
+     .             iml-1, jml, lon_in, lat_in, tmp_var  )
+      ENDIF
+
       CALL gr_int_dyn(tmp_var, psol_dyn, iml-1, jml)
     !
@@ -800,5 +940,6 @@
     !   coming out of the restart file. In case we dor have it we will initialize it.
     !
-        CALL start_init_phys( iml, jml, lon_in, lat_in)
+        CALL start_init_phys( iml, jml, lon_in, lat_in,jml2,lon_in2,
+     .                 lat_in2 , interbar )
       ELSE
         IF ( SIZE(tsol) .NE. SIZE(psol_dyn) ) THEN
@@ -814,5 +955,6 @@
             !     coming out of the restart file. In case we dor have it we will initialize it.
             !
-          CALL start_init_orog( iml, jml, lon_in, lat_in)
+        CALL start_init_orog( iml, jml, lon_in, lat_in, jml2, lon_in2 ,
+     .      lat_in2 , interbar )
             !
       ELSE
@@ -860,5 +1002,5 @@
     !
       SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in,
-     .                           lat_in, pls_in, var3d)
+     .      lat_in, jml2, lon_in2, lat_in2, pls_in, var3d, interbar )
     !
     !    This subroutine gets a variables from a 3D file and does the interpolations needed
@@ -868,15 +1010,21 @@
     !
       CHARACTER*(*) :: varname
-      INTEGER :: iml, jml, lml
+      INTEGER :: iml, jml, lml, jml2
       REAL :: lon_in(iml), lat_in(jml), pls_in(iml, jml, lml)
+      REAL :: lon_in2(iml) , lat_in2(jml2)
       REAL :: var3d(iml, jml, lml)
+      LOGICAL interbar
+      real chmin,chmax
     !
     !  LOCAL
     !
-      INTEGER :: ii, ij, il
+      CHARACTER*25 title
+      INTEGER :: ii, ij, il, jsort,i,j,l
       REAL :: bx, by
       REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
+      REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) , lev_dyn(:)
       REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:)
       REAL, ALLOCATABLE :: ax(:), ay(:), yder(:)
+       REAL, ALLOCATABLE :: varrr(:,:,:)
       INTEGER, ALLOCATABLE :: lind(:)
     !
@@ -886,4 +1034,5 @@
           ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
       ENDIF
+          ALLOCATE(varrr(iml_dyn, jml_dyn, llm_dyn))
     !
     !
@@ -900,16 +1049,28 @@
     !
       ALLOCATE(lon_rad(iml_dyn))
+      ALLOCATE(lon_ini(iml_dyn))
+
       IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
-          lon_rad(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
-      ELSE
-          lon_rad(:) = lon_dyn(:,1) 
-      ENDIF
+          lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
+      ELSE
+          lon_ini(:) = lon_dyn(:,1) 
+      ENDIF
+
       ALLOCATE(lat_rad(jml_dyn))
+      ALLOCATE(lat_ini(jml_dyn))
+
+      ALLOCATE(lev_dyn(llm_dyn))
+
       IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
-          lat_rad(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
-      ELSE
-          lat_rad(:) = lat_dyn(1,:) 
-      ENDIF
-    !
+          lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
+      ELSE
+          lat_ini(:) = lat_dyn(1,:) 
+      ENDIF
+    !
+
+      CALL conf_dat3d ( varname,iml_dyn, jml_dyn, llm_dyn, lon_ini, 
+     . lat_ini, levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d  ,
+     ,  interbar                                                   )
+
       ALLOCATE(var_tmp2d(iml-1, jml))
       ALLOCATE(var_tmp3d(iml, jml, llm_dyn))
@@ -919,8 +1080,22 @@
       ALLOCATE(lind(llm_dyn))
     !
+ 
       DO il=1,llm_dyn
         !
-        CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad, 
-     .var_ana3d(:,:,il), iml-1, jml, lon_in, lat_in, var_tmp2d)
+      IF( interbar )  THEN
+       IF( il.EQ.1 )  THEN
+        WRITE(6,*) '-------------------------------------------------',
+     ,'--------------'
+        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour ', varname
+        WRITE(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 ) 
+      ELSE
+       CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad, 
+     .  var_ana3d(:,:,il), iml-1, jml, lon_in, lat_in, var_tmp2d )
+      ENDIF
         !
         CALL gr_int_dyn(var_tmp2d, var_tmp3d(:,:,il), iml-1, jml)
@@ -928,23 +1103,19 @@
        ENDDO
        !
-       !  IF needed we return the vertical axis. The spline interpolation
-       !  Requires the coordinate to be in increasing order.
-    !
-      IF ( lev_dyn(1) .LT. lev_dyn(llm_dyn)) THEN
-          DO il=1,llm_dyn
-            lind(il) = il
-          ENDDO        
-      ELSE
           DO il=1,llm_dyn
             lind(il) = llm_dyn-il+1
           ENDDO
-      ENDIF
-    !
+    !
+c
+c  ... Pour l'interpolation verticale ,on interpole du haut de l'atmosphere
+c                    vers  le  sol  ...
+c
       DO ij=1,jml
         DO ii=1,iml-1
           !
-          ax(:) = lev_dyn(lind(:)) * 100
+          ax(:) = lev_dyn(lind(:)) 
           ay(:) = var_tmp3d(ii, ij, lind(:))
           !
+         
           CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
           !
@@ -959,4 +1130,18 @@
       ENDDO
 
+      do il=1,lml
+        call minmax(iml*jml,var3d(1,1,il),chmin,chmax)
+      SELECTCASE(varname)
+       CASE('U')
+          WRITE(*,*) ' U  min max l ',il,chmin,chmax
+       CASE('V')
+          WRITE(*,*) ' V  min max l ',il,chmin,chmax
+       CASE('TEMP')
+          WRITE(*,*) ' TEMP  min max l ',il,chmin,chmax
+       CASE('R')
+          WRITE(*,*) ' R  min max l ',il,chmin,chmax
+      END SELECT
+           enddo
+
       DEALLOCATE(lon_rad)
       DEALLOCATE(lat_rad)
