Index: trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F	(revision 134)
+++ 	(revision )
@@ -1,404 +1,0 @@
-!
-! $Id: advtrac.F 1446 2010-10-22 09:27:25Z emillour $
-!
-c
-c
-      SUBROUTINE advtrac(pbaru,pbarv ,
-     *                   p,  masse,q,iapptrac,teta,
-     *                  flxw,
-     *                  pk)
-c     Auteur :  F. Hourdin
-c
-c     Modif. P. Le Van     (20/12/97)
-c            F. Codron     (10/99)
-c            D. Le Croller (07/2001)
-c            M.A Filiberti (04/2002)
-c
-      USE infotrac
-      USE control_mod
- 
-
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comvert.h"
-#include "comdissip.h"
-#include "comgeom2.h"
-#include "logic.h"
-#include "temps.h"
-#include "ener.h"
-#include "description.h"
-#include "iniprint.h"
-
-c-------------------------------------------------------------------
-c     Arguments
-c-------------------------------------------------------------------
-c     Ajout PPM
-c--------------------------------------------------------
-      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
-c--------------------------------------------------------
-      INTEGER iapptrac
-      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
-      REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
-      REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
-      REAL pk(ip1jmp1,llm)
-      REAL flxw(ip1jmp1,llm)
-
-c-------------------------------------------------------------
-c     Variables locales
-c-------------------------------------------------------------
-
-      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
-      REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
-      REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 
-      REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
-      INTEGER iadvtr
-      INTEGER ij,l,iq,iiq
-      REAL zdpmin, zdpmax
-      EXTERNAL  minmax
-      SAVE iadvtr, massem, pbaruc, pbarvc
-      DATA iadvtr/0/
-c----------------------------------------------------------
-c     Rajouts pour PPM
-c----------------------------------------------------------
-      INTEGER indice,n
-      REAL dtbon ! Pas de temps adaptatif pour que CFL<1
-      REAL CFLmaxz,aaa,bbb ! CFL maximum
-      REAL psppm(iim,jjp1) ! pression  au sol
-      REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
-      REAL qppm(iim*jjp1,llm,nqtot)
-      REAL fluxwppm(iim,jjp1,llm)
-      REAL apppm(llmp1), bpppm(llmp1)
-      LOGICAL dum,fill
-      DATA fill/.true./
-      DATA dum/.true./
-
-      integer,save :: countcfl=0
-      real cflx(ip1jmp1,llm)
-      real cfly(ip1jm,llm)
-      real cflz(ip1jmp1,llm)
-      real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
-
-      IF(iadvtr.EQ.0) THEN
-         CALL initial0(ijp1llm,pbaruc)
-         CALL initial0(ijmllm,pbarvc)
-      ENDIF
-
-c   accumulation des flux de masse horizontaux
-      DO l=1,llm
-         DO ij = 1,ip1jmp1
-            pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
-         ENDDO
-         DO ij = 1,ip1jm
-            pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
-         ENDDO
-      ENDDO
-
-c   selection de la masse instantannee des mailles avant le transport.
-      IF(iadvtr.EQ.0) THEN
-
-         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
-ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
-c
-      ENDIF
-
-      iadvtr   = iadvtr+1
-      iapptrac = iadvtr
-
-
-c   Test pour savoir si on advecte a ce pas de temps
-      IF ( iadvtr.EQ.iapp_tracvl ) THEN
-
-cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
-cc
-
-c   traitement des flux de masse avant advection.
-c     1. calcul de w
-c     2. groupement des mailles pres du pole.
-
-        CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
-
-      ! ... Flux de masse diaganostiques traceurs
-      flxw = wg / REAL(iapp_tracvl)
-
-c  test sur l'eventuelle creation de valeurs negatives de la masse
-         DO l=1,llm-1
-            DO ij = iip2+1,ip1jm
-              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
-     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
-     s                  +       wg(ij,l+1)  - wg(ij,l)
-            ENDDO
-            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
-            DO ij = iip2,ip1jm
-               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 
-            ENDDO 
-
-
-            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
-
-            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
-            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
-     s        '   MAX:', zdpmax
-            ENDIF
-
-         ENDDO
-
-
-c-------------------------------------------------------------------
-! Calcul des criteres CFL en X, Y et Z
-c-------------------------------------------------------------------
-
-      if (countcfl == 0. ) then
-          cflxmax(:)=0.
-          cflymax(:)=0.
-          cflzmax(:)=0.
-      endif
-
-      countcfl=countcfl+iapp_tracvl
-      cflx(:,:)=0.
-      cfly(:,:)=0.
-      cflz(:,:)=0.
-      do l=1,llm
-         do ij=iip2,ip1jm-1
-            if (pbarug(ij,l)>=0.) then
-                cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
-            else
-                cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
-            endif 
-         enddo
-      enddo
-      do l=1,llm
-         do ij=iip2,ip1jm-1,iip1
-            cflx(ij+iip1,l)=cflx(ij,l)
-         enddo
-      enddo
-
-      do l=1,llm
-         do ij=1,ip1jm
-            if (pbarvg(ij,l)>=0.) then
-                cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
-            else
-                cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
-            endif 
-         enddo
-      enddo
-
-      do l=2,llm
-         do ij=1,ip1jm
-            if (wg(ij,l)>=0.) then
-                cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
-            else
-                cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
-            endif 
-         enddo
-      enddo
-
-      do l=1,llm
-         cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
-         cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
-         cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
-      enddo
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Par defaut, on sort le diagnostic des CFL tous les jours.
-! Si on veut le sortir a chaque pas d'advection en cas de plantage 
-!     if (countcfl==iapp_tracvl) then
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-      if (countcfl==day_step) then
-         do l=1,llm
-         write(lunout,*) 'L, CFLmax '
-     s   ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
-         enddo
-         countcfl=0
-      endif
-    
-c-------------------------------------------------------------------
-c   Advection proprement dite (Modification Le Croller (07/2001)
-c-------------------------------------------------------------------
-
-c----------------------------------------------------
-c        Calcul des moyennes basées sur la masse
-c----------------------------------------------------
-          call massbar(massem,massebx,masseby)          
-
-c-----------------------------------------------------------
-c     Appel des sous programmes d'advection
-c-----------------------------------------------------------
-      do iq=1,nqtot
-c        call clock(t_initial)
-        if(iadv(iq) == 0) cycle 
-c   ----------------------------------------------------------------
-c   Schema de Van Leer I MUSCL
-c   ----------------------------------------------------------------
-        if(iadv(iq).eq.10) THEN
-            call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
-
-c   ----------------------------------------------------------------
-c   Schema "pseudo amont" + test sur humidite specifique
-C    pour la vapeur d'eau. F. Codron
-c   ----------------------------------------------------------------
-        else if(iadv(iq).eq.14) then
-c
-           CALL vlspltqs( q(1,1,1), 2., massem, wg ,
-     *                 pbarug,pbarvg,dtvr,p,pk,teta )
-c   ----------------------------------------------------------------
-c   Schema de Frederic Hourdin
-c   ----------------------------------------------------------------
-        else if(iadv(iq).eq.12) then
-c            Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
-     s             dtvr,'n=',n
-           endif
-           do indice=1,n
-            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
-           end do
-        else if(iadv(iq).eq.13) then
-c            Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
-     s             dtvr,'n=',n
-           endif
-          do indice=1,n
-            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
-          end do
-c   ----------------------------------------------------------------
-c   Schema de pente SLOPES
-c   ----------------------------------------------------------------
-        else if (iadv(iq).eq.20) then
-            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
-
-c   ----------------------------------------------------------------
-c   Schema de Prather
-c   ----------------------------------------------------------------
-        else if (iadv(iq).eq.30) then
-c            Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
-     s             dtvr,'n=',n
-           endif
-           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
-     s                     n,dtbon)
-
-c   ----------------------------------------------------------------
-c   Schemas PPM Lin et Rood
-c   ----------------------------------------------------------------
-         else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
-     s                     iadv(iq).LE.18)) then
-
-c        Test sur le flux horizontal
-c        Pas de temps adaptatif
-         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-         if (n.GT.1) then
-         write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
-     s             dtvr,'n=',n
-         endif
-c        Test sur le flux vertical
-         CFLmaxz=0.
-         do l=2,llm
-           do ij=iip2,ip1jm
-            aaa=wg(ij,l)*dtvr/massem(ij,l)
-            CFLmaxz=max(CFLmaxz,aaa)
-            bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
-            CFLmaxz=max(CFLmaxz,bbb)
-           enddo
-         enddo
-         if (CFLmaxz.GE.1) then
-            write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
-         endif
-
-c-----------------------------------------------------------
-c        Ss-prg interface LMDZ.4->PPM3d
-c-----------------------------------------------------------
-
-          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
-     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
-     s                 unatppm,vnatppm,psppm)
-
-          do indice=1,n
-c---------------------------------------------------------------------
-c                         VL (version PPM) horiz. et PPM vert.
-c---------------------------------------------------------------------
-                if (iadv(iq).eq.11) then
-c                  Ss-prg PPM3d de Lin
-                  call ppm3d(1,qppm(1,1,iq),
-     s                       psppm,psppm,
-     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
-     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
-     s                       fill,dum,220.)
-
-c----------------------------------------------------------------------
-c                           Monotonic PPM
-c----------------------------------------------------------------------
-               else if (iadv(iq).eq.16) then
-c                  Ss-prg PPM3d de Lin
-                  call ppm3d(1,qppm(1,1,iq),
-     s                       psppm,psppm,
-     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
-     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
-     s                       fill,dum,220.)
-c---------------------------------------------------------------------
-
-c---------------------------------------------------------------------
-c                           Semi Monotonic PPM
-c---------------------------------------------------------------------
-               else if (iadv(iq).eq.17) then
-c                  Ss-prg PPM3d de Lin
-                  call ppm3d(1,qppm(1,1,iq),
-     s                       psppm,psppm,
-     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
-     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
-     s                       fill,dum,220.)
-c---------------------------------------------------------------------
-
-c---------------------------------------------------------------------
-c                         Positive Definite PPM
-c---------------------------------------------------------------------
-                else if (iadv(iq).eq.18) then
-c                  Ss-prg PPM3d de Lin
-                  call ppm3d(1,qppm(1,1,iq),
-     s                       psppm,psppm,
-     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
-     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
-     s                       fill,dum,220.)
-c---------------------------------------------------------------------
-                endif
-            enddo
-c-----------------------------------------------------------------
-c               Ss-prg interface PPM3d-LMDZ.4
-c-----------------------------------------------------------------
-                  call interpost(q(1,1,iq),qppm(1,1,iq))
-            endif
-c----------------------------------------------------------------------
-
-c-----------------------------------------------------------------
-c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
-c et Nord j=1
-c-----------------------------------------------------------------
-
-c                  call traceurpole(q(1,1,iq),massem)
-
-c calcul du temps cpu pour un schema donne
-
-c                  call clock(t_final)
-cym                  tps_cpu=t_final-t_initial
-cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
-
-       end DO
-
-
-c------------------------------------------------------------------
-c   on reinitialise a zero les flux de masse cumules
-c---------------------------------------------------
-          iadvtr=0
-
-       ENDIF ! if iadvtr.EQ.iapp_tracvl
-
-       RETURN
-       END
-
Index: trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90	(revision 270)
+++ trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90	(revision 270)
@@ -0,0 +1,396 @@
+! $Id: advtrac.F90 1549 2011-07-05 08:41:12Z lguez $
+
+SUBROUTINE advtrac(pbaru,pbarv , p,  masse,q,iapptrac,teta, flxw, pk)
+  !     Auteur :  F. Hourdin
+  !
+  !     Modif. P. Le Van     (20/12/97)
+  !            F. Codron     (10/99)
+  !            D. Le Croller (07/2001)
+  !            M.A Filiberti (04/2002)
+  !
+  USE infotrac
+  USE control_mod
+
+
+  IMPLICIT NONE
+  !
+  include "dimensions.h"
+  include "paramet.h"
+  include "comconst.h"
+  include "comvert.h"
+  include "comdissip.h"
+  include "comgeom2.h"
+  include "logic.h"
+  include "temps.h"
+  include "ener.h"
+  include "description.h"
+  include "iniprint.h"
+
+  !-------------------------------------------------------------------
+  !     Arguments
+  !-------------------------------------------------------------------
+  !     Ajout PPM
+  !--------------------------------------------------------
+  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
+  !--------------------------------------------------------
+  INTEGER iapptrac
+  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+  REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
+  REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
+  REAL pk(ip1jmp1,llm)
+  REAL flxw(ip1jmp1,llm)
+
+  !-------------------------------------------------------------
+  !     Variables locales
+  !-------------------------------------------------------------
+
+  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
+  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
+  REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 
+  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
+  INTEGER iadvtr
+  INTEGER ij,l,iq,iiq
+  REAL zdpmin, zdpmax
+  EXTERNAL  minmax
+  SAVE iadvtr, massem, pbaruc, pbarvc
+  DATA iadvtr/0/
+  !----------------------------------------------------------
+  !     Rajouts pour PPM
+  !----------------------------------------------------------
+  INTEGER indice,n
+  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
+  REAL CFLmaxz,aaa,bbb ! CFL maximum
+  REAL psppm(iim,jjp1) ! pression  au sol
+  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
+  REAL qppm(iim*jjp1,llm,nqtot)
+  REAL fluxwppm(iim,jjp1,llm)
+  REAL apppm(llmp1), bpppm(llmp1)
+  LOGICAL dum,fill
+  DATA fill/.true./
+  DATA dum/.true./
+
+  integer,save :: countcfl=0
+  real cflx(ip1jmp1,llm)
+  real cfly(ip1jm,llm)
+  real cflz(ip1jmp1,llm)
+  real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
+
+  IF(iadvtr.EQ.0) THEN
+     CALL initial0(ijp1llm,pbaruc)
+     CALL initial0(ijmllm,pbarvc)
+  ENDIF
+
+  !   accumulation des flux de masse horizontaux
+  DO l=1,llm
+     DO ij = 1,ip1jmp1
+        pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
+     ENDDO
+     DO ij = 1,ip1jm
+        pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
+     ENDDO
+  ENDDO
+
+  !   selection de la masse instantannee des mailles avant le transport.
+  IF(iadvtr.EQ.0) THEN
+
+     CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
+     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
+     !
+  ENDIF
+
+  iadvtr   = iadvtr+1
+  iapptrac = iadvtr
+
+
+  !   Test pour savoir si on advecte a ce pas de temps
+  IF ( iadvtr.EQ.iapp_tracvl ) THEN
+
+     !c   ..  Modif P.Le Van  ( 20/12/97 )  ....
+     !c
+
+     !   traitement des flux de masse avant advection.
+     !     1. calcul de w
+     !     2. groupement des mailles pres du pole.
+
+     CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
+
+     ! ... Flux de masse diaganostiques traceurs
+     flxw = wg / REAL(iapp_tracvl)
+
+     !  test sur l'eventuelle creation de valeurs negatives de la masse
+     DO l=1,llm-1
+        DO ij = iip2+1,ip1jm
+           zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l) &
+                - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
+                +       wg(ij,l+1)  - wg(ij,l)
+        ENDDO
+        CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
+        DO ij = iip2,ip1jm
+           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 
+        ENDDO
+
+
+        CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
+
+        IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
+           PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin, &
+                '   MAX:', zdpmax
+        ENDIF
+
+     ENDDO
+
+
+     !-------------------------------------------------------------------
+     ! Calcul des criteres CFL en X, Y et Z
+     !-------------------------------------------------------------------
+
+     if (countcfl == 0. ) then
+        cflxmax(:)=0.
+        cflymax(:)=0.
+        cflzmax(:)=0.
+     endif
+
+     countcfl=countcfl+iapp_tracvl
+     cflx(:,:)=0.
+     cfly(:,:)=0.
+     cflz(:,:)=0.
+     do l=1,llm
+        do ij=iip2,ip1jm-1
+           if (pbarug(ij,l)>=0.) then
+              cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
+           else
+              cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
+           endif
+        enddo
+     enddo
+     do l=1,llm
+        do ij=iip2,ip1jm-1,iip1
+           cflx(ij+iip1,l)=cflx(ij,l)
+        enddo
+     enddo
+
+     do l=1,llm
+        do ij=1,ip1jm
+           if (pbarvg(ij,l)>=0.) then
+              cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
+           else
+              cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
+           endif
+        enddo
+     enddo
+
+     do l=2,llm
+        do ij=1,ip1jm
+           if (wg(ij,l)>=0.) then
+              cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
+           else
+              cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
+           endif
+        enddo
+     enddo
+
+     do l=1,llm
+        cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
+        cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
+        cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
+     enddo
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     ! Par defaut, on sort le diagnostic des CFL tous les jours.
+     ! Si on veut le sortir a chaque pas d'advection en cas de plantage 
+     !     if (countcfl==iapp_tracvl) then
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     if (countcfl==day_step) then
+        do l=1,llm
+           write(lunout,*) 'L, CFLmax ' &
+                ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
+        enddo
+        countcfl=0
+     endif
+
+     !-------------------------------------------------------------------
+     !   Advection proprement dite (Modification Le Croller (07/2001)
+     !-------------------------------------------------------------------
+
+     !----------------------------------------------------
+     !        Calcul des moyennes basées sur la masse
+     !----------------------------------------------------
+     call massbar(massem,massebx,masseby)          
+
+     !-----------------------------------------------------------
+     !     Appel des sous programmes d'advection
+     !-----------------------------------------------------------
+     do iq=1,nqtot
+        !        call clock(t_initial)
+        if(iadv(iq) == 0) cycle 
+        !   ----------------------------------------------------------------
+        !   Schema de Van Leer I MUSCL
+        !   ----------------------------------------------------------------
+        if(iadv(iq).eq.10) THEN
+           call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
+
+           !   ----------------------------------------------------------------
+           !   Schema "pseudo amont" + test sur humidite specifique
+           !    pour la vapeur d'eau. F. Codron
+           !   ----------------------------------------------------------------
+        else if(iadv(iq).eq.14) then
+           !
+           CALL vlspltqs( q(1,1,1), 2., massem, wg , &
+                pbarug,pbarvg,dtvr,p,pk,teta )
+           !   ----------------------------------------------------------------
+           !   Schema de Frederic Hourdin
+           !   ----------------------------------------------------------------
+        else if(iadv(iq).eq.12) then
+           !            Pas de temps adaptatif
+           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           do indice=1,n
+              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
+           end do
+        else if(iadv(iq).eq.13) then
+           !            Pas de temps adaptatif
+           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           do indice=1,n
+              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
+           end do
+           !   ----------------------------------------------------------------
+           !   Schema de pente SLOPES
+           !   ----------------------------------------------------------------
+        else if (iadv(iq).eq.20) then
+           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
+
+           !   ----------------------------------------------------------------
+           !   Schema de Prather
+           !   ----------------------------------------------------------------
+        else if (iadv(iq).eq.30) then
+           !            Pas de temps adaptatif
+           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
+                n,dtbon)
+
+           !   ----------------------------------------------------------------
+           !   Schemas PPM Lin et Rood
+           !   ----------------------------------------------------------------
+        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
+             iadv(iq).LE.18)) then
+
+           !        Test sur le flux horizontal
+           !        Pas de temps adaptatif
+           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           !        Test sur le flux vertical
+           CFLmaxz=0.
+           do l=2,llm
+              do ij=iip2,ip1jm
+                 aaa=wg(ij,l)*dtvr/massem(ij,l)
+                 CFLmaxz=max(CFLmaxz,aaa)
+                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
+                 CFLmaxz=max(CFLmaxz,bbb)
+              enddo
+           enddo
+           if (CFLmaxz.GE.1) then
+              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
+           endif
+
+           !-----------------------------------------------------------
+           !        Ss-prg interface LMDZ.4->PPM3d
+           !-----------------------------------------------------------
+
+           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
+                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
+                unatppm,vnatppm,psppm)
+
+           do indice=1,n
+              !----------------------------------------------------------------
+              !                         VL (version PPM) horiz. et PPM vert.
+              !----------------------------------------------------------------
+              if (iadv(iq).eq.11) then
+                 !                  Ss-prg PPM3d de Lin
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+
+                 !-------------------------------------------------------------
+                 !                           Monotonic PPM
+                 !-------------------------------------------------------------
+              else if (iadv(iq).eq.16) then
+                 !                  Ss-prg PPM3d de Lin
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+                 !-------------------------------------------------------------
+
+                 !-------------------------------------------------------------
+                 !                           Semi Monotonic PPM
+                 !-------------------------------------------------------------
+              else if (iadv(iq).eq.17) then
+                 !                  Ss-prg PPM3d de Lin
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+                 !-------------------------------------------------------------
+
+                 !-------------------------------------------------------------
+                 !                         Positive Definite PPM
+                 !-------------------------------------------------------------
+              else if (iadv(iq).eq.18) then
+                 !                  Ss-prg PPM3d de Lin
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+                 !-------------------------------------------------------------
+              endif
+           enddo
+           !-----------------------------------------------------------------
+           !               Ss-prg interface PPM3d-LMDZ.4
+           !-----------------------------------------------------------------
+           call interpost(q(1,1,iq),qppm(1,1,iq))
+        endif
+        !----------------------------------------------------------------------
+
+        !-----------------------------------------------------------------
+        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
+        ! et Nord j=1
+        !-----------------------------------------------------------------
+
+        !                  call traceurpole(q(1,1,iq),massem)
+
+        ! calcul du temps cpu pour un schema donne
+
+        !                  call clock(t_final)
+        !ym                  tps_cpu=t_final-t_initial
+        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
+
+     end DO
+
+
+     !------------------------------------------------------------------
+     !   on reinitialise a zero les flux de masse cumules
+     !---------------------------------------------------
+     iadvtr=0
+
+  ENDIF ! if iadvtr.EQ.iapp_tracvl
+
+END SUBROUTINE advtrac
Index: trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90	(revision 270)
@@ -1,4 +1,4 @@
 !
-! $Id: ce0l.F90 1425 2010-09-02 13:45:23Z lguez $
+! $Id: ce0l.F90 1511 2011-04-28 15:21:47Z jghattas $
 !
 !-------------------------------------------------------------------------------
@@ -27,5 +27,4 @@
 #endif
   IMPLICIT NONE
-#include "iniprint.h"
 #ifndef CPP_EARTH
   WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
@@ -37,10 +36,14 @@
 #include "paramet.h"
 #include "indicesol.h"
+#include "iniprint.h"
 #include "temps.h"
 #include "logic.h"
+  INTEGER, PARAMETER            :: longcles=20
+  REAL,    DIMENSION(longcles)  :: clesphy0
   REAL,    DIMENSION(iip1,jjp1) :: masque
   CHARACTER(LEN=15)             :: calnd
+  REAL,    DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol
 !-------------------------------------------------------------------------------
-  CALL conf_gcm( 99, .TRUE. )
+  CALL conf_gcm( 99, .TRUE. , clesphy0 )
 
   CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
@@ -78,5 +81,5 @@
   WRITE(lunout,'(//)')
   WRITE(lunout,*) ' interbar = ',interbar
-  CALL etat0_netcdf(interbar,masque,ok_etat0)
+  CALL etat0_netcdf(interbar,masque,phis,ok_etat0)
 
   IF(ok_limit) THEN
@@ -89,4 +92,12 @@
   END IF
 
+  IF (grilles_gcm_netcdf) THEN
+     WRITE(lunout,'(//)')
+     WRITE(lunout,*) '  ***************************  '
+     WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
+     WRITE(lunout,*) '  ***************************  '
+     WRITE(lunout,'(//)')
+     CALL grilles_gcm_netcdf_sub(masque,phis)
+  END IF
 #endif
 ! of #ifndef CPP_EARTH #else
Index: trunk/LMDZ.COMMON/libf/dyn3d/comdissipn.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/comdissipn.h	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/comdissipn.h	(revision 270)
@@ -2,15 +2,19 @@
 ! $Header$
 !
-c-----------------------------------------------------------------------
-c INCLUDE comdissipn.h
+!  Attention : ce fichier include est compatible format fixe/format libre
+!                 veillez à n'utiliser que des ! pour les commentaires
+!                 et à bien positionner les & des lignes de continuation 
+!                 (les placer en colonne 6 et en colonne 73)
+!-----------------------------------------------------------------------
+! INCLUDE comdissipn.h
 
       REAL  tetaudiv, tetaurot, tetah, cdivu, crot, cdivh
-c
-      COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm)   ,
-     1                        cdivu,      crot,         cdivh
+!
+      COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm)   ,     &
+     &                        cdivu,      crot,         cdivh
 
-c
-c    Les parametres de ce common proviennent des calculs effectues dans 
-c             Inidissip  .
-c
-c-----------------------------------------------------------------------
+!
+!    Les parametres de ce common proviennent des calculs effectues dans 
+!             Inidissip  .
+!
+!-----------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 270)
@@ -230,11 +230,12 @@
        CALL getin('output_grads_dyn',output_grads_dyn)
 
-!Config  Key  = idissip
+!Config  Key  = dissip_period
 !Config  Desc = periode de la dissipation 
-!Config  Def  = 10
+!Config  Def  = 0
 !Config  Help = periode de la dissipation 
-!Config         (en pas) ... a completer !
-       idissip = 10
-       CALL getin('idissip',idissip)
+!Config  dissip_period=0 => la valeur sera calcule dans inidissip       
+!Config  dissip_period>0 => on prend cette valeur
+       dissip_period = 0
+       CALL getin('dissip_period',dissip_period)
 
 ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
@@ -652,5 +653,5 @@
       write(lunout,*)' periodav = ', periodav 
       write(lunout,*)' output_grads_dyn = ', output_grads_dyn
-      write(lunout,*)' idissip = ', idissip
+      write(lunout,*)' dissip_period = ', dissip_period
       write(lunout,*)' lstardis = ', lstardis
       write(lunout,*)' nitergdiv = ', nitergdiv
@@ -876,4 +877,10 @@
       CALL getin('ok_etat0',ok_etat0)
 
+!Config  Key  = grilles_gcm_netcdf
+!Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
+!Config  Def  = n
+      grilles_gcm_netcdf = .FALSE.
+      CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
+
       write(lunout,*)' #########################################'
       write(lunout,*)' Configuration des parametres de cel0'
@@ -890,5 +897,5 @@
       write(lunout,*)' periodav = ', periodav 
       write(lunout,*)' output_grads_dyn = ', output_grads_dyn
-      write(lunout,*)' idissip = ', idissip
+      write(lunout,*)' dissip_period = ', dissip_period
       write(lunout,*)' lstardis = ', lstardis
       write(lunout,*)' nitergdiv = ', nitergdiv
@@ -922,4 +929,5 @@
       write(lunout,*)' ok_limit = ', ok_limit
       write(lunout,*)' ok_etat0 = ', ok_etat0
+      write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
 c
       RETURN
Index: trunk/LMDZ.COMMON/libf/dyn3d/control_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/control_mod.F90	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/control_mod.F90	(revision 270)
@@ -12,5 +12,5 @@
   REAL    :: periodav
   INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys
-  INTEGER :: iconser,iecri,idissip,iphysiq,iecrimoy
+  INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy
   INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn
   LOGICAL :: offline
Index: trunk/LMDZ.COMMON/libf/dyn3d/defrun.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/defrun.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/defrun.F	(revision 270)
@@ -139,7 +139,7 @@
 
       READ (tapedef,9001) ch1,ch4
-      READ (tapedef,*)    idissip
-      WRITE(tapeout,9001) ch1,'idissip'
-      WRITE(tapeout,*)    idissip
+      READ (tapedef,*)    dissip_period
+      WRITE(tapeout,9001) ch1,'dissip_period'
+      WRITE(tapeout,*)    dissip_period
 
 ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
Index: trunk/LMDZ.COMMON/libf/dyn3d/gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 270)
@@ -427,5 +427,5 @@
 c   -------------------------------
 
-      IF (call_iniphys.and.(iflag_phys.eq.1)) THEN
+      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
          latfi(1)=rlatu(1)
          lonfi(1)=0.
@@ -486,4 +486,12 @@
 #endif
 
+#ifdef CPP_EARTH
+! Create start file (startphy.nc) and boundary conditions (limit.nc) 
+! for the Earth verstion
+       if (iflag_phys>=100) then
+          call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
+       endif
+#endif
+
       if (planet_type.eq."mars") then
 ! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem0_mars
Index: trunk/LMDZ.COMMON/libf/dyn3d/grilles_gcm_netcdf_sub.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/grilles_gcm_netcdf_sub.F90	(revision 270)
+++ trunk/LMDZ.COMMON/libf/dyn3d/grilles_gcm_netcdf_sub.F90	(revision 270)
@@ -0,0 +1,237 @@
+!
+! $Header$
+!
+! This subroutine creates the file grilles_gcm.nc containg longitudes and
+! latitudes in degrees for grid u and v. This subroutine is called from
+! ce0l if grilles_gcm_netcdf=TRUE. This subroutine corresponds to the first 
+! part in the program create_fausse_var.
+!
+SUBROUTINE grilles_gcm_netcdf_sub(masque,phis)
+
+  IMPLICIT NONE
+
+  INCLUDE "dimensions.h"
+  INCLUDE "paramet.h"
+  INCLUDE "comconst.h"
+  INCLUDE "comgeom.h"
+  INCLUDE "comvert.h"
+  INCLUDE "netcdf.inc"
+  INCLUDE "serre.h"
+
+
+  REAL,DIMENSION(iip1,jjp1),INTENT(IN)  :: masque ! masque terre/mer
+  REAL,DIMENSION(iip1,jjp1),INTENT(IN)  :: phis   ! geopotentiel au sol
+
+  REAL temp(iim+1,jjm+1)
+  ! Attributs netcdf sortie
+  INTEGER ncid_out,rcode_out
+  INTEGER out_lonuid,out_lonvid,out_latuid,out_latvid,out_levid
+  INTEGER out_varid
+  INTEGER out_lonudim,out_lonvdim
+  INTEGER out_latudim,out_latvdim,out_dim(3)
+  INTEGER out_levdim
+
+  INTEGER, PARAMETER :: longcles = 20
+  REAL  clesphy0(longcles)
+
+  INTEGER start(4),COUNT(4)
+
+  INTEGER status,i,j
+  REAL rlatudeg(jjp1),rlatvdeg(jjm),rlevdeg(llm)
+  REAL rlonudeg(iip1),rlonvdeg(iip1)
+
+  REAL dlon1(iip1),dlon2(iip1),dlat1(jjp1),dlat2(jjp1)
+  REAL acoslat,dxkm,dykm,resol(iip1,jjp1)
+  REAL,DIMENSION(iip1,jjp1)  :: phis_loc
+  INTEGER masque_int(iip1,jjp1)
+  INTEGER :: phis_id
+  INTEGER :: area_id
+  INTEGER :: mask_id
+  
+  rad = 6400000
+  omeg = 7.272205e-05
+  g = 9.8
+  kappa = 0.285716
+  daysec = 86400
+  cpp = 1004.70885
+
+  preff = 101325.
+  pa= 50000.
+
+  CALL conf_gcm( 99, .TRUE. , clesphy0 )
+  CALL iniconst
+  CALL inigeom
+
+  DO j=1,jjp1
+     rlatudeg(j)=rlatu(j)*180./pi
+  ENDDO
+  DO j=1,jjm
+     rlatvdeg(j)=rlatv(j)*180./pi
+  ENDDO
+
+  DO i=1,iip1
+     rlonudeg(i)=rlonu(i)*180./pi + 360.
+     rlonvdeg(i)=rlonv(i)*180./pi + 360.
+  ENDDO
+
+
+  !  2 ----- OUVERTURE DE LA SORTIE NETCDF
+  ! ---------------------------------------------------
+  ! CREATION OUTPUT
+  ! ouverture fichier netcdf de sortie out
+  status=NF_CREATE('grilles_gcm.nc',NF_NOCLOBBER,ncid_out)
+  status=NF_DEF_DIM(ncid_out,'lonu',iim+1,out_lonudim)
+  status=NF_DEF_DIM(ncid_out,'lonv',iim+1,out_lonvdim)
+  status=NF_DEF_DIM(ncid_out,'latu',jjm+1,out_latudim)
+  status=NF_DEF_DIM(ncid_out,'latv',jjm,out_latvdim)
+
+
+  !   Longitudes en u
+  status=NF_DEF_VAR(ncid_out,'lonu',NF_FLOAT,1,out_lonudim, out_lonuid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_lonuid,'units', 12,'degrees_east')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_lonuid,'long_name',9,'Longitude en u')
+
+  !   Longitudes en v
+  status=NF_DEF_VAR(ncid_out,'lonv',NF_FLOAT,1,out_lonvdim, out_lonvid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_lonvid,'units', 12,'degrees_east')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_lonvid,'long_name', 9,'Longitude en v')
+
+  !   Latitude en u
+  status=NF_DEF_VAR(ncid_out,'latu',NF_FLOAT,1,out_latudim, out_latuid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_latuid,'units', 13,'degrees_north')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_latuid,'long_name', 8,'Latitude en u')
+
+  !  Latitude en v
+  status=NF_DEF_VAR(ncid_out,'latv',NF_FLOAT,1,out_latvdim, out_latvid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_latvid,'units', 13,'degrees_north')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_latvid,'long_name', 8,'Latitude en v')
+
+  !   ecriture de la grille u
+  out_dim(1)=out_lonudim
+  out_dim(2)=out_latudim
+  status=NF_DEF_VAR(ncid_out,'grille_u',NF_FLOAT,2,out_dim, out_varid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units', 6,'Kelvin')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name', 16,'Grille aux point u')
+
+  !   ecriture de la grille v
+  out_dim(1)=out_lonvdim
+  out_dim(2)=out_latvdim
+  status=NF_DEF_VAR(ncid_out,'grille_v',NF_FLOAT,2,out_dim, out_varid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units', 6,'Kelvin')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name', 16,'Grille aux point v')
+
+  !   ecriture de la grille u
+  out_dim(1)=out_lonvdim
+  out_dim(2)=out_latudim
+  status=NF_DEF_VAR(ncid_out,'grille_s',NF_FLOAT,2,out_dim, out_varid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units', 6,'Kelvin')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name',16,'Grille aux point u')
+
+  status=NF_ENDDEF(ncid_out)
+  ! 5) ----- FERMETURE DES FICHIERS NETCDF------------------
+  ! --------------------------------------------------------
+  ! 3-b- Ecriture de la grille pour la sortie
+  ! rajoute l'ecriture de la grille
+
+#ifdef NC_DOUBLE
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_lonuid,1,iim+1,rlonudeg)
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_lonvid,1,iim+1,rlonvdeg)
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_latuid,1,jjm+1,rlatudeg)
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_latvid,1,jjm,rlatvdeg)
+#else
+  status=NF_PUT_VARA_REAL(ncid_out,out_lonuid,1,iim+1,rlonudeg)
+  status=NF_PUT_VARA_REAL(ncid_out,out_lonvid,1,iim+1,rlonvdeg)
+  status=NF_PUT_VARA_REAL(ncid_out,out_latuid,1,jjm+1,rlatudeg)
+  status=NF_PUT_VARA_REAL(ncid_out,out_latvid,1,jjm,rlatvdeg)
+#endif
+
+  start(1)=1
+  start(2)=1
+  start(3)=1
+  start(4)=1
+
+  COUNT(1)=iim+1
+  COUNT(2)=jjm+1
+  COUNT(3)=1
+  COUNT(4)=1
+
+  DO j=1,jjm+1
+     DO i=1,iim+1
+        temp(i,j)=MOD(i,2)+MOD(j,2)
+     ENDDO
+  ENDDO
+
+#ifdef NC_DOUBLE
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_varid,start, count,temp)
+#else
+  status=NF_PUT_VARA_REAL(ncid_out,out_varid,start, count,temp)
+#endif
+
+  ! On re-ouvre le fichier pour rajouter 4 nouvelles variables necessaire pour INCA
+! lev - phis - aire - mask
+  rlevdeg(:) = presnivs
+  phis_loc(:,:) = phis(:,:)/g
+
+! niveaux de pression verticaux
+  status = NF_REDEF (ncid_out)
+  status=NF_DEF_DIM(ncid_out,'lev',llm,out_levdim)
+  
+! fields
+  out_dim(1)=out_lonvdim
+  out_dim(2)=out_latudim
+
+  status = nf_def_var(ncid_out,'phis',NF_FLOAT,2,out_dim,phis_id)
+  CALL handle_err(status)
+  status = nf_def_var(ncid_out,'aire',NF_FLOAT,2,out_dim,area_id)
+  CALL handle_err(status)
+  status = nf_def_var(ncid_out,'mask',NF_INT  ,2,out_dim,mask_id)
+  CALL handle_err(status)
+
+  status=NF_ENDDEF(ncid_out)
+
+  ! ecriture des variables
+#ifdef NC_DOUBLE
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_levid,1,llm,rlevdeg)
+#else
+  status=NF_PUT_VARA_REAL(ncid_out,out_levid,1,llm,rlevdeg)
+#endif
+
+  start(1)=1
+  start(2)=1
+  start(3)=1
+  start(4)=0
+  COUNT(1)=iip1
+  COUNT(2)=jjp1
+  COUNT(3)=1
+  COUNT(4)=0
+
+  status = nf_put_vara_double(ncid_out, phis_id,start,count, phis_loc)
+  status = nf_put_vara_double(ncid_out, area_id,start,count, aire)
+  masque_int(:,:) = nINT(masque(:,:))
+  status = nf_put_vara_int(ncid_out, mask_id,start,count,masque_int)
+  CALL handle_err(status)
+  
+  ! fermeture du fichier netcdf
+  CALL ncclos(ncid_out,rcode_out)
+
+END SUBROUTINE grilles_gcm_netcdf_sub
+
+
+
+SUBROUTINE handle_err(status)
+  INCLUDE "netcdf.inc"
+
+  INTEGER status
+  IF (status.NE.nf_noerr) THEN
+     PRINT *,NF_STRERROR(status)
+     CALL abort_gcm('grilles_gcm_netcdf','netcdf error',1)
+  ENDIF
+END SUBROUTINE handle_err
+
Index: trunk/LMDZ.COMMON/libf/dyn3d/ini_paramLMDZ_dyn.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/ini_paramLMDZ_dyn.h	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/ini_paramLMDZ_dyn.h	(revision 270)
@@ -78,5 +78,5 @@
      .                "once", dt_cum,dt_cum)
 c
-         CALL histdef(nid_ctesGCM, "idissip", 
+         CALL histdef(nid_ctesGCM, "dissip_period", 
      .  "periode de la dissipation (en pas) ... a completer",
      .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
Index: trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90	(revision 270)
@@ -1,4 +1,4 @@
 !
-! $Id: iniacademic.F90 1520 2011-05-23 11:37:09Z emillour $
+! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $
 !
 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
@@ -115,5 +115,5 @@
   endif
 
-  academic_case: if (iflag_phys == 2) then
+  academic_case: if (iflag_phys >= 2) then
      ! initializations
 
@@ -208,5 +208,9 @@
      IF (.NOT. read_start) THEN
         ! surface pressure
-        ps(:)=preff
+        if (iflag_phys>2) then
+           ps(:)=108080.  ! Earth aqua/terra planets
+        else
+           ps(:)=preff
+        endif
         ! ground geopotential
         phis(:)=0.
Index: trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F	(revision 270)
@@ -47,5 +47,4 @@
 c-----------------------------------------------------------------------
 
-      dtdiss  = idissip * dtvr
       dtphys  = iphysiq * dtvr
       unsim   = 1./iim
Index: trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F	(revision 134)
+++ 	(revision )
@@ -1,244 +1,0 @@
-!
-! $Id: inidissip.F 1403 2010-07-01 09:02:53Z fairhead $
-!
-      SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,
-     *                       tetagdiv,tetagrot,tetatemp             )
-c=======================================================================
-c   initialisation de la dissipation horizontale
-c=======================================================================
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-      USE control_mod
-
-      IMPLICIT NONE
-#include "dimensions.h"
-#include "paramet.h"
-#include "comdissipn.h"
-#include "comconst.h"
-#include "comvert.h"
-#include "logic.h"
-
-      LOGICAL lstardis
-      INTEGER nitergdiv,nitergrot,niterh
-      REAL    tetagdiv,tetagrot,tetatemp
-      REAL fact,zvert(llm),zz
-      REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
-      REAL ullm,vllm,umin,vmin,zhmin,zhmax
-      REAL zllm,z1llm
-
-      INTEGER l,ij,idum,ii
-      REAL tetamin
-      REAL Pup
-
-      REAL ran1
-
-
-c-----------------------------------------------------------------------
-c
-c   calcul des valeurs propres des operateurs par methode iterrative:
-c   -----------------------------------------------------------------
-
-      crot     = -1.
-      cdivu    = -1.
-      cdivh    = -1.
-
-c   calcul de la valeur propre de divgrad:
-c   --------------------------------------
-      idum = 0
-      DO l = 1, llm
-       DO ij = 1, ip1jmp1
-        deltap(ij,l) = 1.
-       ENDDO
-      ENDDO
-
-      idum  = -1
-      zh(1) = RAN1(idum)-.5
-      idum  = 0
-      DO ij = 2, ip1jmp1
-        zh(ij) = RAN1(idum) -.5
-      ENDDO
-
-      CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
-
-      CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
-
-      IF ( zhmin .GE. zhmax  )     THEN
-         PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
-         STOP'probleme generateur alleatoire dans inidissip'
-      ENDIF
-
-      zllm = ABS( zhmax )
-      DO l = 1,50
-         IF(lstardis) THEN
-            CALL divgrad2(1,zh,deltap,niterh,zh)
-         ELSE
-            CALL divgrad (1,zh,niterh,zh)
-         ENDIF
-
-        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
-
-         zllm  = ABS( zhmax )
-         z1llm = 1./zllm
-         DO ij = 1,ip1jmp1
-            zh(ij) = zh(ij)* z1llm
-         ENDDO
-      ENDDO
-
-      IF(lstardis) THEN
-         cdivh = 1./ zllm
-      ELSE
-         cdivh = zllm ** ( -1./niterh )
-      ENDIF
-
-c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
-c   -----------------------------------------------------------------
-      print*,'calcul des valeurs propres'
-
-      DO  20  ii = 1, 2
-c
-         DO ij = 1, ip1jmp1
-           zu(ij)  = RAN1(idum) -.5
-         ENDDO
-         CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
-         DO ij = 1, ip1jm
-            zv(ij) = RAN1(idum) -.5
-         ENDDO
-         CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
-
-         CALL minmax(iip1*jjp1,zu,umin,ullm )
-         CALL minmax(iip1*jjm, zv,vmin,vllm )
-
-         ullm = ABS ( ullm )
-         vllm = ABS ( vllm )
-
-         DO  5  l = 1, 50
-            IF(ii.EQ.1) THEN
-ccccc             CALL covcont( 1,zu,zv,zu,zv )
-               IF(lstardis) THEN
-                  CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
-               ELSE
-                  CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
-               ENDIF
-            ELSE
-               IF(lstardis) THEN
-                  CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
-               ELSE
-                  CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
-               ENDIF
-            ENDIF
-
-            CALL minmax(iip1*jjp1,zu,umin,ullm )
-            CALL minmax(iip1*jjm, zv,vmin,vllm )
-
-            ullm = ABS  ( ullm )
-            vllm = ABS  ( vllm )
-
-            zllm  = MAX( ullm,vllm )
-            z1llm = 1./ zllm
-            DO ij = 1, ip1jmp1
-              zu(ij) = zu(ij)* z1llm
-            ENDDO
-            DO ij = 1, ip1jm
-               zv(ij) = zv(ij)* z1llm
-            ENDDO
- 5       CONTINUE
-
-         IF ( ii.EQ.1 ) THEN
-            IF(lstardis) THEN
-               cdivu  = 1./zllm
-            ELSE
-               cdivu  = zllm **( -1./nitergdiv )
-            ENDIF
-         ELSE
-            IF(lstardis) THEN
-               crot   = 1./ zllm
-            ELSE
-               crot   = zllm **( -1./nitergrot )
-            ENDIF
-         ENDIF
-
- 20   CONTINUE
-
-c   petit test pour les operateurs non star:
-c   ----------------------------------------
-
-c     IF(.NOT.lstardis) THEN
-         fact    = rad*24./REAL(jjm)
-         fact    = fact*fact
-         PRINT*,'coef u ', fact/cdivu, 1./cdivu
-         PRINT*,'coef r ', fact/crot , 1./crot
-         PRINT*,'coef h ', fact/cdivh, 1./cdivh
-c     ENDIF
-
-c-----------------------------------------------------------------------
-c   variation verticale du coefficient de dissipation:
-c   --------------------------------------------------
-
-c First step: going from 1 to dissip_fac_mid (in gcm.def)
-c============
-      DO l=1,llm
-         zz      = 1. - preff/presnivs(l)
-         zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
-      ENDDO
-
-         write(*,*) 'Dissipation : '
-         write(*,*) 'Multiplication de la dissipation en altitude :'
-         write(*,*) '  dissip_fac_mid =', dissip_fac_mid
-
-c Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
-c==========================
-c Utilisation de la fonction tangente hyperbolique pour augmenter
-c arbitrairement la dissipation et donc la stabilite du modele a 
-c partir d'une certaine altitude.
-
-c   Le facteur multiplicatif de basse atmosphere etant deja pris 
-c   en compte, il faut diviser le facteur multiplicatif de haute 
-c   atmosphere par celui-ci.
-
-      if (ok_strato) then
-
-       Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
-       do l=1,llm
-         zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1)
-     &   *(1-(0.5*(1+tanh(-6./dissip_deltaz*
-     &    (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
-       enddo 
-
-         write(*,*) '  dissip_fac_up =', dissip_fac_up
-         write(*,*) 'Transition mid /up:  Pupstart,delta =',
-     &             dissip_pupstart,'Pa', dissip_deltaz , '(km)'
-
-      endif
-
-
-      PRINT*,'Constantes de temps de la diffusion horizontale'
-
-      tetamin =  1.e+6
-
-      DO l=1,llm
-        tetaudiv(l)   = zvert(l)/tetagdiv
-        tetaurot(l)   = zvert(l)/tetagrot
-        tetah(l)      = zvert(l)/tetatemp
-
-        IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
-        IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
-        IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
-      ENDDO
-
-      PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
-      idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
-      PRINT *,' INIDI tetamin idissip ',tetamin,idissip
-      idissip = MAX(iperiod,idissip)
-      dtdiss  = idissip * dtvr
-      PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
-
-      DO l = 1,llm
-         PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),
-     *                   dtdiss*tetah(l)
-      ENDDO
-
-c
-      RETURN
-      END
Index: trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90	(revision 270)
+++ trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90	(revision 270)
@@ -0,0 +1,250 @@
+!
+! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $
+!
+SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
+     tetagdiv,tetagrot,tetatemp             )
+  !=======================================================================
+  !   initialisation de la dissipation horizontale
+  !=======================================================================
+  !-----------------------------------------------------------------------
+  !   declarations:
+  !   -------------
+
+  USE control_mod, only : dissip_period,iperiod
+
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comdissipn.h"
+  include "comconst.h"
+  include "comvert.h"
+  include "logic.h"
+  include "iniprint.h"
+
+  LOGICAL,INTENT(in) :: lstardis
+  INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
+  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
+
+! Local variables:
+  REAL fact,zvert(llm),zz
+  REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
+  REAL ullm,vllm,umin,vmin,zhmin,zhmax
+  REAL zllm,z1llm
+
+  INTEGER l,ij,idum,ii
+  REAL tetamin
+  REAL Pup
+  character (len=80) :: abort_message
+
+  REAL ran1
+
+
+  !-----------------------------------------------------------------------
+  !
+  !   calcul des valeurs propres des operateurs par methode iterrative:
+  !   -----------------------------------------------------------------
+
+  crot     = -1.
+  cdivu    = -1.
+  cdivh    = -1.
+
+  !   calcul de la valeur propre de divgrad:
+  !   --------------------------------------
+  idum = 0
+  DO l = 1, llm
+     DO ij = 1, ip1jmp1
+        deltap(ij,l) = 1.
+     ENDDO
+  ENDDO
+
+  idum  = -1
+  zh(1) = RAN1(idum)-.5
+  idum  = 0
+  DO ij = 2, ip1jmp1
+     zh(ij) = RAN1(idum) -.5
+  ENDDO
+
+  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
+
+  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
+
+  IF ( zhmin .GE. zhmax  )     THEN
+     write(lunout,*)'  Inidissip  zh min max  ',zhmin,zhmax
+     abort_message='probleme generateur alleatoire dans inidissip'
+     call abort_gcm('inidissip',abort_message,1)
+  ENDIF
+
+  zllm = ABS( zhmax )
+  DO l = 1,50
+     IF(lstardis) THEN
+        CALL divgrad2(1,zh,deltap,niterh,zh)
+     ELSE
+        CALL divgrad (1,zh,niterh,zh)
+     ENDIF
+
+     CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
+
+     zllm  = ABS( zhmax )
+     z1llm = 1./zllm
+     DO ij = 1,ip1jmp1
+        zh(ij) = zh(ij)* z1llm
+     ENDDO
+  ENDDO
+
+  IF(lstardis) THEN
+     cdivh = 1./ zllm
+  ELSE
+     cdivh = zllm ** ( -1./niterh )
+  ENDIF
+
+  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
+  !   -----------------------------------------------------------------
+  write(lunout,*)'inidissip: calcul des valeurs propres'
+
+  DO    ii = 1, 2
+     !
+     DO ij = 1, ip1jmp1
+        zu(ij)  = RAN1(idum) -.5
+     ENDDO
+     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
+     DO ij = 1, ip1jm
+        zv(ij) = RAN1(idum) -.5
+     ENDDO
+     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
+
+     CALL minmax(iip1*jjp1,zu,umin,ullm )
+     CALL minmax(iip1*jjm, zv,vmin,vllm )
+
+     ullm = ABS ( ullm )
+     vllm = ABS ( vllm )
+
+     DO    l = 1, 50
+        IF(ii.EQ.1) THEN
+           !cccc             CALL covcont( 1,zu,zv,zu,zv )
+           IF(lstardis) THEN
+              CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
+           ELSE
+              CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
+           ENDIF
+        ELSE
+           IF(lstardis) THEN
+              CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
+           ELSE
+              CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
+           ENDIF
+        ENDIF
+
+        CALL minmax(iip1*jjp1,zu,umin,ullm )
+        CALL minmax(iip1*jjm, zv,vmin,vllm )
+
+        ullm = ABS  ( ullm )
+        vllm = ABS  ( vllm )
+
+        zllm  = MAX( ullm,vllm )
+        z1llm = 1./ zllm
+        DO ij = 1, ip1jmp1
+           zu(ij) = zu(ij)* z1llm
+        ENDDO
+        DO ij = 1, ip1jm
+           zv(ij) = zv(ij)* z1llm
+        ENDDO
+     end DO
+
+     IF ( ii.EQ.1 ) THEN
+        IF(lstardis) THEN
+           cdivu  = 1./zllm
+        ELSE
+           cdivu  = zllm **( -1./nitergdiv )
+        ENDIF
+     ELSE
+        IF(lstardis) THEN
+           crot   = 1./ zllm
+        ELSE
+           crot   = zllm **( -1./nitergrot )
+        ENDIF
+     ENDIF
+
+  end DO
+
+  !   petit test pour les operateurs non star:
+  !   ----------------------------------------
+
+  !     IF(.NOT.lstardis) THEN
+  fact    = rad*24./REAL(jjm)
+  fact    = fact*fact
+  write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
+  write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
+  write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
+  !     ENDIF
+
+  !-----------------------------------------------------------------------
+  !   variation verticale du coefficient de dissipation:
+  !   --------------------------------------------------
+
+! First step: going from 1 to dissip_fac_mid (in gcm.def)
+!============
+  DO l=1,llm
+     zz      = 1. - preff/presnivs(l)
+     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
+  ENDDO
+
+  write(lunout,*) 'Dissipation : '
+  write(lunout,*) 'Multiplication de la dissipation en altitude :'
+  write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
+
+! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
+!==========================
+! Utilisation de la fonction tangente hyperbolique pour augmenter
+! arbitrairement la dissipation et donc la stabilite du modele a 
+! partir d'une certaine altitude.
+
+!   Le facteur multiplicatif de basse atmosphere etant deja pris 
+!   en compte, il faut diviser le facteur multiplicatif de haute 
+!   atmosphere par celui-ci.
+
+  if (ok_strato) then
+
+    Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
+    do l=1,llm
+      zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
+                *(1-(0.5*(1+tanh(-6./dissip_deltaz*              &
+               (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
+    enddo 
+
+    write(*,*) '  dissip_fac_up =', dissip_fac_up
+    write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
+                   dissip_pupstart,'Pa', dissip_deltaz , '(km)'
+
+  endif
+
+
+  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
+
+  tetamin =  1.e+6
+
+  DO l=1,llm
+     tetaudiv(l)   = zvert(l)/tetagdiv
+     tetaurot(l)   = zvert(l)/tetagrot
+     tetah(l)      = zvert(l)/tetatemp
+
+     IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
+     IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
+     IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
+  ENDDO
+
+  ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
+  IF (dissip_period == 0) THEN
+     dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
+     write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period
+     dissip_period = MAX(iperiod,dissip_period)
+  END IF
+
+  dtdiss  = dissip_period * dtvr
+  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
+
+  DO l = 1,llm
+     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
+          dtdiss*tetah(l)
+  ENDDO
+
+END SUBROUTINE inidissip
Index: trunk/LMDZ.COMMON/libf/dyn3d/integrd.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/integrd.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/integrd.F	(revision 270)
@@ -1,4 +1,4 @@
 !
-! $Id: integrd.F 1446 2010-10-22 09:27:25Z emillour $
+! $Id: integrd.F 1550 2011-07-05 09:44:55Z lguez $
 !
       SUBROUTINE integrd
@@ -89,5 +89,6 @@
         IF( ps(ij).LT.0. ) THEN
          PRINT*,' Au point ij = ',ij, ' , pression sol neg. ', ps(ij)
-         STOP' dans integrd'
+         print *, ' dans integrd'
+         stop 1
         ENDIF
       ENDDO
Index: trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 270)
@@ -152,4 +152,5 @@
       logical ok_sync
       parameter (ok_sync = .true.) 
+      logical physic
 
       data callinigrads/.true./
@@ -222,4 +223,7 @@
 
       itau = 0
+      physic=.true.
+      if (iflag_phys==0.or.iflag_phys==2) physic=.false.
+
 c      iday = day_ini+itau/day_step
 c      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
@@ -322,12 +326,14 @@
       ! Purely Matsuno time stepping
          IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
-         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
+         IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward ) 
+     s        apdiss = .TRUE.
          IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 
-     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
+     s          .and. physic                        ) apphys = .TRUE.
       ELSE
       ! Leapfrog/Matsuno time stepping 
          IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
-         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
-         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
+         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
+     s        apdiss = .TRUE.
+         IF( MOD(itau+1,iphysiq).EQ.0.AND.physic       ) apphys=.TRUE.
       END IF
 
@@ -357,6 +363,7 @@
 c   -------------------------------------------------------------
 
-      IF( forward. OR . leapf )  THEN
-
+!      IF( forward. OR . leapf )  THEN
+      IF((.not.forward).OR. leapf )  THEN
+        ! Ehouarn: gather mass fluxes during backward Matsuno or LF step
          CALL caladvtrac(q,pbaru,pbarv,
      *        p, masse, dq,  teta,
@@ -497,19 +504,11 @@
           ENDDO
         ENDDO ! of DO l=1,llm 
-          call friction(ucov,vcov,dtvr)
     
-       if (planet_type.eq."giant") then
-          ! Intrinsic heat flux
-          ! Aymeric -- for giant planets
-          if (ihf .gt. 1.e-6) then
-          !print *, '**** INTRINSIC HEAT FLUX ****', ihf 
-            DO ij=1,ip1jmp1
-              teta(ij,1) = teta(ij,1)
-     &        + dtvr * aire(ij) * ihf / cpp / masse(ij,1)
-            ENDDO
-          !print *, '**** d teta ' 
-          !print *, dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
-          endif
-       endif
+        if (planet_type.eq."giant") then
+          ! add an intrinsic heat flux at the base of the atmosphere
+          teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
+        endif
+
+        call friction(ucov,vcov,dtvr)
 
     
Index: trunk/LMDZ.COMMON/libf/dyn3d/limit_netcdf.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/limit_netcdf.F90	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/limit_netcdf.F90	(revision 270)
@@ -1,4 +1,4 @@
 !
-! $Id: limit_netcdf.F90 1441 2010-10-13 13:06:56Z emillour $
+! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $
 !-------------------------------------------------------------------------------
 !
@@ -42,5 +42,5 @@
   REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
 #ifndef CPP_EARTH
-  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
+  CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1)
 #else
 !-------------------------------------------------------------------------------
@@ -52,18 +52,15 @@
 #include "indicesol.h"
 
-!--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------
-  LOGICAL, PARAMETER :: fracterre=.TRUE.
-
 !--- INPUT NETCDF FILES NAMES --------------------------------------------------
   CHARACTER(LEN=25) :: icefile, sstfile, dumstr
   CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc        ',        &
                                   famipsic='amipbc_sic_1x1.nc        ',        &
-                                  fclimsst='amipbc_sst_1x1_clim.nc   ',        &
-                                  fclimsic='amipbc_sic_1x1_clim.nc   ',        &
                                   fcpldsst='cpl_atm_sst.nc           ',        &
                                   fcpldsic='cpl_atm_sic.nc           ',        &
+                                  fhistsst='histmth_sst.nc           ',        &
+                                  fhistsic='histmth_sic.nc           ',        &
                                   frugo   ='Rugos.nc                 ',        &
                                   falbe   ='Albedo.nc                '
-
+  CHARACTER(LEN=10) :: varname
 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
   REAL,   DIMENSION(klon)                :: fi_ice, verif
@@ -80,5 +77,4 @@
   INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC
   INTEGER :: NF90_FORMAT
-  LOGICAL :: lCPL                    !--- T: IPCC-IPSL cpl model output files
   INTEGER :: ndays                   !--- Depending on the output calendar
 
@@ -104,110 +100,97 @@
 
 !--- RUGOSITY TREATMENT --------------------------------------------------------
-  WRITE(lunout,*) 'Traitement de la rugosite'
-  CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
+  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite'
+  varname='RUGOS'
+  CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
 
 !--- OCEAN TREATMENT -----------------------------------------------------------
-  PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.
+  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique'
 
 ! Input SIC file selection
-  icefile='fake'
-  IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic)
-  IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic)
-  IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic)
-  SELECT CASE(icefile)
-    CASE(famipsic); dumstr='Amip.'
-    CASE(fclimsic); dumstr='Amip climatologique.'
-    CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
-    CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1)
-  END SELECT
+! Open file only to test if available
+  IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     icefile=TRIM(famipsic)
+     varname='sicbcs'
+  ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     icefile=TRIM(fcpldsic)
+     varname='SIICECOV'
+  ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     icefile=TRIM(fhistsic)
+     varname='pourc_sic'
+  ELSE
+     WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
+     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic)
+     CALL abort_gcm('limit_netcdf','No sea-ice file was found',1)
+  END IF
   ierr=NF90_CLOSE(nid)
-  WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)
-
-  CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL)
+  IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile)
+
+  CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
 
   ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
   DO k=1,ndays
-    fi_ice=phy_ice(:,k)
-    WHERE(fi_ice>=1.0  ) fi_ice=1.0
-    WHERE(fi_ice<EPSFRA) fi_ice=0.0
-    IF(fracterre) THEN
-      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
-      pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
-      IF(lCPL) THEN                               ! SIC=pICE*(1-LIC-TER)
-        pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter))
-      ELSE                                        ! SIC=pICE-LIC
+     fi_ice=phy_ice(:,k)
+     WHERE(fi_ice>=1.0  ) fi_ice=1.0
+     WHERE(fi_ice<EPSFRA) fi_ice=0.0
+     pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
+     pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
+     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
+        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
+     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
+        pctsrf_t(:,is_sic,k)=fi_ice(:)
+     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
         pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
-      END IF
-      WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
-      WHERE(1.0-zmasq<EPSFRA)
+     END IF
+     WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
+     WHERE(1.0-zmasq<EPSFRA)
         pctsrf_t(:,is_sic,k)=0.0
         pctsrf_t(:,is_oce,k)=0.0
-      ELSEWHERE
+     ELSEWHERE
         WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
-          pctsrf_t(:,is_sic,k)=1.0-zmasq
-          pctsrf_t(:,is_oce,k)=0.0
+           pctsrf_t(:,is_sic,k)=1.0-zmasq
+           pctsrf_t(:,is_oce,k)=0.0
         ELSEWHERE
-          pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
-          WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
-            pctsrf_t(:,is_oce,k)=0.0
-            pctsrf_t(:,is_sic,k)=1.0-zmasq
-          END WHERE
+           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
+           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
+              pctsrf_t(:,is_oce,k)=0.0
+              pctsrf_t(:,is_sic,k)=1.0-zmasq
+           END WHERE
         END WHERE
-      END WHERE
-      nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
-      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
-      nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
-      IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
-    ELSE 
-      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)
-      WHERE(NINT(pctsrf(:,is_ter))==1)
-        pctsrf_t(:,is_sic,k)=0.
-        pctsrf_t(:,is_oce,k)=0.                  
-        WHERE(fi_ice>=1.e-5)
-          pctsrf_t(:,is_lic,k)=fi_ice
-          pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k)
-        ELSEWHERE
-          pctsrf_t(:,is_lic,k)=0.0
-        END WHERE
-      ELSEWHERE
-        pctsrf_t(:,is_lic,k) = 0.0 
-        WHERE(fi_ice>=1.e-5)
-          pctsrf_t(:,is_ter,k)=0.0
-          pctsrf_t(:,is_sic,k)=fi_ice
-          pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k)
-        ELSEWHERE
-          pctsrf_t(:,is_sic,k)=0.0
-          pctsrf_t(:,is_oce,k)=1.0
-        END WHERE
-      END WHERE
-      verif=sum(pctsrf_t(:,:,k),dim=2)
-      nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5)
-      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
-    END IF 
+     END WHERE
+     nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
+     IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
+     nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
+     IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
   END DO
   DEALLOCATE(phy_ice)
 
 !--- SST TREATMENT -------------------------------------------------------------
-  WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.
+  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst'
 
 ! Input SST file selection
-  sstfile='fake'
-  IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst)
-  IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst)
-  IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst)
-  SELECT CASE(icefile)
-    CASE(famipsic); dumstr='Amip.'
-    CASE(fclimsic); dumstr='Amip climatologique.'
-    CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
-    CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1)
-  END SELECT
+! Open file only to test if available
+  IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     sstfile=TRIM(famipsst)
+     varname='tosbcs'
+  ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     sstfile=TRIM(fcpldsst)
+     varname='SISUTESW'
+  ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     sstfile=TRIM(fhistsst)
+     varname='tsol_oce'
+  ELSE
+     WRITE(lunout,*) 'ERROR! No sst input file was found.'
+     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst)
+     CALL abort_gcm('limit_netcdf','No sst file was found',1)
+  END IF
   ierr=NF90_CLOSE(nid)
-  WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)
-
-  CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)
+  IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile)
+
+  CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
 
 !--- ALBEDO TREATMENT ----------------------------------------------------------
-  WRITE(lunout,*) 'Traitement de l albedo'
-  CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb)
+  IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo'
+  varname='ALBEDO'
+  CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb)
 
 !--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
@@ -215,5 +198,5 @@
 
 !--- OUTPUT FILE WRITING -------------------------------------------------------
-  WRITE(lunout,*) 'Ecriture du fichier limit'
+  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut'
 
   !--- File creation
@@ -264,4 +247,6 @@
   ierr=NF90_CLOSE(nid)
 
+  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin'
+
   DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
 
@@ -276,5 +261,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
+SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask)
 !
 !-----------------------------------------------------------------------------
@@ -304,11 +289,11 @@
 ! Arguments:
   CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
+  CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
   CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
   LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
   INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
-  REAL,    POINTER,  DIMENSION(:, :) :: champo   ! output field = f(t)
+  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
   LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
   REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
-  LOGICAL, OPTIONAL, INTENT(IN)     :: lCPL     ! Coupled model flag (for ICE)
 !------------------------------------------------------------------------------
 ! Local variables:
@@ -316,5 +301,4 @@
   INTEGER :: ncid, varid                  ! NetCDF identifiers
   CHARACTER(LEN=30)               :: dnam       ! dimension name
-  CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
 !--- dimensions
   INTEGER,           DIMENSION(4) :: dids       ! NetCDF dimensions identifiers
@@ -331,4 +315,5 @@
 !--- input files
   CHARACTER(LEN=20)                 :: cal_in   ! calendar
+  CHARACTER(LEN=20)                 :: unit_sic ! attribute unit in sea-ice file
   INTEGER                           :: ndays_in ! number of days
 !--- misc
@@ -337,26 +322,45 @@
   CHARACTER(LEN=25)                 :: title    ! for messages
   LOGICAL                           :: extrp    ! flag for extrapolation
+  LOGICAL                           :: oldice   ! flag for old way ice computation 
   REAL                              :: chmin, chmax
   INTEGER ierr
   integer n_extrap ! number of extrapolated points
   logical skip
+
 !------------------------------------------------------------------------------
 !---Variables depending on keyword 'mode' -------------------------------------
   NULLIFY(champo)
+
   SELECT CASE(mode)
-    CASE('RUG'); varname='RUGOS';  title='Rugosite'
-    CASE('SIC'); varname='sicbcs'; title='Sea-ice'
-    CASE('SST'); varname='tosbcs'; title='SST'
-    CASE('ALB'); varname='ALBEDO'; title='Albedo'
+  CASE('RUG'); title='Rugosite'
+  CASE('SIC'); title='Sea-ice'
+  CASE('SST'); title='SST'
+  CASE('ALB'); title='Albedo'
   END SELECT
+  
+
   extrp=.FALSE. 
+  oldice=.FALSE.
   IF ( PRESENT(flag) ) THEN 
     IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 
+    IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. 
   END IF
 
 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
+  IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam
   ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
-  ierr=NF90_INQ_VARID(ncid, varname, varid);            CALL ncerr(ierr, fnam)
+  ierr=NF90_INQ_VARID(ncid, trim(varname), varid);            CALL ncerr(ierr, fnam)
   ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
+
+!--- Read unit for sea-ice variable only
+  IF (mode=='SIC') THEN
+     ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic)
+     IF(ierr/=NF90_NOERR) THEN
+        IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value'
+        unit_sic='X'
+     ELSE
+        IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic
+     END IF
+  END IF
 
 !--- Longitude
@@ -365,5 +369,5 @@
   ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
   ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
-  WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
+  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
 
 !--- Latitude
@@ -372,5 +376,5 @@
   ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
   ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
-  WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
+  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
 
 !--- Time (variable is not needed - it is rebuilt - but calendar is)
@@ -385,10 +389,11 @@
       CASE('SIC', 'SST'); cal_in='gregorian'
     END SELECT
-    WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
+    IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
          // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
   END IF
-  WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
+  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
        cal_in
 
+  
 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
   !--- Determining input file number of days, depending on calendar
@@ -398,7 +403,6 @@
 !--- If input records are not monthly, time sampling has to be constant !
   timeyear=mid_months(anneeref, cal_in, lmdep)
-  IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' &
-       // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, &
-       ' enregistrements.'
+  IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
+       ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
 
 !--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
@@ -406,7 +410,6 @@
   IF(extrp) ALLOCATE(work(imdep, jmdep))
 
-  WRITE(lunout, *)
-  WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, &
-       ' CHAMPS.'
+  IF (prt_level>5) WRITE(lunout, *)
+  IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.'
   ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
   DO l=1, lmdep
@@ -419,12 +422,9 @@
          work)
 
-    IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
-      IF(l==1) THEN
-        WRITE(lunout, *)                                                      &
-  '-------------------------------------------------------------------------'
-        WRITE(lunout, *)                                                     &
-  'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
-        WRITE(lunout, *)                                                      &
-  '-------------------------------------------------------------------------'
+    IF(ibar .AND. .NOT.oldice) THEN
+      IF(l==1 .AND. prt_level>5) THEN
+        WRITE(lunout, *) '-------------------------------------------------------------------------'
+        WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$'
+        WRITE(lunout, *) '-------------------------------------------------------------------------'
       END IF
       IF(mode=='RUG') champ=LOG(champ)
@@ -453,8 +453,11 @@
 
 !--- TIME INTERPOLATION ------------------------------------------------------
-  WRITE(lunout, *)
-  WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
-  WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear
-  WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays
+  IF (prt_level>5) THEN
+     WRITE(lunout, *)
+     WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
+     WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
+     WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
+  END IF
+
   ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
   skip = .false.
@@ -471,5 +474,5 @@
   END DO
   if (n_extrap /= 0) then
-     print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap
+     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
   end if
   champan(iip1, :, :)=champan(1, :, :)
@@ -479,10 +482,10 @@
   DO j=1, jjp1
     CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
-    WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j
+    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
   END DO
 
 !--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
   IF(mode=='SST') THEN
-    WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
+    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
     WHERE(champan<271.38) champan=271.38
   END IF
@@ -490,10 +493,20 @@
 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
   IF(mode=='SIC') THEN
-    WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
-    IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100.
+    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
+
+    IF (unit_sic=='1') THEN
+       ! Nothing to be done for sea-ice field is already in fraction of 1
+       ! This is the case for sea-ice in file cpl_atm_sic.nc
+       IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1'
+    ELSE
+       ! Convert sea ice from percentage to fraction of 1
+       IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' 
+       champan(:, :, :)=champan(:, :, :)/100.
+    END IF
+
     champan(iip1, :, :)=champan(1, :, :)
     WHERE(champan>1.0) champan=1.0
     WHERE(champan<0.0) champan=0.0
-  END IF
+ END IF
 
 !--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3d/vlsplt.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/vlsplt.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/vlsplt.F	(revision 270)
@@ -1,6 +1,4 @@
-!
-! $Header$
-!
-c
+c
+c $Id: vlsplt.F 1550 2011-07-05 09:44:55Z lguez $
 c
 
@@ -131,9 +129,10 @@
       IMPLICIT NONE
 c
-#include "dimensions.h"
-#include "paramet.h"
-#include "logic.h"
-#include "comvert.h"
-#include "comconst.h"
+      include "dimensions.h"
+      include "paramet.h"
+      include "logic.h"
+      include "comvert.h"
+      include "comconst.h"
+      include "iniprint.h"
 c
 c
@@ -353,5 +352,6 @@
 
       IF(n0.gt.0) THEN
-      PRINT*,'Nombre de points pour lesquels on advect plus que le'
+      if (prt_level > 2) PRINT *,
+     $        'Nombre de points pour lesquels on advect plus que le'
      &       ,'contenu de la maille : ',n0
 
Index: trunk/LMDZ.COMMON/libf/dyn3d/write_paramLMDZ_dyn.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/write_paramLMDZ_dyn.h	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3d/write_paramLMDZ_dyn.h	(revision 270)
@@ -51,6 +51,6 @@
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
-      zx_tmp_2d(1:iip1,1:jjp1)=REAL(idissip)
-      CALL histwrite(nid_ctesGCM, "idissip", itau_w,
+      zx_tmp_2d(1:iip1,1:jjp1)=REAL(dissip_period)
+      CALL histwrite(nid_ctesGCM, "dissip_period", itau_w,
      .               zx_tmp_2d,iip1*jjp1,ndex2d)
 c
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F	(revision 134)
+++ 	(revision )
@@ -1,505 +1,0 @@
-!
-! $Id: advtrac_p.F 1446 2010-10-22 09:27:25Z emillour $
-!
-c
-c
-      SUBROUTINE advtrac_p(pbaru,pbarv ,
-     *                   p,  masse,q,iapptrac,teta,
-     *                  flxw,
-     *                  pk   )
-
-c     Auteur :  F. Hourdin
-c
-c     Modif. P. Le Van     (20/12/97)
-c            F. Codron     (10/99)
-c            D. Le Croller (07/2001)
-c            M.A Filiberti (04/2002)
-c
-      USE parallel
-      USE Write_Field_p
-      USE Bands
-      USE mod_hallo
-      USE Vampir
-      USE times
-      USE infotrac
-      USE control_mod
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comvert.h"
-#include "comdissip.h"
-#include "comgeom2.h"
-#include "logic.h"
-#include "temps.h"
-#include "ener.h"
-#include "description.h"
-
-c-------------------------------------------------------------------
-c     Arguments
-c-------------------------------------------------------------------
-c     Ajout PPM
-c--------------------------------------------------------
-      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
-c--------------------------------------------------------
-      INTEGER iapptrac
-      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
-      REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
-      REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
-      REAL pk(ip1jmp1,llm)
-      REAL               :: flxw(ip1jmp1,llm)
-
-c-------------------------------------------------------------
-c     Variables locales
-c-------------------------------------------------------------
-
-      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
-      REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
-      REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 
-      REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
-      INTEGER iadvtr
-      INTEGER ij,l,iq,iiq
-      REAL zdpmin, zdpmax
-      SAVE iadvtr, massem, pbaruc, pbarvc
-      DATA iadvtr/0/
-c$OMP THREADPRIVATE(iadvtr)
-c----------------------------------------------------------
-c     Rajouts pour PPM
-c----------------------------------------------------------
-      INTEGER indice,n
-      REAL dtbon ! Pas de temps adaptatif pour que CFL<1
-      REAL CFLmaxz,aaa,bbb ! CFL maximum
-      REAL psppm(iim,jjp1) ! pression  au sol
-      REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
-      REAL qppm(iim*jjp1,llm,nqtot)
-      REAL fluxwppm(iim,jjp1,llm)
-      REAL apppm(llmp1), bpppm(llmp1)
-      LOGICAL dum,fill
-      DATA fill/.true./
-      DATA dum/.true./
-      REAL,SAVE :: finmasse(ip1jmp1,llm)
-      integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
-      type(Request) :: Request_vanleer
-      REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
-      REAL,SAVE :: teta_tmp(ip1jmp1,llm)
-      REAL,SAVE :: pk_tmp(ip1jmp1,llm)
-
-      ijb_u=ij_begin
-      ije_u=ij_end
-      
-      ijb_v=ij_begin-iip1
-      ije_v=ij_end
-      if (pole_nord) ijb_v=ij_begin
-      if (pole_sud)  ije_v=ij_end-iip1
-
-      IF(iadvtr.EQ.0) THEN
-c         CALL initial0(ijp1llm,pbaruc)
-c         CALL initial0(ijmllm,pbarvc)
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)        
-	DO l=1,llm   
-          pbaruc(ijb_u:ije_u,l)=0.
-          pbarvc(ijb_v:ije_v,l)=0.
-        ENDDO
-c$OMP END DO NOWAIT  
-      ENDIF
-
-c   accumulation des flux de masse horizontaux
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-         DO ij = ijb_u,ije_u
-            pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
-         ENDDO
-         DO ij = ijb_v,ije_v
-            pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
-         ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
-
-c   selection de la masse instantannee des mailles avant le transport.
-      IF(iadvtr.EQ.0) THEN
-
-c         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
-          ijb=ij_begin
-          ije=ij_end
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-       DO l=1,llm
-          massem(ijb:ije,l)=masse(ijb:ije,l)
-       ENDDO
-c$OMP END DO NOWAIT
-
-ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
-c
-      ENDIF ! of IF(iadvtr.EQ.0)
-
-      iadvtr   = iadvtr+1
-
-c$OMP MASTER
-      iapptrac = iadvtr
-c$OMP END MASTER
-
-c   Test pour savoir si on advecte a ce pas de temps
-
-      IF ( iadvtr.EQ.iapp_tracvl ) THEN
-c$OMP MASTER
-        call suspend_timer(timer_caldyn)
-c$OMP END MASTER
-      
-      ijb=ij_begin
-      ije=ij_end
-      
-
-cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
-cc
-
-c   traitement des flux de masse avant advection.
-c     1. calcul de w
-c     2. groupement des mailles pres du pole.
-
-        CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
-
-c$OMP BARRIER
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llmp1
-        p_tmp(ijb:ije,l)=p(ijb:ije,l)
-      ENDDO
-c$OMP END DO NOWAIT
-      
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-        pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
-        teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
-      ENDDO
-c$OMP END DO NOWAIT
-
-c$OMP MASTER
-      call VTb(VTHallo)
-c$OMP END MASTER
-
-      call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm,
-     *                             jj_Nb_vanleer,0,0,Request_vanleer)
-      call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm,
-     *                             jj_Nb_vanleer,1,0,Request_vanleer)
-      call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm,
-     *                             jj_Nb_vanleer,0,0,Request_vanleer)
-      call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm,
-     *                             jj_Nb_vanleer,0,0,Request_vanleer)
-      call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm,
-     *                             jj_Nb_vanleer,1,1,Request_vanleer)
-      call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1,
-     *                             jj_Nb_vanleer,1,1,Request_vanleer)
-      call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,
-     *                             jj_Nb_vanleer,1,1,Request_vanleer)
-      do j=1,nqtot
-        call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
-     *                             jj_nb_vanleer,0,0,Request_vanleer)
-      enddo
-
-      call SendRequest(Request_vanleer)
-c$OMP BARRIER
-      call WaitRequest(Request_vanleer)
-
-
-c$OMP BARRIER
-c$OMP MASTER      
-      call SetDistrib(jj_nb_vanleer)
-      call VTe(VTHallo)
-      call VTb(VTadvection)
-      call start_timer(timer_vanleer)
-c$OMP END MASTER
-c$OMP BARRIER
-      
-      ! ... Flux de masse diaganostiques traceurs
-         ijb=ij_begin
-         ije=ij_end
-         flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
-
-c  test sur l'eventuelle creation de valeurs negatives de la masse
-         ijb=ij_begin
-         ije=ij_end
-         if (pole_nord) ijb=ij_begin+iip1
-         if (pole_sud) ije=ij_end-iip1
-         
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)          
-         DO l=1,llm-1
-            DO ij = ijb+1,ije
-              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
-     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
-     s                  +       wg(ij,l+1)  - wg(ij,l)
-            ENDDO
-            
-c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
-c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
-            
-            do ij=ijb,ije-iip1+1,iip1
-              zdp(ij)=zdp(ij+iip1-1)
-            enddo
-            
-            DO ij = ijb,ije
-               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 
-            ENDDO 
-
-
-c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
-c  ym ---> eventuellement a revoir
-            CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
-            
-            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
-            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
-     s        '   MAX:', zdpmax
-            ENDIF
-
-         ENDDO
-c$OMP END DO NOWAIT
-
-c-------------------------------------------------------------------
-c   Advection proprement dite (Modification Le Croller (07/2001)
-c-------------------------------------------------------------------
-
-c----------------------------------------------------
-c        Calcul des moyennes basées sur la masse
-c----------------------------------------------------
-
-cym      ----> Normalement, inutile pour les schémas classiques
-cym      ----> Revérifier lors de la parallélisation des autres schemas
-   
-cym          call massbar_p(massem,massebx,masseby)  
-
-          call vlspltgen_p( q,iadv, 2., massem, wg ,
-     *                    pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
-
-         
-	  GOTO 1234     
-c-----------------------------------------------------------
-c     Appel des sous programmes d'advection
-c-----------------------------------------------------------
-      do iq=1,nqtot
-c        call clock(t_initial)
-        if(iadv(iq) == 0) cycle 
-c   ----------------------------------------------------------------
-c   Schema de Van Leer I MUSCL
-c   ----------------------------------------------------------------
-        if(iadv(iq).eq.10) THEN
-     
-	    call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
-
-c   ----------------------------------------------------------------
-c   Schema "pseudo amont" + test sur humidite specifique
-C    pour la vapeur d'eau. F. Codron
-c   ----------------------------------------------------------------
-        else if(iadv(iq).eq.14) then
-c
-cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
-           CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
-     *                 pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
-c   ----------------------------------------------------------------
-c   Schema de Frederic Hourdin
-c   ----------------------------------------------------------------
-        else if(iadv(iq).eq.12) then
-          stop 'advtrac : schema non parallelise'
-c            Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
-     s             dtvr,'n=',n
-           endif
-           do indice=1,n
-            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
-           end do
-        else if(iadv(iq).eq.13) then
-          stop 'advtrac : schema non parallelise'
-c            Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
-     s             dtvr,'n=',n
-           endif
-          do indice=1,n
-            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
-          end do
-c   ----------------------------------------------------------------
-c   Schema de pente SLOPES
-c   ----------------------------------------------------------------
-        else if (iadv(iq).eq.20) then
-          stop 'advtrac : schema non parallelise'
-
-            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
-
-c   ----------------------------------------------------------------
-c   Schema de Prather
-c   ----------------------------------------------------------------
-        else if (iadv(iq).eq.30) then
-          stop 'advtrac : schema non parallelise'
-c            Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
-     s             dtvr,'n=',n
-           endif
-           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
-     s                     n,dtbon)
-c   ----------------------------------------------------------------
-c   Schemas PPM Lin et Rood
-c   ----------------------------------------------------------------
-       else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
-     s                     iadv(iq).LE.18)) then
-
-           stop 'advtrac : schema non parallelise'
-
-c        Test sur le flux horizontal
-c        Pas de temps adaptatif
-         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-         if (n.GT.1) then
-           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
-     s             dtvr,'n=',n
-         endif
-c        Test sur le flux vertical
-         CFLmaxz=0.
-         do l=2,llm
-           do ij=iip2,ip1jm
-            aaa=wg(ij,l)*dtvr/massem(ij,l)
-            CFLmaxz=max(CFLmaxz,aaa)
-            bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
-            CFLmaxz=max(CFLmaxz,bbb)
-           enddo
-         enddo
-         if (CFLmaxz.GE.1) then
-            write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
-         endif
-
-c-----------------------------------------------------------
-c        Ss-prg interface LMDZ.4->PPM3d
-c-----------------------------------------------------------
-
-          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
-     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
-     s                 unatppm,vnatppm,psppm)
-
-          do indice=1,n
-c---------------------------------------------------------------------
-c                         VL (version PPM) horiz. et PPM vert.
-c---------------------------------------------------------------------
-                if (iadv(iq).eq.11) then
-c                  Ss-prg PPM3d de Lin
-                  call ppm3d(1,qppm(1,1,iq),
-     s                       psppm,psppm,
-     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
-     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
-     s                       fill,dum,220.)
-
-c----------------------------------------------------------------------
-c                           Monotonic PPM
-c----------------------------------------------------------------------
-               else if (iadv(iq).eq.16) then
-c                  Ss-prg PPM3d de Lin
-                  call ppm3d(1,qppm(1,1,iq),
-     s                       psppm,psppm,
-     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
-     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
-     s                       fill,dum,220.)
-c---------------------------------------------------------------------
-
-c---------------------------------------------------------------------
-c                           Semi Monotonic PPM
-c---------------------------------------------------------------------
-               else if (iadv(iq).eq.17) then
-c                  Ss-prg PPM3d de Lin
-                  call ppm3d(1,qppm(1,1,iq),
-     s                       psppm,psppm,
-     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
-     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
-     s                       fill,dum,220.)
-c---------------------------------------------------------------------
-
-c---------------------------------------------------------------------
-c                         Positive Definite PPM
-c---------------------------------------------------------------------
-                else if (iadv(iq).eq.18) then
-c                  Ss-prg PPM3d de Lin
-                  call ppm3d(1,qppm(1,1,iq),
-     s                       psppm,psppm,
-     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
-     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
-     s                       fill,dum,220.)
-c---------------------------------------------------------------------
-                endif
-            enddo
-c-----------------------------------------------------------------
-c               Ss-prg interface PPM3d-LMDZ.4
-c-----------------------------------------------------------------
-                  call interpost(q(1,1,iq),qppm(1,1,iq))
-            endif
-c----------------------------------------------------------------------
-
-c-----------------------------------------------------------------
-c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
-c et Nord j=1
-c-----------------------------------------------------------------
-
-c                  call traceurpole(q(1,1,iq),massem)
-
-c calcul du temps cpu pour un schema donne
-
-c                  call clock(t_final)
-cym                  tps_cpu=t_final-t_initial
-cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
-
-       end DO
-
-1234  CONTINUE
-c$OMP BARRIER
-
-      if (planet_type=="earth") then
-
-        ijb=ij_begin
-        ije=ij_end
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-        DO l = 1, llm
-         DO ij = ijb, ije
-           finmasse(ij,l) =  p(ij,l) - p(ij,l+1) 
-         ENDDO
-        ENDDO
-c$OMP END DO
-
-        CALL qminimum_p( q, 2, finmasse )
-
-c------------------------------------------------------------------
-c   on reinitialise a zero les flux de masse cumules
-c---------------------------------------------------
-c          iadvtr=0
-
-c$OMP MASTER
-	call VTe(VTadvection)
-        call stop_timer(timer_vanleer)
-        call VTb(VThallo)
-c$OMP END MASTER
-
-	do j=1,nqtot
-          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
-     *                             jj_nb_caldyn,0,0,Request_vanleer)
-        enddo
-
-        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
-     *       jj_nb_caldyn,0,0,Request_vanleer)
-
-        call SendRequest(Request_vanleer)
-c$OMP BARRIER
-        call WaitRequest(Request_vanleer)      
-
-c$OMP BARRIER
-c$OMP MASTER
-        call SetDistrib(jj_nb_caldyn)
-	call VTe(VThallo)
-	call resume_timer(timer_caldyn)
-c$OMP END MASTER
-c$OMP BARRIER	
-          iadvtr=0
-      endif ! of if (planet_type=="earth")
-       ENDIF ! if iadvtr.EQ.iapp_tracvl
-
-       RETURN
-       END
-
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90	(revision 270)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90	(revision 270)
@@ -0,0 +1,497 @@
+! $Id: advtrac_p.F90 1549 2011-07-05 08:41:12Z lguez $
+
+SUBROUTINE advtrac_p(pbaru,pbarv , p,  masse,q,iapptrac,teta, flxw, pk)
+
+  !     Auteur :  F. Hourdin
+  !
+  !     Modif. P. Le Van     (20/12/97)
+  !            F. Codron     (10/99)
+  !            D. Le Croller (07/2001)
+  !            M.A Filiberti (04/2002)
+  !
+  USE parallel
+  USE Write_Field_p
+  USE Bands
+  USE mod_hallo
+  USE Vampir
+  USE times
+  USE infotrac
+  USE control_mod
+  IMPLICIT NONE
+  !
+  include "dimensions.h"
+  include "paramet.h"
+  include "comconst.h"
+  include "comvert.h"
+  include "comdissip.h"
+  include "comgeom2.h"
+  include "logic.h"
+  include "temps.h"
+  include "ener.h"
+  include "description.h"
+
+  !-------------------------------------------------------------------
+  !     Arguments
+  !-------------------------------------------------------------------
+  !     Ajout PPM
+  !--------------------------------------------------------
+  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
+  !--------------------------------------------------------
+  INTEGER iapptrac
+  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+  REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
+  REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
+  REAL pk(ip1jmp1,llm)
+  REAL               :: flxw(ip1jmp1,llm)
+
+  !-------------------------------------------------------------
+  !     Variables locales
+  !-------------------------------------------------------------
+
+  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
+  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
+  REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 
+  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
+  INTEGER iadvtr
+  INTEGER ij,l,iq,iiq
+  REAL zdpmin, zdpmax
+  SAVE iadvtr, massem, pbaruc, pbarvc
+  DATA iadvtr/0/
+  !$OMP THREADPRIVATE(iadvtr)
+  !----------------------------------------------------------
+  !     Rajouts pour PPM
+  !----------------------------------------------------------
+  INTEGER indice,n
+  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
+  REAL CFLmaxz,aaa,bbb ! CFL maximum
+  REAL psppm(iim,jjp1) ! pression  au sol
+  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
+  REAL qppm(iim*jjp1,llm,nqtot)
+  REAL fluxwppm(iim,jjp1,llm)
+  REAL apppm(llmp1), bpppm(llmp1)
+  LOGICAL dum,fill
+  DATA fill/.true./
+  DATA dum/.true./
+  REAL,SAVE :: finmasse(ip1jmp1,llm)
+  integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
+  type(Request) :: Request_vanleer
+  REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
+  REAL,SAVE :: teta_tmp(ip1jmp1,llm)
+  REAL,SAVE :: pk_tmp(ip1jmp1,llm)
+
+  ijb_u=ij_begin
+  ije_u=ij_end
+
+  ijb_v=ij_begin-iip1
+  ije_v=ij_end
+  if (pole_nord) ijb_v=ij_begin
+  if (pole_sud)  ije_v=ij_end-iip1
+
+  IF(iadvtr.EQ.0) THEN
+     !         CALL initial0(ijp1llm,pbaruc)
+     !         CALL initial0(ijmllm,pbarvc)
+     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)        
+     DO l=1,llm   
+        pbaruc(ijb_u:ije_u,l)=0.
+        pbarvc(ijb_v:ije_v,l)=0.
+     ENDDO
+     !$OMP END DO NOWAIT  
+  ENDIF
+
+  !   accumulation des flux de masse horizontaux
+  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+  DO l=1,llm
+     DO ij = ijb_u,ije_u
+        pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
+     ENDDO
+     DO ij = ijb_v,ije_v
+        pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
+     ENDDO
+  ENDDO
+  !$OMP END DO NOWAIT
+
+  !   selection de la masse instantannee des mailles avant le transport.
+  IF(iadvtr.EQ.0) THEN
+
+     !         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
+     ijb=ij_begin
+     ije=ij_end
+
+     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+     DO l=1,llm
+        massem(ijb:ije,l)=masse(ijb:ije,l)
+     ENDDO
+     !$OMP END DO NOWAIT
+
+     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
+     !
+  ENDIF ! of IF(iadvtr.EQ.0)
+
+  iadvtr   = iadvtr+1
+
+  !$OMP MASTER
+  iapptrac = iadvtr
+  !$OMP END MASTER
+
+  !   Test pour savoir si on advecte a ce pas de temps
+
+  IF ( iadvtr.EQ.iapp_tracvl ) THEN
+     !$OMP MASTER
+     call suspend_timer(timer_caldyn)
+     !$OMP END MASTER
+
+     ijb=ij_begin
+     ije=ij_end
+
+
+     !c   ..  Modif P.Le Van  ( 20/12/97 )  ....
+     !c
+
+     !   traitement des flux de masse avant advection.
+     !     1. calcul de w
+     !     2. groupement des mailles pres du pole.
+
+     CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
+
+     !$OMP BARRIER
+
+     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+     DO l=1,llmp1
+        p_tmp(ijb:ije,l)=p(ijb:ije,l)
+     ENDDO
+     !$OMP END DO NOWAIT
+
+     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+     DO l=1,llm
+        pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
+        teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
+     ENDDO
+     !$OMP END DO NOWAIT
+
+     !$OMP MASTER
+     call VTb(VTHallo)
+     !$OMP END MASTER
+
+     call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm, &
+          jj_Nb_vanleer,0,0,Request_vanleer)
+     call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm, &
+          jj_Nb_vanleer,1,0,Request_vanleer)
+     call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm, &
+          jj_Nb_vanleer,0,0,Request_vanleer)
+     call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm, &
+          jj_Nb_vanleer,0,0,Request_vanleer)
+     call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm, &
+          jj_Nb_vanleer,1,1,Request_vanleer)
+     call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1, &
+          jj_Nb_vanleer,1,1,Request_vanleer)
+     call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm, &
+          jj_Nb_vanleer,1,1,Request_vanleer)
+     do j=1,nqtot
+        call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
+             jj_nb_vanleer,0,0,Request_vanleer)
+     enddo
+
+     call SendRequest(Request_vanleer)
+     !$OMP BARRIER
+     call WaitRequest(Request_vanleer)
+
+
+     !$OMP BARRIER
+     !$OMP MASTER      
+     call SetDistrib(jj_nb_vanleer)
+     call VTe(VTHallo)
+     call VTb(VTadvection)
+     call start_timer(timer_vanleer)
+     !$OMP END MASTER
+     !$OMP BARRIER
+
+     ! ... Flux de masse diaganostiques traceurs
+     ijb=ij_begin
+     ije=ij_end
+     flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
+
+     !  test sur l'eventuelle creation de valeurs negatives de la masse
+     ijb=ij_begin
+     ije=ij_end
+     if (pole_nord) ijb=ij_begin+iip1
+     if (pole_sud) ije=ij_end-iip1
+
+     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)          
+     DO l=1,llm-1
+        DO ij = ijb+1,ije
+           zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l) &
+                - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
+                +       wg(ij,l+1)  - wg(ij,l)
+        ENDDO
+
+        !            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
+        ! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
+
+        do ij=ijb,ije-iip1+1,iip1
+           zdp(ij)=zdp(ij+iip1-1)
+        enddo
+
+        DO ij = ijb,ije
+           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 
+        ENDDO
+
+
+        !            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
+        !  ym ---> eventuellement a revoir
+        CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
+
+        IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
+           PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin, &
+                '   MAX:', zdpmax
+        ENDIF
+
+     ENDDO
+     !$OMP END DO NOWAIT
+
+     !-------------------------------------------------------------------
+     !   Advection proprement dite (Modification Le Croller (07/2001)
+     !-------------------------------------------------------------------
+
+     !----------------------------------------------------
+     !        Calcul des moyennes basées sur la masse
+     !----------------------------------------------------
+
+     !ym      ----> Normalement, inutile pour les schémas classiques
+     !ym      ----> Revérifier lors de la parallélisation des autres schemas
+
+     !ym          call massbar_p(massem,massebx,masseby)  
+
+     call vlspltgen_p( q,iadv, 2., massem, wg , &
+          pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
+
+
+     GOTO 1234     
+     !-----------------------------------------------------------
+     !     Appel des sous programmes d'advection
+     !-----------------------------------------------------------
+     do iq=1,nqtot
+        !        call clock(t_initial)
+        if(iadv(iq) == 0) cycle 
+        !   ----------------------------------------------------------------
+        !   Schema de Van Leer I MUSCL
+        !   ----------------------------------------------------------------
+        if(iadv(iq).eq.10) THEN
+
+           call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
+
+           !   ----------------------------------------------------------------
+           !   Schema "pseudo amont" + test sur humidite specifique
+           !    pour la vapeur d'eau. F. Codron
+           !   ----------------------------------------------------------------
+        else if(iadv(iq).eq.14) then
+           !
+           !ym        stop 'advtrac : appel à vlspltqs :schema non parallelise'
+           CALL vlspltqs_p( q(1,1,1), 2., massem, wg , &
+                pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
+           !   ----------------------------------------------------------------
+           !   Schema de Frederic Hourdin
+           !   ----------------------------------------------------------------
+        else if(iadv(iq).eq.12) then
+           stop 'advtrac : schema non parallelise'
+           !            Pas de temps adaptatif
+           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           do indice=1,n
+              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
+           end do
+        else if(iadv(iq).eq.13) then
+           stop 'advtrac : schema non parallelise'
+           !            Pas de temps adaptatif
+           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           do indice=1,n
+              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
+           end do
+           !   ----------------------------------------------------------------
+           !   Schema de pente SLOPES
+           !   ----------------------------------------------------------------
+        else if (iadv(iq).eq.20) then
+           stop 'advtrac : schema non parallelise'
+
+           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
+
+           !   ----------------------------------------------------------------
+           !   Schema de Prather
+           !   ----------------------------------------------------------------
+        else if (iadv(iq).eq.30) then
+           stop 'advtrac : schema non parallelise'
+           !            Pas de temps adaptatif
+           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
+                n,dtbon)
+           !   ----------------------------------------------------------------
+           !   Schemas PPM Lin et Rood
+           !   ----------------------------------------------------------------
+        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
+             iadv(iq).LE.18)) then
+
+           stop 'advtrac : schema non parallelise'
+
+           !        Test sur le flux horizontal
+           !        Pas de temps adaptatif
+           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
+           if (n.GT.1) then
+              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
+                   dtvr,'n=',n
+           endif
+           !        Test sur le flux vertical
+           CFLmaxz=0.
+           do l=2,llm
+              do ij=iip2,ip1jm
+                 aaa=wg(ij,l)*dtvr/massem(ij,l)
+                 CFLmaxz=max(CFLmaxz,aaa)
+                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
+                 CFLmaxz=max(CFLmaxz,bbb)
+              enddo
+           enddo
+           if (CFLmaxz.GE.1) then
+              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
+           endif
+
+           !-----------------------------------------------------------
+           !        Ss-prg interface LMDZ.4->PPM3d
+           !-----------------------------------------------------------
+
+           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
+                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
+                unatppm,vnatppm,psppm)
+
+           do indice=1,n
+              !----------------------------------------------------------------
+              !                         VL (version PPM) horiz. et PPM vert.
+              !----------------------------------------------------------------
+              if (iadv(iq).eq.11) then
+                 !                  Ss-prg PPM3d de Lin
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+
+                 !-------------------------------------------------------------
+                 !                           Monotonic PPM
+                 !-------------------------------------------------------------
+              else if (iadv(iq).eq.16) then
+                 !                  Ss-prg PPM3d de Lin
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+                 !-------------------------------------------------------------
+
+                 !-------------------------------------------------------------
+                 !                           Semi Monotonic PPM
+                 !-------------------------------------------------------------
+              else if (iadv(iq).eq.17) then
+                 !                  Ss-prg PPM3d de Lin
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+                 !-------------------------------------------------------------
+
+                 !-------------------------------------------------------------
+                 !                         Positive Definite PPM
+                 !-------------------------------------------------------------
+              else if (iadv(iq).eq.18) then
+                 !                  Ss-prg PPM3d de Lin
+                 call ppm3d(1,qppm(1,1,iq), &
+                      psppm,psppm, &
+                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
+                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
+                      fill,dum,220.)
+                 !-------------------------------------------------------------
+              endif
+           enddo
+           !-----------------------------------------------------------------
+           !               Ss-prg interface PPM3d-LMDZ.4
+           !-----------------------------------------------------------------
+           call interpost(q(1,1,iq),qppm(1,1,iq))
+        endif
+        !----------------------------------------------------------------------
+
+        !-----------------------------------------------------------------
+        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
+        ! et Nord j=1
+        !-----------------------------------------------------------------
+
+        !                  call traceurpole(q(1,1,iq),massem)
+
+        ! calcul du temps cpu pour un schema donne
+
+        !                  call clock(t_final)
+        !ym                  tps_cpu=t_final-t_initial
+        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
+
+     end DO
+
+1234 CONTINUE
+     !$OMP BARRIER
+
+     if (planet_type=="earth") then
+
+        ijb=ij_begin
+        ije=ij_end
+
+        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+        DO l = 1, llm
+           DO ij = ijb, ije
+              finmasse(ij,l) =  p(ij,l) - p(ij,l+1) 
+           ENDDO
+        ENDDO
+        !$OMP END DO
+
+        CALL qminimum_p( q, 2, finmasse )
+
+        !------------------------------------------------------------------
+        !   on reinitialise a zero les flux de masse cumules
+        !---------------------------------------------------
+        !          iadvtr=0
+
+        !$OMP MASTER
+	call VTe(VTadvection)
+        call stop_timer(timer_vanleer)
+        call VTb(VThallo)
+        !$OMP END MASTER
+
+	do j=1,nqtot
+           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
+                jj_nb_caldyn,0,0,Request_vanleer)
+        enddo
+
+        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, &
+             jj_nb_caldyn,0,0,Request_vanleer)
+
+        call SendRequest(Request_vanleer)
+        !$OMP BARRIER
+        call WaitRequest(Request_vanleer)      
+
+        !$OMP BARRIER
+        !$OMP MASTER
+        call SetDistrib(jj_nb_caldyn)
+	call VTe(VThallo)
+	call resume_timer(timer_caldyn)
+ !$OMP END MASTER
+ !$OMP BARRIER	
+        iadvtr=0
+     endif ! of if (planet_type=="earth")
+  ENDIF ! if iadvtr.EQ.iapp_tracvl
+
+END SUBROUTINE advtrac_p
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/bilan_dyn_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/bilan_dyn_p.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/bilan_dyn_p.F	(revision 270)
@@ -389,5 +389,4 @@
          Q_cum(:,jjb:jje,:,:)=0.
          flux_uQ_cum(:,jjb:jje,:,:)=0.
-         flux_v_cum(:,jjb:jje,:)=0.
          if (pole_sud) jje=jj_end-1
          flux_v_cum(:,jjb:jje,:)=0.
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90	(revision 270)
@@ -1,4 +1,4 @@
 !
-! $Id: ce0l.F90 1482 2011-02-09 15:03:09Z jghattas $
+! $Id: ce0l.F90 1511 2011-04-28 15:21:47Z jghattas $
 !
 !-------------------------------------------------------------------------------
@@ -19,8 +19,5 @@
   USE dimphy
   USE comgeomphy
-  USE mod_phys_lmdz_para
-  USE mod_const_mpi
   USE infotrac
-  USE parallel, ONLY: finalize_parallel
 
 #ifdef CPP_IOIPSL
@@ -30,5 +27,4 @@
 #endif
   IMPLICIT NONE
-#include "iniprint.h"
 #ifndef CPP_EARTH
   WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
@@ -40,21 +36,17 @@
 #include "paramet.h"
 #include "indicesol.h"
+#include "iniprint.h"
 #include "temps.h"
 #include "logic.h"
+  INTEGER, PARAMETER            :: longcles=20
+  REAL,    DIMENSION(longcles)  :: clesphy0
   REAL,    DIMENSION(iip1,jjp1) :: masque
   CHARACTER(LEN=15)             :: calnd
+  REAL,    DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol
 !-------------------------------------------------------------------------------
-  CALL conf_gcm( 99, .TRUE. )
-
-  CALL init_mpi
+  CALL conf_gcm( 99, .TRUE. , clesphy0 )
 
   CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
   WRITE(lunout,*)'---> klon=',klon
-  IF (mpi_size>1 .OR. omp_size>1) THEN
-       CALL abort_gcm('ce0l','In parallel mode,                         &
- &                 ce0l must be called only                             &
- &                 for 1 process and 1 task',1)
-  ENDIF
-
   CALL InitComgeomphy
 
@@ -89,5 +81,5 @@
   WRITE(lunout,'(//)')
   WRITE(lunout,*) ' interbar = ',interbar
-  CALL etat0_netcdf(interbar,masque,ok_etat0)
+  CALL etat0_netcdf(interbar,masque,phis,ok_etat0)
 
   IF(ok_limit) THEN
@@ -100,8 +92,12 @@
   END IF
 
-!$OMP MASTER
-  CALL finalize_parallel
-!$OMP END MASTER
-
+  IF (grilles_gcm_netcdf) THEN
+     WRITE(lunout,'(//)')
+     WRITE(lunout,*) '  ***************************  '
+     WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
+     WRITE(lunout,*) '  ***************************  '
+     WRITE(lunout,'(//)')
+     CALL grilles_gcm_netcdf_sub(masque,phis)
+  END IF
 #endif
 ! of #ifndef CPP_EARTH #else
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/comdissipn.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/comdissipn.h	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/comdissipn.h	(revision 270)
@@ -2,15 +2,19 @@
 ! $Header$
 !
-c-----------------------------------------------------------------------
-c INCLUDE comdissipn.h
+!  Attention : ce fichier include est compatible format fixe/format libre
+!                 veillez à n'utiliser que des ! pour les commentaires
+!                 et à bien positionner les & des lignes de continuation 
+!                 (les placer en colonne 6 et en colonne 73)
+!-----------------------------------------------------------------------
+! INCLUDE comdissipn.h
 
       REAL  tetaudiv, tetaurot, tetah, cdivu, crot, cdivh
-c
-      COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm)   ,
-     1                        cdivu,      crot,         cdivh
+!
+      COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm)   ,     &
+     &                        cdivu,      crot,         cdivh
 
-c
-c    Les parametres de ce common proviennent des calculs effectues dans 
-c             Inidissip  .
-c
-c-----------------------------------------------------------------------
+!
+!    Les parametres de ce common proviennent des calculs effectues dans 
+!             Inidissip  .
+!
+!-----------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 270)
@@ -1,4 +1,4 @@
 !
-! $Id: conf_gcm.F 1438 2010-10-08 10:19:34Z jghattas $
+! $Id: conf_gcm.F 1418 2010-07-19 15:11:24Z jghattas $
 !
 c
@@ -6,4 +6,5 @@
       SUBROUTINE conf_gcm( tapedef, etatinit )
 c
+      USE control_mod
 #ifdef CPP_IOIPSL
       use IOIPSL
@@ -12,9 +13,4 @@
       use ioipsl_getincom
 #endif
-      use misc_mod
-      use mod_filtre_fft, ONLY : use_filtre_fft
-      use mod_hallo, ONLY : use_mpi_alloc
-      use parallel, ONLY : omp_chunk
-      USE control_mod
       IMPLICIT NONE
 c-----------------------------------------------------------------------
@@ -37,10 +33,10 @@
 #include "serre.h"
 #include "comdissnew.h"
-!#include "clesphys.h"
-#include "iniprint.h"
 #include "temps.h"
 #include "comconst.h"
 
 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
+! #include "clesphys.h"
+#include "iniprint.h"
 c
 c
@@ -53,5 +49,5 @@
       LOGICAL  fxyhypbb, ysinuss
       INTEGER i
-      
+      LOGICAL use_filtre_fft
 c
 c  -------------------------------------------------------------------
@@ -82,11 +78,5 @@
 c   initialisations:
 c   ----------------
-      adjust=.false.
-      call getin('adjust',adjust)
-      
-      itaumax=0
-      call getin('itaumax',itaumax);
-      if (itaumax<=0) itaumax=HUGE(itaumax)
-      
+
 !Config  Key  = lunout
 !Config  Desc = unite de fichier pour les impressions
@@ -190,8 +180,7 @@
 
 !Config  Key  = nsplit_phys
-!Config  Desc = nombre d'iteration de la physique
-!Config  Def  = 240 
-!Config  Help = nombre d'itration de la physique
-!
+!Config  Desc = nombre de subdivisions par pas physique
+!Config  Def  = 1 
+!Config  Help = nombre de subdivisions par pas physique
        nsplit_phys = 1 
        CALL getin('nsplit_phys',nsplit_phys)
@@ -241,11 +230,12 @@
        CALL getin('output_grads_dyn',output_grads_dyn)
 
-!Config  Key  = idissip
+!Config  Key  = dissip_period
 !Config  Desc = periode de la dissipation 
-!Config  Def  = 10
+!Config  Def  = 0
 !Config  Help = periode de la dissipation 
-!Config         (en pas) ... a completer !
-       idissip = 10
-       CALL getin('idissip',idissip)
+!Config  dissip_period=0 => la valeur sera calcule dans inidissip       
+!Config  dissip_period>0 => on prend cette valeur
+       dissip_period = 0
+       CALL getin('dissip_period',dissip_period)
 
 ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
@@ -340,5 +330,4 @@
        CALL getin('tau_top_bound',tau_top_bound)
 
-!
 !Config  Key  = coefdis
 !Config  Desc = coefficient pour gamdissip
@@ -407,6 +396,4 @@
        ip_ebil_dyn = 0
        CALL getin('ip_ebil_dyn',ip_ebil_dyn)
-!
-
 
 ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
@@ -613,13 +600,5 @@
        offline = .FALSE.
        CALL getin('offline',offline)
-       IF (offline .AND. adjust) THEN
-          WRITE(lunout,*) 
-     &         'WARNING : option offline does not work with adjust=y :'
-          WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ', 
-     &         'and fluxstokev.nc will not be created'
-          WRITE(lunout,*) 
-     &         'only the file phystoke.nc will still be created ' 
-       END IF
-       
+
 !Config  Key  = config_inca
 !Config  Desc = Choix de configuration de INCA
@@ -655,4 +634,5 @@
       ok_dyn_ave = .FALSE. 
       CALL getin('ok_dyn_ave',ok_dyn_ave) 
+
 
       write(lunout,*)' #########################################'
@@ -669,10 +649,9 @@
       write(lunout,*)' day_step = ', day_step
       write(lunout,*)' iperiod = ', iperiod
-      write(lunout,*)' nsplit_phys = ', nsplit_phys
       write(lunout,*)' iconser = ', iconser
       write(lunout,*)' iecri = ', iecri
       write(lunout,*)' periodav = ', periodav 
       write(lunout,*)' output_grads_dyn = ', output_grads_dyn
-      write(lunout,*)' idissip = ', idissip
+      write(lunout,*)' dissip_period = ', dissip_period
       write(lunout,*)' lstardis = ', lstardis
       write(lunout,*)' nitergdiv = ', nitergdiv
@@ -816,12 +795,4 @@
        offline = .FALSE.
        CALL getin('offline',offline)
-       IF (offline .AND. adjust) THEN
-          WRITE(lunout,*) 
-     &         'WARNING : option offline does not work with adjust=y :'
-          WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ', 
-     &         'and fluxstokev.nc will not be created'
-          WRITE(lunout,*) 
-     &         'only the file phystoke.nc will still be created ' 
-       END IF
 
 !Config  Key  = config_inca
@@ -836,10 +807,10 @@
 
 !Config  Key  = ok_dynzon 
-!Config  Desc = calcul et sortie des transports 
+!Config  Desc = sortie des transports zonaux dans la dynamique
 !Config  Def  = n 
-!Config  Help = Permet de mettre en route le calcul des transports 
+!Config  Help =  
 !Config          
-      ok_dynzon = .FALSE. 
-      CALL getin('ok_dynzon',ok_dynzon) 
+       ok_dynzon = .FALSE.
+       CALL getin('ok_dynzon',ok_dynzon) 
 
 !Config  Key  = ok_dyn_ins
@@ -863,38 +834,15 @@
 !Config  Def  = false
 !Config  Help = permet d'activer l'utilisation des FFT pour effectuer
-!Config         le filtrage aux poles. 
+!Config         le filtrage aux poles.
+! Le filtre fft n'est pas implemente dans dyn3d 
       use_filtre_fft=.FALSE.
       CALL getin('use_filtre_fft',use_filtre_fft)
 
-      IF (use_filtre_fft .AND. grossismx /= 1.0) THEN
-        write(lunout,*)'WARNING !!! '
-        write(lunout,*)"Le zoom en longitude est incompatible",
-     &                 " avec l'utilisation du filtre FFT ",
-     &                 "---> filtre FFT dÃ©sactivÃ© "
-       use_filtre_fft=.FALSE.
+      IF (use_filtre_fft) THEN
+        write(lunout,*)'STOP !!!'
+        write(lunout,*)'use_filtre_fft n est pas implemente dans dyn3d'
+        STOP
       ENDIF
       
- 
-      
-!Config  Key  = use_mpi_alloc
-!Config  Desc = Utilise un buffer MPI en mï¿½moire globale
-!Config  Def  = false
-!Config  Help = permet d'activer l'utilisation d'un buffer MPI
-!Config         en mï¿½moire globale a l'aide de la fonction MPI_ALLOC.
-!Config         Cela peut amï¿½liorer la bande passante des transferts MPI
-!Config         d'un facteur 2  
-      use_mpi_alloc=.FALSE.
-      CALL getin('use_mpi_alloc',use_mpi_alloc)
-
-!Config  Key  = omp_chunk
-!Config  Desc = taille des blocs openmp
-!Config  Def  = 1
-!Config  Help = defini la taille des packets d'itï¿½ration openmp
-!Config         distribuï¿½e ï¿½ chaque tï¿½che lors de l'entrï¿½e dans une
-!Config         boucle parallï¿½lisï¿½e
-  
-      omp_chunk=1
-      CALL getin('omp_chunk',omp_chunk)
-
 !Config key = ok_strato
 !Config  Desc = activation de la version strato
@@ -902,5 +850,5 @@
 !Config  Help = active la version stratosphérique de LMDZ de F. Lott
 
-      ok_strato=.FALSE.
+      ok_strato=.TRUE.
       CALL getin('ok_strato',ok_strato)
 
@@ -928,4 +876,10 @@
       ok_etat0 = .TRUE.
       CALL getin('ok_etat0',ok_etat0)
+
+!Config  Key  = grilles_gcm_netcdf
+!Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
+!Config  Def  = n
+      grilles_gcm_netcdf = .FALSE.
+      CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
 
       write(lunout,*)' #########################################'
@@ -943,5 +897,5 @@
       write(lunout,*)' periodav = ', periodav 
       write(lunout,*)' output_grads_dyn = ', output_grads_dyn
-      write(lunout,*)' idissip = ', idissip
+      write(lunout,*)' dissip_period = ', dissip_period
       write(lunout,*)' lstardis = ', lstardis
       write(lunout,*)' nitergdiv = ', nitergdiv
@@ -968,14 +922,12 @@
       write(lunout,*)' offline = ', offline
       write(lunout,*)' config_inca = ', config_inca
-      write(lunout,*)' ok_dynzon = ', ok_dynzon 
+      write(lunout,*)' ok_dynzon = ', ok_dynzon
       write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
       write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave 
-      write(lunout,*)' use_filtre_fft = ', use_filtre_fft
-      write(lunout,*)' use_mpi_alloc = ', use_mpi_alloc
-      write(lunout,*)' omp_chunk = ', omp_chunk
       write(lunout,*)' ok_strato = ', ok_strato
       write(lunout,*)' ok_gradsfile = ', ok_gradsfile
       write(lunout,*)' ok_limit = ', ok_limit
       write(lunout,*)' ok_etat0 = ', ok_etat0
+      write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
 c
       RETURN
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/control_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/control_mod.F90	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/control_mod.F90	(revision 270)
@@ -12,5 +12,5 @@
   REAL    :: periodav
   INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys
-  INTEGER :: iconser,iecri,idissip,iphysiq,iecrimoy
+  INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy
   INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn
   LOGICAL :: offline
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/defrun.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/defrun.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/defrun.F	(revision 270)
@@ -6,5 +6,12 @@
       SUBROUTINE defrun( tapedef, etatinit, clesphy0 )
 c
+! ========================== ATTENTION =============================
+! COMMENTAIRE SL : 
+! NE SERT PLUS APPAREMMENT
+! DONC PAS MIS A JOUR POUR L'UTILISATION AVEC LES PLANETES
+! ==================================================================
+
       USE control_mod
+ 
       IMPLICIT NONE
 c-----------------------------------------------------------------------
@@ -132,7 +139,7 @@
 
       READ (tapedef,9001) ch1,ch4
-      READ (tapedef,*)    idissip
-      WRITE(tapeout,9001) ch1,'idissip'
-      WRITE(tapeout,*)    idissip
+      READ (tapedef,*)    dissip_period
+      WRITE(tapeout,9001) ch1,'dissip_period'
+      WRITE(tapeout,*)    dissip_period
 
 ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
@@ -196,6 +203,4 @@
 
 
-ccc   .... P.Le Van, ajout le 03/01/96 pour l'ecriture phys ...
-c
       READ (tapedef,9001) ch1,ch4
       READ (tapedef,*)    cycle_diurne
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F	(revision 270)
@@ -125,5 +125,5 @@
       endif
 !$OMP END MASTER
-
+!$OMP BARRIER
         jjb=jj_begin
         jje=jj_end
@@ -171,4 +171,5 @@
       endif
 c$OMP END MASTER
+c$OMP BARRIER
 c
 c
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F	(revision 270)
@@ -121,5 +121,5 @@
       endif
 !$OMP END MASTER
-
+!$OMP BARRIER
         jjb=jj_begin
         jje=jj_end
@@ -169,4 +169,5 @@
       endif
 c$OMP END MASTER
+c$OMP BARRIER
 c
 c
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/grilles_gcm_netcdf_sub.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/grilles_gcm_netcdf_sub.F90	(revision 270)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/grilles_gcm_netcdf_sub.F90	(revision 270)
@@ -0,0 +1,237 @@
+!
+! $Header$
+!
+! This subroutine creates the file grilles_gcm.nc containg longitudes and
+! latitudes in degrees for grid u and v. This subroutine is called from
+! ce0l if grilles_gcm_netcdf=TRUE. This subroutine corresponds to the first 
+! part in the program create_fausse_var.
+!
+SUBROUTINE grilles_gcm_netcdf_sub(masque,phis)
+
+  IMPLICIT NONE
+
+  INCLUDE "dimensions.h"
+  INCLUDE "paramet.h"
+  INCLUDE "comconst.h"
+  INCLUDE "comgeom.h"
+  INCLUDE "comvert.h"
+  INCLUDE "netcdf.inc"
+  INCLUDE "serre.h"
+
+
+  REAL,DIMENSION(iip1,jjp1),INTENT(IN)  :: masque ! masque terre/mer
+  REAL,DIMENSION(iip1,jjp1),INTENT(IN)  :: phis   ! geopotentiel au sol
+
+  REAL temp(iim+1,jjm+1)
+  ! Attributs netcdf sortie
+  INTEGER ncid_out,rcode_out
+  INTEGER out_lonuid,out_lonvid,out_latuid,out_latvid,out_levid
+  INTEGER out_varid
+  INTEGER out_lonudim,out_lonvdim
+  INTEGER out_latudim,out_latvdim,out_dim(3)
+  INTEGER out_levdim
+
+  INTEGER, PARAMETER :: longcles = 20
+  REAL  clesphy0(longcles)
+
+  INTEGER start(4),COUNT(4)
+
+  INTEGER status,i,j
+  REAL rlatudeg(jjp1),rlatvdeg(jjm),rlevdeg(llm)
+  REAL rlonudeg(iip1),rlonvdeg(iip1)
+
+  REAL dlon1(iip1),dlon2(iip1),dlat1(jjp1),dlat2(jjp1)
+  REAL acoslat,dxkm,dykm,resol(iip1,jjp1)
+  REAL,DIMENSION(iip1,jjp1)  :: phis_loc
+  INTEGER masque_int(iip1,jjp1)
+  INTEGER :: phis_id
+  INTEGER :: area_id
+  INTEGER :: mask_id
+  
+  rad = 6400000
+  omeg = 7.272205e-05
+  g = 9.8
+  kappa = 0.285716
+  daysec = 86400
+  cpp = 1004.70885
+
+  preff = 101325.
+  pa= 50000.
+
+  CALL conf_gcm( 99, .TRUE. , clesphy0 )
+  CALL iniconst
+  CALL inigeom
+
+  DO j=1,jjp1
+     rlatudeg(j)=rlatu(j)*180./pi
+  ENDDO
+  DO j=1,jjm
+     rlatvdeg(j)=rlatv(j)*180./pi
+  ENDDO
+
+  DO i=1,iip1
+     rlonudeg(i)=rlonu(i)*180./pi + 360.
+     rlonvdeg(i)=rlonv(i)*180./pi + 360.
+  ENDDO
+
+
+  !  2 ----- OUVERTURE DE LA SORTIE NETCDF
+  ! ---------------------------------------------------
+  ! CREATION OUTPUT
+  ! ouverture fichier netcdf de sortie out
+  status=NF_CREATE('grilles_gcm.nc',NF_NOCLOBBER,ncid_out)
+  status=NF_DEF_DIM(ncid_out,'lonu',iim+1,out_lonudim)
+  status=NF_DEF_DIM(ncid_out,'lonv',iim+1,out_lonvdim)
+  status=NF_DEF_DIM(ncid_out,'latu',jjm+1,out_latudim)
+  status=NF_DEF_DIM(ncid_out,'latv',jjm,out_latvdim)
+
+
+  !   Longitudes en u
+  status=NF_DEF_VAR(ncid_out,'lonu',NF_FLOAT,1,out_lonudim, out_lonuid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_lonuid,'units', 12,'degrees_east')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_lonuid,'long_name',9,'Longitude en u')
+
+  !   Longitudes en v
+  status=NF_DEF_VAR(ncid_out,'lonv',NF_FLOAT,1,out_lonvdim, out_lonvid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_lonvid,'units', 12,'degrees_east')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_lonvid,'long_name', 9,'Longitude en v')
+
+  !   Latitude en u
+  status=NF_DEF_VAR(ncid_out,'latu',NF_FLOAT,1,out_latudim, out_latuid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_latuid,'units', 13,'degrees_north')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_latuid,'long_name', 8,'Latitude en u')
+
+  !  Latitude en v
+  status=NF_DEF_VAR(ncid_out,'latv',NF_FLOAT,1,out_latvdim, out_latvid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_latvid,'units', 13,'degrees_north')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_latvid,'long_name', 8,'Latitude en v')
+
+  !   ecriture de la grille u
+  out_dim(1)=out_lonudim
+  out_dim(2)=out_latudim
+  status=NF_DEF_VAR(ncid_out,'grille_u',NF_FLOAT,2,out_dim, out_varid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units', 6,'Kelvin')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name', 16,'Grille aux point u')
+
+  !   ecriture de la grille v
+  out_dim(1)=out_lonvdim
+  out_dim(2)=out_latvdim
+  status=NF_DEF_VAR(ncid_out,'grille_v',NF_FLOAT,2,out_dim, out_varid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units', 6,'Kelvin')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name', 16,'Grille aux point v')
+
+  !   ecriture de la grille u
+  out_dim(1)=out_lonvdim
+  out_dim(2)=out_latudim
+  status=NF_DEF_VAR(ncid_out,'grille_s',NF_FLOAT,2,out_dim, out_varid)
+  CALL handle_err(status)
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units', 6,'Kelvin')
+  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name',16,'Grille aux point u')
+
+  status=NF_ENDDEF(ncid_out)
+  ! 5) ----- FERMETURE DES FICHIERS NETCDF------------------
+  ! --------------------------------------------------------
+  ! 3-b- Ecriture de la grille pour la sortie
+  ! rajoute l'ecriture de la grille
+
+#ifdef NC_DOUBLE
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_lonuid,1,iim+1,rlonudeg)
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_lonvid,1,iim+1,rlonvdeg)
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_latuid,1,jjm+1,rlatudeg)
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_latvid,1,jjm,rlatvdeg)
+#else
+  status=NF_PUT_VARA_REAL(ncid_out,out_lonuid,1,iim+1,rlonudeg)
+  status=NF_PUT_VARA_REAL(ncid_out,out_lonvid,1,iim+1,rlonvdeg)
+  status=NF_PUT_VARA_REAL(ncid_out,out_latuid,1,jjm+1,rlatudeg)
+  status=NF_PUT_VARA_REAL(ncid_out,out_latvid,1,jjm,rlatvdeg)
+#endif
+
+  start(1)=1
+  start(2)=1
+  start(3)=1
+  start(4)=1
+
+  COUNT(1)=iim+1
+  COUNT(2)=jjm+1
+  COUNT(3)=1
+  COUNT(4)=1
+
+  DO j=1,jjm+1
+     DO i=1,iim+1
+        temp(i,j)=MOD(i,2)+MOD(j,2)
+     ENDDO
+  ENDDO
+
+#ifdef NC_DOUBLE
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_varid,start, count,temp)
+#else
+  status=NF_PUT_VARA_REAL(ncid_out,out_varid,start, count,temp)
+#endif
+
+  ! On re-ouvre le fichier pour rajouter 4 nouvelles variables necessaire pour INCA
+! lev - phis - aire - mask
+  rlevdeg(:) = presnivs
+  phis_loc(:,:) = phis(:,:)/g
+
+! niveaux de pression verticaux
+  status = NF_REDEF (ncid_out)
+  status=NF_DEF_DIM(ncid_out,'lev',llm,out_levdim)
+  
+! fields
+  out_dim(1)=out_lonvdim
+  out_dim(2)=out_latudim
+
+  status = nf_def_var(ncid_out,'phis',NF_FLOAT,2,out_dim,phis_id)
+  CALL handle_err(status)
+  status = nf_def_var(ncid_out,'aire',NF_FLOAT,2,out_dim,area_id)
+  CALL handle_err(status)
+  status = nf_def_var(ncid_out,'mask',NF_INT  ,2,out_dim,mask_id)
+  CALL handle_err(status)
+
+  status=NF_ENDDEF(ncid_out)
+
+  ! ecriture des variables
+#ifdef NC_DOUBLE
+  status=NF_PUT_VARA_DOUBLE(ncid_out,out_levid,1,llm,rlevdeg)
+#else
+  status=NF_PUT_VARA_REAL(ncid_out,out_levid,1,llm,rlevdeg)
+#endif
+
+  start(1)=1
+  start(2)=1
+  start(3)=1
+  start(4)=0
+  COUNT(1)=iip1
+  COUNT(2)=jjp1
+  COUNT(3)=1
+  COUNT(4)=0
+
+  status = nf_put_vara_double(ncid_out, phis_id,start,count, phis_loc)
+  status = nf_put_vara_double(ncid_out, area_id,start,count, aire)
+  masque_int(:,:) = nINT(masque(:,:))
+  status = nf_put_vara_int(ncid_out, mask_id,start,count,masque_int)
+  CALL handle_err(status)
+  
+  ! fermeture du fichier netcdf
+  CALL ncclos(ncid_out,rcode_out)
+
+END SUBROUTINE grilles_gcm_netcdf_sub
+
+
+
+SUBROUTINE handle_err(status)
+  INCLUDE "netcdf.inc"
+
+  INTEGER status
+  IF (status.NE.nf_noerr) THEN
+     PRINT *,NF_STRERROR(status)
+     CALL abort_gcm('grilles_gcm_netcdf','netcdf error',1)
+  ENDIF
+END SUBROUTINE handle_err
+
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/iniacademic.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/iniacademic.F90	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/iniacademic.F90	(revision 270)
@@ -1,4 +1,4 @@
 !
-! $Id: iniacademic.F90 1520 2011-05-23 11:37:09Z emillour $
+! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $
 !
 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
@@ -115,5 +115,5 @@
   endif
 
-  academic_case: if (iflag_phys == 2) then
+  academic_case: if (iflag_phys >= 2) then
      ! initializations
 
@@ -208,5 +208,9 @@
      IF (.NOT. read_start) THEN
         ! surface pressure
-        ps(:)=preff
+        if (iflag_phys>2) then
+           ps(:)=108080.  ! Earth aqua/terra planets
+        else
+           ps(:)=preff
+        endif
         ! ground geopotential
         phis(:)=0.
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/iniconst.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/iniconst.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/iniconst.F	(revision 270)
@@ -47,5 +47,4 @@
 c-----------------------------------------------------------------------
 
-      dtdiss  = idissip * dtvr
       dtphys  = iphysiq * dtvr
       unsim   = 1./iim
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F	(revision 134)
+++ 	(revision )
@@ -1,244 +1,0 @@
-!
-! $Id: inidissip.F 1403 2010-07-01 09:02:53Z fairhead $
-!
-      SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,
-     *                       tetagdiv,tetagrot,tetatemp             )
-c=======================================================================
-c   initialisation de la dissipation horizontale
-c=======================================================================
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-      USE control_mod
-
-      IMPLICIT NONE
-#include "dimensions.h"
-#include "paramet.h"
-#include "comdissipn.h"
-#include "comconst.h"
-#include "comvert.h"
-#include "logic.h"
-
-      LOGICAL lstardis
-      INTEGER nitergdiv,nitergrot,niterh
-      REAL    tetagdiv,tetagrot,tetatemp
-      REAL fact,zvert(llm),zz
-      REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
-      REAL ullm,vllm,umin,vmin,zhmin,zhmax
-      REAL zllm,z1llm
-
-      INTEGER l,ij,idum,ii
-      REAL tetamin
-      REAL Pup
-
-      REAL ran1
-
-
-c-----------------------------------------------------------------------
-c
-c   calcul des valeurs propres des operateurs par methode iterrative:
-c   -----------------------------------------------------------------
-
-      crot     = -1.
-      cdivu    = -1.
-      cdivh    = -1.
-
-c   calcul de la valeur propre de divgrad:
-c   --------------------------------------
-      idum = 0
-      DO l = 1, llm
-       DO ij = 1, ip1jmp1
-        deltap(ij,l) = 1.
-       ENDDO
-      ENDDO
-
-      idum  = -1
-      zh(1) = RAN1(idum)-.5
-      idum  = 0
-      DO ij = 2, ip1jmp1
-        zh(ij) = RAN1(idum) -.5
-      ENDDO
-
-      CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
-
-      CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
-
-      IF ( zhmin .GE. zhmax  )     THEN
-         PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
-         STOP'probleme generateur alleatoire dans inidissip'
-      ENDIF
-
-      zllm = ABS( zhmax )
-      DO l = 1,50
-         IF(lstardis) THEN
-            CALL divgrad2(1,zh,deltap,niterh,zh)
-         ELSE
-            CALL divgrad (1,zh,niterh,zh)
-         ENDIF
-
-        CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
-
-         zllm  = ABS( zhmax )
-         z1llm = 1./zllm
-         DO ij = 1,ip1jmp1
-            zh(ij) = zh(ij)* z1llm
-         ENDDO
-      ENDDO
-
-      IF(lstardis) THEN
-         cdivh = 1./ zllm
-      ELSE
-         cdivh = zllm ** ( -1./niterh )
-      ENDIF
-
-c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
-c   -----------------------------------------------------------------
-      print*,'calcul des valeurs propres'
-
-      DO  20  ii = 1, 2
-c
-         DO ij = 1, ip1jmp1
-           zu(ij)  = RAN1(idum) -.5
-         ENDDO
-         CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
-         DO ij = 1, ip1jm
-            zv(ij) = RAN1(idum) -.5
-         ENDDO
-         CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
-
-         CALL minmax(iip1*jjp1,zu,umin,ullm )
-         CALL minmax(iip1*jjm, zv,vmin,vllm )
-
-         ullm = ABS ( ullm )
-         vllm = ABS ( vllm )
-
-         DO  5  l = 1, 50
-            IF(ii.EQ.1) THEN
-ccccc             CALL covcont( 1,zu,zv,zu,zv )
-               IF(lstardis) THEN
-                  CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
-               ELSE
-                  CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
-               ENDIF
-            ELSE
-               IF(lstardis) THEN
-                  CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
-               ELSE
-                  CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
-               ENDIF
-            ENDIF
-
-            CALL minmax(iip1*jjp1,zu,umin,ullm )
-            CALL minmax(iip1*jjm, zv,vmin,vllm )
-
-            ullm = ABS  ( ullm )
-            vllm = ABS  ( vllm )
-
-            zllm  = MAX( ullm,vllm )
-            z1llm = 1./ zllm
-            DO ij = 1, ip1jmp1
-              zu(ij) = zu(ij)* z1llm
-            ENDDO
-            DO ij = 1, ip1jm
-               zv(ij) = zv(ij)* z1llm
-            ENDDO
- 5       CONTINUE
-
-         IF ( ii.EQ.1 ) THEN
-            IF(lstardis) THEN
-               cdivu  = 1./zllm
-            ELSE
-               cdivu  = zllm **( -1./nitergdiv )
-            ENDIF
-         ELSE
-            IF(lstardis) THEN
-               crot   = 1./ zllm
-            ELSE
-               crot   = zllm **( -1./nitergrot )
-            ENDIF
-         ENDIF
-
- 20   CONTINUE
-
-c   petit test pour les operateurs non star:
-c   ----------------------------------------
-
-c     IF(.NOT.lstardis) THEN
-         fact    = rad*24./REAL(jjm)
-         fact    = fact*fact
-         PRINT*,'coef u ', fact/cdivu, 1./cdivu
-         PRINT*,'coef r ', fact/crot , 1./crot
-         PRINT*,'coef h ', fact/cdivh, 1./cdivh
-c     ENDIF
-
-c-----------------------------------------------------------------------
-c   variation verticale du coefficient de dissipation:
-c   --------------------------------------------------
-
-c First step: going from 1 to dissip_fac_mid (in gcm.def)
-c============
-      DO l=1,llm
-         zz      = 1. - preff/presnivs(l)
-         zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
-      ENDDO
-
-         write(*,*) 'Dissipation : '
-         write(*,*) 'Multiplication de la dissipation en altitude :'
-         write(*,*) '  dissip_fac_mid =', dissip_fac_mid
-
-c Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
-c==========================
-c Utilisation de la fonction tangente hyperbolique pour augmenter
-c arbitrairement la dissipation et donc la stabilite du modele a 
-c partir d'une certaine altitude.
-
-c   Le facteur multiplicatif de basse atmosphere etant deja pris 
-c   en compte, il faut diviser le facteur multiplicatif de haute 
-c   atmosphere par celui-ci.
-
-      if (ok_strato) then
-
-       Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
-       do l=1,llm
-         zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1)
-     &   *(1-(0.5*(1+tanh(-6./dissip_deltaz*
-     &    (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
-       enddo 
-
-         write(*,*) '  dissip_fac_up =', dissip_fac_up
-         write(*,*) 'Transition mid /up:  Pupstart,delta =',
-     &             dissip_pupstart,'Pa', dissip_deltaz , '(km)'
-
-      endif
-
-
-      PRINT*,'Constantes de temps de la diffusion horizontale'
-
-      tetamin =  1.e+6
-
-      DO l=1,llm
-        tetaudiv(l)   = zvert(l)/tetagdiv
-        tetaurot(l)   = zvert(l)/tetagrot
-        tetah(l)      = zvert(l)/tetatemp
-
-        IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
-        IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
-        IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
-      ENDDO
-
-      PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
-      idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
-      PRINT *,' INIDI tetamin idissip ',tetamin,idissip
-      idissip = MAX(iperiod,idissip)
-      dtdiss  = idissip * dtvr
-      PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
-
-      DO l = 1,llm
-         PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),
-     *                   dtdiss*tetah(l)
-      ENDDO
-
-c
-      RETURN
-      END
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90	(revision 270)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90	(revision 270)
@@ -0,0 +1,250 @@
+!
+! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $
+!
+SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
+     tetagdiv,tetagrot,tetatemp             )
+  !=======================================================================
+  !   initialisation de la dissipation horizontale
+  !=======================================================================
+  !-----------------------------------------------------------------------
+  !   declarations:
+  !   -------------
+
+  USE control_mod, only : dissip_period,iperiod
+
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comdissipn.h"
+  include "comconst.h"
+  include "comvert.h"
+  include "logic.h"
+  include "iniprint.h"
+
+  LOGICAL,INTENT(in) :: lstardis
+  INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
+  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
+
+! Local variables:
+  REAL fact,zvert(llm),zz
+  REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
+  REAL ullm,vllm,umin,vmin,zhmin,zhmax
+  REAL zllm,z1llm
+
+  INTEGER l,ij,idum,ii
+  REAL tetamin
+  REAL Pup
+  character (len=80) :: abort_message
+
+  REAL ran1
+
+
+  !-----------------------------------------------------------------------
+  !
+  !   calcul des valeurs propres des operateurs par methode iterrative:
+  !   -----------------------------------------------------------------
+
+  crot     = -1.
+  cdivu    = -1.
+  cdivh    = -1.
+
+  !   calcul de la valeur propre de divgrad:
+  !   --------------------------------------
+  idum = 0
+  DO l = 1, llm
+     DO ij = 1, ip1jmp1
+        deltap(ij,l) = 1.
+     ENDDO
+  ENDDO
+
+  idum  = -1
+  zh(1) = RAN1(idum)-.5
+  idum  = 0
+  DO ij = 2, ip1jmp1
+     zh(ij) = RAN1(idum) -.5
+  ENDDO
+
+  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
+
+  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
+
+  IF ( zhmin .GE. zhmax  )     THEN
+     write(lunout,*)'  Inidissip  zh min max  ',zhmin,zhmax
+     abort_message='probleme generateur alleatoire dans inidissip'
+     call abort_gcm('inidissip',abort_message,1)
+  ENDIF
+
+  zllm = ABS( zhmax )
+  DO l = 1,50
+     IF(lstardis) THEN
+        CALL divgrad2(1,zh,deltap,niterh,zh)
+     ELSE
+        CALL divgrad (1,zh,niterh,zh)
+     ENDIF
+
+     CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
+
+     zllm  = ABS( zhmax )
+     z1llm = 1./zllm
+     DO ij = 1,ip1jmp1
+        zh(ij) = zh(ij)* z1llm
+     ENDDO
+  ENDDO
+
+  IF(lstardis) THEN
+     cdivh = 1./ zllm
+  ELSE
+     cdivh = zllm ** ( -1./niterh )
+  ENDIF
+
+  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
+  !   -----------------------------------------------------------------
+  write(lunout,*)'inidissip: calcul des valeurs propres'
+
+  DO    ii = 1, 2
+     !
+     DO ij = 1, ip1jmp1
+        zu(ij)  = RAN1(idum) -.5
+     ENDDO
+     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
+     DO ij = 1, ip1jm
+        zv(ij) = RAN1(idum) -.5
+     ENDDO
+     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
+
+     CALL minmax(iip1*jjp1,zu,umin,ullm )
+     CALL minmax(iip1*jjm, zv,vmin,vllm )
+
+     ullm = ABS ( ullm )
+     vllm = ABS ( vllm )
+
+     DO    l = 1, 50
+        IF(ii.EQ.1) THEN
+           !cccc             CALL covcont( 1,zu,zv,zu,zv )
+           IF(lstardis) THEN
+              CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
+           ELSE
+              CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
+           ENDIF
+        ELSE
+           IF(lstardis) THEN
+              CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
+           ELSE
+              CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
+           ENDIF
+        ENDIF
+
+        CALL minmax(iip1*jjp1,zu,umin,ullm )
+        CALL minmax(iip1*jjm, zv,vmin,vllm )
+
+        ullm = ABS  ( ullm )
+        vllm = ABS  ( vllm )
+
+        zllm  = MAX( ullm,vllm )
+        z1llm = 1./ zllm
+        DO ij = 1, ip1jmp1
+           zu(ij) = zu(ij)* z1llm
+        ENDDO
+        DO ij = 1, ip1jm
+           zv(ij) = zv(ij)* z1llm
+        ENDDO
+     end DO
+
+     IF ( ii.EQ.1 ) THEN
+        IF(lstardis) THEN
+           cdivu  = 1./zllm
+        ELSE
+           cdivu  = zllm **( -1./nitergdiv )
+        ENDIF
+     ELSE
+        IF(lstardis) THEN
+           crot   = 1./ zllm
+        ELSE
+           crot   = zllm **( -1./nitergrot )
+        ENDIF
+     ENDIF
+
+  end DO
+
+  !   petit test pour les operateurs non star:
+  !   ----------------------------------------
+
+  !     IF(.NOT.lstardis) THEN
+  fact    = rad*24./REAL(jjm)
+  fact    = fact*fact
+  write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
+  write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
+  write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
+  !     ENDIF
+
+  !-----------------------------------------------------------------------
+  !   variation verticale du coefficient de dissipation:
+  !   --------------------------------------------------
+
+! First step: going from 1 to dissip_fac_mid (in gcm.def)
+!============
+  DO l=1,llm
+     zz      = 1. - preff/presnivs(l)
+     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
+  ENDDO
+
+  write(lunout,*) 'Dissipation : '
+  write(lunout,*) 'Multiplication de la dissipation en altitude :'
+  write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
+
+! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
+!==========================
+! Utilisation de la fonction tangente hyperbolique pour augmenter
+! arbitrairement la dissipation et donc la stabilite du modele a 
+! partir d'une certaine altitude.
+
+!   Le facteur multiplicatif de basse atmosphere etant deja pris 
+!   en compte, il faut diviser le facteur multiplicatif de haute 
+!   atmosphere par celui-ci.
+
+  if (ok_strato) then
+
+    Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
+    do l=1,llm
+      zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
+                *(1-(0.5*(1+tanh(-6./dissip_deltaz*              &
+               (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
+    enddo 
+
+    write(*,*) '  dissip_fac_up =', dissip_fac_up
+    write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
+                   dissip_pupstart,'Pa', dissip_deltaz , '(km)'
+
+  endif
+
+
+  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
+
+  tetamin =  1.e+6
+
+  DO l=1,llm
+     tetaudiv(l)   = zvert(l)/tetagdiv
+     tetaurot(l)   = zvert(l)/tetagrot
+     tetah(l)      = zvert(l)/tetatemp
+
+     IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
+     IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
+     IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
+  ENDDO
+
+  ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
+  IF (dissip_period == 0) THEN
+     dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
+     write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period
+     dissip_period = MAX(iperiod,dissip_period)
+  END IF
+
+  dtdiss  = dissip_period * dtvr
+  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
+
+  DO l = 1,llm
+     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
+          dtdiss*tetah(l)
+  ENDDO
+
+END SUBROUTINE inidissip
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F	(revision 270)
@@ -1,4 +1,4 @@
 !
-! $Id: integrd_p.F 1446 2010-10-22 09:27:25Z emillour $
+! $Id: integrd_p.F 1550 2011-07-05 09:44:55Z lguez $
 !
       SUBROUTINE integrd_p
@@ -128,5 +128,6 @@
          PRINT*,' Au point ij = ',stop_it, ' , pression sol neg. '
      &         , ps(stop_it)
-         STOP' dans integrd'
+         print *, ' dans integrd'
+         stop 1
         ENDIF
 
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 270)
@@ -168,5 +168,5 @@
 
       character*80 dynhist_file, dynhistave_file
-      character*20 modname
+      character(len=*),parameter :: modname="leapfrog"
       character*80 abort_message
 
@@ -226,5 +226,4 @@
       endif
       itaufinp1 = itaufin +1
-      modname="leapfrog_p"
 
       itau = 0
@@ -392,5 +391,6 @@
       ! Purely Matsuno time stepping
          IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
-         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
+         IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward ) 
+     s        apdiss = .TRUE.
          IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 
      s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
@@ -398,5 +398,6 @@
       ! Leapfrog/Matsuno time stepping 
          IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
-         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
+         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
+     s        apdiss = .TRUE.
          IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
       END IF
@@ -650,5 +651,7 @@
 c   -------------------------------------------------------------
 
-      IF( forward. OR . leapf )  THEN
+!      IF( forward. OR . leapf )  THEN
+      IF((.not.forward).OR. leapf )  THEN
+        ! Ehouarn: gather mass fluxes during backward Matsuno or LF step
 cc$OMP PARALLEL DEFAULT(SHARED) 
 c
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/limit_netcdf.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/limit_netcdf.F90	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/limit_netcdf.F90	(revision 270)
@@ -1,4 +1,4 @@
 !
-! $Id: limit_netcdf.F90 1441 2010-10-13 13:06:56Z emillour $
+! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $
 !-------------------------------------------------------------------------------
 !
@@ -42,5 +42,5 @@
   REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
 #ifndef CPP_EARTH
-  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
+  CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1)
 #else
 !-------------------------------------------------------------------------------
@@ -52,18 +52,15 @@
 #include "indicesol.h"
 
-!--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------
-  LOGICAL, PARAMETER :: fracterre=.TRUE.
-
 !--- INPUT NETCDF FILES NAMES --------------------------------------------------
   CHARACTER(LEN=25) :: icefile, sstfile, dumstr
   CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc        ',        &
                                   famipsic='amipbc_sic_1x1.nc        ',        &
-                                  fclimsst='amipbc_sst_1x1_clim.nc   ',        &
-                                  fclimsic='amipbc_sic_1x1_clim.nc   ',        &
                                   fcpldsst='cpl_atm_sst.nc           ',        &
                                   fcpldsic='cpl_atm_sic.nc           ',        &
+                                  fhistsst='histmth_sst.nc           ',        &
+                                  fhistsic='histmth_sic.nc           ',        &
                                   frugo   ='Rugos.nc                 ',        &
                                   falbe   ='Albedo.nc                '
-
+  CHARACTER(LEN=10) :: varname
 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
   REAL,   DIMENSION(klon)                :: fi_ice, verif
@@ -80,5 +77,4 @@
   INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC
   INTEGER :: NF90_FORMAT
-  LOGICAL :: lCPL                    !--- T: IPCC-IPSL cpl model output files
   INTEGER :: ndays                   !--- Depending on the output calendar
 
@@ -104,110 +100,97 @@
 
 !--- RUGOSITY TREATMENT --------------------------------------------------------
-  WRITE(lunout,*) 'Traitement de la rugosite'
-  CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
+  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite'
+  varname='RUGOS'
+  CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
 
 !--- OCEAN TREATMENT -----------------------------------------------------------
-  PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.
+  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique'
 
 ! Input SIC file selection
-  icefile='fake'
-  IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic)
-  IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic)
-  IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic)
-  SELECT CASE(icefile)
-    CASE(famipsic); dumstr='Amip.'
-    CASE(fclimsic); dumstr='Amip climatologique.'
-    CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
-    CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1)
-  END SELECT
+! Open file only to test if available
+  IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     icefile=TRIM(famipsic)
+     varname='sicbcs'
+  ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     icefile=TRIM(fcpldsic)
+     varname='SIICECOV'
+  ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     icefile=TRIM(fhistsic)
+     varname='pourc_sic'
+  ELSE
+     WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
+     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic)
+     CALL abort_gcm('limit_netcdf','No sea-ice file was found',1)
+  END IF
   ierr=NF90_CLOSE(nid)
-  WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)
-
-  CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL)
+  IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile)
+
+  CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
 
   ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
   DO k=1,ndays
-    fi_ice=phy_ice(:,k)
-    WHERE(fi_ice>=1.0  ) fi_ice=1.0
-    WHERE(fi_ice<EPSFRA) fi_ice=0.0
-    IF(fracterre) THEN
-      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
-      pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
-      IF(lCPL) THEN                               ! SIC=pICE*(1-LIC-TER)
-        pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter))
-      ELSE                                        ! SIC=pICE-LIC
+     fi_ice=phy_ice(:,k)
+     WHERE(fi_ice>=1.0  ) fi_ice=1.0
+     WHERE(fi_ice<EPSFRA) fi_ice=0.0
+     pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
+     pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
+     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
+        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
+     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
+        pctsrf_t(:,is_sic,k)=fi_ice(:)
+     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
         pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
-      END IF
-      WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
-      WHERE(1.0-zmasq<EPSFRA)
+     END IF
+     WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
+     WHERE(1.0-zmasq<EPSFRA)
         pctsrf_t(:,is_sic,k)=0.0
         pctsrf_t(:,is_oce,k)=0.0
-      ELSEWHERE
+     ELSEWHERE
         WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
-          pctsrf_t(:,is_sic,k)=1.0-zmasq
-          pctsrf_t(:,is_oce,k)=0.0
+           pctsrf_t(:,is_sic,k)=1.0-zmasq
+           pctsrf_t(:,is_oce,k)=0.0
         ELSEWHERE
-          pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
-          WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
-            pctsrf_t(:,is_oce,k)=0.0
-            pctsrf_t(:,is_sic,k)=1.0-zmasq
-          END WHERE
+           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
+           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
+              pctsrf_t(:,is_oce,k)=0.0
+              pctsrf_t(:,is_sic,k)=1.0-zmasq
+           END WHERE
         END WHERE
-      END WHERE
-      nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
-      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
-      nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
-      IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
-    ELSE 
-      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)
-      WHERE(NINT(pctsrf(:,is_ter))==1)
-        pctsrf_t(:,is_sic,k)=0.
-        pctsrf_t(:,is_oce,k)=0.                  
-        WHERE(fi_ice>=1.e-5)
-          pctsrf_t(:,is_lic,k)=fi_ice
-          pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k)
-        ELSEWHERE
-          pctsrf_t(:,is_lic,k)=0.0
-        END WHERE
-      ELSEWHERE
-        pctsrf_t(:,is_lic,k) = 0.0 
-        WHERE(fi_ice>=1.e-5)
-          pctsrf_t(:,is_ter,k)=0.0
-          pctsrf_t(:,is_sic,k)=fi_ice
-          pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k)
-        ELSEWHERE
-          pctsrf_t(:,is_sic,k)=0.0
-          pctsrf_t(:,is_oce,k)=1.0
-        END WHERE
-      END WHERE
-      verif=sum(pctsrf_t(:,:,k),dim=2)
-      nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5)
-      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
-    END IF 
+     END WHERE
+     nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
+     IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
+     nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
+     IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
   END DO
   DEALLOCATE(phy_ice)
 
 !--- SST TREATMENT -------------------------------------------------------------
-  WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.
+  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst'
 
 ! Input SST file selection
-  sstfile='fake'
-  IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst)
-  IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst)
-  IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst)
-  SELECT CASE(icefile)
-    CASE(famipsic); dumstr='Amip.'
-    CASE(fclimsic); dumstr='Amip climatologique.'
-    CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
-    CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1)
-  END SELECT
+! Open file only to test if available
+  IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     sstfile=TRIM(famipsst)
+     varname='tosbcs'
+  ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     sstfile=TRIM(fcpldsst)
+     varname='SISUTESW'
+  ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+     sstfile=TRIM(fhistsst)
+     varname='tsol_oce'
+  ELSE
+     WRITE(lunout,*) 'ERROR! No sst input file was found.'
+     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst)
+     CALL abort_gcm('limit_netcdf','No sst file was found',1)
+  END IF
   ierr=NF90_CLOSE(nid)
-  WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)
-
-  CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)
+  IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile)
+
+  CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
 
 !--- ALBEDO TREATMENT ----------------------------------------------------------
-  WRITE(lunout,*) 'Traitement de l albedo'
-  CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb)
+  IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo'
+  varname='ALBEDO'
+  CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb)
 
 !--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
@@ -215,5 +198,5 @@
 
 !--- OUTPUT FILE WRITING -------------------------------------------------------
-  WRITE(lunout,*) 'Ecriture du fichier limit'
+  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut'
 
   !--- File creation
@@ -264,4 +247,6 @@
   ierr=NF90_CLOSE(nid)
 
+  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin'
+
   DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
 
@@ -276,5 +261,5 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
+SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask)
 !
 !-----------------------------------------------------------------------------
@@ -304,11 +289,11 @@
 ! Arguments:
   CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
+  CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
   CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
   LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
   INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
-  REAL,    POINTER,  DIMENSION(:, :) :: champo   ! output field = f(t)
+  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
   LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
   REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
-  LOGICAL, OPTIONAL, INTENT(IN)     :: lCPL     ! Coupled model flag (for ICE)
 !------------------------------------------------------------------------------
 ! Local variables:
@@ -316,5 +301,4 @@
   INTEGER :: ncid, varid                  ! NetCDF identifiers
   CHARACTER(LEN=30)               :: dnam       ! dimension name
-  CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
 !--- dimensions
   INTEGER,           DIMENSION(4) :: dids       ! NetCDF dimensions identifiers
@@ -331,4 +315,5 @@
 !--- input files
   CHARACTER(LEN=20)                 :: cal_in   ! calendar
+  CHARACTER(LEN=20)                 :: unit_sic ! attribute unit in sea-ice file
   INTEGER                           :: ndays_in ! number of days
 !--- misc
@@ -337,26 +322,45 @@
   CHARACTER(LEN=25)                 :: title    ! for messages
   LOGICAL                           :: extrp    ! flag for extrapolation
+  LOGICAL                           :: oldice   ! flag for old way ice computation 
   REAL                              :: chmin, chmax
   INTEGER ierr
   integer n_extrap ! number of extrapolated points
   logical skip
+
 !------------------------------------------------------------------------------
 !---Variables depending on keyword 'mode' -------------------------------------
   NULLIFY(champo)
+
   SELECT CASE(mode)
-    CASE('RUG'); varname='RUGOS';  title='Rugosite'
-    CASE('SIC'); varname='sicbcs'; title='Sea-ice'
-    CASE('SST'); varname='tosbcs'; title='SST'
-    CASE('ALB'); varname='ALBEDO'; title='Albedo'
+  CASE('RUG'); title='Rugosite'
+  CASE('SIC'); title='Sea-ice'
+  CASE('SST'); title='SST'
+  CASE('ALB'); title='Albedo'
   END SELECT
+  
+
   extrp=.FALSE. 
+  oldice=.FALSE.
   IF ( PRESENT(flag) ) THEN 
     IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 
+    IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. 
   END IF
 
 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
+  IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam
   ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
-  ierr=NF90_INQ_VARID(ncid, varname, varid);            CALL ncerr(ierr, fnam)
+  ierr=NF90_INQ_VARID(ncid, trim(varname), varid);            CALL ncerr(ierr, fnam)
   ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
+
+!--- Read unit for sea-ice variable only
+  IF (mode=='SIC') THEN
+     ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic)
+     IF(ierr/=NF90_NOERR) THEN
+        IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value'
+        unit_sic='X'
+     ELSE
+        IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic
+     END IF
+  END IF
 
 !--- Longitude
@@ -365,5 +369,5 @@
   ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
   ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
-  WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
+  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
 
 !--- Latitude
@@ -372,5 +376,5 @@
   ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
   ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
-  WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
+  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
 
 !--- Time (variable is not needed - it is rebuilt - but calendar is)
@@ -385,10 +389,11 @@
       CASE('SIC', 'SST'); cal_in='gregorian'
     END SELECT
-    WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
+    IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
          // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
   END IF
-  WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
+  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
        cal_in
 
+  
 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
   !--- Determining input file number of days, depending on calendar
@@ -398,7 +403,6 @@
 !--- If input records are not monthly, time sampling has to be constant !
   timeyear=mid_months(anneeref, cal_in, lmdep)
-  IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' &
-       // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, &
-       ' enregistrements.'
+  IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
+       ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
 
 !--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
@@ -406,7 +410,6 @@
   IF(extrp) ALLOCATE(work(imdep, jmdep))
 
-  WRITE(lunout, *)
-  WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, &
-       ' CHAMPS.'
+  IF (prt_level>5) WRITE(lunout, *)
+  IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.'
   ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
   DO l=1, lmdep
@@ -419,12 +422,9 @@
          work)
 
-    IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
-      IF(l==1) THEN
-        WRITE(lunout, *)                                                      &
-  '-------------------------------------------------------------------------'
-        WRITE(lunout, *)                                                     &
-  'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
-        WRITE(lunout, *)                                                      &
-  '-------------------------------------------------------------------------'
+    IF(ibar .AND. .NOT.oldice) THEN
+      IF(l==1 .AND. prt_level>5) THEN
+        WRITE(lunout, *) '-------------------------------------------------------------------------'
+        WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$'
+        WRITE(lunout, *) '-------------------------------------------------------------------------'
       END IF
       IF(mode=='RUG') champ=LOG(champ)
@@ -453,8 +453,11 @@
 
 !--- TIME INTERPOLATION ------------------------------------------------------
-  WRITE(lunout, *)
-  WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
-  WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear
-  WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays
+  IF (prt_level>5) THEN
+     WRITE(lunout, *)
+     WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
+     WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
+     WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
+  END IF
+
   ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
   skip = .false.
@@ -471,5 +474,5 @@
   END DO
   if (n_extrap /= 0) then
-     print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap
+     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
   end if
   champan(iip1, :, :)=champan(1, :, :)
@@ -479,10 +482,10 @@
   DO j=1, jjp1
     CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
-    WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j
+    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
   END DO
 
 !--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
   IF(mode=='SST') THEN
-    WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
+    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
     WHERE(champan<271.38) champan=271.38
   END IF
@@ -490,10 +493,20 @@
 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
   IF(mode=='SIC') THEN
-    WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
-    IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100.
+    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
+
+    IF (unit_sic=='1') THEN
+       ! Nothing to be done for sea-ice field is already in fraction of 1
+       ! This is the case for sea-ice in file cpl_atm_sic.nc
+       IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1'
+    ELSE
+       ! Convert sea ice from percentage to fraction of 1
+       IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' 
+       champan(:, :, :)=champan(:, :, :)/100.
+    END IF
+
     champan(iip1, :, :)=champan(1, :, :)
     WHERE(champan>1.0) champan=1.0
     WHERE(champan<0.0) champan=0.0
-  END IF
+ END IF
 
 !--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
