Index: /LMDZ5/trunk/arch/arch-IA64_PLATINE.fcm
===================================================================
--- /LMDZ5/trunk/arch/arch-IA64_PLATINE.fcm	(revision 1453)
+++ /LMDZ5/trunk/arch/arch-IA64_PLATINE.fcm	(revision 1454)
@@ -5,5 +5,5 @@
 %FPP_FLAGS           -P -traditional
 %FPP_DEF             NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM FFT_MKL
-%BASE_FFLAGS         -i4 -r8 -automatic -align all -I/applications/intel/mkl/10.0.1.014/include
+%BASE_FFLAGS         -i4 -r8 -automatic -align all -I${MKLROOT}/include
 %PROD_FFLAGS         -O3
 %DEV_FFLAGS          -p -g -O3 -traceback
@@ -11,5 +11,5 @@
 %MPI_FFLAGS
 %OMP_FFLAGS          -openmp
-%BASE_LD             -p -i4 -r8 -automatic -L/applications/intel/mkl/10.0.1.014/lib/64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lguide
+%BASE_LD             -p -i4 -r8 -automatic $MKL_LIBS
 %MPI_LD
 %OMP_LD              -openmp
Index: /LMDZ5/trunk/arch/arch-X64_TITANE.fcm
===================================================================
--- /LMDZ5/trunk/arch/arch-X64_TITANE.fcm	(revision 1453)
+++ /LMDZ5/trunk/arch/arch-X64_TITANE.fcm	(revision 1454)
@@ -11,5 +11,5 @@
 %MPI_FFLAGS
 %OMP_FFLAGS          -openmp
-%BASE_LD             -p -i4 -r8 -automatic -L/applications/intel/mkl/10.1.1.019/lib/em64t -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lguide
+%BASE_LD             -p -i4 -r8 -automatic $MKL_LIBS
 %MPI_LD
 %OMP_LD              -openmp
Index: /LMDZ5/trunk/libf/cosp/icarus.F
===================================================================
--- /LMDZ5/trunk/libf/cosp/icarus.F	(revision 1453)
+++ /LMDZ5/trunk/libf/cosp/icarus.F	(revision 1454)
@@ -450,5 +450,4 @@
               write (6,*) 'rangevec:'
               write (6,*) rangevec
-              call flush(6)
               STOP
         endif
Index: /LMDZ5/trunk/libf/dyn3d/academic.h
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/academic.h	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/academic.h	(revision 1454)
@@ -1,5 +1,9 @@
 !
-! $Header$
+! $Id$
 !
-      real tetarappel(ip1jmp1,llm),taurappel
-      common/academic/tetarappel,taurappel
+      common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4
+      real :: tetarappel(ip1jmp1,llm)
+      real :: knewt_t(llm)
+      real :: kfrict(llm)
+      real :: knewt_g
+      real :: clat4(ip1jmp1)
Index: /LMDZ5/trunk/libf/dyn3d/addfi.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/addfi.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/addfi.F	(revision 1454)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
       SUBROUTINE addfi(pdt, leapf, forward,
@@ -7,4 +7,5 @@
 
       USE infotrac, ONLY : nqtot
+      USE control_mod, ONLY : planet_type
       IMPLICIT NONE
 c
@@ -116,5 +117,7 @@
       ENDDO
  
-      DO iq = 1, 2
+      if (planet_type=="earth") then
+      ! earth case, special treatment for first 2 tracers (water)
+       DO iq = 1, 2
          DO k = 1,llm
             DO j = 1,ip1jmp1
@@ -123,7 +126,7 @@
             ENDDO
          ENDDO
-      ENDDO
+       ENDDO
 
-      DO iq = 3, nqtot
+       DO iq = 3, nqtot
          DO k = 1,llm
             DO j = 1,ip1jmp1
@@ -132,5 +135,16 @@
             ENDDO
          ENDDO
-      ENDDO
+       ENDDO
+      else
+      ! general case, treat all tracers equally)
+       DO iq = 1, nqtot
+         DO k = 1,llm
+            DO j = 1,ip1jmp1
+               pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt
+               pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt )
+            ENDDO
+         ENDDO
+       ENDDO
+      endif ! of if (planet_type=="earth")
 
 
Index: /LMDZ5/trunk/libf/dyn3d/advtrac.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/advtrac.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/advtrac.F	(revision 1454)
@@ -236,5 +236,4 @@
             call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
 
-
 c   ----------------------------------------------------------------
 c   Schema "pseudo amont" + test sur humidite specifique
Index: /LMDZ5/trunk/libf/dyn3d/caladvtrac.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/caladvtrac.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/caladvtrac.F	(revision 1454)
@@ -8,6 +8,6 @@
      *                   flxw, pk)
 c
-      USE infotrac
-      USE control_mod
+      USE infotrac, ONLY : nqtot
+      USE control_mod, ONLY : iapp_tracvl,planet_type
  
       IMPLICIT NONE
@@ -30,5 +30,6 @@
 c   ----------
       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm),masse(ip1jmp1,llm)
-      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot),dq( ip1jmp1,llm,2 )
+      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot)
+      real :: dq(ip1jmp1,llm,nqtot)
       REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm)
       REAL               :: flxw(ip1jmp1,llm)
@@ -49,13 +50,13 @@
 cc
 c
+! Earth-specific stuff for the first 2 tracers (water)
+      if (planet_type.eq."earth") then
 C initialisation
-        dq = 0.
-
-        CALL SCOPY( 2 * ijp1llm, q, 1, dq, 1 )
-
+        dq(:,:,1:2)=q(:,:,1:2)
+       
 c  test des valeurs minmax
 cc        CALL minmaxq(q(1,1,1),1.e33,-1.e33,'Eau vapeur (a) ')
 cc        CALL minmaxq(q(1,1,2),1.e33,-1.e33,'Eau liquide(a) ')
-
+      endif ! of if (planet_type.eq."earth")
 c   advection
 
@@ -63,7 +64,10 @@
      *       p,  masse,q,iapptrac, teta,
      .       flxw, pk)
+
 c
 
-         IF( iapptrac.EQ.iapp_tracvl ) THEN
+      IF( iapptrac.EQ.iapp_tracvl ) THEN
+        if (planet_type.eq."earth") then
+! Earth-specific treatment for the first 2 tracers (water)
 c
 cc          CALL minmaxq(q(1,1,1),1.e33,-1.e33,'Eau vapeur     ')
@@ -78,8 +82,5 @@
           ENDDO
           
-	  if (planet_type.eq."earth") then
-! Earth-specific treatment of first 2 tracers (water)
-            CALL qminimum( q, 2, finmasse )
-	  endif
+          CALL qminimum( q, 2, finmasse )
 
           CALL SCOPY   ( ip1jmp1*llm, masse, 1, finmasse,       1 )
@@ -100,22 +101,12 @@
            ENDDO
 c
-         ELSE
-           DO iq = 1 , 2
-           DO l  = 1, llm
-             DO ij = 1,ip1jmp1
-              dq(ij,l,iq)  = 0.
-             ENDDO
-           ENDDO
-           ENDDO
+        endif ! of if (planet_type.eq."earth")
+      ELSE
+        if (planet_type.eq."earth") then
+! Earth-specific treatment for the first 2 tracers (water)
+          dq(:,:,1:2)=0.
+        endif ! of if (planet_type.eq."earth")
+      ENDIF ! of IF( iapptrac.EQ.iapp_tracvl )
 
-
-         ENDIF
-
-c
-
-c  ... On appelle  qminimum uniquement  pour l'eau vapeur et liquide  ..
-
- 
-      RETURN
       END
 
Index: /LMDZ5/trunk/libf/dyn3d/comconst.h
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/comconst.h	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/comconst.h	(revision 1454)
@@ -5,18 +5,33 @@
 ! INCLUDE comconst.h
 
-      COMMON/comconst/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,           &
-     & dtvr,daysec,                                                     &
+      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
+     &                 iflag_top_bound
+      COMMON/comconstr/dtvr,daysec,                                     &
      & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
      &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
-     &                   ,iflag_top_bound,tau_top_bound
+     &                   ,tau_top_bound,                                &
+     & daylen,year_day,molmass
 
 
       INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
-      REAL dtvr,daysec
-      REAL pi,dtphys,dtdiss,rad,r,cpp,kappa
-      REAL cotot,unsim,g,omeg
+      REAL dtvr ! dynamical time step (in s)
+      REAL daysec !length (in s) of a standard day
+      REAL pi    ! something like 3.14159....
+      REAL dtphys ! (s) time step for the physics
+      REAL dtdiss ! (s) time step for the dissipation
+      REAL rad ! (m) radius of the planet
+      REAL r ! Gas constant R=8.31 J.K-1.mol-1
+      REAL cpp   ! Cp
+      REAL kappa ! kappa=R/Cp 
+      REAL cotot
+      REAL unsim ! = 1./iim
+      REAL g ! (m/s2) gravity
+      REAL omeg ! (rad/s) rotation rate of the planet
       REAL dissip_factz,dissip_deltaz,dissip_zref
       INTEGER iflag_top_bound
       REAL tau_top_bound
+      REAL daylen ! length of solar day, in 'standard' day length
+      REAL year_day ! Number of standard days in a year
+      REAL molmass ! (g/mol) molar mass of the atmosphere
 
 
Index: /LMDZ5/trunk/libf/dyn3d/conf_planete.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/conf_planete.F90	(revision 1454)
+++ /LMDZ5/trunk/libf/dyn3d/conf_planete.F90	(revision 1454)
@@ -0,0 +1,66 @@
+!
+! $Id$
+!
+SUBROUTINE conf_planete
+!
+#ifdef CPP_IOIPSL
+USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+USE ioipsl_getincom
+#endif
+IMPLICIT NONE
+!
+!
+!   Declarations :
+!   --------------
+#include "dimensions.h"
+#include "comconst.h"
+#include "comvert.h"
+!
+!   local:
+!   ------
+
+! ---------------------------------------------
+! Initialisations de constantes de la dynamique
+! ---------------------------------------------
+! Pi
+pi=2.*asin(1.)
+
+!Reference surface pressure (Pa)
+preff=101325.
+CALL getin('preff', preff)
+! Reference pressure at which hybrid coord. become purely pressure
+! pa=50000.
+pa=preff/2.
+CALL getin('pa', pa)
+! Gravity
+g=9.80665
+CALL getin('g',g)
+! Molar mass of the atmosphere
+molmass = 28.9644
+CALL getin('molmass',molmass)
+! kappa=R/Cp et Cp      
+kappa = 2./7.
+CALL getin('kappa',kappa)
+cpp=8.3145/molmass/kappa*1000.
+CALL getin('cpp',cpp)
+! Radius of the planet
+rad = 6371229. 
+CALL getin('radius',rad)
+! Length of a standard day (s)
+daysec=86400.
+CALL getin('daysec',daysec)
+! Rotation rate of the planet:
+! Length of a solar day, in standard days
+daylen = 1.
+CALL getin('daylen',daylen)
+! Number of days (standard) per year:
+year_day = 365.25
+CALL getin('year_day',year_day)
+! Omega
+! omeg=2.*pi/86400.
+omeg=2.*pi/daysec*(1./daylen+1./year_day)
+CALL getin('omeg',omeg)
+
+END SUBROUTINE conf_planete
Index: /LMDZ5/trunk/libf/dyn3d/ener.h
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/ener.h	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/ener.h	(revision 1454)
@@ -1,14 +1,18 @@
 !
-! $Header$
+! $Id$
 !
-c-----------------------------------------------------------------------
-c INCLUDE 'ener.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 'ener.h'
 
-      COMMON/ener/ang0,etot0,ptot0,ztot0,stot0,
-     *            ang,etot,ptot,ztot,stot,rmsdpdt ,
-     *            rmsv,gtot(llmm1)
+      COMMON/ener/ang0,etot0,ptot0,ztot0,stot0,                         &
+     &            ang,etot,ptot,ztot,stot,rmsdpdt ,                     &
+     &            rmsv,gtot(llmm1)
 
-      REAL ang0,etot0,ptot0,ztot0,stot0,
-     s     ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot
+      REAL ang0,etot0,ptot0,ztot0,stot0,                                &
+     &     ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot
 
-c-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
Index: /LMDZ5/trunk/libf/dyn3d/friction.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/friction.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/friction.F	(revision 1454)
@@ -6,18 +6,24 @@
 
       USE control_mod
- 
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      USE ioipsl_getincom
+#endif
+      
       IMPLICIT NONE
 
-c=======================================================================
-c
-c
-c   Objet:
-c   ------
-c
-c  ***********
-c    Friction
-c  ***********
-c
-c=======================================================================
+!=======================================================================
+!
+!   Friction for the Newtonian case:
+!   --------------------------------
+!    2 possibilities (depending on flag 'friction_type'
+!     friction_type=0 : A friction that is only applied to the lowermost
+!                       atmospheric layer
+!     friction_type=1 : Friction applied on all atmospheric layer (but
+!       (default)       with stronger magnitude near the surface; see
+!                       iniacademic.F)
+!=======================================================================
 
 #include "dimensions.h"
@@ -25,15 +31,36 @@
 #include "comgeom2.h"
 #include "comconst.h"
+#include "iniprint.h"
+#include "academic.h"
 
-      REAL pdt
+! arguments:
+      REAL,INTENT(out) :: ucov( iip1,jjp1,llm )
+      REAL,INTENT(out) :: vcov( iip1,jjm,llm )
+      REAL,INTENT(in) :: pdt ! time step
+
+! local variables:
+
       REAL modv(iip1,jjp1),zco,zsi
       REAL vpn,vps,upoln,upols,vpols,vpoln
       REAL u2(iip1,jjp1),v2(iip1,jjm)
-      REAL ucov( iip1,jjp1,llm ),vcov( iip1,jjm,llm )
-      INTEGER  i,j
-      REAL cfric
-      parameter (cfric=1.e-5)
+      INTEGER  i,j,l
+      REAL,PARAMETER :: cfric=1.e-5
+      LOGICAL,SAVE :: firstcall=.true.
+      INTEGER,SAVE :: friction_type=1
+      CHARACTER(len=20) :: modname="friction"
+      CHARACTER(len=80) :: abort_message
+      
+      IF (firstcall) THEN
+        ! set friction type
+        call getin("friction_type",friction_type)
+        if ((friction_type.lt.0).or.(friction_type.gt.1)) then
+          abort_message="wrong friction type"
+          write(lunout,*)'Friction: wrong friction type',friction_type
+          call abort_gcm(modname,abort_message,42)
+        endif
+        firstcall=.false.
+      ENDIF
 
-
+      if (friction_type.eq.0) then
 c   calcul des composantes au carre du vent naturel
       do j=1,jjp1
@@ -96,5 +123,13 @@
          vcov(iip1,j,1)=vcov(1,j,1)
       enddo
+      endif ! of if (friction_type.eq.0)
 
+      if (friction_type.eq.1) then
+        do l=1,llm
+          ucov(:,:,l)=ucov(:,:,l)*(1.-pdt*kfrict(l))
+          vcov(:,:,l)=vcov(:,:,l)*(1.-pdt*kfrict(l))
+        enddo
+      endif
+      
       RETURN
       END
Index: /LMDZ5/trunk/libf/dyn3d/gcm.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/gcm.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/gcm.F	(revision 1454)
@@ -249,18 +249,9 @@
         endif
 
-        if (planet_type.eq."earth") then
-#ifdef CPP_EARTH
+!        if (planet_type.eq."earth") then
 ! Load an Earth-format start file
          CALL dynetat0("start.nc",vcov,ucov,
      &              teta,q,masse,ps,phis, time_0)
-#else
-        ! SW model also has Earth-format start files
-        ! (but can be used without the CPP_EARTH directive)
-          if (iflag_phys.eq.0) then
-            CALL dynetat0("start.nc",vcov,ucov,
-     &              teta,q,masse,ps,phis, time_0)
-          endif
-#endif
-        endif ! of if (planet_type.eq."earth")
+!        endif ! of if (planet_type.eq."earth")
         
 c       write(73,*) 'ucov',ucov
@@ -468,7 +459,8 @@
 #endif
 
-      if (planet_type.eq."earth") then
+!      if (planet_type.eq."earth") then
+! Write an Earth-format restart file
         CALL dynredem0("restart.nc", day_end, phis)
-      endif
+!      endif
 
       ecripar = .TRUE.
Index: /LMDZ5/trunk/libf/dyn3d/grid_noro.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/grid_noro.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/grid_noro.F	(revision 1454)
@@ -458,11 +458,8 @@
 C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
 
-      PARAMETER (ISMo=400,JSMo=200)
-      REAL X(IMAR,JMAR),XF(ISMo,JSMo)
+      REAL X(IMAR,JMAR),XF(IMAR,JMAR)
       real WEIGHTpb(-1:1,-1:1)
 
-      if(imar.gt.ismo) stop'surdimensionner ismo dans mva9 (grid_noro)'
-      if(jmar.gt.jsmo) stop'surdimensionner jsmo dans mva9 (grid_noro)'
-      
+
       SUM=0.
       DO IS=-1,1
Index: /LMDZ5/trunk/libf/dyn3d/infotrac.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/infotrac.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/infotrac.F90	(revision 1454)
@@ -31,7 +31,5 @@
 
   SUBROUTINE infotrac_init
-
     USE control_mod
- 
     IMPLICIT NONE
 !=======================================================================
@@ -63,10 +61,10 @@
     CHARACTER(len=1), DIMENSION(3)  :: txts
     CHARACTER(len=2), DIMENSION(9)  :: txtp
-    CHARACTER(len=13)               :: str1,str2
+    CHARACTER(len=23)               :: str1,str2
   
     INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
     INTEGER :: iq, new_iq, iiq, jq, ierr
-    INTEGER, EXTERNAL :: lnblnk
- 
+
+    character(len=*),parameter :: modname="infotrac_init"
 !-----------------------------------------------------------------------
 ! Initialization :
@@ -102,13 +100,23 @@
        OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
        IF(ierr.EQ.0) THEN
-          WRITE(lunout,*) 'Open traceur.def : ok'
+          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
           READ(90,*) nqtrue
        ELSE 
-          WRITE(lunout,*) 'Problem in opening traceur.def'
-          WRITE(lunout,*) 'ATTENTION using defaut values'
-          nqtrue=4 ! Defaut value
-       END IF
-       ! Attention! Only for planet_type=='earth'
-       nbtr=nqtrue-2
+          WRITE(lunout,*) trim(modname),': Problem in opening traceur.def'
+          WRITE(lunout,*) trim(modname),': WARNING using defaut values'
+          if (planet_type=='earth') then
+            nqtrue=4 ! Default value for Earth
+          else
+            nqtrue=1 ! Default value for other planets
+          endif
+       END IF
+       if ( planet_type=='earth') then
+         ! For Earth, water vapour & liquid tracers are not in the physics
+         nbtr=nqtrue-2
+       else
+         ! Other planets (for now); we have the same number of tracers
+         ! in the dynamics than in the physics
+         nbtr=nqtrue
+       endif
     ELSE
        ! nbtr has been read from INCA by init_cont_lmdz() in gcm.F 
@@ -116,6 +124,6 @@
     END IF
 
-    IF (nqtrue < 2) THEN
-       WRITE(lunout,*) 'nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum'
+    IF ((planet_type=="earth").and.(nqtrue < 2)) THEN
+       WRITE(lunout,*) trim(modname),': nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum'
        CALL abort_gcm('infotrac_init','Not enough tracers',1)
     END IF
@@ -158,8 +166,10 @@
           ! Continue to read tracer.def
           DO iq=1,nqtrue
-             READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
+             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
           END DO
           CLOSE(90)  
-       ELSE ! Without tracer.def
+       ELSE ! Without tracer.def, set default values 
+         if (planet_type=="earth") then
+          ! for Earth, default is to have 4 tracers
           hadv(1) = 14
           vadv(1) = 14
@@ -174,8 +184,13 @@
           vadv(4) = 10
           tnom_0(4) = 'PB'
+         else ! default for other planets
+          hadv(1) = 10
+          vadv(1) = 10
+          tnom_0(1) = 'dummy'
+         endif ! of if (planet_type=="earth")
        END IF
        
-       WRITE(lunout,*) 'Valeur de traceur.def :'
-       WRITE(lunout,*) 'nombre de traceurs ',nqtrue
+       WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
+       WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue
        DO iq=1,nqtrue
           WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq)
@@ -219,5 +234,5 @@
           new_iq=new_iq+10 ! 9 tracers added
        ELSE
-          WRITE(lunout,*) 'This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
+          WRITE(lunout,*) trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
           CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1)
        END IF
@@ -229,8 +244,8 @@
        nqtot = new_iq
 
-       WRITE(lunout,*) 'The choice of advection schema for one or more tracers'
+       WRITE(lunout,*) trim(modname),': The choice of advection schema for one or more tracers'
        WRITE(lunout,*) 'makes it necessary to add tracers'
-       WRITE(lunout,*) nqtrue,' is the number of true tracers'
-       WRITE(lunout,*) nqtot, ' is the total number of tracers needed'
+       WRITE(lunout,*) trim(modname)//': ',nqtrue,' is the number of true tracers'
+       WRITE(lunout,*) trim(modname)//': ',nqtot, ' is the total number of tracers needed'
 
     ELSE
@@ -260,5 +275,6 @@
           iadv(new_iq)=11
        ELSE
-          WRITE(lunout,*)'This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
+          WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
+
           CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1)
        END IF
@@ -267,7 +283,7 @@
        tname(new_iq)= tnom_0(iq)
        IF (iadv(new_iq)==0) THEN
-          ttext(new_iq)=str1(1:lnblnk(str1))
+          ttext(new_iq)=trim(str1)
        ELSE
-          ttext(new_iq)=str1(1:lnblnk(str1))//descrq(iadv(new_iq))
+          ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
        END IF
 
@@ -278,6 +294,6 @@
              new_iq=new_iq+1
              iadv(new_iq)=-20
-             ttext(new_iq)=str2(1:lnblnk(str2))//txts(jq)
-             tname(new_iq)=str1(1:lnblnk(str1))//txts(jq)
+             ttext(new_iq)=trim(str2)//txts(jq)
+             tname(new_iq)=trim(str1)//txts(jq)
           END DO
        ELSE IF (iadv(new_iq)==30) THEN
@@ -285,6 +301,6 @@
              new_iq=new_iq+1
              iadv(new_iq)=-30
-             ttext(new_iq)=str2(1:lnblnk(str2))//txtp(jq)
-             tname(new_iq)=str1(1:lnblnk(str1))//txtp(jq)
+             ttext(new_iq)=trim(str2)//txtp(jq)
+             tname(new_iq)=trim(str1)//txtp(jq)
           END DO
        END IF
@@ -305,8 +321,9 @@
 
 
-    WRITE(lunout,*) 'Information stored in infotrac :'
-    WRITE(lunout,*) 'iadv  niadv tname  ttext :'
+    WRITE(lunout,*) trim(modname),': Information stored in infotrac :'
+    WRITE(lunout,*) trim(modname),': iadv  niadv tname  ttext :'
     DO iq=1,nqtot
-       WRITE(lunout,*) iadv(iq),niadv(iq), tname(iq), ttext(iq)
+       WRITE(lunout,*) iadv(iq),niadv(iq),&
+       ' ',trim(tname(iq)),' ',trim(ttext(iq))
     END DO
 
@@ -317,8 +334,8 @@
     DO iq=1,nqtot
        IF (iadv(iq)/=10 .AND. iadv(iq)/=14 .AND. iadv(iq)/=0) THEN
-          WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
+          WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
           CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1)
        ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN
-          WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
+          WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
           CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1)
        END IF
@@ -331,6 +348,4 @@
     DEALLOCATE(tracnam)
 
-999 FORMAT (i2,1x,i2,1x,a15)
-
   END SUBROUTINE infotrac_init
 
Index: /LMDZ5/trunk/libf/dyn3d/iniacademic.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/iniacademic.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/iniacademic.F	(revision 1454)
@@ -8,5 +8,12 @@
       USE filtreg_mod
       USE infotrac, ONLY : nqtot
-      USE control_mod
+      USE control_mod, ONLY: day_step,planet_type
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      USE ioipsl_getincom
+#endif
+      USE Write_Field
  
 
@@ -70,8 +77,14 @@
       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL ddsin,tetarappelj,tetarappell,zsig
+      REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
       real tetajl(jjp1,llm)
       INTEGER i,j,l,lsup,ij
 
+      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
+      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
+      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
+      LOGICAL ok_pv                ! Polar Vortex
+      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
+      
       real zz,ran1
       integer idum
@@ -82,23 +95,17 @@
 ! 1. Initializations for Earth-like case
 ! --------------------------------------
-      if (planet_type=="earth") then
-c
+c
+        ! initialize planet radius, rotation rate,...
+        call conf_planete
+
         time_0=0.
-        day_ref=0
+        day_ref=1
         annee_ref=0
 
         im         = iim
         jm         = jjm
-        day_ini    = 0
-        omeg       = 4.*asin(1.)/86400.
-        rad    = 6371229.
-        g      = 9.8
-        daysec = 86400.
+        day_ini    = 1
         dtvr    = daysec/REAL(day_step)
         zdtvr=dtvr
-        kappa  = 0.2857143
-        cpp    = 1004.70885
-        preff     = 101325.
-        pa        =  50000.
         etot0      = 0.
         ptot0      = 0.
@@ -120,7 +127,5 @@
           if (.not.read_start) then
             phis(:)=0.
-            q(:,:,1)=1.e-10
-            q(:,:,2)=1.e-15
-            q(:,:,3:nqtot)=0.
+            q(:,:,:)=0
             CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
           endif
@@ -129,30 +134,78 @@
         if (iflag_phys.eq.2) then
           ! initializations for the academic case
-          ps(:)=1.e5
-          phis(:)=0.
-c---------------------------------------------------------------------
-
-          taurappel=10.*daysec
-
-c---------------------------------------------------------------------
-c   Calcul de la temperature potentielle :
-c   --------------------------------------
-
+          
+!         if (planet_type=="earth") then
+
+          ! 1. local parameters
+          ! by convention, winter is in the southern hemisphere
+          ! Geostrophic wind or no wind?
+          ok_geost=.TRUE.
+          CALL getin('ok_geost',ok_geost)
+          ! Constants for Newtonian relaxation and friction
+          k_f=1.                !friction 
+          CALL getin('k_j',k_f)
+          k_f=1./(daysec*k_f)
+          k_c_s=4.  !cooling surface
+          CALL getin('k_c_s',k_c_s)
+          k_c_s=1./(daysec*k_c_s)
+          k_c_a=40. !cooling free atm
+          CALL getin('k_c_a',k_c_a)
+          k_c_a=1./(daysec*k_c_a)
+          ! Constants for Teta equilibrium profile
+          teta0=315.     ! mean Teta (S.H. 315K)
+          CALL getin('teta0',teta0)
+          ttp=200.       ! Tropopause temperature (S.H. 200K)
+          CALL getin('ttp',ttp)
+          eps=0.         ! Deviation to N-S symmetry(~0-20K)
+          CALL getin('eps',eps)
+          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
+          CALL getin('delt_y',delt_y)
+          delt_z=10.     ! Vertical Gradient (S.H. 10K)
+          CALL getin('delt_z',delt_z)
+          ! Polar vortex
+          ok_pv=.false.
+          CALL getin('ok_pv',ok_pv)
+          phi_pv=-50.            ! Latitude of edge of vortex
+          CALL getin('phi_pv',phi_pv)
+          phi_pv=phi_pv*pi/180.
+          dphi_pv=5.             ! Width of the edge
+          CALL getin('dphi_pv',dphi_pv)
+          dphi_pv=dphi_pv*pi/180.
+          gam_pv=4.              ! -dT/dz vortex (in K/km)
+          CALL getin('gam_pv',gam_pv)
+          
+          ! 2. Initialize fields towards which to relax
+          ! Friction
+          knewt_g=k_c_a
           DO l=1,llm
-            zsig=ap(l)/preff+bp(l)
-            if (zsig.gt.0.3) then
-             lsup=l
-             tetarappell=1./8.*(-log(zsig)-.5)
-             DO j=1,jjp1
-             ddsin=sin(rlatu(j))-sin(pi/20.)
-             tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell)
-             ENDDO
-            else
-c   Choix isotherme au-dessus de 300 mbar
-             do j=1,jjp1
-               tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
-             enddo
-            endif ! of if (zsig.gt.0.3)
+           zsig=presnivs(l)/preff
+           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
+           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
+          ENDDO
+          DO j=1,jjp1
+            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
+          ENDDO
+          
+          ! Potential temperature 
+          DO l=1,llm
+           zsig=presnivs(l)/preff
+           tetastrat=ttp*zsig**(-kappa)
+           tetapv=tetastrat
+           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
+               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
+           ENDIF
+           DO j=1,jjp1
+             ! Troposphere
+             ddsin=sin(rlatu(j))
+             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
+     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
+             ! Profil stratospherique isotherme (+vortex)
+             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
+             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
+             tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
+           ENDDO
           ENDDO ! of DO l=1,llm
+ 
+!          CALL writefield('theta_eq',tetajl)
 
           do l=1,llm
@@ -165,56 +218,68 @@
           enddo
 
-c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
-
-          CALL pression ( ip1jmp1, ap, bp, ps, p       )
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-          CALL massdair(p,masse)
-
-c  intialisation du vent et de la temperature
-          teta(:,:)=tetarappel(:,:)
-          CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-          call ugeostr(phi,ucov)
-          vcov=0.
-          q(:,:,1   )=1.e-10
-          q(:,:,2   )=1.e-15
-          q(:,:,3:nqtot)=0.
-
-
-c   perturbation aleatoire sur la temperature
-          idum  = -1
-          zz = ran1(idum)
-          idum  = 0
-          do l=1,llm
-            do ij=iip2,ip1jm
-              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+
+!         else
+!          write(lunout,*)"iniacademic: planet types other than earth",
+!     &                   " not implemented (yet)."
+!          stop
+!         endif ! of if (planet_type=="earth")
+
+          ! 3. Initialize fields (if necessary)
+          IF (.NOT. read_start) THEN
+            ! surface pressure
+            ps(:)=preff
+            ! ground geopotential
+            phis(:)=0.
+            
+            CALL pression ( ip1jmp1, ap, bp, ps, p       )
+            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+            CALL massdair(p,masse)
+
+            ! bulk initialization of temperature
+            teta(:,:)=tetarappel(:,:)
+            
+            ! geopotential
+            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
+            
+            ! winds
+            if (ok_geost) then
+              call ugeostr(phi,ucov)
+            else
+              ucov(:,:)=0.
+            endif
+            vcov(:,:)=0.
+            
+            ! bulk initialization of tracers
+            if (planet_type=="earth") then
+              ! Earth: first two tracers will be water
+              do i=1,nqtot
+                if (i.eq.1) q(:,:,i)=1.e-10
+                if (i.eq.2) q(:,:,i)=1.e-15
+                if (i.gt.2) q(:,:,i)=0.
+              enddo
+            else
+              q(:,:,:)=0
+            endif ! of if (planet_type=="earth")
+
+            ! add random perturbation to temperature
+            idum  = -1
+            zz = ran1(idum)
+            idum  = 0
+            do l=1,llm
+              do ij=iip2,ip1jm
+                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+              enddo
             enddo
-          enddo
-
-          do l=1,llm
-            do ij=1,ip1jmp1,iip1
-              teta(ij+iim,l)=teta(ij,l)
+
+            ! maintain periodicity in longitude
+            do l=1,llm
+              do ij=1,ip1jmp1,iip1
+                teta(ij+iim,l)=teta(ij,l)
+              enddo
             enddo
-          enddo
-
-
-
-c     PRINT *,' Appel test_period avec tetarappel '
-c     CALL  test_period ( ucov,vcov,tetarappel,q,p,phis )
-c     PRINT *,' Appel test_period avec teta '
-c     CALL  test_period ( ucov,vcov,teta,q,p,phis )
-
-c   initialisation d'un traceur sur une colonne
-          j=jjp1*3/4
-          i=iip1/2
-          ij=(j-1)*iip1+i
-          q(ij,:,3)=1.
+
+          ENDIF ! of IF (.NOT. read_start)
         endif ! of if (iflag_phys.eq.2)
         
-      else
-        write(lunout,*)"iniacademic: planet types other than earth",
-     &                 " not implemented (yet)."
-        stop
-      endif ! of if (planet_type=="earth")
-      return
       END
 c-----------------------------------------------------------------------
Index: /LMDZ5/trunk/libf/dyn3d/integrd.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/integrd.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/integrd.F	(revision 1454)
@@ -6,5 +6,5 @@
      $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis,finvmaold )
 
-      USE control_mod
+      use control_mod, only : planet_type
 
       IMPLICIT NONE
@@ -81,8 +81,8 @@
       CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
 
-      DO 2 ij = 1,ip1jmp1
+      DO ij = 1,ip1jmp1
        pscr (ij)    = ps(ij)
        ps (ij)      = psm1(ij) + dt * dp(ij)
-   2  CONTINUE
+      ENDDO
 c
       DO ij = 1,ip1jmp1
@@ -115,45 +115,45 @@
 c    ............   integration  de  ucov, vcov,  h     ..............
 
-      DO 10 l = 1,llm
-
-      DO 4 ij = iip2,ip1jm
-      uscr( ij )   =  ucov( ij,l )
-      ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
-   4  CONTINUE
-
-      DO 5 ij = 1,ip1jm
-      vscr( ij )   =  vcov( ij,l )
-      vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
-   5  CONTINUE
-
-      DO 6 ij = 1,ip1jmp1
-      hscr( ij )    =  teta(ij,l)
-      teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l) 
-     $                + dt * dteta(ij,l) / masse(ij,l)
-   6  CONTINUE
+      DO l = 1,llm
+
+       DO ij = iip2,ip1jm
+        uscr( ij )   =  ucov( ij,l )
+        ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
+       ENDDO
+
+       DO ij = 1,ip1jm
+        vscr( ij )   =  vcov( ij,l )
+        vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
+       ENDDO
+
+       DO ij = 1,ip1jmp1
+        hscr( ij )    =  teta(ij,l)
+        teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l) 
+     &                + dt * dteta(ij,l) / masse(ij,l)
+       ENDDO
 
 c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
 c
 c
-      DO  ij   = 1, iim
+       DO  ij   = 1, iim
         tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
         tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
-      ENDDO
+       ENDDO
         tpn      = SSUM(iim,tppn,1)/apoln
         tps      = SSUM(iim,tpps,1)/apols
 
-      DO ij   = 1, iip1
+       DO ij   = 1, iip1
         teta(   ij   ,l)  = tpn
         teta(ij+ip1jm,l)  = tps
-      ENDDO
-c
-
-      IF(leapf)  THEN
+       ENDDO
+c
+
+       IF(leapf)  THEN
          CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
          CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
          CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
-      END IF
-
-  10  CONTINUE
+       END IF
+
+      ENDDO ! of DO l = 1,llm
 
 
@@ -185,14 +185,13 @@
 c$$$      ENDIF
 
-	 if (planet_type.eq."earth") then
+      if (planet_type.eq."earth") then
 ! Earth-specific treatment of first 2 tracers (water)
-          DO l = 1, llm
-           DO ij = 1, ip1jmp1
+        DO l = 1, llm
+          DO ij = 1, ip1jmp1
             deltap(ij,l) =  p(ij,l) - p(ij,l+1) 
-           ENDDO
           ENDDO
-
-          CALL qminimum( q, nq, deltap )
-	 endif ! of if (planet_type.eq."earth")
+        ENDDO
+
+        CALL qminimum( q, nq, deltap )
 
 c
@@ -200,5 +199,5 @@
 c
 
-      DO iq = 1, nq
+       DO iq = 1, nq
         DO l = 1, llm
 
@@ -216,13 +215,13 @@
 
         ENDDO
-      ENDDO
-
-
-         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
+       ENDDO
+
+
+      CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
+
+      endif ! of if (planet_type.eq."earth")
 c
 c
 c     .....   FIN  de l'integration  de   q    .......
-
-15    continue
 
 c    .................................................................
Index: /LMDZ5/trunk/libf/dyn3d/leapfrog.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/leapfrog.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/leapfrog.F	(revision 1454)
@@ -209,5 +209,5 @@
 c   --------------------------------------------------
 
-      dq=0.
+      dq(:,:,:)=0.
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
@@ -253,7 +253,4 @@
       CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
       CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
-
-! Ehouarn: what is this for? zqmin & zqmax are not used anyway ...
-!      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
 
    2  CONTINUE
@@ -430,10 +427,27 @@
 
       IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
-c   Calcul academique de la physique = Rappel Newtonien + friction 
-c   --------------------------------------------------------------
-       teta(:,:)=teta(:,:)
-     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
-       call friction(ucov,vcov,iphysiq*dtvr)
-      ENDIF
+!   Academic case : Simple friction and Newtonan relaxation 
+!   -------------------------------------------------------
+        DO l=1,llm   
+          DO ij=1,ip1jmp1
+           teta(ij,l)=teta(ij,l)-dtvr*
+     &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
+          ENDDO
+        ENDDO ! of DO l=1,llm 
+          call friction(ucov,vcov,dtvr)
+        
+        ! Sponge layer (if any)
+        IF (ok_strato) THEN
+          dufi(:,:)=0.
+          dvfi(:,:)=0.
+          dtetafi(:,:)=0.
+          dqfi(:,:,:)=0.
+          dpfi(:)=0.
+          CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
+          CALL addfi( dtvr, leapf, forward   ,
+     $                  ucov, vcov, teta , q   ,ps ,
+     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
+        ENDIF ! of IF (ok_strato) 
+      ENDIF ! of IF (iflag_phys.EQ.2)
 
 
@@ -615,9 +629,9 @@
 
 
-              if (planet_type.eq."earth") then
+!              if (planet_type.eq."earth") then
 ! Write an Earth-format restart file
                 CALL dynredem1("restart.nc",0.0,
      &                         vcov,ucov,teta,q,masse,ps)
-              endif ! of if (planet_type.eq."earth")
+!              endif ! of if (planet_type.eq."earth")
 
               CLOSE(99)
@@ -727,8 +741,8 @@
 
               IF(itau.EQ.itaufin) THEN
-                if (planet_type.eq."earth") then
+!                if (planet_type.eq."earth") then
                   CALL dynredem1("restart.nc",0.0,
      &                           vcov,ucov,teta,q,masse,ps)
-                endif ! of if (planet_type.eq."earth")
+!                endif ! of if (planet_type.eq."earth")
               ENDIF ! of IF(itau.EQ.itaufin)
 
Index: /LMDZ5/trunk/libf/dyn3d/limit_netcdf.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/limit_netcdf.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3d/limit_netcdf.F90	(revision 1454)
@@ -97,5 +97,5 @@
   kappa = 0.2857143
   cpp   = 1004.70885
-  dtvr  = daysec/FLOAT(day_step)
+  dtvr  = daysec/REAL(day_step)
   CALL inigeom
 
@@ -265,6 +265,4 @@
 
   DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
-#endif
-! of #ifdef CPP_EARTH
 
 
@@ -592,7 +590,7 @@
 
 !--- Mid-months times
-  mid_months(1)=0.5*FLOAT(mnth(1))
+  mid_months(1)=0.5*REAL(mnth(1))
   DO k=2,nm
-    mid_months(k)=mid_months(k-1)+0.5*FLOAT(mnth(k-1)+mnth(k))
+    mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
   END DO
 
@@ -626,4 +624,6 @@
 !-------------------------------------------------------------------------------
 
+#endif
+! of #ifdef CPP_EARTH
 
 END SUBROUTINE limit_netcdf
Index: /LMDZ5/trunk/libf/dyn3dpar/academic.h
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/academic.h	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/academic.h	(revision 1454)
@@ -1,5 +1,9 @@
 !
-! $Header$
+! $Id$
 !
-      real tetarappel(ip1jmp1,llm),taurappel
-      common/academic/tetarappel,taurappel
+      common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4
+      real :: tetarappel(ip1jmp1,llm)
+      real :: knewt_t(llm)
+      real :: kfrict(llm)
+      real :: knewt_g
+      real :: clat4(ip1jmp1)
Index: /LMDZ5/trunk/libf/dyn3dpar/addfi_p.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/addfi_p.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/addfi_p.F	(revision 1454)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
       SUBROUTINE addfi_p(pdt, leapf, forward,
@@ -7,4 +7,5 @@
       USE parallel
       USE infotrac, ONLY : nqtot
+      USE control_mod, ONLY : planet_type
       IMPLICIT NONE
 c
@@ -154,5 +155,7 @@
 c$OMP END MASTER
  
-      DO iq = 1, 2
+      if (planet_type=="earth") then
+      ! earth case, special treatment for first 2 tracers (water)
+       DO iq = 1, 2
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO k = 1,llm
@@ -163,7 +166,7 @@
          ENDDO
 c$OMP END DO NOWAIT
-      ENDDO
-
-      DO iq = 3, nqtot
+       ENDDO
+
+       DO iq = 3, nqtot
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO k = 1,llm
@@ -174,5 +177,18 @@
          ENDDO
 c$OMP END DO NOWAIT
-      ENDDO
+       ENDDO
+      else
+      ! general case, treat all tracers equally)
+       DO iq = 1, nqtot
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO k = 1,llm
+            DO j = ijb,ije
+               pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt
+               pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt )
+            ENDDO
+         ENDDO
+c$OMP END DO NOWAIT
+       ENDDO
+      endif ! of if (planet_type=="earth")
 
 c$OMP MASTER
Index: /LMDZ5/trunk/libf/dyn3dpar/advtrac_p.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/advtrac_p.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/advtrac_p.F	(revision 1454)
@@ -132,5 +132,5 @@
 ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
 c
-      ENDIF
+      ENDIF ! of IF(iadvtr.EQ.0)
 
       iadvtr   = iadvtr+1
@@ -266,5 +266,5 @@
 cym      ----> Revérifier lors de la parallélisation des autres schemas
    
-cym          call massbar_p(massem,massebx,masseby)          
+cym          call massbar_p(massem,massebx,masseby)  
 
           call vlspltgen_p( q,iadv, 2., massem, wg ,
@@ -452,16 +452,18 @@
 c$OMP BARRIER
 
-      ijb=ij_begin
-      ije=ij_end
+      if (planet_type=="earth") then
+
+        ijb=ij_begin
+        ije=ij_end
 
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-       DO l = 1, llm
+        DO l = 1, llm
          DO ij = ijb, ije
            finmasse(ij,l) =  p(ij,l) - p(ij,l+1) 
          ENDDO
-       ENDDO
+        ENDDO
 c$OMP END DO
 
-       CALL qminimum_p( q, 2, finmasse )
+        CALL qminimum_p( q, 2, finmasse )
 
 c------------------------------------------------------------------
@@ -496,4 +498,5 @@
 c$OMP BARRIER	
           iadvtr=0
+      endif ! of if (planet_type=="earth")
        ENDIF ! if iadvtr.EQ.iapp_tracvl
 
Index: /LMDZ5/trunk/libf/dyn3dpar/caladvtrac_p.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/caladvtrac_p.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/caladvtrac_p.F	(revision 1454)
@@ -8,6 +8,6 @@
      *                   flxw, pk, iapptrac)
       USE parallel
-      USE infotrac
-      USE control_mod
+      USE infotrac, ONLY : nqtot
+      USE control_mod, ONLY : iapp_tracvl,planet_type
 c
       IMPLICIT NONE
@@ -30,5 +30,6 @@
 c   ----------
       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm),masse(ip1jmp1,llm)
-      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot),dq( ip1jmp1,llm,2 )
+      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot)
+      real :: dq( ip1jmp1,llm,nqtot)
       REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm)
       REAL               :: flxw(ip1jmp1,llm)
Index: /LMDZ5/trunk/libf/dyn3dpar/comconst.h
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/comconst.h	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/comconst.h	(revision 1454)
@@ -5,18 +5,33 @@
 ! INCLUDE comconst.h
 
-      COMMON/comconst/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,           &
-     & dtvr,daysec,                                                     &
+      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
+     &                 iflag_top_bound
+      COMMON/comconstr/dtvr,daysec,                                     &
      & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
      &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
-     &                   ,iflag_top_bound,tau_top_bound
+     &                   ,tau_top_bound,                                &
+     & daylen,year_day,molmass
 
 
       INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
-      REAL dtvr,daysec
-      REAL pi,dtphys,dtdiss,rad,r,cpp,kappa
-      REAL cotot,unsim,g,omeg
+      REAL dtvr ! dynamical time step (in s)
+      REAL daysec !length (in s) of a standard day
+      REAL pi    ! something like 3.14159....
+      REAL dtphys ! (s) time step for the physics
+      REAL dtdiss ! (s) time step for the dissipation
+      REAL rad ! (m) radius of the planet
+      REAL r ! Gas constant R=8.31 J.K-1.mol-1
+      REAL cpp   ! Cp
+      REAL kappa ! kappa=R/Cp 
+      REAL cotot
+      REAL unsim ! = 1./iim
+      REAL g ! (m/s2) gravity
+      REAL omeg ! (rad/s) rotation rate of the planet
       REAL dissip_factz,dissip_deltaz,dissip_zref
       INTEGER iflag_top_bound
       REAL tau_top_bound
+      REAL daylen ! length of solar day, in 'standard' day length
+      REAL year_day ! Number of standard days in a year
+      REAL molmass ! (g/mol) molar mass of the atmosphere
 
 
Index: /LMDZ5/trunk/libf/dyn3dpar/conf_gcm.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/conf_gcm.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/conf_gcm.F	(revision 1454)
@@ -578,5 +578,13 @@
        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
@@ -768,4 +776,12 @@
        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
Index: /LMDZ5/trunk/libf/dyn3dpar/conf_planete.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/conf_planete.F90	(revision 1454)
+++ /LMDZ5/trunk/libf/dyn3dpar/conf_planete.F90	(revision 1454)
@@ -0,0 +1,66 @@
+!
+! $Id$
+!
+SUBROUTINE conf_planete
+!
+#ifdef CPP_IOIPSL
+USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+USE ioipsl_getincom
+#endif
+IMPLICIT NONE
+!
+!
+!   Declarations :
+!   --------------
+#include "dimensions.h"
+#include "comconst.h"
+#include "comvert.h"
+!
+!   local:
+!   ------
+
+! ---------------------------------------------
+! Initialisations de constantes de la dynamique
+! ---------------------------------------------
+! Pi
+pi=2.*asin(1.)
+
+!Reference surface pressure (Pa)
+preff=101325.
+CALL getin('preff', preff)
+! Reference pressure at which hybrid coord. become purely pressure
+! pa=50000.
+pa=preff/2.
+CALL getin('pa', pa)
+! Gravity
+g=9.80665
+CALL getin('g',g)
+! Molar mass of the atmosphere
+molmass = 28.9644
+CALL getin('molmass',molmass)
+! kappa=R/Cp et Cp      
+kappa = 2./7.
+CALL getin('kappa',kappa)
+cpp=8.3145/molmass/kappa*1000.
+CALL getin('cpp',cpp)
+! Radius of the planet
+rad = 6371229. 
+CALL getin('radius',rad)
+! Length of a standard day (s)
+daysec=86400.
+CALL getin('daysec',daysec)
+! Rotation rate of the planet:
+! Length of a solar day, in standard days
+daylen = 1.
+CALL getin('daylen',daylen)
+! Number of days (standard) per year:
+year_day = 365.25
+CALL getin('year_day',year_day)
+! Omega
+! omeg=2.*pi/86400.
+omeg=2.*pi/daysec*(1./daylen+1./year_day)
+CALL getin('omeg',omeg)
+
+END SUBROUTINE conf_planete
Index: /LMDZ5/trunk/libf/dyn3dpar/ener.h
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/ener.h	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/ener.h	(revision 1454)
@@ -1,14 +1,18 @@
 !
-! $Header$
+! $Id$
 !
-c-----------------------------------------------------------------------
-c INCLUDE 'ener.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 'ener.h'
 
-      COMMON/ener/ang0,etot0,ptot0,ztot0,stot0,
-     *            ang,etot,ptot,ztot,stot,rmsdpdt ,
-     *            rmsv,gtot(llmm1)
+      COMMON/ener/ang0,etot0,ptot0,ztot0,stot0,                         &
+     &            ang,etot,ptot,ztot,stot,rmsdpdt ,                     &
+     &            rmsv,gtot(llmm1)
 
-      REAL ang0,etot0,ptot0,ztot0,stot0,
-     s     ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot
+      REAL ang0,etot0,ptot0,ztot0,stot0,                                &
+     &     ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot
 
-c-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
Index: /LMDZ5/trunk/libf/dyn3dpar/fluxstokenc_p.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/fluxstokenc_p.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/fluxstokenc_p.F	(revision 1454)
@@ -57,5 +57,5 @@
 
 c AC initialisations
-cym      pbarug(:,:)   = 0.
+      pbarug(:,:)   = 0.
 cym      pbarvg(:,:,:) = 0.
 cym      wg(:,:)       = 0.
Index: /LMDZ5/trunk/libf/dyn3dpar/friction_p.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/friction_p.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/friction_p.F	(revision 1454)
@@ -6,17 +6,23 @@
       USE parallel
       USE control_mod
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      USE ioipsl_getincom
+#endif
       IMPLICIT NONE
 
-c=======================================================================
-c
-c
-c   Objet:
-c   ------
-c
-c  ***********
-c    Friction
-c  ***********
-c
-c=======================================================================
+!=======================================================================
+!
+!   Friction for the Newtonian case:
+!   --------------------------------
+!    2 possibilities (depending on flag 'friction_type'
+!     friction_type=0 : A friction that is only applied to the lowermost
+!                       atmospheric layer
+!     friction_type=1 : Friction applied on all atmospheric layer (but
+!       (default)       with stronger magnitude near the surface; see
+!                       iniacademic.F)
+!=======================================================================
 
 #include "dimensions.h"
@@ -24,16 +30,40 @@
 #include "comgeom2.h"
 #include "comconst.h"
-
-      REAL pdt
+#include "iniprint.h"
+#include "academic.h"
+
+! arguments:
+      REAL,INTENT(out) :: ucov( iip1,jjp1,llm )
+      REAL,INTENT(out) :: vcov( iip1,jjm,llm )
+      REAL,INTENT(in) :: pdt ! time step
+
+! local variables:
       REAL modv(iip1,jjp1),zco,zsi
       REAL vpn,vps,upoln,upols,vpols,vpoln
       REAL u2(iip1,jjp1),v2(iip1,jjm)
-      REAL ucov( iip1,jjp1,llm ),vcov( iip1,jjm,llm )
-      INTEGER  i,j
-      REAL cfric
-      parameter (cfric=1.e-5)
+      INTEGER  i,j,l
+      REAL,PARAMETER :: cfric=1.e-5
+      LOGICAL,SAVE :: firstcall=.true.
+      INTEGER,SAVE :: friction_type=1
+      CHARACTER(len=20) :: modname="friction_p"
+      CHARACTER(len=80) :: abort_message
+!$OMP THREADPRIVATE(firstcall,friction_type)
       integer :: jjb,jje
 
-
+!$OMP SINGLE
+      IF (firstcall) THEN
+        ! set friction type
+        call getin("friction_type",friction_type)
+        if ((friction_type.lt.0).or.(friction_type.gt.1)) then
+          abort_message="wrong friction type"
+          write(lunout,*)'Friction: wrong friction type',friction_type
+          call abort_gcm(modname,abort_message,42)
+        endif
+        firstcall=.false.
+      ENDIF
+!$OMP END SINGLE COPYPRIVATE(friction_type,firstcall)
+
+      if (friction_type.eq.0) then ! friction on first layer only
+!$OMP SINGLE
 c   calcul des composantes au carre du vent naturel
       jjb=jj_begin
@@ -138,4 +168,33 @@
          vcov(iip1,j,1)=vcov(1,j,1)
       enddo
+!$OMP END SINGLE
+      endif ! of if (friction_type.eq.0)
+
+      if (friction_type.eq.1) then
+       ! for ucov() 
+        jjb=jj_begin
+        jje=jj_end
+        if (pole_nord) jjb=jj_begin+1
+        if (pole_sud) jje=jj_end-1
+
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,llm
+          ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)*
+     &                                  (1.-pdt*kfrict(l))
+        enddo
+!$OMP END DO NOWAIT
+        
+       ! for vcoc()
+        jjb=jj_begin
+        jje=jj_end
+        if (pole_sud) jje=jj_end-1
+        
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,llm
+          vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)*
+     &                                  (1.-pdt*kfrict(l))
+        enddo
+!$OMP END DO
+      endif ! of if (friction_type.eq.1)
 
       RETURN
Index: /LMDZ5/trunk/libf/dyn3dpar/gcm.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/gcm.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/gcm.F	(revision 1454)
@@ -276,18 +276,10 @@
         endif
 
-        if (planet_type.eq."earth") then
-#ifdef CPP_EARTH
+!        if (planet_type.eq."earth") then
 ! Load an Earth-format start file
          CALL dynetat0("start.nc",vcov,ucov,
      &              teta,q,masse,ps,phis, time_0)
-#else
-        ! SW model also has Earth-format start files
-        ! (but can be used without the CPP_EARTH directive)
-          if (iflag_phys.eq.0) then
-            CALL dynetat0("start.nc",vcov,ucov,
-     &              teta,q,masse,ps,phis, time_0)
-          endif
-#endif
-        endif ! of if (planet_type.eq."earth")
+!        endif ! of if (planet_type.eq."earth")
+
 c       write(73,*) 'ucov',ucov
 c       write(74,*) 'vcov',vcov
@@ -494,7 +486,8 @@
 #endif
 
-      if (planet_type.eq."earth") then
+!      if (planet_type.eq."earth") then
+! Write an Earth-format restart file
         CALL dynredem0_p("restart.nc", day_end, phis)
-      endif
+!      endif
 
       ecripar = .TRUE.
Index: /LMDZ5/trunk/libf/dyn3dpar/grid_noro.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/grid_noro.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/grid_noro.F	(revision 1454)
@@ -458,10 +458,7 @@
 C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
 
-      PARAMETER (ISMo=300,JSMo=200)
-      REAL X(IMAR,JMAR),XF(ISMo,JSMo)
+      REAL X(IMAR,JMAR),XF(IMAR,JMAR)
       real WEIGHTpb(-1:1,-1:1)
 
-      if(imar.gt.ismo) stop'surdimensionner ismo dans mva9 (grid_noro)'
-      if(jmar.gt.jsmo) stop'surdimensionner jsmo dans mva9 (grid_noro)'
       
       SUM=0.
Index: /LMDZ5/trunk/libf/dyn3dpar/infotrac.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/infotrac.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/infotrac.F90	(revision 1454)
@@ -61,10 +61,10 @@
     CHARACTER(len=1), DIMENSION(3)  :: txts
     CHARACTER(len=2), DIMENSION(9)  :: txtp
-    CHARACTER(len=13)               :: str1,str2
+    CHARACTER(len=23)               :: str1,str2
   
     INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
     INTEGER :: iq, new_iq, iiq, jq, ierr
-    INTEGER, EXTERNAL :: lnblnk
- 
+
+    character(len=*),parameter :: modname="infotrac_init"
 !-----------------------------------------------------------------------
 ! Initialization :
@@ -100,13 +100,23 @@
        OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
        IF(ierr.EQ.0) THEN
-          WRITE(lunout,*) 'Open traceur.def : ok'
+          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
           READ(90,*) nqtrue
        ELSE 
-          WRITE(lunout,*) 'Problem in opening traceur.def'
-          WRITE(lunout,*) 'ATTENTION using defaut values'
-          nqtrue=4 ! Defaut value
-       END IF
-       ! Attention! Only for planet_type=='earth'
-       nbtr=nqtrue-2
+          WRITE(lunout,*) trim(modname),': Problem in opening traceur.def'
+          WRITE(lunout,*) trim(modname),': WARNING using defaut values'
+          if (planet_type=='earth') then
+            nqtrue=4 ! Default value for Earth
+          else
+            nqtrue=1 ! Default value for other planets
+          endif
+       END IF
+       if ( planet_type=='earth') then
+         ! For Earth, water vapour & liquid tracers are not in the physics
+         nbtr=nqtrue-2
+       else
+         ! Other planets (for now); we have the same number of tracers
+         ! in the dynamics than in the physics
+         nbtr=nqtrue
+       endif
     ELSE
        ! nbtr has been read from INCA by init_cont_lmdz() in gcm.F 
@@ -114,6 +124,6 @@
     END IF
 
-    IF (nqtrue < 2) THEN
-       WRITE(lunout,*) 'nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum'
+    IF ((planet_type=="earth").and.(nqtrue < 2)) THEN
+       WRITE(lunout,*) trim(modname),': nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum'
        CALL abort_gcm('infotrac_init','Not enough tracers',1)
     END IF
@@ -156,8 +166,10 @@
           ! Continue to read tracer.def
           DO iq=1,nqtrue
-             READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
+             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
           END DO
           CLOSE(90)  
-       ELSE ! Without tracer.def
+       ELSE ! Without tracer.def, set default values 
+         if (planet_type=="earth") then
+          ! for Earth, default is to have 4 tracers
           hadv(1) = 14
           vadv(1) = 14
@@ -172,8 +184,13 @@
           vadv(4) = 10
           tnom_0(4) = 'PB'
+         else ! default for other planets
+          hadv(1) = 10
+          vadv(1) = 10
+          tnom_0(1) = 'dummy'
+         endif ! of if (planet_type=="earth")
        END IF
        
-       WRITE(lunout,*) 'Valeur de traceur.def :'
-       WRITE(lunout,*) 'nombre de traceurs ',nqtrue
+       WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
+       WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue
        DO iq=1,nqtrue
           WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq)
@@ -217,5 +234,5 @@
           new_iq=new_iq+10 ! 9 tracers added
        ELSE
-          WRITE(lunout,*) 'This choice of advection schema is not available'
+          WRITE(lunout,*) trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
           CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1)
        END IF
@@ -227,8 +244,8 @@
        nqtot = new_iq
 
-       WRITE(lunout,*) 'The choice of advection schema for one or more tracers'
+       WRITE(lunout,*) trim(modname),': The choice of advection schema for one or more tracers'
        WRITE(lunout,*) 'makes it necessary to add tracers'
-       WRITE(lunout,*) nqtrue,' is the number of true tracers'
-       WRITE(lunout,*) nqtot, ' is the total number of tracers needed'
+       WRITE(lunout,*) trim(modname)//': ',nqtrue,' is the number of true tracers'
+       WRITE(lunout,*) trim(modname)//': ',nqtot, ' is the total number of tracers needed'
 
     ELSE
@@ -258,5 +275,6 @@
           iadv(new_iq)=11
        ELSE
-          WRITE(lunout,*)'This choice of advection schema is not available'
+          WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
+
           CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1)
        END IF
@@ -265,7 +283,7 @@
        tname(new_iq)= tnom_0(iq)
        IF (iadv(new_iq)==0) THEN
-          ttext(new_iq)=str1(1:lnblnk(str1))
+          ttext(new_iq)=trim(str1)
        ELSE
-          ttext(new_iq)=str1(1:lnblnk(str1))//descrq(iadv(new_iq))
+          ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
        END IF
 
@@ -276,6 +294,6 @@
              new_iq=new_iq+1
              iadv(new_iq)=-20
-             ttext(new_iq)=str2(1:lnblnk(str2))//txts(jq)
-             tname(new_iq)=str1(1:lnblnk(str1))//txts(jq)
+             ttext(new_iq)=trim(str2)//txts(jq)
+             tname(new_iq)=trim(str1)//txts(jq)
           END DO
        ELSE IF (iadv(new_iq)==30) THEN
@@ -283,6 +301,6 @@
              new_iq=new_iq+1
              iadv(new_iq)=-30
-             ttext(new_iq)=str2(1:lnblnk(str2))//txtp(jq)
-             tname(new_iq)=str1(1:lnblnk(str1))//txtp(jq)
+             ttext(new_iq)=trim(str2)//txtp(jq)
+             tname(new_iq)=trim(str1)//txtp(jq)
           END DO
        END IF
@@ -303,8 +321,9 @@
 
 
-    WRITE(lunout,*) 'Information stored in infotrac :'
-    WRITE(lunout,*) 'iadv  niadv tname  ttext :'
+    WRITE(lunout,*) trim(modname),': Information stored in infotrac :'
+    WRITE(lunout,*) trim(modname),': iadv  niadv tname  ttext :'
     DO iq=1,nqtot
-       WRITE(lunout,*) iadv(iq),niadv(iq), tname(iq), ttext(iq)
+       WRITE(lunout,*) iadv(iq),niadv(iq),&
+       ' ',trim(tname(iq)),' ',trim(ttext(iq))
     END DO
 
@@ -315,8 +334,8 @@
     DO iq=1,nqtot
        IF (iadv(iq)/=10 .AND. iadv(iq)/=14 .AND. iadv(iq)/=0) THEN
-          WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
+          WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
           CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1)
        ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN
-          WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
+          WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
           CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1)
        END IF
@@ -329,6 +348,4 @@
     DEALLOCATE(tracnam)
 
-999 FORMAT (i2,1x,i2,1x,a15)
-
   END SUBROUTINE infotrac_init
 
Index: /LMDZ5/trunk/libf/dyn3dpar/iniacademic.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/iniacademic.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/iniacademic.F	(revision 1454)
@@ -8,5 +8,12 @@
       USE filtreg_mod
       USE infotrac, ONLY : nqtot
-      USE control_mod
+      USE control_mod, ONLY: day_step,planet_type
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      USE ioipsl_getincom
+#endif
+      USE Write_Field
  
 
@@ -70,8 +77,14 @@
       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL ddsin,tetarappelj,tetarappell,zsig
+      REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
       real tetajl(jjp1,llm)
       INTEGER i,j,l,lsup,ij
 
+      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
+      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
+      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
+      LOGICAL ok_pv                ! Polar Vortex
+      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
+      
       real zz,ran1
       integer idum
@@ -82,23 +95,17 @@
 ! 1. Initializations for Earth-like case
 ! --------------------------------------
-      if (planet_type=="earth") then
-c
+c
+        ! initialize planet radius, rotation rate,...
+        call conf_planete
+
         time_0=0.
-        day_ref=0
+        day_ref=1
         annee_ref=0
 
         im         = iim
         jm         = jjm
-        day_ini    = 0
-        omeg       = 4.*asin(1.)/86400.
-        rad    = 6371229.
-        g      = 9.8
-        daysec = 86400.
+        day_ini    = 1
         dtvr    = daysec/REAL(day_step)
         zdtvr=dtvr
-        kappa  = 0.2857143
-        cpp    = 1004.70885
-        preff     = 101325.
-        pa        =  50000.
         etot0      = 0.
         ptot0      = 0.
@@ -120,7 +127,5 @@
           if (.not.read_start) then
             phis(:)=0.
-            q(:,:,1)=1.e-10
-            q(:,:,2)=1.e-15
-            q(:,:,3:nqtot)=0.
+            q(:,:,:)=0
             CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
           endif
@@ -129,30 +134,78 @@
         if (iflag_phys.eq.2) then
           ! initializations for the academic case
-          ps(:)=1.e5
-          phis(:)=0.
-c---------------------------------------------------------------------
-
-          taurappel=10.*daysec
-
-c---------------------------------------------------------------------
-c   Calcul de la temperature potentielle :
-c   --------------------------------------
-
+          
+!         if (planet_type=="earth") then
+
+          ! 1. local parameters
+          ! by convention, winter is in the southern hemisphere
+          ! Geostrophic wind or no wind?
+          ok_geost=.TRUE.
+          CALL getin('ok_geost',ok_geost)
+          ! Constants for Newtonian relaxation and friction
+          k_f=1.                !friction 
+          CALL getin('k_j',k_f)
+          k_f=1./(daysec*k_f)
+          k_c_s=4.  !cooling surface
+          CALL getin('k_c_s',k_c_s)
+          k_c_s=1./(daysec*k_c_s)
+          k_c_a=40. !cooling free atm
+          CALL getin('k_c_a',k_c_a)
+          k_c_a=1./(daysec*k_c_a)
+          ! Constants for Teta equilibrium profile
+          teta0=315.     ! mean Teta (S.H. 315K)
+          CALL getin('teta0',teta0)
+          ttp=200.       ! Tropopause temperature (S.H. 200K)
+          CALL getin('ttp',ttp)
+          eps=0.         ! Deviation to N-S symmetry(~0-20K)
+          CALL getin('eps',eps)
+          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
+          CALL getin('delt_y',delt_y)
+          delt_z=10.     ! Vertical Gradient (S.H. 10K)
+          CALL getin('delt_z',delt_z)
+          ! Polar vortex
+          ok_pv=.false.
+          CALL getin('ok_pv',ok_pv)
+          phi_pv=-50.            ! Latitude of edge of vortex
+          CALL getin('phi_pv',phi_pv)
+          phi_pv=phi_pv*pi/180.
+          dphi_pv=5.             ! Width of the edge
+          CALL getin('dphi_pv',dphi_pv)
+          dphi_pv=dphi_pv*pi/180.
+          gam_pv=4.              ! -dT/dz vortex (in K/km)
+          CALL getin('gam_pv',gam_pv)
+          
+          ! 2. Initialize fields towards which to relax
+          ! Friction
+          knewt_g=k_c_a
           DO l=1,llm
-            zsig=ap(l)/preff+bp(l)
-            if (zsig.gt.0.3) then
-             lsup=l
-             tetarappell=1./8.*(-log(zsig)-.5)
-             DO j=1,jjp1
-             ddsin=sin(rlatu(j))-sin(pi/20.)
-             tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell)
-             ENDDO
-            else
-c   Choix isotherme au-dessus de 300 mbar
-             do j=1,jjp1
-               tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
-             enddo
-            endif ! of if (zsig.gt.0.3)
+           zsig=presnivs(l)/preff
+           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
+           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
+          ENDDO
+          DO j=1,jjp1
+            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
+          ENDDO
+          
+          ! Potential temperature 
+          DO l=1,llm
+           zsig=presnivs(l)/preff
+           tetastrat=ttp*zsig**(-kappa)
+           tetapv=tetastrat
+           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
+               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
+           ENDIF
+           DO j=1,jjp1
+             ! Troposphere
+             ddsin=sin(rlatu(j))
+             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
+     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
+             ! Profil stratospherique isotherme (+vortex)
+             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
+             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
+             tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
+           ENDDO
           ENDDO ! of DO l=1,llm
+ 
+!          CALL writefield('theta_eq',tetajl)
 
           do l=1,llm
@@ -165,56 +218,68 @@
           enddo
 
-c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
-
-          CALL pression ( ip1jmp1, ap, bp, ps, p       )
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-          CALL massdair(p,masse)
-
-c  intialisation du vent et de la temperature
-          teta(:,:)=tetarappel(:,:)
-          CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-          call ugeostr(phi,ucov)
-          vcov=0.
-          q(:,:,1   )=1.e-10
-          q(:,:,2   )=1.e-15
-          q(:,:,3:nqtot)=0.
-
-
-c   perturbation aleatoire sur la temperature
-          idum  = -1
-          zz = ran1(idum)
-          idum  = 0
-          do l=1,llm
-            do ij=iip2,ip1jm
-              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+
+!         else
+!          write(lunout,*)"iniacademic: planet types other than earth",
+!     &                   " not implemented (yet)."
+!          stop
+!         endif ! of if (planet_type=="earth")
+
+          ! 3. Initialize fields (if necessary)
+          IF (.NOT. read_start) THEN
+            ! surface pressure
+            ps(:)=preff
+            ! ground geopotential
+            phis(:)=0.
+            
+            CALL pression ( ip1jmp1, ap, bp, ps, p       )
+            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+            CALL massdair(p,masse)
+
+            ! bulk initialization of temperature
+            teta(:,:)=tetarappel(:,:)
+            
+            ! geopotential
+            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
+            
+            ! winds
+            if (ok_geost) then
+              call ugeostr(phi,ucov)
+            else
+              ucov(:,:)=0.
+            endif
+            vcov(:,:)=0.
+            
+            ! bulk initialization of tracers
+            if (planet_type=="earth") then
+              ! Earth: first two tracers will be water
+              do i=1,nqtot
+                if (i.eq.1) q(:,:,i)=1.e-10
+                if (i.eq.2) q(:,:,i)=1.e-15
+                if (i.gt.2) q(:,:,i)=0.
+              enddo
+            else
+              q(:,:,:)=0
+            endif ! of if (planet_type=="earth")
+
+            ! add random perturbation to temperature
+            idum  = -1
+            zz = ran1(idum)
+            idum  = 0
+            do l=1,llm
+              do ij=iip2,ip1jm
+                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+              enddo
             enddo
-          enddo
-
-          do l=1,llm
-            do ij=1,ip1jmp1,iip1
-              teta(ij+iim,l)=teta(ij,l)
+
+            ! maintain periodicity in longitude
+            do l=1,llm
+              do ij=1,ip1jmp1,iip1
+                teta(ij+iim,l)=teta(ij,l)
+              enddo
             enddo
-          enddo
-
-
-
-c     PRINT *,' Appel test_period avec tetarappel '
-c     CALL  test_period ( ucov,vcov,tetarappel,q,p,phis )
-c     PRINT *,' Appel test_period avec teta '
-c     CALL  test_period ( ucov,vcov,teta,q,p,phis )
-
-c   initialisation d'un traceur sur une colonne
-          j=jjp1*3/4
-          i=iip1/2
-          ij=(j-1)*iip1+i
-          q(ij,:,3)=1.
+
+          ENDIF ! of IF (.NOT. read_start)
         endif ! of if (iflag_phys.eq.2)
         
-      else
-        write(lunout,*)"iniacademic: planet types other than earth",
-     &                 " not implemented (yet)."
-        stop
-      endif ! of if (planet_type=="earth")
-      return
       END
 c-----------------------------------------------------------------------
Index: /LMDZ5/trunk/libf/dyn3dpar/initfluxsto_p.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/initfluxsto_p.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/initfluxsto_p.F	(revision 1454)
@@ -203,9 +203,10 @@
      .              llm, nivsigs, zvertiid)
 c pour le fichier def
-      nivd(1) = 1
-      call histvert(filedid, 'sig_s', 'Niveaux sigma',
-     .  'sigma_level',
-     .              1, nivd, dvertiid)
-
+      if (mpi_rank==0) then
+         nivd(1) = 1
+         call histvert(filedid, 'sig_s', 'Niveaux sigma',
+     .        'sigma_level',
+     .        1, nivd, dvertiid)
+      endif
 C
 C  Appels a histdef pour la definition des variables a sauvegarder
@@ -282,9 +283,9 @@
       call histend(fileid)
       call histend(filevid)
-      call histend(filedid)
+      if (mpi_rank==0) call histend(filedid)
       if (ok_sync) then
         call histsync(fileid)
         call histsync(filevid)
-        call histsync(filedid)
+        if (mpi_rank==0) call histsync(filedid)
       endif
 	
Index: /LMDZ5/trunk/libf/dyn3dpar/integrd_p.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/integrd_p.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/integrd_p.F	(revision 1454)
@@ -6,5 +6,5 @@
      $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis,finvmaold)
       USE parallel
-      USE control_mod
+      USE control_mod, only : planet_type
       IMPLICIT NONE
 
@@ -279,5 +279,4 @@
 
           CALL qminimum_p( q, nq, deltap )
-	 endif ! of if (planet_type.eq."earth")
 c
 c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
@@ -337,4 +336,7 @@
       ENDDO
 c$OMP END DO NOWAIT
+
+      endif ! of if (planet_type.eq."earth")
+
 c
 c
Index: /LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 1454)
@@ -234,5 +234,5 @@
 
 c$OMP MASTER
-      dq=0.
+      dq(:,:,:)=0.
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
@@ -596,15 +596,10 @@
      .        flxw,pk, iapptrac)
 
-       IF (offline) THEN
-Cmaf stokage du flux de masse pour traceurs OFF-LINE
-
-#ifdef CPP_IOIPSL
-           CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
-     .   dtvr, itau)
-#endif
-
-
-         ENDIF ! of IF (offline)
-c
+C        Stokage du flux de masse pour traceurs OFF-LINE
+         IF (offline .AND. .NOT. adjust) THEN
+            CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
+     .           dtvr, itau)
+         ENDIF
+
       ENDIF ! of IF( forward. OR . leapf )
 cc$OMP END PARALLEL
@@ -968,15 +963,14 @@
 
       IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
-c   Calcul academique de la physique = Rappel Newtonien + fritcion 
-c   --------------------------------------------------------------
-cym       teta(:,:)=teta(:,:)
-cym     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
+!   Academic case : Simple friction and Newtonan relaxation 
+!   -------------------------------------------------------
        ijb=ij_begin
        ije=ij_end
 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
        do l=1,llm
-       teta(ijb:ije,l)=teta(ijb:ije,l)
-     &  -iphysiq*dtvr*(teta(ijb:ije,l)-tetarappel(ijb:ije,l))/taurappel
-       enddo
+        teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
+     &         (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
+     &                  (knewt_g+knewt_t(l)*clat4(ijb:ije))
+       enddo ! of do l=1,llm
 !$OMP END DO
 
@@ -987,8 +981,37 @@
        call WaitRequest(Request_Physic)     
 c$OMP BARRIER
-!$OMP MASTER
-       call friction_p(ucov,vcov,iphysiq*dtvr)
-!$OMP END MASTER
+       call friction_p(ucov,vcov,dtvr)
 !$OMP BARRIER
+
+        ! Sponge layer (if any)
+        IF (ok_strato) THEN
+          ! set dufi,dvfi,... to zero
+          ijb=ij_begin
+          ije=ij_end
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+          do l=1,llm
+            dufi(ijb:ije,l)=0
+            dtetafi(ijb:ije,l)=0
+            dqfi(ijb:ije,l,1:nqtot)=0
+          enddo
+!$OMP END DO
+!$OMP SINGLE
+          dpfi(ijb:ije)=0
+!$OMP END SINGLE
+          ijb=ij_begin
+          ije=ij_end
+          if (pole_sud) ije=ije-iip1
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+          do l=1,llm
+            dvfi(ijb:ije,l)=0
+          enddo
+!$OMP END DO
+
+          CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
+          CALL addfi_p( dtvr, leapf, forward   ,
+     $                  ucov, vcov, teta , q   ,ps ,
+     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
+!$OMP BARRIER
+        ENDIF ! of IF (ok_strato) 
       ENDIF ! of IF(iflag_phys.EQ.2)
 
@@ -1465,9 +1488,9 @@
 c$OMP MASTER
 
-              if (planet_type.eq."earth") then
+!              if (planet_type.eq."earth") then
 ! Write an Earth-format restart file
                 CALL dynredem1_p("restart.nc",0.0,
      &                           vcov,ucov,teta,q,masse,ps)
-              endif ! of if (planet_type.eq."earth")
+!              endif ! of if (planet_type.eq."earth")
 
 !              CLOSE(99)
@@ -1658,10 +1681,10 @@
 
               IF(itau.EQ.itaufin) THEN
-                if (planet_type.eq."earth") then
+!                if (planet_type.eq."earth") then
 c$OMP MASTER
                    CALL dynredem1_p("restart.nc",0.0,
      .                               vcov,ucov,teta,q,masse,ps)
 c$OMP END MASTER
-                endif ! of if (planet_type.eq."earth")
+!                endif ! of if (planet_type.eq."earth")
               ENDIF ! of IF(itau.EQ.itaufin)
 
Index: /LMDZ5/trunk/libf/dyn3dpar/limit_netcdf.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/limit_netcdf.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/dyn3dpar/limit_netcdf.F90	(revision 1454)
@@ -97,5 +97,5 @@
   kappa = 0.2857143
   cpp   = 1004.70885
-  dtvr  = daysec/FLOAT(day_step)
+  dtvr  = daysec/REAL(day_step)
   CALL inigeom
 
@@ -265,6 +265,4 @@
 
   DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
-#endif
-! of #ifdef CPP_EARTH
 
 
@@ -592,7 +590,7 @@
 
 !--- Mid-months times
-  mid_months(1)=0.5*FLOAT(mnth(1))
+  mid_months(1)=0.5*REAL(mnth(1))
   DO k=2,nm
-    mid_months(k)=mid_months(k-1)+0.5*FLOAT(mnth(k-1)+mnth(k))
+    mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
   END DO
 
@@ -626,4 +624,6 @@
 !-------------------------------------------------------------------------------
 
+#endif
+! of #ifdef CPP_EARTH
 
 END SUBROUTINE limit_netcdf
Index: /LMDZ5/trunk/libf/phylmd/calcul_divers.h
===================================================================
--- /LMDZ5/trunk/libf/phylmd/calcul_divers.h	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/calcul_divers.h	(revision 1454)
@@ -2,28 +2,23 @@
 c $Header$
 c
-c
-c initialisations diverses au "debut" du mois
-c
-      IF(MOD(itap,INT(ecrit_mth/dtime)).EQ.1) THEN
+
+c     Initialisations diverses au "debut" du mois
+      IF(debut) THEN
+         nday_rain(:)=0.
+
+c        surface terre
+         paire_ter(:)=0.
          DO i=1, klon
-          nday_rain(i)=0.
-         ENDDO !i
-c
-c surface terre
-       DO i=1, klon
-         IF(pctsrf(i,is_ter).GT.0.) THEN
-            paire_ter(i)=airephy(i)*pctsrf(i,is_ter)
-         ENDIF 
-       ENDDO
-c
-      ENDIF !MOD(itap,INT(ecrit_mth)).EQ.1
-c
-      IF(MOD(itap,INT(ecrit_day/dtime)).EQ.0) THEN
-c
-cIM calcul total_rain, nday_rain
-c
-       DO i = 1, klon
-        total_rain(i)=rain_fall(i)+snow_fall(i)  
-        IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.
-       ENDDO
-      ENDIF !itap.EQ.ecrit_mth
+            IF(pctsrf(i,is_ter).GT.0.) THEN
+               paire_ter(i)=airephy(i)*pctsrf(i,is_ter)
+            ENDIF 
+         ENDDO
+      ENDIF
+
+cIM   Calcul une fois par jour : total_rain, nday_rain
+      IF(MOD(itap,INT(un_jour/dtime)).EQ.0) THEN
+         DO i = 1, klon
+            total_rain(i)=rain_fall(i)+snow_fall(i)  
+            IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.
+         ENDDO
+      ENDIF
Index: /LMDZ5/trunk/libf/phylmd/carbon_cycle_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/carbon_cycle_mod.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/carbon_cycle_mod.F90	(revision 1454)
@@ -1,4 +1,9 @@
 MODULE carbon_cycle_mod
-
+! Controle module for the carbon CO2 tracers :
+!   - Identification
+!   - Get concentrations comming from coupled model or read from file to tracers
+!   - Calculate new RCO2 for radiation scheme
+!   - Calculate new carbon flux for sending to coupled models (PISCES and ORCHIDEE)
+!
 ! Author : Josefine GHATTAS, Patricia CADULE
 
@@ -13,5 +18,10 @@
   LOGICAL, PUBLIC :: carbon_cycle_cpl       ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/OCEAN(PISCES) 
 !$OMP THREADPRIVATE(carbon_cycle_cpl)
+
   LOGICAL :: carbon_cycle_emis_comp=.FALSE. ! Calculation of emission compatible
+!$OMP THREADPRIVATE(carbon_cycle_emis_comp)
+
+  LOGICAL :: RCO2_inter  ! RCO2 interactive : if true calculate new value RCO2 for the radiation scheme
+!$OMP THREADPRIVATE(RCO2_inter)
 
 ! Scalare values when no transport, from physiq.def
@@ -21,16 +31,14 @@
 !$OMP THREADPRIVATE(emis_land_s)
 
-  INTEGER :: ntr_co2                ! Number of tracers concerning the carbon cycle
-  INTEGER :: id_fco2_tot            ! Tracer index
-  INTEGER :: id_fco2_ocn            !  - " -
-  INTEGER :: id_fco2_land           !  - " -
-  INTEGER :: id_fco2_land_use       !  - " -
-  INTEGER :: id_fco2_fos_fuel       !  - " -
-!$OMP THREADPRIVATE(ntr_co2, id_fco2_tot, id_fco2_ocn, id_fco2_land, id_fco2_land_use, id_fco2_fos_fuel)  
-
-  REAL, DIMENSION(:), ALLOCATABLE :: fos_fuel        ! CO2 fossil fuel emission from file [gC/m2/d]
-!$OMP THREADPRIVATE(fos_fuel)
-  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day ! flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]
+  REAL :: airetot     ! Total area of the earth surface
+!$OMP THREADPRIVATE(airetot)
+
+  INTEGER :: ntr_co2  ! Number of tracers concerning the carbon cycle
+!$OMP THREADPRIVATE(ntr_co2)
+
+! fco2_ocn_day : flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]. Allocation and initalization done in cpl_mod
+  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day 
 !$OMP THREADPRIVATE(fco2_ocn_day)
+
   REAL, DIMENSION(:), ALLOCATABLE :: fco2_land_day   ! flux CO2 from land for 1 day (cumulated)  [gC/m2/d]
 !$OMP THREADPRIVATE(fco2_land_day)
@@ -38,5 +46,8 @@
 !$OMP THREADPRIVATE(fco2_lu_day)
 
-! Following 2 fields will be initialized in surf_land_orchidee at each time step
+  REAL, DIMENSION(:,:), ALLOCATABLE :: dtr_add       ! Tracer concentration to be injected 
+!$OMP THREADPRIVATE(dtr_add)
+
+! Following 2 fields will be allocated and initialized in surf_land_orchidee
   REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_inst  ! flux CO2 from land at one time step
 !$OMP THREADPRIVATE(fco2_land_inst)
@@ -45,14 +56,37 @@
 
 ! Calculated co2 field to be send to the ocean via the coupler and to ORCHIDEE 
-  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send 
+  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send ! Field allocated in phyetat0
 !$OMP THREADPRIVATE(co2_send)
+
+
+  TYPE, PUBLIC ::   co2_trac_type
+     CHARACTER(len = 8) :: name       ! Tracer name in tracer.def
+     INTEGER            :: id         ! Index in total tracer list, tr_seri
+     CHARACTER(len=30)  :: file       ! File name
+     LOGICAL            :: cpl        ! True if this tracers is coupled from ORCHIDEE or PISCES. 
+                                      ! False if read from file.
+     INTEGER            :: updatefreq ! Frequence to inject in second
+     INTEGER            :: readstep   ! Actual time step to read in file
+     LOGICAL            :: updatenow  ! True if this tracer should be updated this time step
+  END TYPE co2_trac_type
+  INTEGER,PARAMETER :: maxco2trac=5  ! Maximum number of different CO2 fluxes
+  TYPE(co2_trac_type), DIMENSION(maxco2trac) :: co2trac
 
 CONTAINS
   
-  SUBROUTINE carbon_cycle_init(tr_seri, aerosol, radio)
+  SUBROUTINE carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
+! This subroutine is called from traclmdz_init, only at first timestep.
+! - Read controle parameters from .def input file
+! - Search for carbon tracers and set default values
+! - Allocate variables
+! - Test for compatibility
+
     USE dimphy
+    USE comgeomphy
+    USE mod_phys_lmdz_transfert_para
     USE infotrac
     USE IOIPSL
     USE surface_data, ONLY : ok_veget, type_ocean
+    USE phys_cal_mod, ONLY : mth_len
 
     IMPLICIT NONE
@@ -62,4 +96,5 @@
 ! Input argument
     REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! Concentration Traceur [U/KgA]  
+    REAL,INTENT(IN)                           :: pdtphys ! length of time step in physiq (sec)
 
 ! InOutput arguments
@@ -68,78 +103,13 @@
 
 ! Local variables
-    INTEGER               :: ierr, it, iiq
-    REAL, DIMENSION(klon) :: tr_seri_sum
-
-
-! 0) Test for compatibility
-    IF (carbon_cycle_cpl .AND. type_ocean/='couple') &
-         CALL abort_gcm('carbon_cycle_init', 'Coupling with ocean model is needed for carbon_cycle_cpl',1)
-    IF (carbon_cycle_cpl .AND..NOT. ok_veget) &
-         CALL abort_gcm('carbon_cycle_init', 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl',1)
-
-
-! 1) Check if transport of one tracer flux CO2 or 4 separated tracers
-    IF (carbon_cycle_tr) THEN
-       id_fco2_tot=0
-       id_fco2_ocn=0
-       id_fco2_land=0
-       id_fco2_land_use=0
-       id_fco2_fos_fuel=0
-       
-       ! Search in tracer list
-       DO it=1,nbtr
-          iiq=niadv(it+2)
-          IF (tname(iiq) == "fCO2" ) THEN
-             id_fco2_tot=it
-          ELSE IF (tname(iiq) == "fCO2_ocn" ) THEN
-             id_fco2_ocn=it
-          ELSE IF (tname(iiq) == "fCO2_land" ) THEN
-             id_fco2_land=it
-          ELSE IF (tname(iiq) == "fCO2_land_use" ) THEN
-             id_fco2_land_use=it
-          ELSE IF (tname(iiq) == "fCO2_fos_fuel" ) THEN
-             id_fco2_fos_fuel=it
-          END IF
-       END DO
-
-       ! Count tracers found
-       IF (id_fco2_tot /= 0 .AND. &
-            id_fco2_ocn==0 .AND. id_fco2_land==0 .AND. id_fco2_land_use==0 .AND. id_fco2_fos_fuel==0) THEN
-          
-          ! transport  1 tracer flux CO2
-          ntr_co2 = 1
-          
-       ELSE IF (id_fco2_tot==0 .AND. &
-            id_fco2_ocn /=0 .AND. id_fco2_land/=0 .AND. id_fco2_land_use/=0 .AND. id_fco2_fos_fuel/=0) THEN
-          
-          ! transport 4 tracers seperatively
-          ntr_co2 = 4
-          
-       ELSE
-          CALL abort_gcm('carbon_cycle_init', 'error in coherence between traceur.def and gcm.def',1)
-       END IF
-
-       ! Definition of control varaiables for the tracers
-       DO it=1,nbtr
-          IF (it==id_fco2_tot .OR. it==id_fco2_ocn .OR. it==id_fco2_land .OR. &
-               it==id_fco2_land_use .OR. it==id_fco2_fos_fuel) THEN
-             aerosol(it) = .FALSE.
-             radio(it)   = .FALSE.
-          END IF
-       END DO
-
-    ELSE 
-       ! No transport of CO2
-       ntr_co2 = 0
-    END IF ! carbon_cycle_tr
-
-
-! 2) Allocate variable for CO2 fossil fuel emission
-    IF (carbon_cycle_tr) THEN
-       ! Allocate 2D variable
-       ALLOCATE(fos_fuel(klon), stat=ierr)
-       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 1',1)
-    ELSE
-       ! No transport : read value from .def
+    INTEGER               :: ierr, it, iiq, itc
+    INTEGER               :: teststop
+
+
+
+! 1) Read controle parameters from .def input file
+! ------------------------------------------------
+    ! Read fosil fuel value if no transport
+    IF (.NOT. carbon_cycle_tr) THEN
        fos_fuel_s = 0.
        CALL getin ('carbon_cycle_fos_fuel',fos_fuel_s)
@@ -148,68 +118,177 @@
 
 
-! 3) Allocate and initialize fluxes
-    IF (carbon_cycle_cpl) THEN
-       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
-       ALLOCATE(fco2_land_day(klon), stat=ierr)
-       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
-       ALLOCATE(fco2_lu_day(klon), stat=ierr)
-       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 4',1)
-
-       fco2_land_day(:) = 0.  ! JG : Doit prend valeur de restart
-       fco2_lu_day(:)   = 0.  ! JG : Doit prend valeur de restart
-
-       ! fco2_ocn_day is allocated in cpl_mod
-       ! fco2_land_inst and fco2_lu_inst are allocated in surf_land_orchidee
-       
-       ALLOCATE(co2_send(klon), stat=ierr)
-       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 7',1)
-       
-       ! Calculate using restart tracer values
-       IF (carbon_cycle_tr) THEN
-          IF (ntr_co2==1) THEN
-             co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
-          ELSE ! ntr_co2==4
-             ! Calculate the delta CO2 flux
-             tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
-                  tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
-             co2_send(:) = tr_seri_sum(:) + co2_ppm0
-          END IF
-       ELSE
-          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
-          co2_send(:) = co2_ppm
-       END IF
-
-
-    ELSE 
-       IF (carbon_cycle_tr) THEN
-          ! No coupling of CO2 fields : 
-          ! corresponding fields will instead be read from files
-          ALLOCATE(fco2_ocn_day(klon), stat=ierr)
-          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 8',1)
-          ALLOCATE(fco2_land_day(klon), stat=ierr)
-          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 9',1)
-          ALLOCATE(fco2_lu_day(klon), stat=ierr)
-          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 10',1)       
-       END IF
-    END IF
-
-! 4) Read parmeter for calculation of emission compatible
+    ! Read parmeter for calculation compatible emission
     IF (.NOT. carbon_cycle_tr) THEN
        carbon_cycle_emis_comp=.FALSE.
        CALL getin('carbon_cycle_emis_comp',carbon_cycle_emis_comp)
        WRITE(lunout,*) 'carbon_cycle_emis_comp = ',carbon_cycle_emis_comp
-    END IF
+       IF (carbon_cycle_emis_comp) THEN
+          CALL abort_gcm('carbon_cycle_init', 'carbon_cycle_emis_comp option not yet implemented!!',1)
+       END IF
+    END IF
+
+    ! Read parameter for interactive calculation of the CO2 value for the radiation scheme
+    RCO2_inter=.FALSE.
+    CALL getin('RCO2_inter',RCO2_inter)
+    WRITE(lunout,*) 'RCO2_inter = ', RCO2_inter
+    IF (RCO2_inter) THEN
+       WRITE(lunout,*) 'RCO2 will be recalculated once a day'
+       WRITE(lunout,*) 'RCO2 initial = ', RCO2
+    END IF
+
+
+! 2) Search for carbon tracers and set default values
+! ---------------------------------------------------
+    itc=0
+    DO it=1,nbtr
+       iiq=niadv(it+2)
+       
+       SELECT CASE(tname(iiq))
+       CASE("fCO2_ocn")
+          itc = itc + 1
+          co2trac(itc)%name='fCO2_ocn'
+          co2trac(itc)%id=it
+          co2trac(itc)%file='fl_co2_ocean.nc'
+          IF (carbon_cycle_cpl .AND. type_ocean=='couple') THEN 
+             co2trac(itc)%cpl=.TRUE.
+             co2trac(itc)%updatefreq = 86400 ! Once a day as the coupling with OASIS/PISCES
+          ELSE
+             co2trac(itc)%cpl=.FALSE.
+             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
+          END IF
+       CASE("fCO2_land")
+          itc = itc + 1
+          co2trac(itc)%name='fCO2_land'
+          co2trac(itc)%id=it
+          co2trac(itc)%file='fl_co2_land.nc'
+          IF (carbon_cycle_cpl .AND. ok_veget) THEN 
+             co2trac(itc)%cpl=.TRUE.
+             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
+          ELSE
+             co2trac(itc)%cpl=.FALSE.
+!             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
+             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
+          END IF
+       CASE("fCO2_land_use")
+          itc = itc + 1
+          co2trac(itc)%name='fCO2_land_use'
+          co2trac(itc)%id=it
+          co2trac(itc)%file='fl_co2_land_use.nc'
+          IF (carbon_cycle_cpl .AND. ok_veget) THEN 
+             co2trac(it)%cpl=.TRUE.
+             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
+          ELSE
+             co2trac(itc)%cpl=.FALSE.
+             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
+          END IF
+       CASE("fCO2_fos_fuel")
+          itc = itc + 1
+          co2trac(itc)%name='fCO2_fos_fuel'
+          co2trac(itc)%id=it
+          co2trac(itc)%file='fossil_fuel.nc'
+          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
+!         co2trac(itc)%updatefreq = 86400  ! 86400sec = 24H Cadule case
+          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
+       CASE("fCO2_bbg")
+          itc = itc + 1
+          co2trac(itc)%name='fCO2_bbg'
+          co2trac(itc)%id=it
+          co2trac(itc)%file='fl_co2_bbg.nc'
+          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
+          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
+       CASE("fCO2")
+          ! fCO2 : One tracer transporting the total CO2 flux
+          itc = itc + 1
+          co2trac(itc)%name='fCO2'
+          co2trac(itc)%id=it
+          co2trac(itc)%file='fl_co2.nc'
+          IF (carbon_cycle_cpl) THEN 
+             co2trac(itc)%cpl=.TRUE.
+          ELSE
+             co2trac(itc)%cpl=.FALSE.
+          END IF
+          co2trac(itc)%updatefreq = 86400
+          ! DOES THIS WORK ???? Problematic due to implementation of the coupled fluxes...
+          CALL abort_gcm('carbon_cycle_init','transport of total CO2 has to be implemented and tested',1)
+       END SELECT
+    END DO
+
+    ! Total number of carbon CO2 tracers
+    ntr_co2 = itc 
+    
+    ! Definition of control varaiables for the tracers
+    DO it=1,ntr_co2
+       aerosol(co2trac(it)%id) = .FALSE.
+       radio(co2trac(it)%id)   = .FALSE.
+    END DO
+    
+    ! Vector indicating which timestep to read for each tracer
+    ! Always start read in the beginning of the file
+    co2trac(:)%readstep = 0
+   
+
+! 3) Allocate variables
+! ---------------------
+    ! Allocate vector for storing fluxes to inject
+    ALLOCATE(dtr_add(klon,maxco2trac), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 11',1)       
+    
+    ! Allocate variables for cumulating fluxes from ORCHIDEE
+    IF (RCO2_inter) THEN
+       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
+          ALLOCATE(fco2_land_day(klon), stat=ierr)
+          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
+          fco2_land_day(1:klon) = 0.
+          
+          ALLOCATE(fco2_lu_day(klon), stat=ierr)
+          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
+          fco2_lu_day(1:klon)   = 0.
+       END IF
+    END IF
+
+
+! 4) Test for compatibility
+! -------------------------
+!    IF (carbon_cycle_cpl .AND. type_ocean/='couple') THEN
+!       WRITE(lunout,*) 'Coupling with ocean model is needed for carbon_cycle_cpl'
+!       CALL abort_gcm('carbon_cycle_init', 'coupled ocean is needed for carbon_cycle_cpl',1)
+!    END IF
+!
+!    IF (carbon_cycle_cpl .AND..NOT. ok_veget) THEN
+!       WRITE(lunout,*) 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl'
+!       CALL abort_gcm('carbon_cycle_init', 'ok_veget is needed for carbon_cycle_cpl',1)
+!    END IF
+
+    ! Compiler test : following should never happen
+    teststop=0
+    DO it=1,teststop
+       CALL abort_gcm('carbon_cycle_init', 'Entering loop from 1 to 0',1)
+    END DO
+
+    IF (ntr_co2==0) THEN
+       ! No carbon tracers found in tracer.def. It is not possible to do carbon cycle 
+       WRITE(lunout,*) 'No carbon tracers found in tracer.def. Not ok with carbon_cycle_tr and/or carbon_cycle_cp'
+       CALL abort_gcm('carbon_cycle_init', 'No carbon tracers found in tracer.def',1)
+    END IF
+    
+! 5) Calculate total area of the earth surface
+! --------------------------------------------
+    CALL reduce_sum(SUM(airephy),airetot)
+    CALL bcast(airetot)
 
   END SUBROUTINE carbon_cycle_init
 
+
+  SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
+! Subroutine for injection of co2 in the tracers
 !
-!
-!
-
-  SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
-    
+! - Find out if it is time to update
+! - Get tracer from coupled model or from file
+! - Calculate new RCO2 value for the radiation scheme
+! - Calculate CO2 flux to send to ocean and land models (PISCES and ORCHIDEE)
+
     USE infotrac
     USE dimphy
-    USE mod_phys_lmdz_transfert_para, ONLY : reduce_sum
+    USE mod_phys_lmdz_transfert_para
     USE phys_cal_mod, ONLY : mth_cur, mth_len
     USE phys_cal_mod, ONLY : day_cur
@@ -220,12 +299,16 @@
     INCLUDE "clesphys.h"
     INCLUDE "indicesol.h"
+    INCLUDE "iniprint.h"
+    INCLUDE "YOMCST.h"
 
 ! In/Output arguments
     INTEGER,INTENT(IN) :: nstep      ! time step in physiq
     REAL,INTENT(IN)    :: pdtphys    ! length of time step in physiq (sec)
-    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
-    REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri
+    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf            ! Surface fraction
+    REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri ! All tracers
+    REAL, DIMENSION(klon,nbtr), INTENT(INOUT)       :: source  ! Source for all tracers
 
 ! Local variables
+    INTEGER :: it
     LOGICAL :: newmonth ! indicates if a new month just started
     LOGICAL :: newday   ! indicates if a new day just started
@@ -233,10 +316,11 @@
 
     REAL, PARAMETER :: fact=1.E-15/2.12  ! transformation factor from gC/m2/day => ppm/m2/day
-    REAL, DIMENSION(klon) :: fco2_tmp, tr_seri_sum
+    REAL, DIMENSION(klon) :: fco2_tmp
     REAL :: sumtmp
-    REAL :: airetot     ! Total area the earth
     REAL :: delta_co2_ppm
     
-! -) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
+
+! 1) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
+! -------------------------------------------------------------------------------------------------------
 
     newday = .FALSE.; endday = .FALSE.; newmonth = .FALSE.
@@ -245,152 +329,134 @@
     IF (MOD(nstep,INT(86400./pdtphys))==0) endday=.TRUE.
     IF (newday .AND. day_cur==1) newmonth=.TRUE.
-    
-! -) Read new maps if new month started
-    IF (newmonth .AND. carbon_cycle_tr) THEN
-       CALL read_map2D('fossil_fuel.nc','fos_fuel', mth_cur, .FALSE., fos_fuel)
-       
-       ! division by month lenght to get dayly value
-       fos_fuel(:) = fos_fuel(:)/mth_len
-       
-       IF (.NOT. carbon_cycle_cpl) THEN
-          ! Get dayly values from monthly fluxes
-          CALL read_map2D('fl_co2_ocean.nc','CO2_OCN',mth_cur,.FALSE.,fco2_ocn_day)
-          CALL read_map2D('fl_co2_land.nc','CO2_LAND', mth_cur,.FALSE.,fco2_land_day)
-          CALL read_map2D('fl_co2_land_use.nc','CO2_LAND_USE',mth_cur,.FALSE.,fco2_lu_day)
-       END IF
-    END IF
-    
-
-
-! -) Update tracers at beginning of a new day. Beginning of a new day correspond to a new coupling period in cpl_mod.
-    IF (newday) THEN
+
+! 2)  For each carbon tracer find out if it is time to inject (update)
+! --------------------------------------------------------------------
+    DO it = 1, ntr_co2
+       IF ( MOD(nstep,INT(co2trac(it)%updatefreq/pdtphys)) == 1 ) THEN
+          co2trac(it)%updatenow = .TRUE.
+       ELSE
+          co2trac(it)%updatenow = .FALSE.
+       END IF
+    END DO
+
+! 3) Get tracer update
+! --------------------------------------
+    DO it = 1, ntr_co2
+       IF ( co2trac(it)%updatenow ) THEN
+          IF ( co2trac(it)%cpl ) THEN
+             ! Get tracer from coupled model
+             SELECT CASE(co2trac(it)%name)
+             CASE('fCO2_land')     ! from ORCHIDEE
+                dtr_add(:,it) = fco2_land_inst(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
+             CASE('fCO2_land_use') ! from ORCHIDEE
+                dtr_add(:,it) = fco2_lu_inst(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
+             CASE('fCO2_ocn')      ! from PISCES
+                dtr_add(:,it) = fco2_ocn_day(:)  *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
+             CASE DEFAULT
+                WRITE(lunout,*) 'Error with tracer ',co2trac(it)%name
+                CALL abort_gcm('carbon_cycle', 'No coupling implemented for this tracer',1)
+             END SELECT
+          ELSE
+             ! Read tracer from file
+             co2trac(it)%readstep = co2trac(it)%readstep + 1 ! increment time step in file
+! Patricia   CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.FALSE.,dtr_add(:,it))
+             CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.TRUE.,dtr_add(:,it))
+
+             ! Converte from kgC/m2/h to kgC/m2/s
+             dtr_add(:,it) = dtr_add(:,it)/3600
+             ! Add individual treatment of values read from file
+             SELECT CASE(co2trac(it)%name)
+             CASE('fCO2_land')
+                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
+             CASE('fCO2_land_use')
+                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
+             CASE('fCO2_ocn')
+                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_oce)
+! Patricia :
+!             CASE('fCO2_fos_fuel')
+!                dtr_add(:,it) = dtr_add(:,it)/mth_len
+!                co2trac(it)%readstep = 0 ! Always read same value for fossil fuel(Cadule case)
+             END SELECT
+          END IF
+       END IF
+    END DO
+
+! 4) Update co2 tracers : 
+!    Loop over all carbon tracers and add source
+! ------------------------------------------------------------------
+    IF (carbon_cycle_tr) THEN
+       DO it = 1, ntr_co2
+          IF (.FALSE.) THEN
+             tr_seri(1:klon,1,co2trac(it)%id) = tr_seri(1:klon,1,co2trac(it)%id) + dtr_add(1:klon,it)
+             source(1:klon,co2trac(it)%id) = 0.
+          ELSE
+             source(1:klon,co2trac(it)%id) = dtr_add(1:klon,it)
+          END IF
+       END DO
+    END IF
+
+
+! 5) Calculations for new CO2 value for the radiation scheme(instead of reading value from .def)
+! ----------------------------------------------------------------------------------------------
+    IF (RCO2_inter) THEN
+       ! Cumulate fluxes from ORCHIDEE at each timestep
+       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
+          IF (newday) THEN ! Reset cumulative variables once a day 
+             fco2_land_day(1:klon) = 0.
+             fco2_lu_day(1:klon)   = 0.
+          END IF
+          fco2_land_day(1:klon) = fco2_land_day(1:klon) + fco2_land_inst(1:klon) ![gC/m2/day]
+          fco2_lu_day(1:klon)   = fco2_lu_day(1:klon)   + fco2_lu_inst(1:klon)   ![gC/m2/day]
+       END IF
+
+       ! At the end of a new day, calculate a mean scalare value of CO2
+       ! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ? 
+       IF (endday) THEN
+
+          IF (carbon_cycle_tr) THEN
+             ! Sum all co2 tracers to get the total delta CO2 flux
+             fco2_tmp(:) = 0.
+             DO it = 1, ntr_co2
+                fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
+             END DO
+             
+          ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
+             ! Sum co2 fluxes comming from coupled models and parameter for fossil fuel
+             fco2_tmp(1:klon) = fos_fuel_s + ((fco2_lu_day(1:klon) + fco2_land_day(1:klon))*pctsrf(1:klon,is_ter) &
+                  + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
+          END IF
+
+          ! Calculate a global mean value of delta CO2 flux
+          fco2_tmp(1:klon) = fco2_tmp(1:klon) * airephy(1:klon)
+          CALL reduce_sum(SUM(fco2_tmp),sumtmp)
+          CALL bcast(sumtmp)
+          delta_co2_ppm = sumtmp/airetot
+          
+          ! Add initial value for co2_ppm and delta value
+          co2_ppm = co2_ppm0 + delta_co2_ppm
+          
+          ! Transformation of atmospheric CO2 concentration for the radiation code
+          RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97 
+          
+          WRITE(lunout,*) 'RCO2 is now updated! RCO2 = ', RCO2
+       END IF ! endday
+
+    END IF ! RCO2_inter
+
+
+! 6) Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
+! ----------------------------------------------------------------------------
+    IF (carbon_cycle_cpl) THEN
 
        IF (carbon_cycle_tr) THEN
-
-          ! Update tracers
-          IF (ntr_co2 == 1) THEN
-             ! Calculate the new flux CO2
-             tr_seri(:,1,id_fco2_tot) = tr_seri(:,1,id_fco2_tot) + &
-                  (fos_fuel(:) + &
-                  fco2_lu_day(:)  * pctsrf(:,is_ter) + &
-                  fco2_land_day(:)* pctsrf(:,is_ter) + &
-                  fco2_ocn_day(:) * pctsrf(:,is_oce)) * fact
-
-          ELSE ! ntr_co2 == 4
-             tr_seri(:,1,id_fco2_fos_fuel)  = tr_seri(:,1,id_fco2_fos_fuel) + fos_fuel(:) * fact ! [ppm/m2/day]
-
-             tr_seri(:,1,id_fco2_land_use)  = tr_seri(:,1,id_fco2_land_use) + &
-                  fco2_lu_day(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
-
-             tr_seri(:,1,id_fco2_land)      = tr_seri(:,1,id_fco2_land)     + &
-                  fco2_land_day(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
-
-             tr_seri(:,1,id_fco2_ocn)       = tr_seri(:,1,id_fco2_ocn)      + &
-                  fco2_ocn_day(:) *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
-
-          END IF
-
-       ELSE ! no transport
-          IF (carbon_cycle_cpl) THEN
-             IF (carbon_cycle_emis_comp) THEN
-                ! Calcul emission compatible a partir des champs 2D et co2_ppm
-                !!! TO DO!!
-                CALL abort_gcm('carbon_cycle', ' Option carbon_cycle_emis_comp not yet implemented',1)
-             END IF
-          END IF
-
-       END IF ! carbon_cycle_tr
-    
-       ! Reset cumluative variables
-       IF (carbon_cycle_cpl) THEN
-          fco2_land_day(:) = 0.
-          fco2_lu_day(:)   = 0.
-       END IF
-    
-    END IF ! newday
-       
-
-
-! -) Cumulate fluxes from ORCHIDEE at each timestep
-    IF (carbon_cycle_cpl) THEN
-       fco2_land_day(:) = fco2_land_day(:) + fco2_land_inst(:)
-       fco2_lu_day(:)   = fco2_lu_day(:)   + fco2_lu_inst(:)
-    END IF
-
-
-
-! -)  At the end of a new day, calculate a mean scalare value of CO2 to be used by 
-!     the radiation scheme (instead of reading value from .def)
-
-! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ? 
-    IF (endday) THEN 
-       ! Calculte total area of the earth surface
-       CALL reduce_sum(SUM(airephy),airetot)
-       
-
-       IF (carbon_cycle_tr) THEN
-          IF (ntr_co2 == 1) THEN
-          
-             ! Calculate mean value of tracer CO2 to get an scalare value to be used in the 
-             ! radiation scheme (instead of reading value from .def)
-             ! Mean value weighted with the grid cell area
-             
-             ! Calculate mean value
-             fco2_tmp(:) = tr_seri(:,1,id_fco2_tot) * airephy(:) 
-             CALL reduce_sum(SUM(fco2_tmp),sumtmp)
-             co2_ppm = sumtmp/airetot + co2_ppm0
-             
-          ELSE ! ntr_co2 == 4
-             
-             ! Calculate the delta CO2 flux
-             tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
-                  tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
-             
-             ! Calculate mean value of delta CO2 flux
-             fco2_tmp(:) = tr_seri_sum(:) * airephy(:)
-             CALL reduce_sum(SUM(fco2_tmp),sumtmp)
-             delta_co2_ppm = sumtmp/airetot
-             
-             ! Add initial value for co2_ppm to delta value
-             co2_ppm = delta_co2_ppm + co2_ppm0
-          END IF
-
-       ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
-
-          ! Calculate the total CO2 flux and integrate to get scalare value for the radiation scheme
-          fco2_tmp(:) = (fos_fuel(:) + (fco2_lu_day(:) + fco2_land_day(:))*pctsrf(:,is_ter) &
-               + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
-          
-          ! Calculate mean value
-          fco2_tmp(:) = fco2_tmp(:) * airephy(:)
-          CALL reduce_sum(SUM(fco2_tmp),sumtmp)
-          delta_co2_ppm = sumtmp/airetot
-
-          ! Update current value of the atmospheric co2_ppm
-          co2_ppm = co2_ppm + delta_co2_ppm
-          
-       END IF ! carbon_cycle_tr
-
-       ! transformation of the atmospheric CO2 concentration for the radiation code
-       RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97 
-
-    END IF
-
-    ! Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
-    IF (endday .AND. carbon_cycle_cpl) THEN
-       
-       IF (carbon_cycle_tr) THEN
-          IF (ntr_co2==1) THEN
-
-             co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
-
-          ELSE ! ntr_co2 == 4
-
-             co2_send(:) = tr_seri_sum(:) + co2_ppm0
-
-          END IF
+          ! Sum all co2 tracers to get the total delta CO2 flux at first model layer
+          fco2_tmp(:) = 0.
+          DO it = 1, ntr_co2
+             fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
+          END DO
+          co2_send(1:klon) = fco2_tmp(1:klon) + co2_ppm0
        ELSE
           ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
-          co2_send(:) = co2_ppm
+          co2_send(1:klon) = co2_ppm
        END IF
 
Index: /LMDZ5/trunk/libf/phylmd/change_srf_frac_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/change_srf_frac_mod.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/change_srf_frac_mod.F90	(revision 1454)
@@ -9,5 +9,6 @@
 ! 
 ! Change Surface Fractions
-!
+! Author J Ghattas 2008
+
   SUBROUTINE change_srf_frac(itime, dtime, jour, &
        pctsrf, alb1, alb2, tsurf, u10m, v10m, pbl_tke)
@@ -76,5 +77,5 @@
     END SELECT
 
-    IF (is_modified) THEN
+
 !****************************************************************************************
 ! 2) 
@@ -84,4 +85,5 @@
 !
 !****************************************************************************************
+    IF (is_modified) THEN
   
 ! Test and exit if a fraction is negative
@@ -150,4 +152,8 @@
        CALL pbl_surface_newfrac(itime, pctsrf, pctsrf_old, tsurf, alb1, alb2, u10m, v10m, pbl_tke)
 
+    ELSE
+       ! No modifcation should be done
+       pctsrf(:,:) = pctsrf_old(:,:)
+
     END IF ! is_modified
 
Index: /LMDZ5/trunk/libf/phylmd/cpl_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/cpl_mod.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/cpl_mod.F90	(revision 1454)
@@ -100,4 +100,5 @@
   SUBROUTINE cpl_init(dtime, rlon, rlat)
     USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, fco2_ocn_day
+    USE surface_data
 
     INCLUDE "dimensions.h"
@@ -270,4 +271,14 @@
     ENDIF    ! is_sequential
     
+
+!*************************************************************************************
+! compatibility test
+!
+!*************************************************************************************
+    IF (carbon_cycle_cpl .AND. version_ocean=='opa8') THEN
+       abort_message='carbon_cycle_cpl does not work with opa8'
+       CALL abort_gcm(modname,abort_message,1)
+    END IF
+
   END SUBROUTINE cpl_init
   
Index: DZ5/trunk/libf/phylmd/initphysto.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/initphysto.F	(revision 1453)
+++ 	(revision )
@@ -1,305 +1,0 @@
-!
-! $Id$
-!
-C
-C
-      subroutine initphysto
-     .  (infile,
-     .  rlon, rlat, tstep,t_ops,t_wrt,nq,fileid)
-       
-       USE dimphy
-       USE mod_phys_lmdz_para
-       USE IOIPSL
-       USE iophy
-       USE control_mod
-
-      implicit none
-
-C
-C   Routine d'initialisation des ecritures des fichiers histoires LMDZ
-C   au format IOIPSL
-C
-C   Appels succesifs des routines: histbeg
-C                                  histhori
-C                                  histver
-C                                  histdef
-C                                  histend
-C
-C   Entree:
-C
-C      infile: nom du fichier histoire a creer
-C      day0,anne0: date de reference
-C      tstep: duree du pas de temps en seconde
-C      t_ops: frequence de l'operation pour IOIPSL
-C      t_wrt: frequence d'ecriture sur le fichier
-C      nq: nombre de traceurs
-C
-C   Sortie:
-C      fileid: ID du fichier netcdf cree
-C      filevid:ID du fichier netcdf pour la grille v
-C
-C   L. Fairhead, LMD, 03/99
-C
-C =====================================================================
-C
-C   Declarations
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "temps.h"
-#include "ener.h"
-#include "logic.h"
-#include "description.h"
-#include "serre.h"
-#include "indicesol.h"
-cym#include "dimphy.h"
-
-C   Arguments
-      character*(*) infile
-      integer nhoriid, i
-      real tstep, t_ops, t_wrt
-      integer fileid, filevid
-      integer nq,l
-      real nivsigs(llm)
-
-C   Variables locales
-C
-      integer tau0
-      real zjulian
-      character*3 str
-      character*10 ctrac
-      integer iq
-      integer uhoriid, vhoriid, thoriid, zvertiid
-      integer ii,jj
-      integer zan, idayref
-      logical ok_sync
-      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
-C
-      REAL rlon(klon), rlat(klon)
-
-C  Initialisations
-C
-      pi = 4. * atan (1.)
-      str='q  '
-      ctrac = 'traceur   '
-      ok_sync= .true.
-C
-C  Appel a histbeg: creation du fichier netcdf et initialisations diverses
-C         
-
-      zan = annee_ref
-      idayref = day_ref
-      CALL ymds2ju(zan, 1, idayref, 0.0, zjulian)
-      tau0 = 0
-	
-cym	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
-cym         DO i = 1, iim
-cym            zx_lon(i,1) = rlon(i+1)
-cym            zx_lon(i,jjm+1) = rlon(i+1)
-cym         ENDDO
-cym         CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
-
-
-      call histbeg_phy(infile,tau0, zjulian, tstep,
-     .                 nhoriid, fileid)
-
-c$OMP MASTER	
-C  Appel a histvert pour la grille verticale
-C
-	DO l=1,llm
-            nivsigs(l)=REAL(l)
-         ENDDO
-
-	write(*,*) 'avant histvert ds initphysto'
-
-      call histvert(fileid, 'sig_s', 'Niveaux sigma',
-     . 'sigma_level',
-     .              llm, nivsigs, zvertiid)
-C
-C  Appels a histdef pour la definition des variables a sauvegarder
-C
-	write(*,*) 'apres histvert ds initphysto'
-
-       CALL histdef(fileid, "phis", "Surface geop. height", "-",
-     .                iim,jj_nb,nhoriid, 1,1,1, -99, 32,
-     .                "once", t_ops, t_wrt)
-c
-	write(*,*) 'apres phis ds initphysto'
-
-         CALL histdef(fileid, "aire", "Grid area", "-",
-     .                iim,jj_nb,nhoriid, 1,1,1, -99, 32,
-     .                "once", t_ops, t_wrt)
-         write(*,*) 'apres aire ds initphysto'
-
-cym     Attention dtime et istphy ne sont pas �rit ---> a �iminer ?
-         CALL histdef(fileid, "dtime", "tps phys ", "s",
-     .                1,1,nhoriid, 1,1,1, -99, 32,
-     .                "once", t_ops, t_wrt)
-	
-	 CALL histdef(fileid, "istphy", "tps stock", "s",
-     .                1,1,nhoriid, 1,1,1, -99, 32,
-     .                "once", t_ops, t_wrt)
-
-C T 
-C
-      call histdef(fileid, 't', 'Temperature', 'K',
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
-        write(*,*) 'apres t ds initphysto'
-C mfu 
-C
-      call histdef(fileid, 'mfu', 'flx m. pan. mt', 'kg m/s',
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
-        write(*,*) 'apres mfu ds initphysto'
-C
-C mfd 
-C
-      call histdef(fileid, 'mfd', 'flx m. pan. des', 'kg m/s',
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
-
-C
-C en_u 
-C
-      call histdef(fileid, 'en_u', 'flx ent pan mt', 'kg m/s',
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
-               write(*,*) 'apres en_u ds initphysto'
-C
-C de_u 
-C
-      call histdef(fileid, 'de_u', 'flx det pan mt', 'kg m/s',
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
-
-C
-C en_d 
-C
-      call histdef(fileid, 'en_d', 'flx ent pan dt', 'kg m/s',
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
-C
-
-C
-C de_d 
-C
-      call histdef(fileid, 'de_d', 'flx det pan dt', 'kg m/s',
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
-
-c coefh frac_impa,frac_nucl
-	
-	call histdef(fileid, "coefh", " ", " ",
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, "inst(X)", t_ops, t_wrt)
-
-c abderrahmane le 16 09 02
-        call histdef(fileid, "fm_th", " ", " ",
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, "inst(X)", t_ops, t_wrt)
-
-        call histdef(fileid, "en_th", " ", " ",
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, "inst(X)", t_ops, t_wrt)
-c fin aj
-	
-	write(*,*) 'apres coefh ds initphysto'	
-
-	call histdef(fileid, 'frac_impa', ' ', ' ',
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
-	
-	call histdef(fileid, 'frac_nucl', ' ', ' ',
-     .             iim, jj_nb, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
-
-c
-c pyu1
-c
-      CALL histdef(fileid, "pyu1", " ", " ",
-     .                iim,jj_nb,nhoriid, 1,1,1, -99, 32,
-     .                "inst(X)", t_ops, t_wrt)
-
-c
-c pyv1
-c
-	CALL histdef(fileid, "pyv1", " ", " ",
-     .                iim,jj_nb,nhoriid, 1,1,1, -99, 32,
-     .                "inst(X)", t_ops, t_wrt)
-	
-	write(*,*) 'apres pyv1 ds initphysto'
-c
-c ftsol1
-c
-	call histdef(fileid, "ftsol1", " ", " ",
-     .             iim, jj_nb, nhoriid, 1, 1,1, -99,32,
-     .             "inst(X)", t_ops, t_wrt)
-
-c
-c ftsol2
-c
-        call histdef(fileid, "ftsol2", " ", " ",
-     .             iim, jj_nb, nhoriid, 1, 1,1, -99,32,
-     .             "inst(X)", t_ops, t_wrt)
-
-c
-c ftsol3
-c
-        call histdef(fileid, "ftsol3", " ", " ",
-     .             iim, jj_nb, nhoriid, 1, 1,1, -99,
-     .             32, "inst(X)", t_ops, t_wrt)
-
-c
-c ftsol4
-c
-        call histdef(fileid, "ftsol4", " ", " ",
-     .             iim, jj_nb, nhoriid, 1, 1,1, -99,
-     .             32, "inst(X)", t_ops, t_wrt)
-	
-c
-c rain
-c
-        call histdef(fileid, "rain", " ", " ",
-     .             iim, jj_nb, nhoriid, 1, 1,1, -99,
-     .             32, "inst(X)", t_ops, t_wrt)
-
-c
-c psrf1
-c
-	call histdef(fileid, "psrf1", " ", " ",
-     .             iim, jj_nb, nhoriid, 1, 1, 1, -99,
-     .             32, "inst(X)", t_ops, t_wrt)
-	
-c
-c psrf2
-c
-        call histdef(fileid, "psrf2", " ", " ",
-     .             iim, jj_nb, nhoriid, 1, 1, 1, -99,
-     .             32, "inst(X)", t_ops, t_wrt)
-
-c
-c psrf3
-c
-        call histdef(fileid, "psrf3", " ", " ",
-     .             iim, jj_nb, nhoriid, 1, 1, 1, -99,
-     .             32, "inst(X)", t_ops, t_wrt)
-
-c
-c psrf4
-c
-        call histdef(fileid, "psrf4", " ", " ",
-     .             iim, jj_nb, nhoriid, 1, 1, 1, -99,
-     .             32, "inst(X)", t_ops, t_wrt)
-	
-	write(*,*) 'avant histend ds initphysto'	
-
-      call histend(fileid)
-c     if (ok_sync) call histsync(fileid)
-      if (ok_sync) call histsync
-c$OMP END MASTER
-	
-
-      return
-      end
Index: /LMDZ5/trunk/libf/phylmd/initphysto.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/initphysto.F90	(revision 1454)
+++ /LMDZ5/trunk/libf/phylmd/initphysto.F90	(revision 1454)
@@ -0,0 +1,216 @@
+!
+! $Id$
+!
+SUBROUTINE initphysto(infile,tstep,t_ops,t_wrt,fileid)
+  
+  USE dimphy
+  USE mod_phys_lmdz_para
+  USE IOIPSL
+  USE iophy
+  USE control_mod
+  
+  IMPLICIT NONE
+
+!
+!   Routine d'initialisation des ecritures des fichiers histoires LMDZ
+!   au format IOIPSL
+!
+!   Appels succesifs des routines: histbeg
+!                                  histhori
+!                                  histver
+!                                  histdef
+!                                  histend
+!
+!   Entree:
+!
+!      infile: nom du fichier histoire a creer
+!      day0,anne0: date de reference
+!      tstep: duree du pas de temps en seconde
+!      t_ops: frequence de l'operation pour IOIPSL
+!      t_wrt: frequence d'ecriture sur le fichier
+!
+!   Sortie:
+!      fileid: ID du fichier netcdf cree
+!
+!   L. Fairhead, LMD, 03/99
+!
+! =====================================================================
+!
+!   Declarations
+  INCLUDE "dimensions.h"
+  INCLUDE "paramet.h"
+  INCLUDE "comconst.h"
+  INCLUDE "comgeom.h"
+  INCLUDE "temps.h"
+  INCLUDE "logic.h"
+  INCLUDE "description.h"
+  INCLUDE "serre.h"
+  INCLUDE "indicesol.h"
+
+!   Arguments
+  CHARACTER(len=*), INTENT(IN) :: infile
+  REAL, INTENT(IN)             :: tstep
+  REAL, INTENT(IN)             :: t_ops
+  REAL, INTENT(IN)             :: t_wrt
+  INTEGER, INTENT(OUT)         :: fileid
+
+! Variables locales
+  INTEGER nhoriid, i
+  INTEGER l,k
+  REAL nivsigs(llm)
+  INTEGER tau0
+  REAL zjulian
+  INTEGER iq
+  INTEGER uhoriid, vhoriid, thoriid, zvertiid
+  INTEGER ii,jj
+  INTEGER zan, idayref
+  LOGICAL ok_sync
+  REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
+  CHARACTER(len=12) :: nvar
+
+!  Initialisations
+!
+  pi = 4. * ATAN (1.)
+  ok_sync= .TRUE.
+!
+!  Appel a histbeg: creation du fichier netcdf et initialisations diverses
+!         
+
+  zan = annee_ref
+  idayref = day_ref
+  CALL ymds2ju(zan, 1, idayref, 0.0, zjulian)
+  tau0 = 0
+  
+  CALL histbeg_phy(infile,tau0, zjulian, tstep, &
+       nhoriid, fileid)
+
+!$OMP MASTER	
+!  Appel a histvert pour la grille verticale
+!
+  DO l=1,llm
+     nivsigs(l)=REAL(l)
+  ENDDO
+  
+  CALL histvert(fileid, 'sig_s', 'Niveaux sigma', &
+       'sigma_level', &
+       llm, nivsigs, zvertiid)
+!
+!  Appels a histdef pour la definition des variables a sauvegarder
+!
+  CALL histdef(fileid, "phis", "Surface geop. height", "-", &
+       iim,jj_nb,nhoriid, 1,1,1, -99, 32, &
+       "once", t_ops, t_wrt)
+  
+  CALL histdef(fileid, "aire", "Grid area", "-", &
+       iim,jj_nb,nhoriid, 1,1,1, -99, 32, &
+       "once", t_ops, t_wrt)
+
+  CALL histdef(fileid, "longitudes", "longitudes", "-", &
+       iim,jj_nb,nhoriid, 1,1,1, -99, 32, &
+       "once", t_ops, t_wrt)
+
+  CALL histdef(fileid, "latitudes", "latitudes", "-", &
+       iim,jj_nb,nhoriid, 1,1,1, -99, 32, &
+       "once", t_ops, t_wrt)
+! T 
+  CALL histdef(fileid, 't', 'Temperature', 'K', iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid, 32, 'inst(X)', t_ops, t_wrt)
+! mfu 
+  CALL histdef(fileid, 'mfu', 'flx m. pan. mt', 'kg m/s',iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid,32, 'inst(X)', t_ops, t_wrt)
+! mfd 
+  CALL histdef(fileid, 'mfd', 'flx m. pan. des', 'kg m/s',iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid, 32, 'inst(X)', t_ops, t_wrt)
+! en_u 
+  CALL histdef(fileid, 'en_u', 'flx ent pan mt', 'kg m/s', iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid,32, 'inst(X)', t_ops, t_wrt)
+! de_u 
+  CALL histdef(fileid, 'de_u', 'flx det pan mt', 'kg m/s',iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid,32, 'inst(X)', t_ops, t_wrt)
+! en_d 
+  CALL histdef(fileid, 'en_d', 'flx ent pan dt', 'kg m/s', iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid,32, 'inst(X)', t_ops, t_wrt)
+! de_d 
+  CALL histdef(fileid, 'de_d', 'flx det pan dt', 'kg m/s', iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid, 32, 'inst(X)', t_ops, t_wrt)
+! coefh
+  CALL histdef(fileid, "coefh", " ", " ", iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid,32, "inst(X)", t_ops, t_wrt)
+! fm_th
+  CALL histdef(fileid, "fm_th", " ", " ",iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid,32, "inst(X)", t_ops, t_wrt)
+! en_th
+  CALL histdef(fileid, "en_th", " ", " ",iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid,32, "inst(X)", t_ops, t_wrt)
+! frac_impa
+  CALL histdef(fileid, 'frac_impa', ' ', ' ',iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid,32, 'inst(X)', t_ops, t_wrt)
+! frac_nucl
+  CALL histdef(fileid, 'frac_nucl', ' ', ' ',iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid,32, 'inst(X)', t_ops, t_wrt)
+! pyu1
+  CALL histdef(fileid, "pyu1", " ", " ", iim,jj_nb,nhoriid, &
+       1,1,1, -99, 32, "inst(X)", t_ops, t_wrt)
+! pyv1
+  CALL histdef(fileid, "pyv1", " ", " ", iim,jj_nb,nhoriid, &
+       1,1,1, -99, 32,"inst(X)", t_ops, t_wrt)    
+! ftsol1
+  CALL histdef(fileid, "ftsol1", " ", " ",iim, jj_nb, nhoriid, &
+       1, 1,1, -99,32, "inst(X)", t_ops, t_wrt)
+! ftsol2
+  CALL histdef(fileid, "ftsol2", " ", " ",iim, jj_nb, nhoriid, &
+       1, 1,1, -99,32, "inst(X)", t_ops, t_wrt)
+! ftsol3
+  CALL histdef(fileid, "ftsol3", " ", " ", iim, jj_nb, nhoriid, &
+       1, 1,1, -99,32, "inst(X)", t_ops, t_wrt)
+! ftsol4
+  CALL histdef(fileid, "ftsol4", " ", " ",iim, jj_nb, nhoriid, &
+       1, 1,1, -99, 32, "inst(X)", t_ops, t_wrt)
+! psrf1
+  CALL histdef(fileid, "psrf1", " ", " ",iim, jj_nb, nhoriid, &
+       1, 1, 1, -99,32, "inst(X)", t_ops, t_wrt)
+! psrf2
+  CALL histdef(fileid, "psrf2", " ", " ",iim, jj_nb, nhoriid, &
+       1, 1, 1, -99, 32, "inst(X)", t_ops, t_wrt)
+! psrf3
+  CALL histdef(fileid, "psrf3", " ", " ",iim, jj_nb, nhoriid, &
+       1, 1, 1, -99, 32, "inst(X)", t_ops, t_wrt)
+! psrf4
+  CALL histdef(fileid, "psrf4", " ", " ", iim, jj_nb, nhoriid, &
+       1, 1, 1, -99,32, "inst(X)", t_ops, t_wrt)
+! sh
+  CALL histdef(fileid, 'sh', '', '', iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid, 32, 'inst(X)', t_ops, t_wrt)
+! da
+  CALL histdef(fileid, 'da', '', '', iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid, 32, 'inst(X)', t_ops, t_wrt)
+! mp
+  CALL histdef(fileid, 'mp', '', '', iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid, 32, 'inst(X)', t_ops, t_wrt)
+! upwd
+  CALL histdef(fileid, 'upwd', '', '', iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid, 32, 'inst(X)', t_ops, t_wrt)
+! dnwd
+  CALL histdef(fileid, 'dnwd', '', '', iim, jj_nb, nhoriid, &
+       llm, 1, llm, zvertiid, 32, 'inst(X)', t_ops, t_wrt)
+
+! phi
+  DO k=1,llm
+     IF (k<10) THEN
+        WRITE(nvar,'(i1)') k
+     ELSE IF (k<100) THEN
+        WRITE(nvar,'(i2)') k
+     ELSE
+        WRITE(nvar,'(i3)') k
+     END IF
+     nvar='phi_lev'//trim(nvar)
+     
+     CALL histdef(fileid, nvar, '', '', iim, jj_nb, nhoriid, &
+          llm, 1, llm, zvertiid, 32, 'inst(X)', t_ops, t_wrt)
+  END DO
+
+  CALL histend(fileid)
+  IF (ok_sync) CALL histsync
+!$OMP END MASTER
+	
+END SUBROUTINE initphysto
Index: /LMDZ5/trunk/libf/phylmd/phyetat0.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phyetat0.F	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/phyetat0.F	(revision 1454)
@@ -21,5 +21,6 @@
       USE infotrac
       USE traclmdz_mod,    ONLY : traclmdz_from_restart
-      USE carbon_cycle_mod,ONLY : carbon_cycle_tr, carbon_cycle_cpl
+      USE carbon_cycle_mod,ONLY :
+     &     carbon_cycle_tr, carbon_cycle_cpl, co2_send
 
       IMPLICIT none
@@ -133,37 +134,12 @@
 
        
-
-         IF( clesphy0(1).NE.tab_cntrl( 5 ) )  THEN
-             clesphy0(1)=tab_cntrl( 5 )
-         ENDIF
-
-         IF( clesphy0(2).NE.tab_cntrl( 6 ) )  THEN
-             clesphy0(2)=tab_cntrl( 6 )
-         ENDIF
-
-         IF( clesphy0(3).NE.tab_cntrl( 7 ) )  THEN
-             clesphy0(3)=tab_cntrl( 7 )
-         ENDIF
-
-         IF( clesphy0(4).NE.tab_cntrl( 8 ) )  THEN
-             clesphy0(4)=tab_cntrl( 8 )
-         ENDIF
-
-         IF( clesphy0(5).NE.tab_cntrl( 9 ) )  THEN
-             clesphy0(5)=tab_cntrl( 9 )
-         ENDIF
-
-         IF( clesphy0(6).NE.tab_cntrl( 10 ) )  THEN
-             clesphy0(6)=tab_cntrl( 10 )
-         ENDIF
-
-         IF( clesphy0(7).NE.tab_cntrl( 11 ) )  THEN
-             clesphy0(7)=tab_cntrl( 11 )
-         ENDIF
-
-         IF( clesphy0(8).NE.tab_cntrl( 12 ) )  THEN
-             clesphy0(8)=tab_cntrl( 12 )
-         ENDIF
-
+      clesphy0(1)=tab_cntrl( 5 )
+      clesphy0(2)=tab_cntrl( 6 )
+      clesphy0(3)=tab_cntrl( 7 )
+      clesphy0(4)=tab_cntrl( 8 )
+      clesphy0(5)=tab_cntrl( 9 )
+      clesphy0(6)=tab_cntrl( 10 )
+      clesphy0(7)=tab_cntrl( 11 )
+      clesphy0(8)=tab_cntrl( 12 )
 
 c
@@ -1078,6 +1054,17 @@
 
          END DO
-         
          CALL traclmdz_from_restart(trs)
+
+         IF (carbon_cycle_cpl) THEN
+            ALLOCATE(co2_send(klon), stat=ierr)
+            IF (ierr /= 0) CALL abort_gcm
+     &           ('phyetat0','pb allocation co2_send',1)
+            CALL get_field("co2_send",co2_send,found)
+            IF (.NOT. found) THEN
+               PRINT*,"phyetat0: Le champ <co2_send> est absent"
+               PRINT*,"Initialisation uniforme a co2_ppm=",co2_ppm
+               co2_send(:) = co2_ppm
+            END IF
+         END IF
       END IF
 
Index: /LMDZ5/trunk/libf/phylmd/phyredem.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phyredem.F	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/phyredem.F	(revision 1454)
@@ -15,5 +15,5 @@
       USE infotrac
       USE control_mod
-
+      USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
 
       IMPLICIT none
@@ -337,4 +337,12 @@
             CALL put_field("trs_"//tname(iiq),"",trs(:,it))
          END DO
+         IF (carbon_cycle_cpl) THEN
+            IF (.NOT. ALLOCATED(co2_send)) THEN
+               ! This is the case of create_etat0_limit, ce0l
+               ALLOCATE(co2_send(klon))
+               co2_send(:) = co2_ppm0
+            END IF
+            CALL put_field("co2_send","co2_ppm for coupling",co2_send)
+         END IF
       END IF
 
Index: /LMDZ5/trunk/libf/phylmd/physiq.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1454)
@@ -1250,5 +1250,4 @@
 cym Attention pbase pas initialise dans concvl !!!!
           pbase=0
-          paire_ter(:)=0.    
 cIM 180608
 c         pmflxr=0.
@@ -3376,5 +3375,6 @@
      I                   cdragh,coefh,u1,v1,ftsol,pctsrf,
      I                   frac_impa, frac_nucl,
-     I                   pphis,airephy,dtime,itap)
+     I                   pphis,airephy,dtime,itap,
+     I                   rlon,rlat,qx(:,:,ivap),da,phi,mp,upwd,dnwd)
 
 
Index: DZ5/trunk/libf/phylmd/phystokenc.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phystokenc.F	(revision 1453)
+++ 	(revision )
@@ -1,435 +1,0 @@
-!
-c
-c
-      SUBROUTINE phystokenc (
-     I                   nlon,nlev,pdtphys,rlon,rlat,
-     I                   pt,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
-     I                   pfm_therm,pentr_therm,
-     I                   cdragh, pcoefh,yu1,yv1,ftsol,pctsrf,
-     I                   frac_impa,frac_nucl,
-     I                   pphis,paire,dtime,itap)
-      USE ioipsl
-      USE dimphy
-      USE infotrac, ONLY : nqtot
-      USE iophy
-      USE control_mod
-
-      IMPLICIT none
-
-c======================================================================
-c Auteur(s) FH
-c Objet: Moniteur general des tendances traceurs
-c
-
-c======================================================================
-#include "dimensions.h"
-#include "tracstoke.h"
-#include "indicesol.h"
-c======================================================================
-
-c Arguments:
-c
-c   EN ENTREE:
-c   ==========
-c
-c   divers:
-c   -------
-c
-      integer nlon ! nombre de points horizontaux
-      integer nlev ! nombre de couches verticales
-      real pdtphys ! pas d'integration pour la physique (seconde)
-c
-      integer physid, itap
-      save physid
-c$OMP THREADPRIVATE(physid)
-      integer ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
-
-c   convection:
-c   -----------
-c
-      REAL pmfu(klon,klev)  ! flux de masse dans le panache montant
-      REAL pmfd(klon,klev)  ! flux de masse dans le panache descendant
-      REAL pen_u(klon,klev) ! flux entraine dans le panache montant
-      REAL pde_u(klon,klev) ! flux detraine dans le panache montant
-      REAL pen_d(klon,klev) ! flux entraine dans le panache descendant
-      REAL pde_d(klon,klev) ! flux detraine dans le panache descendant
-      real pt(klon,klev)
-      REAL,allocatable,save :: t(:,:)
-c$OMP THREADPRIVATE(t)
-c
-      REAL rlon(klon), rlat(klon), dtime
-      REAL zx_tmp_3d(iim,jjm+1,klev),zx_tmp_2d(iim,jjm+1)
-
-c   Couche limite:
-c   --------------
-c
-      REAL cdragh(klon)          ! cdrag
-      REAL pcoefh(klon,klev)     ! coeff melange CL
-      REAL pcoefh_buf(klon,klev) ! coeff melange CL + cdrag
-      REAL yv1(klon)
-      REAL yu1(klon),pphis(klon),paire(klon)
-
-c   Les Thermiques : (Abderr 25 11 02)
-c   ---------------
-      REAL pfm_therm(klon,klev+1)
-      real fm_therm1(klon,klev)
-      REAL pentr_therm(klon,klev)
-    
-      REAL,allocatable,save :: entr_therm(:,:)
-      REAL,allocatable,save :: fm_therm(:,:)
-c$OMP THREADPRIVATE(entr_therm)
-c$OMP THREADPRIVATE(fm_therm)
-c
-c   Lessivage:
-c   ----------
-c
-      REAL frac_impa(klon,klev)
-      REAL frac_nucl(klon,klev)
-c
-c Arguments necessaires pour les sources et puits de traceur
-C
-      real ftsol(klon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
-      real pctsrf(klon,nbsrf) ! Pourcentage de sol f(nature du sol)
-c======================================================================
-c
-      INTEGER i, k
-c
-      REAL,allocatable,save :: mfu(:,:)  ! flux de masse dans le panache montant
-      REAL,allocatable,save :: mfd(:,:)  ! flux de masse dans le panache descendant
-      REAL,allocatable,save :: en_u(:,:) ! flux entraine dans le panache montant
-      REAL,allocatable,save :: de_u(:,:) ! flux detraine dans le panache montant
-      REAL,allocatable,save :: en_d(:,:) ! flux entraine dans le panache descendant
-      REAL,allocatable,save :: de_d(:,:) ! flux detraine dans le panache descendant
-      REAL,allocatable,save :: coefh(:,:) ! flux detraine dans le panache descendant
-
-      REAL,allocatable,save :: pyu1(:)
-      REAL,allocatable,save :: pyv1(:)
-      REAL,allocatable,save :: pftsol(:,:)
-      REAL,allocatable,save :: ppsrf(:,:)
-c$OMP THREADPRIVATE(mfu,mfd,en_u,de_u,en_d,de_d,coefh)
-c$OMP THREADPRIVATE(pyu1,pyv1,pftsol,ppsrf)
-      real pftsol1(klon),pftsol2(klon),pftsol3(klon),pftsol4(klon)
-      real ppsrf1(klon),ppsrf2(klon),ppsrf3(klon),ppsrf4(klon)
-
-      REAL dtcum
-
-      integer iadvtr,irec
-      real zmin,zmax
-      logical ok_sync
- 
-      save dtcum
-      save iadvtr,irec
-c$OMP THREADPRIVATE(dtcum,iadvtr,irec)
-      data iadvtr,irec/0,1/
-      logical,save :: first=.true.
-c$OMP THREADPRIVATE(first)
-c
-c   Couche limite:
-c======================================================================
-
-c Dans le meme vecteur on recombine le drag et les coeff d'echange
-      pcoefh_buf(:,1)      = cdragh(:)
-      pcoefh_buf(:,2:klev) = pcoefh(:,2:klev)
-
-      ok_sync = .true.
-	print*,'Dans phystokenc.F'
-      print*,'iadvtr= ',iadvtr
-      print*,'istphy= ',istphy
-      print*,'istdyn= ',istdyn
-
-      if (first) then
-      
-        allocate( t(klon,klev))
-        allocate( mfu(klon,klev))  
-        allocate( mfd(klon,klev))  
-        allocate( en_u(klon,klev)) 
-        allocate( de_u(klon,klev)) 
-        allocate( en_d(klon,klev)) 
-        allocate( de_d(klon,klev)) 
-        allocate( coefh(klon,klev)) 
-        allocate( entr_therm(klon,klev))
-        allocate( fm_therm(klon,klev))
-        allocate( pyu1(klon))
-        allocate( pyv1(klon))
-        allocate( pftsol(klon,nbsrf))
-        allocate( ppsrf(klon,nbsrf))
-  
-        first=.false.
-      endif
-      
-      IF (iadvtr.eq.0) THEN
-	
-	CALL initphysto('phystoke',
-     . rlon,rlat,dtime, dtime*istphy,dtime*istphy,nqtot,physid)
-  	
-	write(*,*) 'apres initphysto ds phystokenc' 
-
-	
-      ENDIF
-c
-      ndex2d = 0
-      ndex3d = 0
-      i=itap 
-cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
-      CALL histwrite_phy(physid,"phis",i,pphis)
-c
-      i=itap
-cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
-      CALL histwrite_phy(physid,"aire",i,paire)
-
-      iadvtr=iadvtr+1
-c
-      if (mod(iadvtr,istphy).eq.1.or.istphy.eq.1) then
-	print*,'reinitialisation des champs cumules 
-     s          a iadvtr=',iadvtr
-         do k=1,klev
-            do i=1,klon
-               mfu(i,k)=0.
-               mfd(i,k)=0.
-               en_u(i,k)=0.
-               de_u(i,k)=0.
-               en_d(i,k)=0.
-               de_d(i,k)=0.
-               coefh(i,k)=0.
-                t(i,k)=0.
-		fm_therm(i,k)=0.
-               entr_therm(i,k)=0.
-            enddo
-         enddo
-         do i=1,klon
-            pyv1(i)=0.
-            pyu1(i)=0.
-         end do
-         do k=1,nbsrf
-             do i=1,klon
-               pftsol(i,k)=0.
-               ppsrf(i,k)=0.
-            enddo
-         enddo
-
-         dtcum=0.
-      endif
-
-      do k=1,klev
-         do i=1,klon
-            mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
-            mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
-            en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
-            de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
-            en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
-            de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
-            coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
-                t(i,k)=t(i,k)+pt(i,k)*pdtphys
-       fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
-       entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
-         enddo
-      enddo
-         do i=1,klon
-            pyv1(i)=pyv1(i)+yv1(i)*pdtphys
-            pyu1(i)=pyu1(i)+yu1(i)*pdtphys
-         end do
-         do k=1,nbsrf
-             do i=1,klon
-               pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
-               ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
-            enddo
-         enddo
-
-      dtcum=dtcum+pdtphys
-
-      IF(mod(iadvtr,istphy).eq.0) THEN 
-c
-c   normalisation par le temps cumule
-         do k=1,klev
-            do i=1,klon
-               mfu(i,k)=mfu(i,k)/dtcum
-               mfd(i,k)=mfd(i,k)/dtcum
-               en_u(i,k)=en_u(i,k)/dtcum
-               de_u(i,k)=de_u(i,k)/dtcum
-               en_d(i,k)=en_d(i,k)/dtcum
-               de_d(i,k)=de_d(i,k)/dtcum
-               coefh(i,k)=coefh(i,k)/dtcum
-c Unitel a enlever
-	      t(i,k)=t(i,k)/dtcum	
-               fm_therm(i,k)=fm_therm(i,k)/dtcum
-	       entr_therm(i,k)=entr_therm(i,k)/dtcum
-            enddo
-         enddo
-         do i=1,klon
-            pyv1(i)=pyv1(i)/dtcum
-            pyu1(i)=pyu1(i)/dtcum
-         end do
-         do k=1,nbsrf
-             do i=1,klon
-               pftsol(i,k)=pftsol(i,k)/dtcum
-               pftsol1(i) = pftsol(i,1)
-               pftsol2(i) = pftsol(i,2)
-               pftsol3(i) = pftsol(i,3)
-               pftsol4(i) = pftsol(i,4)
-
-               ppsrf(i,k)=ppsrf(i,k)/dtcum
-               ppsrf1(i) = ppsrf(i,1)
-               ppsrf2(i) = ppsrf(i,2)
-               ppsrf3(i) = ppsrf(i,3)
-               ppsrf4(i) = ppsrf(i,4)
-
-            enddo
-         enddo
-c
-c   ecriture des champs
-c
-         irec=irec+1
-
-ccccc
-cym         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, t, zx_tmp_3d)
-         CALL histwrite_phy(physid,"t",itap,t)
-
-cym         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfu, zx_tmp_3d)
-      CALL histwrite_phy(physid,"mfu",itap,mfu)
-cym	CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfd, zx_tmp_3d)
-      CALL histwrite_phy(physid,"mfd",itap,mfd)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_u, zx_tmp_3d)
-      CALL histwrite_phy(physid,"en_u",itap,en_u)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_u, zx_tmp_3d)
-      CALL histwrite_phy(physid,"de_u",itap,de_u)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_d, zx_tmp_3d)
-      CALL histwrite_phy(physid,"en_d",itap,en_d)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_d, zx_tmp_3d)       
-      CALL histwrite_phy(physid,"de_d",itap,de_d)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, coefh, zx_tmp_3d)         
-      CALL histwrite_phy(physid,"coefh",itap,coefh)	
-
-c ajou...
-	do k=1,klev
-           do i=1,klon
-	 fm_therm1(i,k)=fm_therm(i,k)	
-	   enddo
-	enddo
-
-cym      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, fm_therm1, zx_tmp_3d)
-      CALL histwrite_phy(physid,"fm_th",itap,fm_therm1)
-c
-cym      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, entr_therm, zx_tmp_3d)
-      CALL histwrite_phy(physid,"en_th",itap,entr_therm)
-cccc
-cym       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_impa,zx_tmp_3d)
-        CALL histwrite_phy(physid,"frac_impa",itap,frac_impa)
-
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_nucl,zx_tmp_3d)
-        CALL histwrite_phy(physid,"frac_nucl",itap,frac_nucl)
- 
-cym        CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyu1,zx_tmp_2d)
-      CALL histwrite_phy(physid,"pyu1",itap,pyu1)
-	
-cym	CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyv1,zx_tmp_2d)
-      CALL histwrite_phy(physid,"pyv1",itap,pyv1)
-	
-cym	CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol1, zx_tmp_2d)
-      CALL histwrite_phy(physid,"ftsol1",itap,pftsol1)
-cym         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol2, zx_tmp_2d)
-      CALL histwrite_phy(physid,"ftsol2",itap,pftsol2)
-cym          CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol3, zx_tmp_2d)
-      CALL histwrite_phy(physid,"ftsol3",itap,pftsol3)
-cym         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol4, zx_tmp_2d)
-      CALL histwrite_phy(physid,"ftsol4",itap,pftsol4)
-
-cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf1, zx_tmp_2d)
-      CALL histwrite_phy(physid,"psrf1",itap,ppsrf1)
-cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf2, zx_tmp_2d)
-      CALL histwrite_phy(physid,"psrf2",itap,ppsrf2)
-cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf3, zx_tmp_2d)
-      CALL histwrite_phy(physid,"psrf3",itap,ppsrf3)
-cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf4, zx_tmp_2d)
-      CALL histwrite_phy(physid,"psrf4",itap,ppsrf4)
-
-c$OMP MASTER
-      if (ok_sync) call histsync(physid)
-c$OMP END MASTER
-c     if (ok_sync) call histsync
-	
-c
-cAA Test sur la valeur des coefficients de lessivage 
-c
-         zmin=1e33
-         zmax=-1e33
-         do k=1,klev
-            do i=1,klon
-                  zmax=max(zmax,frac_nucl(i,k))
-                  zmin=min(zmin,frac_nucl(i,k))
-            enddo
-         enddo
-         Print*,'------ coefs de lessivage (min et max) --------'
-         Print*,'facteur de nucleation ',zmin,zmax
-         zmin=1e33
-         zmax=-1e33
-         do k=1,klev
-            do i=1,klon
-                  zmax=max(zmax,frac_impa(i,k))
-                  zmin=min(zmin,frac_impa(i,k))
-            enddo
-         enddo
-         Print*,'facteur d impaction ',zmin,zmax
-
-      ENDIF 
-
-c   reinitialisation des champs cumules
-	go to 768
-      if (mod(iadvtr,istphy).eq.1) then
-         do k=1,klev
-            do i=1,klon
-               mfu(i,k)=0.
-               mfd(i,k)=0.
-               en_u(i,k)=0.
-               de_u(i,k)=0.
-               en_d(i,k)=0.
-               de_d(i,k)=0.
-               coefh(i,k)=0.
-	       t(i,k)=0.
-               fm_therm(i,k)=0.
-	       entr_therm(i,k)=0.
-            enddo
-         enddo
-         do i=1,klon
-            pyv1(i)=0.
-            pyu1(i)=0.
-         end do
-         do k=1,nbsrf
-             do i=1,klon
-               pftsol(i,k)=0.
-               ppsrf(i,k)=0.
-            enddo
-         enddo
-
-         dtcum=0.
-      endif
-
-      do k=1,klev
-         do i=1,klon
-            mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
-            mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
-            en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
-            de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
-            en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
-            de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
-            coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
-		t(i,k)=t(i,k)+pt(i,k)*pdtphys
-       fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
-       entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
-         enddo
-      enddo
-         do i=1,klon
-            pyv1(i)=pyv1(i)+yv1(i)*pdtphys
-            pyu1(i)=pyu1(i)+yu1(i)*pdtphys
-         end do
-         do k=1,nbsrf
-             do i=1,klon
-               pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
-               ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
-            enddo
-         enddo
-
-      dtcum=dtcum+pdtphys
-768   continue
-
-      RETURN
-      END
Index: /LMDZ5/trunk/libf/phylmd/phystokenc.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phystokenc.F90	(revision 1454)
+++ /LMDZ5/trunk/libf/phylmd/phystokenc.F90	(revision 1454)
@@ -0,0 +1,362 @@
+SUBROUTINE phystokenc (nlon,nlev,pdtphys,rlon,rlat, &
+     pt,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
+     pfm_therm,pentr_therm, &
+     cdragh, pcoefh,yu1,yv1,ftsol,pctsrf, &
+     frac_impa,frac_nucl, &
+     pphis,paire,dtime,itap, &
+     psh, pda, pphi, pmp, pupwd, pdnwd)
+  
+  USE ioipsl
+  USE dimphy
+  USE infotrac, ONLY : nqtot
+  USE iophy
+  USE control_mod
+  
+  IMPLICIT NONE
+  
+!======================================================================
+! Auteur(s) FH
+! Objet: Ecriture des variables pour transport offline
+!
+!======================================================================
+  INCLUDE "dimensions.h"
+  INCLUDE "tracstoke.h"
+  INCLUDE "indicesol.h"
+  INCLUDE "iniprint.h"
+!======================================================================
+
+! Arguments:
+!
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: psh   ! humidite specifique
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pda
+  REAL,DIMENSION(klon,klev,klev), INTENT(IN):: pphi
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pmp
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pupwd ! saturated updraft mass flux
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pdnwd ! saturated downdraft mass flux
+
+!   EN ENTREE:
+!   ==========
+!
+!   divers:
+!   -------
+!
+  INTEGER nlon ! nombre de points horizontaux
+  INTEGER nlev ! nombre de couches verticales
+  REAL pdtphys ! pas d'integration pour la physique (seconde)
+  INTEGER itap
+  INTEGER, SAVE :: physid
+!$OMP THREADPRIVATE(physid)
+
+!   convection:
+!   -----------
+!
+  REAL pmfu(klon,klev)  ! flux de masse dans le panache montant
+  REAL pmfd(klon,klev)  ! flux de masse dans le panache descendant
+  REAL pen_u(klon,klev) ! flux entraine dans le panache montant
+  REAL pde_u(klon,klev) ! flux detraine dans le panache montant
+  REAL pen_d(klon,klev) ! flux entraine dans le panache descendant
+  REAL pde_d(klon,klev) ! flux detraine dans le panache descendant
+  REAL pt(klon,klev)
+  REAL,ALLOCATABLE,SAVE :: t(:,:)
+!$OMP THREADPRIVATE(t)
+!
+  REAL rlon(klon), rlat(klon), dtime
+  REAL zx_tmp_3d(iim,jjm+1,klev),zx_tmp_2d(iim,jjm+1)
+
+!   Couche limite:
+!   --------------
+!
+  REAL cdragh(klon)          ! cdrag
+  REAL pcoefh(klon,klev)     ! coeff melange CL
+  REAL pcoefh_buf(klon,klev) ! coeff melange CL + cdrag
+  REAL yv1(klon)
+  REAL yu1(klon),pphis(klon),paire(klon)
+
+!   Les Thermiques : (Abderr 25 11 02)
+!   ---------------
+  REAL, INTENT(IN) ::  pfm_therm(klon,klev+1)
+  REAL pentr_therm(klon,klev)
+  
+  REAL,ALLOCATABLE,SAVE :: entr_therm(:,:)
+  REAL,ALLOCATABLE,SAVE :: fm_therm(:,:)
+!$OMP THREADPRIVATE(entr_therm)
+!$OMP THREADPRIVATE(fm_therm)
+!
+!   Lessivage:
+!   ----------
+!
+  REAL frac_impa(klon,klev)
+  REAL frac_nucl(klon,klev)
+!
+! Arguments necessaires pour les sources et puits de traceur
+!
+  REAL ftsol(klon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
+  REAL pctsrf(klon,nbsrf) ! Pourcentage de sol f(nature du sol)
+!======================================================================
+!
+  INTEGER i, k, kk
+  REAL,ALLOCATABLE,SAVE :: mfu(:,:)  ! flux de masse dans le panache montant
+  REAL,ALLOCATABLE,SAVE :: mfd(:,:)  ! flux de masse dans le panache descendant
+  REAL,ALLOCATABLE,SAVE :: en_u(:,:) ! flux entraine dans le panache montant
+  REAL,ALLOCATABLE,SAVE :: de_u(:,:) ! flux detraine dans le panache montant
+  REAL,ALLOCATABLE,SAVE :: en_d(:,:) ! flux entraine dans le panache descendant
+  REAL,ALLOCATABLE,SAVE :: de_d(:,:) ! flux detraine dans le panache descendant
+  REAL,ALLOCATABLE,SAVE :: coefh(:,:) ! flux detraine dans le panache descendant
+  
+  REAL,ALLOCATABLE,SAVE :: pyu1(:)
+  REAL,ALLOCATABLE,SAVE :: pyv1(:)
+  REAL,ALLOCATABLE,SAVE :: pftsol(:,:)
+  REAL,ALLOCATABLE,SAVE :: ppsrf(:,:)
+!$OMP THREADPRIVATE(mfu,mfd,en_u,de_u,en_d,de_d,coefh)
+!$OMP THREADPRIVATE(pyu1,pyv1,pftsol,ppsrf)
+
+
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: sh  
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: da
+  REAL,DIMENSION(:,:,:), ALLOCATABLE,SAVE   :: phi
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: mp
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: upwd
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: dnwd
+  
+  REAL, SAVE :: dtcum
+  INTEGER, SAVE:: iadvtr=0
+!$OMP THREADPRIVATE(dtcum,iadvtr)
+  REAL zmin,zmax
+  LOGICAL ok_sync
+  CHARACTER(len=12) :: nvar
+!
+!======================================================================
+
+  iadvtr=iadvtr+1
+
+! Dans le meme vecteur on recombine le drag et les coeff d'echange
+  pcoefh_buf(:,1)      = cdragh(:)
+  pcoefh_buf(:,2:klev) = pcoefh(:,2:klev)
+  
+  ok_sync = .TRUE.
+
+! Initialization done only once
+!======================================================================
+  IF (iadvtr==1) THEN
+     ALLOCATE( t(klon,klev))
+     ALLOCATE( mfu(klon,klev))  
+     ALLOCATE( mfd(klon,klev))  
+     ALLOCATE( en_u(klon,klev)) 
+     ALLOCATE( de_u(klon,klev)) 
+     ALLOCATE( en_d(klon,klev)) 
+     ALLOCATE( de_d(klon,klev)) 
+     ALLOCATE( coefh(klon,klev)) 
+     ALLOCATE( entr_therm(klon,klev))
+     ALLOCATE( fm_therm(klon,klev))
+     ALLOCATE( pyu1(klon))
+     ALLOCATE( pyv1(klon))
+     ALLOCATE( pftsol(klon,nbsrf))
+     ALLOCATE( ppsrf(klon,nbsrf))
+     
+     ALLOCATE(sh(klon,klev)) 
+     ALLOCATE(da(klon,klev)) 
+     ALLOCATE(phi(klon,klev,klev)) 
+     ALLOCATE(mp(klon,klev)) 
+     ALLOCATE(upwd(klon,klev)) 
+     ALLOCATE(dnwd(klon,klev)) 
+
+     CALL initphysto('phystoke', dtime, dtime*istphy,dtime*istphy,physid)
+     
+     ! Write field phis and aire only once
+     CALL histwrite_phy(physid,"phis",itap,pphis)
+     CALL histwrite_phy(physid,"aire",itap,paire)
+     CALL histwrite_phy(physid,"longitudes",itap,rlon)
+     CALL histwrite_phy(physid,"latitudes",itap,rlat)
+
+  END IF
+  
+  
+! Set to zero cumulating fields
+!======================================================================
+  IF (MOD(iadvtr,istphy)==1.OR.istphy==1) THEN
+     WRITE(lunout,*)'reinitialisation des champs cumules a iadvtr=',iadvtr
+     mfu(:,:)=0.
+     mfd(:,:)=0.
+     en_u(:,:)=0.
+     de_u(:,:)=0.
+     en_d(:,:)=0.
+     de_d(:,:)=0.
+     coefh(:,:)=0.
+     t(:,:)=0.
+     fm_therm(:,:)=0.
+     entr_therm(:,:)=0.
+     pyv1(:)=0.
+     pyu1(:)=0.
+     pftsol(:,:)=0.
+     ppsrf(:,:)=0.
+     sh(:,:)=0.
+     da(:,:)=0.
+     phi(:,:,:)=0.
+     mp(:,:)=0.
+     upwd(:,:)=0.
+     dnwd(:,:)=0.
+     dtcum=0.
+  ENDIF
+  
+
+! Cumulate fields at each time step
+!======================================================================
+  DO k=1,klev
+     DO i=1,klon
+        mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
+        mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
+        en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
+        de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
+        en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
+        de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
+        coefh(i,k)=coefh(i,k)+pcoefh_buf(i,k)*pdtphys
+        t(i,k)=t(i,k)+pt(i,k)*pdtphys
+        fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
+        entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
+        sh(i,k) = sh(i,k) + psh(i,k)*pdtphys
+        da(i,k) = da(i,k) + pda(i,k)*pdtphys
+        mp(i,k) = mp(i,k) + pmp(i,k)*pdtphys
+        upwd(i,k) = upwd(i,k) + pupwd(i,k)*pdtphys
+        dnwd(i,k) = dnwd(i,k) + pdnwd(i,k)*pdtphys
+     ENDDO
+  ENDDO
+
+  DO kk=1,klev
+     DO k=1,klev
+        DO i=1,klon
+           phi(i,k,kk) = phi(i,k,kk) + pphi(i,k,kk)*pdtphys
+        END DO
+     END DO
+  END DO
+
+  DO i=1,klon
+     pyv1(i)=pyv1(i)+yv1(i)*pdtphys
+     pyu1(i)=pyu1(i)+yu1(i)*pdtphys
+  END DO
+  DO k=1,nbsrf
+     DO i=1,klon
+        pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
+        ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
+     ENDDO
+  ENDDO
+  
+! Add time step to cumulated time
+  dtcum=dtcum+pdtphys
+  
+
+! Write fields to file, if it is time to do so
+!======================================================================
+  IF(MOD(iadvtr,istphy)==0) THEN 
+
+     ! normalize with time period
+     DO k=1,klev
+        DO i=1,klon
+           mfu(i,k)=mfu(i,k)/dtcum
+           mfd(i,k)=mfd(i,k)/dtcum
+           en_u(i,k)=en_u(i,k)/dtcum
+           de_u(i,k)=de_u(i,k)/dtcum
+           en_d(i,k)=en_d(i,k)/dtcum
+           de_d(i,k)=de_d(i,k)/dtcum
+           coefh(i,k)=coefh(i,k)/dtcum
+           t(i,k)=t(i,k)/dtcum	
+           fm_therm(i,k)=fm_therm(i,k)/dtcum
+           entr_therm(i,k)=entr_therm(i,k)/dtcum
+           sh(i,k)=sh(i,k)/dtcum
+           da(i,k)=da(i,k)/dtcum
+           mp(i,k)=mp(i,k)/dtcum
+           upwd(i,k)=upwd(i,k)/dtcum
+           dnwd(i,k)=dnwd(i,k)/dtcum
+        ENDDO
+     ENDDO
+     DO kk=1,klev
+        DO k=1,klev
+           DO i=1,klon
+              phi(i,k,kk) = phi(i,k,kk)/dtcum
+           END DO
+        END DO
+     END DO
+     DO i=1,klon
+        pyv1(i)=pyv1(i)/dtcum
+        pyu1(i)=pyu1(i)/dtcum
+     END DO
+     DO k=1,nbsrf
+        DO i=1,klon
+           pftsol(i,k)=pftsol(i,k)/dtcum
+           ppsrf(i,k)=ppsrf(i,k)/dtcum
+        ENDDO
+     ENDDO
+
+     ! write fields
+     CALL histwrite_phy(physid,"t",itap,t)
+     CALL histwrite_phy(physid,"mfu",itap,mfu)
+     CALL histwrite_phy(physid,"mfd",itap,mfd)
+     CALL histwrite_phy(physid,"en_u",itap,en_u)
+     CALL histwrite_phy(physid,"de_u",itap,de_u)
+     CALL histwrite_phy(physid,"en_d",itap,en_d)
+     CALL histwrite_phy(physid,"de_d",itap,de_d)
+     CALL histwrite_phy(physid,"coefh",itap,coefh)	
+     CALL histwrite_phy(physid,"fm_th",itap,fm_therm)
+     CALL histwrite_phy(physid,"en_th",itap,entr_therm)
+     CALL histwrite_phy(physid,"frac_impa",itap,frac_impa)
+     CALL histwrite_phy(physid,"frac_nucl",itap,frac_nucl)
+     CALL histwrite_phy(physid,"pyu1",itap,pyu1)
+     CALL histwrite_phy(physid,"pyv1",itap,pyv1)
+     CALL histwrite_phy(physid,"ftsol1",itap,pftsol(:,1))
+     CALL histwrite_phy(physid,"ftsol2",itap,pftsol(:,2))
+     CALL histwrite_phy(physid,"ftsol3",itap,pftsol(:,3))
+     CALL histwrite_phy(physid,"ftsol4",itap,pftsol(:,4))
+     CALL histwrite_phy(physid,"psrf1",itap,ppsrf(:,1))
+     CALL histwrite_phy(physid,"psrf2",itap,ppsrf(:,2))
+     CALL histwrite_phy(physid,"psrf3",itap,ppsrf(:,3))
+     CALL histwrite_phy(physid,"psrf4",itap,ppsrf(:,4))
+     CALL histwrite_phy(physid,"sh",itap,sh)
+     CALL histwrite_phy(physid,"da",itap,da)
+     CALL histwrite_phy(physid,"mp",itap,mp)
+     CALL histwrite_phy(physid,"upwd",itap,upwd)
+     CALL histwrite_phy(physid,"dnwd",itap,dnwd)
+
+
+! phi
+     DO k=1,klev
+        IF (k<10) THEN
+           WRITE(nvar,'(i1)') k
+        ELSE IF (k<100) THEN
+           WRITE(nvar,'(i2)') k
+        ELSE
+           WRITE(nvar,'(i3)') k
+        END IF
+        nvar='phi_lev'//trim(nvar)
+        
+        CALL histwrite_phy(physid,nvar,itap,phi(:,:,k))
+     END DO
+     
+     ! Syncronize file
+!$OMP MASTER
+     IF (ok_sync) CALL histsync(physid)
+!$OMP END MASTER
+     
+     
+     ! Calculate min and max values for some fields (coefficients de lessivage)
+     zmin=1e33
+     zmax=-1e33
+     DO k=1,klev
+        DO i=1,klon
+           zmax=MAX(zmax,frac_nucl(i,k))
+           zmin=MIN(zmin,frac_nucl(i,k))
+        ENDDO
+     ENDDO
+     WRITE(lunout,*)'------ coefs de lessivage (min et max) --------'
+     WRITE(lunout,*)'facteur de nucleation ',zmin,zmax
+     zmin=1e33
+     zmax=-1e33
+     DO k=1,klev
+        DO i=1,klon
+           zmax=MAX(zmax,frac_impa(i,k))
+           zmin=MIN(zmin,frac_impa(i,k))
+        ENDDO
+     ENDDO
+     WRITE(lunout,*)'facteur d impaction ',zmin,zmax
+     
+  ENDIF ! IF(MOD(iadvtr,istphy)==0)
+
+END SUBROUTINE phystokenc
Index: /LMDZ5/trunk/libf/phylmd/phytrac.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phytrac.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/phytrac.F90	(revision 1454)
@@ -66,6 +66,6 @@
 !--------
   REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
-  REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! 
-  REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! 
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! variable not used
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! variable not used
   REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
   REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
@@ -118,5 +118,5 @@
 !--------------
 !
-  REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh ! coeff drag pour T et Q
+  REAL,DIMENSION(klon),INTENT(IN)      :: cdragh ! coeff drag pour T et Q
   REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh  ! coeff melange CL (m**2/s)
   REAL,DIMENSION(klon),INTENT(IN)      :: yu1    ! vents au premier niveau
@@ -213,6 +213,5 @@
      SELECT CASE(type_trac)
      CASE('lmdz')
-!IM ajout t_seri, pplay, sh    CALL traclmdz_init(pctsrf, ftsol, tr_seri, aerosol, lessivage)
-        CALL traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, aerosol, lessivage)
+        CALL traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
      CASE('inca')
         source(:,:)=0.
Index: /LMDZ5/trunk/libf/phylmd/read_map2D.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/read_map2D.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/read_map2D.F90	(revision 1454)
@@ -11,6 +11,6 @@
 
 ! Input arguments
-  CHARACTER(len=20), INTENT(IN) :: filename     ! name of file to read
-  CHARACTER(len=20), INTENT(IN) :: varname      ! name of variable in file
+  CHARACTER(len=*), INTENT(IN)  :: filename     ! name of file to read
+  CHARACTER(len=*), INTENT(IN)  :: varname      ! name of variable in file
   INTEGER, INTENT(IN)           :: timestep     ! actual timestep
   LOGICAL, INTENT(IN)           :: inverse      ! TRUE if latitude needs to be inversed
@@ -27,18 +27,21 @@
   REAL, DIMENSION(nbp_lon,nbp_lat) :: var_glo2D_tmp ! 2D global
   REAL, DIMENSION(klon_glo)        :: var_glo1D     ! 1D global
-
+  INCLUDE "iniprint.h"
 
 ! Read variable from file. Done by master process MPI and master thread OpenMP
   IF (is_mpi_root .AND. is_omp_root) THEN
-     ierr = NF90_OPEN (filename, NF90_NOWRITE, nid)
-     IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Problem in opening file '//filename,1)
+     ierr = NF90_OPEN(trim(filename), NF90_NOWRITE, nid)
+     IF (ierr /= NF90_NOERR) CALL write_err_mess('Problem in opening file')
 
-     ierr = NF90_INQ_VARID(nid, varname, nvarid)
-     IF (ierr /= NF90_NOERR) CALL abort_gcm(modname, 'The variable '//varname//' is absent in file',1)
+     ierr = NF90_INQ_VARID(nid, trim(varname), nvarid)
+     IF (ierr /= NF90_NOERR) CALL write_err_mess('The variable is absent in file')
      
      start=(/1,1,timestep/)
      count=(/nbp_lon,nbp_lat,1/)
      ierr = NF90_GET_VAR(nid, nvarid, var_glo2D,start,count)
-     IF (ierr /= NF90_NOERR) CALL abort_gcm(modname, 'Problem in reading varaiable '//varname,1)
+     IF (ierr /= NF90_NOERR) CALL write_err_mess('Problem in reading varaiable')
+
+     ierr = NF90_CLOSE(nid)
+     IF (ierr /= NF90_NOERR) CALL write_err_mess('Problem in closing file')
 
      ! Inverse latitude order
@@ -53,4 +56,7 @@
      CALL grid2Dto1D_glo(var_glo2D,var_glo1D)
 
+     WRITE(lunout,*) 'in read_map2D, filename = ', trim(filename)
+     WRITE(lunout,*) 'in read_map2D, varname  = ', trim(varname)
+     WRITE(lunout,*) 'in read_map2D, timestep = ', timestep
   ENDIF
 
@@ -58,3 +64,17 @@
   CALL scatter(var_glo1D, varout)
 
+  CONTAINS
+    SUBROUTINE write_err_mess(err_mess)
+
+      CHARACTER(len=*), INTENT(IN) :: err_mess
+      INCLUDE "iniprint.h"
+      
+      WRITE(lunout,*) 'Error in read_map2D, filename = ', trim(filename)
+      WRITE(lunout,*) 'Error in read_map2D, varname  = ', trim(varname)
+      WRITE(lunout,*) 'Error in read_map2D, timestep = ', timestep
+
+      CALL abort_gcm(modname, err_mess, 1)
+
+    END SUBROUTINE write_err_mess
+
 END SUBROUTINE read_map2D
Index: /LMDZ5/trunk/libf/phylmd/surf_land_orchidee_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/surf_land_orchidee_mod.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/surf_land_orchidee_mod.F90	(revision 1454)
@@ -43,6 +43,5 @@
     USE mod_surf_para
     USE mod_synchro_omp
-   
-USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, fco2_land_inst, fco2_lu_inst
+    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl
 
 !    
@@ -138,6 +137,4 @@
     INTEGER                                   :: error
     REAL, DIMENSION(klon)                     :: swdown_vrai
-    REAL, DIMENSION(klon)                     :: fco2_land_comp  ! sur grille compresse
-    REAL, DIMENSION(klon)                     :: fco2_lu_comp    ! sur grille compresse
     CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
     CHARACTER (len = 80)                      :: abort_message
@@ -348,14 +345,9 @@
        ENDIF
 !
-! Allocate variables needed for carbon_cycle_mod
+! carbon_cycle_cpl not possible with this interface and version of ORHCHIDEE
 !
        IF (carbon_cycle_cpl) THEN
-          IF (.NOT. ALLOCATED(fco2_land_inst)) THEN
-             ALLOCATE(fco2_land_inst(klon),stat=error)
-             IF (error /= 0)  CALL abort_gcm(modname,'Pb in allocation fco2_land_inst',1)
-             
-             ALLOCATE(fco2_lu_inst(klon),stat=error)
-             IF(error /=0) CALL abort_gcm(modname,'Pb in allocation fco2_lu_inst',1)
-          END IF
+          abort_message='carbon_cycle_cpl not yet possible with this interface of ORCHIDEE'
+          CALL abort_gcm(modname,abort_message,1)
        END IF
        
@@ -464,22 +456,4 @@
 
     
-! JG : TEMPORAIRE!!!! Les variables fco2_land_comp et fco2_lu_comp seront plus tard en sortie d'ORCHIDEE
-!      ici mis a valeur quelquonque pour test. Ces variables sont sur la grille compresse avec uniquement des points de terres
-
-    fco2_land_comp(:) = 1.
-    fco2_lu_comp(:)   = 10.
-
-! Decompress variables for the module carbon_cycle_mod
-    IF (carbon_cycle_cpl) THEN
-       fco2_land_inst(:)=0.
-       fco2_lu_inst(:)=0.
-       
-       DO igrid = 1, knon
-          ireal = knindex(igrid)
-          fco2_land_inst(ireal) = fco2_land_comp(igrid)
-          fco2_lu_inst(ireal)   = fco2_lu_comp(igrid)
-       END DO
-    END IF
-
   END SUBROUTINE surf_land_orchidee
 !
Index: /LMDZ5/trunk/libf/phylmd/surf_land_orchidee_noopenmp_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/surf_land_orchidee_noopenmp_mod.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/surf_land_orchidee_noopenmp_mod.F90	(revision 1454)
@@ -138,7 +138,7 @@
     INTEGER                                   :: ij, jj, igrid, ireal, index
     INTEGER                                   :: error
+    INTEGER, SAVE                             :: nb_fields_cpl ! number of fields for the climate-carbon coupling (between ATM and ORCHIDEE). 
+    REAL, SAVE, ALLOCATABLE, DIMENSION(:,:)   :: fields_cpl    ! Fluxes for the climate-carbon coupling
     REAL, DIMENSION(klon)                     :: swdown_vrai
-    REAL, DIMENSION(klon)                     :: fco2_land_comp  ! sur grille compresse
-    REAL, DIMENSION(klon)                     :: fco2_lu_comp    ! sur grille compresse
     CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
     CHARACTER (len = 80)                      :: abort_message
@@ -210,4 +210,12 @@
   
     IF (debut) THEN
+! Test de coherence
+#ifndef ORCH_NEW
+       ! Compilation avec orchidee nouvelle version necessaire avec carbon_cycle_cpl=y
+       IF (carbon_cycle_cpl) THEN
+          abort_message='You must define preprossing key ORCH_NEW when running carbon_cycle_cpl=y'
+          CALL abort_gcm(modname,abort_message,1)
+       END IF
+#endif
        ALLOCATE(ktindex(knon))
        IF ( .NOT. ALLOCATED(albedo_keep)) THEN
@@ -342,14 +350,21 @@
 !
 ! Allocate variables needed for carbon_cycle_mod
-!
+       IF ( carbon_cycle_cpl ) THEN
+          nb_fields_cpl=2
+       ELSE
+          nb_fields_cpl=1
+       END IF
+
+
        IF (carbon_cycle_cpl) THEN
-          IF (.NOT. ALLOCATED(fco2_land_inst)) THEN
-             ALLOCATE(fco2_land_inst(klon),stat=error)
-             IF (error /= 0)  CALL abort_gcm(modname,'Pb in allocation fco2_land_inst',1)
-             
-             ALLOCATE(fco2_lu_inst(klon),stat=error)
-             IF(error /=0) CALL abort_gcm(modname,'Pb in allocation fco2_lu_inst',1)
-          END IF
+          ALLOCATE(fco2_land_inst(klon),stat=error)
+          IF (error /= 0)  CALL abort_gcm(modname,'Pb in allocation fco2_land_inst',1)
+          
+          ALLOCATE(fco2_lu_inst(klon),stat=error)
+          IF(error /=0) CALL abort_gcm(modname,'Pb in allocation fco2_lu_inst',1)
        END IF
+
+       ALLOCATE(fields_cpl(klon,nb_fields_cpl), stat = error)
+       IF (error /= 0) CALL abort_gcm(modname,'Pb in allocation fields_cpl',1)
 
     ENDIF                          ! (fin debut) 
@@ -406,8 +421,13 @@
                evap, fluxsens, fluxlat, coastalflow, riverflow, &
                tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
-               lon_scat, lat_scat, q2m, t2m)
+               lon_scat, lat_scat, q2m, t2m &
+#ifdef ORCH_NEW
+               , nb_fields_cpl, fields_cpl)
+#else
+               )
+#endif
 
 #else          
-          ! Interface for ORCHIDEE version 1.9 or later(1.9.2, 1.9.3, 1.9.4) compiled in parallel mode(with preprocessing flag CPP_MPI)
+          ! Interface for ORCHIDEE version 1.9 or later(1.9.2, 1.9.3, 1.9.4, 1.9.5) compiled in parallel mode(with preprocessing flag CPP_MPI)
           CALL intersurf_main (itime+itau_phy-1, iim, jjm+1, offset, knon, ktindex, & 
                orch_comm, dtime, lrestart_read, lrestart_write, lalo, &
@@ -418,5 +438,10 @@
                evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
                tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
-               lon_scat, lat_scat, q2m, t2m)
+               lon_scat, lat_scat, q2m, t2m &
+#ifdef ORCH_NEW
+               , nb_fields_cpl, fields_cpl(1:knon,:))
+#else
+               )
+#endif
 #endif
           
@@ -431,5 +456,4 @@
 
     IF (knon /=0) THEN
-    
 #ifndef CPP_MPI
        ! Interface for ORCHIDEE compiled in sequential mode(without preprocessing flag CPP_MPI)
@@ -442,6 +466,10 @@
             evap, fluxsens, fluxlat, coastalflow, riverflow, &
             tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
-            lon_scat, lat_scat, q2m, t2m)
-       
+            lon_scat, lat_scat, q2m, t2m &
+#ifdef ORCH_NEW
+            , nb_fields_cpl, fields_cpl)
+#else
+            )
+#endif
 #else
        ! Interface for ORCHIDEE version 1.9 or later compiled in parallel mode(with preprocessing flag CPP_MPI)
@@ -454,7 +482,11 @@
             evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
             tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
-            lon_scat, lat_scat, q2m, t2m)
-#endif
-       
+            lon_scat, lat_scat, q2m, t2m &
+#ifdef ORCH_NEW
+            , nb_fields_cpl, fields_cpl(1:knon,:))
+#else
+            )
+#endif
+#endif
     ENDIF
 
@@ -478,11 +510,4 @@
 
     IF (debut) lrestart_read = .FALSE.
-
-
-! JG : TEMPORAIRE!!!! Les variables fco2_land_comp et fco2_lu_comp seront plus tard en sortie d'ORCHIDEE
-!      ici mis a valeur quelquonque pour test. Ces variables sont sur la grille compresse avec uniquement des points de terres
-
-    fco2_land_comp(:) = 1.
-    fco2_lu_comp(:)   = 10.
 
 ! Decompress variables for the module carbon_cycle_mod
@@ -493,6 +518,6 @@
        DO igrid = 1, knon
           ireal = knindex(igrid)
-          fco2_land_inst(ireal) = fco2_land_comp(igrid)
-          fco2_lu_inst(ireal)   = fco2_lu_comp(igrid)
+          fco2_land_inst(ireal) = fields_cpl(igrid,1)
+          fco2_lu_inst(ireal)   = fields_cpl(igrid,2)
        END DO
     END IF
Index: /LMDZ5/trunk/libf/phylmd/traclmdz_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/traclmdz_mod.F90	(revision 1453)
+++ /LMDZ5/trunk/libf/phylmd/traclmdz_mod.F90	(revision 1454)
@@ -84,5 +84,5 @@
 
 
-  SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, aerosol, lessivage)
+  SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
     ! This subroutine allocates and initialize module variables and control variables.
     ! Initialization of the tracers should be done here only for those not found in the restart file.
@@ -104,4 +104,5 @@
     REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
     REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
+    REAL,INTENT(IN)                        :: pdtphys ! Pas d'integration pour la physique (seconde)  
 
 ! Output variables
@@ -226,5 +227,5 @@
 ! ----------------------------------------------
     IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
-       CALL carbon_cycle_init(tr_seri, aerosol, radio)
+       CALL carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
     END IF
 
@@ -312,5 +313,5 @@
 !--------------
 !
-    REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh     ! coeff drag pour T et Q
+    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
     REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
     REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
@@ -546,5 +547,5 @@
 !======================================================================
     IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
-       CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
+       CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
     END IF
 
