Index: LMDZ5/trunk/libf/dyn3d/academic.h
===================================================================
--- LMDZ5/trunk/libf/dyn3d/academic.h	(revision 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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 1434)
+++ 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
