Index: LMDZ4/trunk/libf/dyn3dpar/addfi_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/addfi_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/addfi_p.F	(revision 1146)
@@ -2,8 +2,9 @@
 ! $Header$
 !
-      SUBROUTINE addfi_p(nq, pdt, leapf, forward,
+      SUBROUTINE addfi_p(pdt, leapf, forward,
      S          pucov, pvcov, pteta, pq   , pps ,
      S          pdufi, pdvfi, pdhfi,pdqfi, pdpfi  )
       USE parallel
+      USE infotrac, ONLY : nqtot
       IMPLICIT NONE
 c
@@ -53,13 +54,11 @@
 c    -----------
 c
-      INTEGER nq
-
       REAL pdt
 c
       REAL pvcov(ip1jm,llm),pucov(ip1jmp1,llm)
-      REAL pteta(ip1jmp1,llm),pq(ip1jmp1,llm,nq),pps(ip1jmp1)
+      REAL pteta(ip1jmp1,llm),pq(ip1jmp1,llm,nqtot),pps(ip1jmp1)
 c
       REAL pdvfi(ip1jm,llm),pdufi(ip1jmp1,llm)
-      REAL pdqfi(ip1jmp1,llm,nq),pdhfi(ip1jmp1,llm),pdpfi(ip1jmp1)
+      REAL pdqfi(ip1jmp1,llm,nqtot),pdhfi(ip1jmp1,llm),pdpfi(ip1jmp1)
 c
       LOGICAL leapf,forward
@@ -166,5 +165,5 @@
       ENDDO
 
-      DO iq = 3, nq
+      DO iq = 3, nqtot
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO k = 1,llm
@@ -208,5 +207,5 @@
 
       if (pole_nord) then
-        DO iq = 1, nq
+        DO iq = 1, nqtot
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
           DO  k    = 1, llm
@@ -225,5 +224,5 @@
       
       if (pole_sud) then
-        DO iq = 1, nq
+        DO iq = 1, nqtot
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
           DO  k    = 1, llm
Index: LMDZ4/trunk/libf/dyn3dpar/advtrac.h
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/advtrac.h	(revision 1122)
+++ 	(revision )
@@ -1,34 +1,0 @@
-!
-! $Header$
-!
-!-----------------------------------------------------------------------
-! INCLUDE 'advtrac.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)
-!
-
-      COMMON/advtr/iadv, hadv, vadv, tnom, tname, ttext, niadv,         &
-     &   nbtrac, nprath, mmt_adj, hadv_flg, vadv_flg, conv_flg,         &
-     &   pbl_flg, tracnam
-      INTEGER iadv(nqmx) ! indice schema de transport 
-      INTEGER hadv(nqmx) ! indice schema transport horizontal 
-      INTEGER vadv(nqmx) ! indice schema transport vertical 
-      INTEGER niadv(nqmx) ! equivalent dyn / physique
-      character (len=8) :: tnom(nqmx) ! nom court du traceur
-      character (len=10) :: tname(nqmx) ! nom du traceur pour restart
-      character (len=13) :: ttext(nqmx) ! nom long du traceur pour sorties
-
-
-      integer nbtrac
-      integer nprath
-      real    mmt_adj(iim+1,jjm+1,llm, 1) 
-      integer hadv_flg(nqmx)
-      integer vadv_flg(nqmx)
-      integer conv_flg(nqmx-2) 
-      integer pbl_flg(nqmx-2)
-      character (len=8) :: tracnam(nqmx-2) 
-!-----------------------------------------------------------------------
Index: LMDZ4/trunk/libf/dyn3dpar/advtrac_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/advtrac_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/advtrac_p.F	(revision 1146)
@@ -22,4 +22,5 @@
       USE Vampir
       USE times
+      USE infotrac
       IMPLICIT NONE
 c
@@ -35,5 +36,4 @@
 #include "ener.h"
 #include "description.h"
-#include "advtrac.h"
 
 c-------------------------------------------------------------------
@@ -46,5 +46,5 @@
       INTEGER iapptrac
       REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
-      REAL q(ip1jmp1,llm,nqmx),masse(ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
       REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
       REAL pk(ip1jmp1,llm)
@@ -59,7 +59,4 @@
       REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 
       REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
-      real cpuadv(nqmx)
-      common/cpuadv/cpuadv
-
       INTEGER iadvtr
       INTEGER ij,l,iq,iiq
@@ -76,5 +73,5 @@
       REAL psppm(iim,jjp1) ! pression  au sol
       REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
-      REAL qppm(iim*jjp1,llm,nqmx)
+      REAL qppm(iim*jjp1,llm,nqtot)
       REAL fluxwppm(iim,jjp1,llm)
       REAL apppm(llmp1), bpppm(llmp1)
@@ -88,5 +85,5 @@
       REAL,SAVE :: teta_tmp(ip1jmp1,llm)
       REAL,SAVE :: pk_tmp(ip1jmp1,llm)
-      
+
       ijb_u=ij_begin
       ije_u=ij_end
@@ -196,5 +193,5 @@
       call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,
      *                             jj_Nb_vanleer,1,1,Request_vanleer)
-      do j=1,nqmx
+      do j=1,nqtot
         call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
      *                             jj_nb_vanleer,0,0,Request_vanleer)
@@ -279,5 +276,5 @@
 c     Appel des sous programmes d'advection
 c-----------------------------------------------------------
-      do iq=1,nqmx
+      do iq=1,nqtot
 c        call clock(t_initial)
         if(iadv(iq) == 0) cycle 
@@ -479,5 +476,5 @@
 c$OMP END MASTER
 
-	do j=1,nqmx
+	do j=1,nqtot
           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
      *                             jj_nb_caldyn,0,0,Request_vanleer)
Index: LMDZ4/trunk/libf/dyn3dpar/caladvtrac_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/caladvtrac_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/caladvtrac_p.F	(revision 1146)
@@ -8,4 +8,5 @@
      *                   flxw, pk, iapptrac)
       USE parallel
+      USE infotrac
 c
       IMPLICIT NONE
@@ -25,10 +26,9 @@
 #include "comconst.h"
 #include "control.h"
-#include "advtrac.h"
 
 c   Arguments:
 c   ----------
       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm),masse(ip1jmp1,llm)
-      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqmx),dq( ip1jmp1,llm,2 )
+      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot),dq( ip1jmp1,llm,2 )
       REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm)
       REAL               :: flxw(ip1jmp1,llm)
@@ -48,5 +48,4 @@
       REAL finmasse(ip1jmp1,llm), dtvrtrac
       
-
 cc
 c
Index: LMDZ4/trunk/libf/dyn3dpar/calfis_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/calfis_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/calfis_p.F	(revision 1146)
@@ -4,6 +4,5 @@
 C
 C
-      SUBROUTINE calfis_p(nq,
-     $                  lafin,
+      SUBROUTINE calfis_p(lafin,
      $                  rdayvrai,
      $                  heure,
@@ -40,4 +39,6 @@
       USE Times
       USE IOPHY
+      USE infotrac
+
       IMPLICIT NONE
 c=======================================================================
@@ -98,7 +99,6 @@
 #include "paramet.h"
 #include "temps.h"
-#include "advtrac.h"
-
-      INTEGER ngridmx,nq
+
+      INTEGER ngridmx
       PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
 
@@ -119,5 +119,5 @@
       REAL pteta(iip1,jjp1,llm)
       REAL pmasse(iip1,jjp1,llm)
-      REAL pq(iip1,jjp1,llm,nqmx)
+      REAL pq(iip1,jjp1,llm,nqtot)
       REAL pphis(iip1,jjp1)
       REAL pphi(iip1,jjp1,llm)
@@ -126,5 +126,5 @@
       REAL pducov(iip1,jjp1,llm)
       REAL pdteta(iip1,jjp1,llm)
-      REAL pdq(iip1,jjp1,llm,nqmx)
+      REAL pdq(iip1,jjp1,llm,nqtot)
 c
       REAL pps(iip1,jjp1)
@@ -135,5 +135,5 @@
       REAL pdufi(iip1,jjp1,llm)
       REAL pdhfi(iip1,jjp1,llm)
-      REAL pdqfi(iip1,jjp1,llm,nqmx)
+      REAL pdqfi(iip1,jjp1,llm,nqtot)
       REAL pdpsfi(iip1,jjp1)
 
@@ -253,10 +253,10 @@
       ALLOCATE(zphi(klon,llm),zphis(klon))
       ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
-      ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqmx))
+      ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
       ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
       ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
 c      ALLOCATE(pvervel(klon,llm))
       ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
-      ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqmx))
+      ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
       ALLOCATE(zdpsrf(klon))
       ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
@@ -335,5 +335,5 @@
 c
 
-      DO iq=1,nq
+      DO iq=1,nqtot
          iiq=niadv(iq)
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
@@ -369,7 +369,6 @@
       CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
 
-c$OMP MASTER
       CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
-c$OMP END MASTER
+
 c$OMP BARRIER
 
@@ -527,5 +526,5 @@
 cc$OMP  PARALLEL DEFAULT(NONE) 
 cc$OMP+ PRIVATE(i,l,offset,iq)
-cc$OMP+ SHARED(klon_omp_nb,nq,klon_omp_begin,
+cc$OMP+ SHARED(klon_omp_nb,nqtot,klon_omp_begin,
 cc$OMP+        debut,lafin,rdayvrai,heure,dtphys,zplev,zplay,
 cc$OMP+        zphi,zphis,presnivs,clesphy0,zufi,zvfi,ztfi,
@@ -549,10 +548,10 @@
         allocate(zvfi_omp(klon,llm))
         allocate(ztfi_omp(klon,llm))
-        allocate(zqfi_omp(klon,llm,nq))
+        allocate(zqfi_omp(klon,llm,nqtot))
 c        allocate(pvervel_omp(klon,llm))
         allocate(zdufi_omp(klon,llm))
         allocate(zdvfi_omp(klon,llm))
         allocate(zdtfi_omp(klon,llm))
-        allocate(zdqfi_omp(klon,llm,nq))
+        allocate(zdqfi_omp(klon,llm,nqtot))
         allocate(zdpsrf_omp(klon))
         allocate(flxwfi_omp(klon,llm))
@@ -609,5 +608,5 @@
       enddo
 	
-      do iq=1,nq
+      do iq=1,nqtot
         do l=1,llm
           do i=1,klon
@@ -641,5 +640,5 @@
       enddo
 	
-      do iq=1,nq
+      do iq=1,nqtot
         do l=1,llm
           do i=1,klon
@@ -664,5 +663,4 @@
       CALL physiq (klon,
      .             llm,
-     .             nq,
      .             debut,
      .             lafin,
@@ -743,5 +741,5 @@
       enddo
 	
-      do iq=1,nq
+      do iq=1,nqtot
         do l=1,llm
           do i=1,klon
@@ -775,5 +773,5 @@
       enddo
 	
-      do iq=1,nq
+      do iq=1,nqtot
         do l=1,llm
           do i=1,klon
@@ -896,7 +894,5 @@
 c  tendance sur la pression :
 c  -----------------------------------
-c$OMP MASTER
       CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
-c$OMP END MASTER
 c
 c   62. enthalpie potentielle
@@ -937,5 +933,5 @@
 c   ---------------------
 
-      DO iq=1,nqmx
+      DO iq=1,nqtot
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l=1,llm
@@ -976,5 +972,5 @@
 C
 
-      DO iq=1,nq
+      DO iq=1,nqtot
          iiq=niadv(iq)
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
Index: LMDZ4/trunk/libf/dyn3dpar/comdiss.h
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/comdiss.h	(revision 1122)
+++ 	(revision )
@@ -1,15 +1,0 @@
-!
-! $Header$
-!
-c-----------------------------------------------------------------------
-c INCLUDE 'comdiss.h'
-
-      INTEGER niter
-      REAL tetau(llm),tetah(llm),cdivu,cdivh,crot
-      LOGICAL dissipst
-
-      COMMON/comdissi/niter
-      COMMON/comdissr/tetau,tetah,cdivu,cdivh,crot
-      COMMON/comdissl/dissipst
-
-c-----------------------------------------------------------------------
Index: LMDZ4/trunk/libf/dyn3dpar/comdissip.h
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/comdissip.h	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/comdissip.h	(revision 1146)
@@ -2,15 +2,14 @@
 ! $Header$
 !
-c-----------------------------------------------------------------------
-c INCLUDE dissip.h
+!-----------------------------------------------------------------------
+! INCLUDE comdissip.h
 
-      COMMON/comdissip/
-     $    lstardis,niterdis,coefdis,tetavel,tetatemp,gamdissip
+      COMMON/comdissip/                                                 &
+     &    niterdis,coefdis,tetavel,tetatemp,gamdissip
 
 
-      LOGICAL lstardis
       INTEGER niterdis
 
       REAL tetavel,tetatemp,coefdis,gamdissip
 
-c-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
Index: LMDZ4/trunk/libf/dyn3dpar/comgeom.h
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/comgeom.h	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/comgeom.h	(revision 1146)
@@ -2,32 +2,32 @@
 ! $Header$
 !
-*CDK comgeom
-      COMMON/comgeom/
-     1 cu(ip1jmp1),cv(ip1jm),unscu2(ip1jmp1),unscv2(ip1jm),
-     2 aire(ip1jmp1),airesurg(ip1jmp1),aireu(ip1jmp1),
-     3 airev(ip1jm),unsaire(ip1jmp1),apoln,apols,
-     4 unsairez(ip1jm),airuscv2(ip1jm),airvscu2(ip1jm),
-     5 aireij1(ip1jmp1),aireij2(ip1jmp1),aireij3(ip1jmp1),
-     6 aireij4(ip1jmp1),alpha1(ip1jmp1),alpha2(ip1jmp1),
-     7 alpha3(ip1jmp1),alpha4(ip1jmp1),alpha1p2(ip1jmp1),
-     8 alpha1p4(ip1jmp1),alpha2p3(ip1jmp1),alpha3p4(ip1jmp1),
-     9 fext(ip1jm),constang(ip1jmp1),rlatu(jjp1),rlatv(jjm),
-     1 rlonu(iip1),rlonv(iip1),cuvsurcv(ip1jm),cvsurcuv(ip1jm),
-     1 cvusurcu(ip1jmp1),cusurcvu(ip1jmp1),cuvscvgam1(ip1jm),
-     2 cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1),
-     3 cvuscugam2(ip1jmp1),cvscuvgam(ip1jm),cuscvugam(ip1jmp1),
-     4 unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2,
-     5 unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1),unsairz_gam(ip1jm),
-     6 aivscu2gam(ip1jm),aiuscv2gam(ip1jm),xprimu(iip1),xprimv(iip1)
+!CDK comgeom
+      COMMON/comgeom/                                                   &
+     & cu(ip1jmp1),cv(ip1jm),unscu2(ip1jmp1),unscv2(ip1jm),             &
+     & aire(ip1jmp1),airesurg(ip1jmp1),aireu(ip1jmp1),                  &
+     & airev(ip1jm),unsaire(ip1jmp1),apoln,apols,                       &
+     & unsairez(ip1jm),airuscv2(ip1jm),airvscu2(ip1jm),                 &
+     & aireij1(ip1jmp1),aireij2(ip1jmp1),aireij3(ip1jmp1),              &
+     & aireij4(ip1jmp1),alpha1(ip1jmp1),alpha2(ip1jmp1),                &
+     & alpha3(ip1jmp1),alpha4(ip1jmp1),alpha1p2(ip1jmp1),               &
+     & alpha1p4(ip1jmp1),alpha2p3(ip1jmp1),alpha3p4(ip1jmp1),           &
+     & fext(ip1jm),constang(ip1jmp1),rlatu(jjp1),rlatv(jjm),            &
+     & rlonu(iip1),rlonv(iip1),cuvsurcv(ip1jm),cvsurcuv(ip1jm),         &
+     & cvusurcu(ip1jmp1),cusurcvu(ip1jmp1),cuvscvgam1(ip1jm),           &
+     & cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1),                           &
+     & cvuscugam2(ip1jmp1),cvscuvgam(ip1jm),cuscvugam(ip1jmp1),         &
+     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2,                 &
+     & unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1),unsairz_gam(ip1jm),    &
+     & aivscu2gam(ip1jm),aiuscv2gam(ip1jm),xprimu(iip1),xprimv(iip1)
 
-c
-        REAL 
-     1 cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,unsaire,apoln     ,
-     2 apols,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4,
-     3 alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 ,
-     4 fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     ,
-     5 cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1,unsapolnga2 
-     6 ,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2,unsairz_gam    ,
-     7 aivscu2gam ,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu,cusurcvu,xprimu
-     8 , xprimv
-c
+!
+        REAL                                                            &
+     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,unsaire,apoln     ,&
+     & apols,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4,&
+     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 ,&
+     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     ,&
+     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1,unsapolnga2&
+     & ,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2,unsairz_gam    ,&
+     & aivscu2gam ,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu,cusurcvu,xprimu&
+     & , xprimv
+!
Index: LMDZ4/trunk/libf/dyn3dpar/conf_gcm.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/conf_gcm.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/conf_gcm.F	(revision 1146)
@@ -6,5 +6,10 @@
       SUBROUTINE conf_gcm( tapedef, etatinit, clesphy0 )
 c
+#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 misc_mod
       use mod_filtre_fft, ONLY : use_filtre_fft
@@ -109,4 +114,10 @@
 c  Parametres de controle du run:
 c-----------------------------------------------------------------------
+!Config  Key  = planet_type
+!Config  Desc = planet type ("earth", "mars", "venus", ...)
+!Config  Def  = earth
+!Config  Help = this flag sets the type of atymosphere that is considered
+      planet_type="earth"
+      CALL getin('planet_type',planet_type)
 
 !Config  Key  = dayref
@@ -189,4 +200,11 @@
        CALL getin('periodav',periodav)
 
+!Config  Key  = output_grads_dyn
+!Config  Desc = output dynamics diagnostics in 'dyn.dat' file
+!Config  Def  = n
+!Config  Help = output dynamics diagnostics in Grads-readable 'dyn.dat' file
+       output_grads_dyn=.false.
+       CALL getin('output_grads_dyn',output_grads_dyn)
+
 !Config  Key  = idissip
 !Config  Desc = periode de la dissipation 
@@ -284,4 +302,12 @@
 c    ...............................................................
 
+!Config  Key  =  read_start
+!Config  Desc = Initialize model using a 'start.nc' file
+!Config  Def  = y
+!Config  Help = y: intialize dynamical fields using a 'start.nc' file
+!               n: fields are initialized by 'iniacademic' routine
+       read_start= .true.
+       CALL getin('read_start',read_start)
+
 !Config  Key  = iflag_phys
 !Config  Desc = Avec ls physique 
@@ -341,6 +367,6 @@
 c
       IF( ABS(clat - clatt).GE. 0.001 )  THEN
-        PRINT *,' La valeur de clat passee par run.def est differente de 
-     *  celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de clat passee par run.def',
+     &    ' est differente de celle lue sur le fichier  start '
         STOP
       ENDIF
@@ -356,6 +382,6 @@
 
       IF( ABS(grossismx - grossismxx).GE. 0.001 )  THEN
-        PRINT *,' La valeur de grossismx passee par run.def est differente  
-     *  de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de grossismx passee par ',
+     &  'run.def est differente de celle lue sur le fichier  start '
         STOP
       ENDIF
@@ -370,11 +396,12 @@
 
       IF( ABS(grossismy - grossismyy).GE. 0.001 )  THEN
-        PRINT *,' La valeur de grossismy passee par run.def est differen
-     * te de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de grossismy passee par ',
+     & 'run.def est differente de celle lue sur le fichier  start '
         STOP
       ENDIF
       
       IF( grossismx.LT.1. )  THEN
-        PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
+        write(lunout,*)
+     &       'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
          STOP
       ELSE
@@ -384,5 +411,6 @@
 
       IF( grossismy.LT.1. )  THEN
-        PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
+        write(lunout,*)
+     &       'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
          STOP
       ELSE
@@ -390,5 +418,5 @@
       ENDIF
 
-      PRINT *,' alphax alphay defrun ',alphax,alphay
+      write(lunout,*)'conf_gcm: alphax alphay',alphax,alphay
 c
 c    alphax et alphay sont les anciennes formulat. des grossissements
@@ -405,17 +433,17 @@
 
       IF( .NOT.fxyhypb )  THEN
-           IF( fxyhypbb )     THEN
-              PRINT *,' ********  PBS DANS  DEFRUN  ******** '
-              PRINT *,' *** fxyhypb lu sur le fichier start est F ',
-     *       'alors  qu il est  T  sur  run.def  ***'
+         IF( fxyhypbb )     THEN
+            write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
+            write(lunout,*)' *** fxyhypb lu sur le fichier start est ',
+     *       'F alors  qu il est  T  sur  run.def  ***'
               STOP
-           ENDIF
+         ENDIF
       ELSE
-           IF( .NOT.fxyhypbb )   THEN
-              PRINT *,' ********  PBS DANS  DEFRUN  ******** '
-              PRINT *,' ***  fxyhypb lu sur le fichier start est T ',
-     *        'alors  qu il est  F  sur  run.def  ****  '
+         IF( .NOT.fxyhypbb )   THEN
+            write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
+            write(lunout,*)' ***  fxyhypb lu sur le fichier start est ',
+     *        'T alors  qu il est  F  sur  run.def  ****  '
               STOP
-           ENDIF
+         ENDIF
       ENDIF
 c
@@ -430,6 +458,6 @@
       IF( fxyhypb )  THEN
        IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
-        PRINT *,' La valeur de dzoomx passee par run.def est differente 
-     *  de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de dzoomx passee par ',
+     *  'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -446,6 +474,6 @@
       IF( fxyhypb )  THEN
        IF( ABS(dzoomy - dzoomyy).GE. 0.001 )  THEN
-        PRINT *,' La valeur de dzoomy passee par run.def est differente 
-     * de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de dzoomy passee par ',
+     * 'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -461,6 +489,6 @@
       IF( fxyhypb )  THEN
        IF( ABS(taux - tauxx).GE. 0.001 )  THEN
-        PRINT *,' La valeur de taux passee par run.def est differente
-     * de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de taux passee par ',
+     * 'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -476,6 +504,6 @@
       IF( fxyhypb )  THEN
        IF( ABS(tauy - tauyy).GE. 0.001 )  THEN
-        PRINT *,' La valeur de tauy passee par run.def est differente
-     * de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de tauy passee par ',
+     * 'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -495,19 +523,19 @@
 
         IF( .NOT.ysinus )  THEN
-           IF( ysinuss )     THEN
-              PRINT *,' ********  PBS DANS  DEFRUN  ******** '
-              PRINT *,' *** ysinus lu sur le fichier start est F ',
-     *       'alors  qu il est  T  sur  run.def  ***'
+          IF( ysinuss )     THEN
+            write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
+            write(lunout,*)' *** ysinus lu sur le fichier start est F',
+     *       ' alors  qu il est  T  sur  run.def  ***'
+            STOP
+          ENDIF
+        ELSE
+          IF( .NOT.ysinuss )   THEN
+            write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
+            write(lunout,*)' *** ysinus lu sur le fichier start est T',
+     *        ' alors  qu il est  F  sur  run.def  ****  '
               STOP
-           ENDIF
-        ELSE
-           IF( .NOT.ysinuss )   THEN
-              PRINT *,' ********  PBS DANS  DEFRUN  ******** '
-              PRINT *,' ***  ysinus lu sur le fichier start est T ',
-     *        'alors  qu il est  F  sur  run.def  ****  '
-              STOP
-           ENDIF
+          ENDIF
         ENDIF
-      ENDIF
+      ENDIF ! of IF( .NOT.fxyhypb  )
 c
 !Config  Key  = offline
@@ -529,7 +557,16 @@
       CALL getin('config_inca',config_inca)
 
+!Config  Key  = ok_dynzon 
+!Config  Desc = calcul et sortie des transports 
+!Config  Def  = n 
+!Config  Help = Permet de mettre en route le calcul des transports 
+!Config          
+      ok_dynzon = .FALSE. 
+      CALL getin('ok_dynzon',ok_dynzon) 
+
 
       write(lunout,*)' #########################################'
       write(lunout,*)' Configuration des parametres du gcm: '
+      write(lunout,*)' planet_type = ', planet_type
       write(lunout,*)' dayref = ', dayref
       write(lunout,*)' anneeref = ', anneeref
@@ -540,4 +577,5 @@
       write(lunout,*)' iecri = ', iecri
       write(lunout,*)' periodav = ', periodav 
+      write(lunout,*)' output_grads_dyn = ', output_grads_dyn
       write(lunout,*)' idissip = ', idissip
       write(lunout,*)' lstardis = ', lstardis
@@ -550,4 +588,5 @@
       write(lunout,*)' coefdis = ', coefdis
       write(lunout,*)' purmats = ', purmats
+      write(lunout,*)' read_start = ', read_start
       write(lunout,*)' iflag_phys = ', iflag_phys
       write(lunout,*)' clonn = ', clonn 
@@ -562,4 +601,5 @@
       write(lunout,*)' offline = ', offline
       write(lunout,*)' config_inca = ', config_inca
+      write(lunout,*)' ok_dynzon = ', ok_dynzon 
 
       RETURN
@@ -600,5 +640,6 @@
 
       IF( grossismx.LT.1. )  THEN
-        PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
+        write(lunout,*)
+     &   'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
          STOP
       ELSE
@@ -608,5 +649,6 @@
 
       IF( grossismy.LT.1. )  THEN
-        PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
+        write(lunout,*)
+     &  'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
          STOP
       ELSE
@@ -614,5 +656,5 @@
       ENDIF
 
-      PRINT *,' alphax alphay defrun ',alphax,alphay
+      write(lunout,*)'conf_gcm: alphax alphay ',alphax,alphay
 c
 c    alphax et alphay sont les anciennes formulat. des grossissements
@@ -685,4 +727,12 @@
       CALL getin('config_inca',config_inca)
 
+!Config  Key  = ok_dynzon 
+!Config  Desc = calcul et sortie des transports 
+!Config  Def  = n 
+!Config  Help = Permet de mettre en route le calcul des transports 
+!Config          
+      ok_dynzon = .FALSE. 
+      CALL getin('ok_dynzon',ok_dynzon) 
+
 !Config  Key  = use_filtre_fft
 !Config  Desc = flag d'activation des FFT pour le filtre
@@ -692,11 +742,21 @@
       use_filtre_fft=.FALSE.
       CALL getin('use_filtre_fft',use_filtre_fft)
+
+      IF (use_filtre_fft .AND. grossismx /= 1.0) THEN
+        write(lunout,*)'WARNING !!! '
+        write(lunout,*)"Le zoom en longitude est incompatible",
+     &                 " avec l'utilisation du filtre FFT ",
+     &                 "---> filtre FFT dÃ©sactivÃ© "
+       use_filtre_fft=.FALSE.
+      ENDIF
+      
+ 
       
 !Config  Key  = use_mpi_alloc
-!Config  Desc = Utilise un buffer MPI en mémoire globale
+!Config  Desc = Utilise un buffer MPI en mï¿½moire globale
 !Config  Def  = false
 !Config  Help = permet d'activer l'utilisation d'un buffer MPI
-!Config         en mémoire globale a l'aide de la fonction MPI_ALLOC.
-!Config         Cela peut améliorer la bande passante des transferts MPI
+!Config         en mï¿½moire globale a l'aide de la fonction MPI_ALLOC.
+!Config         Cela peut amï¿½liorer la bande passante des transferts MPI
 !Config         d'un facteur 2  
       use_mpi_alloc=.FALSE.
@@ -706,7 +766,7 @@
 !Config  Desc = taille des blocs openmp
 !Config  Def  = 1
-!Config  Help = defini la taille des packets d'itération openmp
-!Config         distribuée à chaque tâche lors de l'entrée dans une
-!Config         boucle parallélisée
+!Config  Help = defini la taille des packets d'itï¿½ration openmp
+!Config         distribuï¿½e ï¿½ chaque tï¿½che lors de l'entrï¿½e dans une
+!Config         boucle parallï¿½lisï¿½e
   
       omp_chunk=1
@@ -716,5 +776,5 @@
 !Config  Desc = activation de la version strato
 !Config  Def  = .FALSE.
-!Config  Help = active la version stratosphérique de LMDZ de F. Lott
+!Config  Help = active la version stratosphï¿½rique de LMDZ de F. Lott
 
       ok_strato=.FALSE.
@@ -731,4 +791,5 @@
       write(lunout,*)' #########################################'
       write(lunout,*)' Configuration des parametres du gcm: '
+      write(lunout,*)' planet_type = ', planet_type
       write(lunout,*)' dayref = ', dayref
       write(lunout,*)' anneeref = ', anneeref
@@ -739,4 +800,5 @@
       write(lunout,*)' iecri = ', iecri
       write(lunout,*)' periodav = ', periodav 
+      write(lunout,*)' output_grads_dyn = ', output_grads_dyn
       write(lunout,*)' idissip = ', idissip
       write(lunout,*)' lstardis = ', lstardis
@@ -749,4 +811,5 @@
       write(lunout,*)' coefdis = ', coefdis
       write(lunout,*)' purmats = ', purmats
+      write(lunout,*)' read_start = ', read_start
       write(lunout,*)' iflag_phys = ', iflag_phys
       write(lunout,*)' clon = ', clon 
@@ -754,5 +817,5 @@
       write(lunout,*)' grossismx = ', grossismx
       write(lunout,*)' grossismy = ', grossismy
-      write(lunout,*)' fxyhypbb = ', fxyhypbb
+      write(lunout,*)' fxyhypb = ', fxyhypb
       write(lunout,*)' dzoomx = ', dzoomx
       write(lunout,*)' dzoomy = ', dzoomy
@@ -761,4 +824,5 @@
       write(lunout,*)' offline = ', offline
       write(lunout,*)' config_inca = ', config_inca
+      write(lunout,*)' ok_dynzon = ', ok_dynzon 
       write(lunout,*)' use_filtre_fft = ', use_filtre_fft
       write(lunout,*)' use_mpi_alloc = ', use_mpi_alloc
Index: LMDZ4/trunk/libf/dyn3dpar/control.h
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/control.h	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/control.h	(revision 1146)
@@ -14,5 +14,6 @@
      &              iperiod,iapp_tracvl,iconser,iecri,idissip,iphysiq , &
      &              periodav,iecrimoy,dayref,anneeref,                  &
-     &              raz_date,offline,ip_ebil_dyn,config_inca
+     &              raz_date,offline,ip_ebil_dyn,config_inca,           &
+     &              planet_type,output_grads_dyn,ok_dynzon
 
       INTEGER   nday,day_step,iperiod,iapp_tracvl,iconser,iecri,        &
@@ -21,4 +22,8 @@
       REAL periodav
       logical offline
-      CHARACTER*4 config_inca
+      CHARACTER (len=4) :: config_inca
+      CHARACTER(len=10) :: planet_type ! planet type ('earth','mars',...)
+      LOGICAL :: output_grads_dyn ! output dynamics diagnostics in
+                                  ! binary grads file 'dyn.dat' (y/n)
+      LOGICAL :: ok_dynzon 
 !-----------------------------------------------------------------------
Index: LMDZ4/trunk/libf/dyn3dpar/covnat_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/covnat_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/covnat_p.F	(revision 1146)
@@ -61,6 +61,7 @@
       END DO
 
-      ijb=ij_begin
+      ijb=ij_begin-iip1
       ije=ij_end
+      if (pole_nord) ijb=ij_begin
       if (pole_sud)  ije=ij_end-iip1
      
Index: LMDZ4/trunk/libf/dyn3dpar/create_etat0_limit.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/create_etat0_limit.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/create_etat0_limit.F	(revision 1146)
@@ -8,4 +8,5 @@
        USE mod_const_mpi
        USE phys_state_var_mod      
+       USE infotrac
        IMPLICIT NONE
 c
@@ -31,5 +32,4 @@
 #include "paramet.h"
 #include "indicesol.h"
-#include "advtrac.h"
 #include  "control.h"
 #include "clesphys.h"
@@ -37,21 +37,11 @@
 !      REAL :: pctsrf(iim*(jjm-1)+2, nbsrf)
 
-c initialisation traceurs
-      hadv_flg(:) = 0.
-      vadv_flg(:) = 0.
-      conv_flg(:) = 0.
-      pbl_flg(:)  = 0. 
-      tracnam(:)  = '        ' 
-      nprath = 1
-      nbtrac = 0
-      mmt_adj(:,:,:,:) = 1
-
       IF (config_inca /= 'none') THEN
 #ifdef INCA
          call init_const_lmdz(
-     $        nbtrac,anneeref,dayref,
+     $        nbtr,anneeref,dayref,
      $        iphysiq,day_step,nday)
 #endif
-         print *, 'nbtrac =' , nbtrac 
+         print *, 'nbtr =' , nbtr 
       END IF
 
@@ -59,5 +49,5 @@
 
 
-      CALL Init_Phys_lmdz(iim,jjp1,llm,nqmx-2,1,(/(jjm-1)*iim+2/))
+      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
       PRINT *,'---> klon=',klon
 
Index: LMDZ4/trunk/libf/dyn3dpar/diagedyn.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/diagedyn.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/diagedyn.F	(revision 1146)
@@ -58,6 +58,7 @@
 #include "paramet.h"
 #include "comgeom.h"
-
-#ifdef CPP_PHYS
+#include "iniprint.h"
+
+#ifdef CPP_EARTH
 #include "../phylmd/YOMCST.h"
 #include "../phylmd/YOETHF.h"
@@ -139,5 +140,5 @@
 
 
-#ifdef CPP_PHYS
+#ifdef CPP_EARTH
 c======================================================================
 C     Compute Kinetic enrgy
@@ -314,6 +315,7 @@
 C
 #else
-      print*,'Pour l instant diagedyn a besoin de la physique'
+      write(lunout,*),'diagedyn: Needs Earth physics to function'
 #endif
+! #endif of #ifdef CPP_EARTH 
       RETURN 
       END 
Index: LMDZ4/trunk/libf/dyn3dpar/disvert0.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/disvert0.F	(revision 1122)
+++ 	(revision )
@@ -1,183 +1,0 @@
-!
-! $Header$
-!
-c
-c
-      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
-
-c    Auteur :  P. Le Van .
-c
-      IMPLICIT NONE
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "iniprint.h"
-#include "logic.h"
-c
-c=======================================================================
-c
-c
-c    s = sigma ** kappa   :  coordonnee  verticale
-c    dsig(l)            : epaisseur de la couche l ds la coord.  s
-c    sig(l)             : sigma a l'interface des couches l et l-1
-c    ds(l)              : distance entre les couches l et l-1 en coord.s
-c
-c=======================================================================
-c
-      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
-      REAL pa,preff,presnivs(llm)
-c
-c   declarations:
-c   -------------
-c
-      REAL sig(llm+1),dsig(llm)
-       real zzz(1:llm+1)
-       real dzz(1:llm)
-      real zk,zkm1,dzk1,dzk2,k0,k1
-c
-      INTEGER l
-      REAL snorm
-      REAL alpha,beta,gama,delta,deltaz,h
-      INTEGER np,ierr
-      REAL pi,x
-
-      REAL SSUM
-
-      pa        =  50 000.
-      preff     = 101 325.
-
-c
-c-----------------------------------------------------------------------
-c
-      pi=2.*ASIN(1.)
-
-      OPEN(99,file='sigma.def',status='old',form='formatted',
-     s   iostat=ierr)
-
-c-----------------------------------------------------------------------
-c   cas 1 on lit les options dans sigma.def:
-c   ----------------------------------------
-
-      IF (ierr.eq.0) THEN
-
-      READ(99,*) h           ! hauteur d'echelle 8.
-      READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
-      READ(99,*) beta        ! facteur d'acroissement en haut 1.3
-      READ(99,*) k0          ! nombre de couches dans la transition surf
-      READ(99,*) k1          ! nombre de couches dans la transition haute
-      CLOSE(99)
-      alpha=deltaz/(llm*h)
-      WRITE(lunout,*)'h,alpha,k0,k1,beta'
-
-c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
-
-      alpha=deltaz/tanh(1./k0)*2.
-      zkm1=0.
-      sig(1)=1.
-      do l=1,llm
-        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
-     + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
-        zk=-h*log(sig(l+1))
-
-        dzk1=alpha*tanh(l/k0)
-        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
-        WRITE(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
-        zkm1=zk
-      enddo
-
-      sig(llm+1)=0.
-
-c
-       DO 2  l = 1, llm
-       dsig(l) = sig(l)-sig(l+1)
-   2   CONTINUE
-c
-
-      ELSE
-c-----------------------------------------------------------------------
-c   cas 2 ancienne discretisation (LMD5...):
-c   ----------------------------------------
-
-      write(lunout,*)'WARNING!!! Ancienne discretisation verticale'
-
-      h=7.
-      snorm  = 0.
-      DO l = 1, llm
-         x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
-
-         IF (ok_strato) THEN
-           dsig(l) =(1.0 + 7.0 * SIN(x)**2)
-     &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2        
-         ELSE
-           dsig(l) = 1.0 + 7.0 * SIN(x)**2
-         ENDIF
-
-         snorm = snorm + dsig(l)
-      ENDDO
-      snorm = 1./snorm
-      DO l = 1, llm
-         dsig(l) = dsig(l)*snorm
-      ENDDO
-      sig(llm+1) = 0.
-      DO l = llm, 1, -1
-         sig(l) = sig(l+1) + dsig(l)
-      ENDDO
-
-      ENDIF
-
-
-      DO l=1,llm
-        nivsigs(l) = FLOAT(l)
-      ENDDO
-
-      DO l=1,llmp1
-        nivsig(l)= FLOAT(l)
-      ENDDO
-
-c
-c    ....  Calculs  de ap(l) et de bp(l)  ....
-c    .........................................
-c
-c
-c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
-c
-
-      bp(llmp1) =   0.
-
-      DO l = 1, llm
-cc
-ccc    ap(l) = 0.
-ccc    bp(l) = sig(l)
-
-      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
-      ap(l) = pa * ( sig(l) - bp(l) )
-c
-      ENDDO
-      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
-
-      write(lunout,*)' BP '
-      write(lunout,*)  bp
-      write(lunout,*)' AP '
-      write(lunout,*)  ap
-
-      write(lunout,*)
-     .'Niveaux de pressions approximatifs aux centres des'
-      write(lunout,*)
-     .'couches calcules pour une pression de surface =',preff
-      write(lunout,*)
-     .'et altitudes equivalentes pour une hauteur d echelle de'
-      write(lunout,*)'8km'
-      DO l = 1, llm
-       dpres(l) = bp(l) - bp(l+1)
-       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
-       write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
-     .        log(preff/presnivs(l))*8.
-     .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
-     .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
-      ENDDO
-
-      write(lunout,*)' PRESNIVS '
-      write(lunout,*)presnivs
-
-      RETURN
-      END
Index: LMDZ4/trunk/libf/dyn3dpar/dynetat0.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/dynetat0.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/dynetat0.F	(revision 1146)
@@ -2,6 +2,7 @@
 ! $Header$
 !
-      SUBROUTINE dynetat0(fichnom,nq,vcov,ucov,
+      SUBROUTINE dynetat0(fichnom,vcov,ucov,
      .                    teta,q,masse,ps,phis,time)
+      USE infotrac
       IMPLICIT NONE
 
@@ -32,5 +33,4 @@
 #include "serre.h"
 #include "logic.h"
-#include "advtrac.h"
 
 c   Arguments:
@@ -38,7 +38,6 @@
 
       CHARACTER*(*) fichnom
-      INTEGER nq
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
-      REAL q(ip1jmp1,llm,nq),masse(ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
       REAL ps(ip1jmp1),phis(ip1jmp1)
 
@@ -53,5 +52,4 @@
 
 c-----------------------------------------------------------------------
-
 c  Ouverture NetCDF du fichier etat initial
 
@@ -315,6 +313,5 @@
 
 
-      IF(nq.GE.1) THEN
-      DO iq=1,nq
+      DO iq=1,nqtot
         ierr =  NF_INQ_VARID (nid, tname(iq), nvarid)
         IF (ierr .NE. NF_NOERR) THEN
@@ -334,5 +331,4 @@
         ENDIF
       ENDDO
-      ENDIF
 
       ierr = NF_INQ_VARID (nid, "masse", nvarid)
Index: LMDZ4/trunk/libf/dyn3dpar/dynredem.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/dynredem.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/dynredem.F	(revision 1146)
@@ -3,6 +3,7 @@
 !
 c
-      SUBROUTINE dynredem0(fichnom,iday_end,phis,nq)
+      SUBROUTINE dynredem0(fichnom,iday_end,phis)
       USE IOIPSL
+      USE infotrac
       IMPLICIT NONE
 c=======================================================================
@@ -22,5 +23,4 @@
 #include "description.h"
 #include "serre.h"
-#include "advtrac.h"
 
 c   Arguments:
@@ -29,5 +29,4 @@
       REAL phis(ip1jmp1)
       CHARACTER*(*) fichnom
-      INTEGER nq
 
 c   Local:
@@ -458,6 +457,6 @@
       dims4(3) = idim_s
       dims4(4) = idim_tim
-      IF(nq.GE.1) THEN
-      DO iq=1,nq
+
+      DO iq=1,nqtot
 cIM 220306 BEG
 #ifdef NC_DOUBLE
@@ -469,5 +468,4 @@
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12,ttext(iq))
       ENDDO
-      ENDIF
 c
       dims4(1) = idim_rlonv
@@ -508,5 +506,6 @@
       END
       SUBROUTINE dynredem1(fichnom,time,
-     .                     vcov,ucov,teta,q,nq,masse,ps)
+     .                     vcov,ucov,teta,q,masse,ps)
+      USE infotrac
       IMPLICIT NONE
 c=================================================================
@@ -519,13 +518,12 @@
 #include "comvert.h"
 #include "comgeom.h"
-#include "advtrac.h"
 #include "temps.h"
 #include "control.h"
 
-      INTEGER nq, l
+      INTEGER l
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) 
       REAL teta(ip1jmp1,llm)                   
       REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
-      REAL q(ip1jmp1,llm,nq)
+      REAL q(ip1jmp1,llm,nqtot)
       CHARACTER*(*) fichnom
      
@@ -633,6 +631,5 @@
       END IF
 
-      IF(nq.GE.1) THEN
-      do iq=1,nq 
+      do iq=1,nqtot 
 
          IF (config_inca == 'none') THEN
@@ -704,5 +701,4 @@
       
       ENDDO
-      ENDIF
 c
       ierr = NF_INQ_VARID(nid, "masse", nvarid)
Index: LMDZ4/trunk/libf/dyn3dpar/dynredem_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/dynredem_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/dynredem_p.F	(revision 1146)
@@ -3,7 +3,8 @@
 !
 c
-      SUBROUTINE dynredem0_p(fichnom,iday_end,phis,nq)
+      SUBROUTINE dynredem0_p(fichnom,iday_end,phis)
       USE IOIPSL
       USE parallel
+      USE infotrac
       IMPLICIT NONE
 c=======================================================================
@@ -23,5 +24,4 @@
 #include "description.h"
 #include "serre.h"
-#include "advtrac.h"
 
 c   Arguments:
@@ -30,5 +30,4 @@
       REAL phis(ip1jmp1)
       CHARACTER*(*) fichnom
-      INTEGER nq
 
 c   Local:
@@ -54,5 +53,4 @@
       INTEGER yyears0,jjour0, mmois0
       character*30 unites
-
 
 c-----------------------------------------------------------------------
@@ -461,6 +459,6 @@
       dims4(3) = idim_s
       dims4(4) = idim_tim
-      IF(nq.GE.1) THEN
-      DO iq=1,nq
+
+      DO iq=1,nqtot
 cIM 220306 BEG
 #ifdef NC_DOUBLE
@@ -472,5 +470,4 @@
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12,ttext(iq))
       ENDDO
-      ENDIF
 c
       dims4(1) = idim_rlonv
@@ -513,6 +510,7 @@
       END
       SUBROUTINE dynredem1_p(fichnom,time,
-     .                     vcov,ucov,teta,q,nq,masse,ps)
+     .                     vcov,ucov,teta,q,masse,ps)
       USE parallel
+      USE infotrac
       IMPLICIT NONE
 c=================================================================
@@ -525,13 +523,12 @@
 #include "comvert.h"
 #include "comgeom.h"
-#include "advtrac.h"
 #include "temps.h"
 #include "control.h"
 
-      INTEGER nq, l
+      INTEGER l
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) 
       REAL teta(ip1jmp1,llm)                   
       REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
-      REAL q(ip1jmp1,llm,nq)
+      REAL q(ip1jmp1,llm,nqtot)
       CHARACTER*(*) fichnom
      
@@ -559,5 +556,5 @@
       call Gather_Field(ps,ip1jmp1,1,0)
       
-      do iq=1,nq
+      do iq=1,nqtot
         call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
       enddo
@@ -660,6 +657,5 @@
       END IF
 
-      IF(nq.GE.1) THEN
-      do iq=1,nq 
+      do iq=1,nqtot 
 
          IF (config_inca == 'none') THEN
@@ -731,6 +727,4 @@
       
       ENDDO
-      ENDIF
-
 
 
Index: LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F	(revision 1146)
@@ -5,5 +5,5 @@
 c
       SUBROUTINE etat0_netcdf (interbar, masque)
-    
+#ifdef CPP_EARTH        
       USE startvar
       USE ioipsl
@@ -12,4 +12,8 @@
       USE pbl_surface_mod
       USE phys_state_var_mod
+      USE filtreg_mod
+      USE infotrac
+#endif
+!#endif of #ifdef CPP_EARTH
       !
       IMPLICIT NONE
@@ -23,4 +27,5 @@
 !     .KLON=KFDIA-KIDIA+1,KLEV=llm
       !
+#ifdef CPP_EARTH    
 #include "comgeom2.h"
 #include "comvert.h"
@@ -29,8 +34,14 @@
 #include "dimsoil.h"
 #include "temps.h"
-      !
+#endif
+!#endif of #ifdef CPP_EARTH
+      ! arguments:
       LOGICAL interbar
+      REAL :: masque(iip1,jjp1)
+
+#ifdef CPP_EARTH
+      ! local variables:
       REAL :: latfi(klon), lonfi(klon)
-      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
+      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1),
      . psol(iip1, jjp1), phis(iip1, jjp1)
       REAL :: p3d(iip1, jjp1, llm+1)
@@ -38,5 +49,6 @@
       REAL :: vvent(iip1, jjm, llm)
       REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
-      REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
+      REAL :: qsat(iip1, jjp1, llm)
+      REAL,ALLOCATABLE :: q3d(:, :, :,:)
       REAL :: tsol(klon), qsol(klon), sn(klon)
       REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 
@@ -141,4 +153,5 @@
       !
       preff     = 101325.
+      pa        =  50000.
       unskap = 1./kappa
       !
@@ -164,8 +177,11 @@
       print*,'dtvr',dtvr
 
-      CALL inicons0()
+      CALL iniconst()
       CALL inigeom()
       !
       CALL inifilr()
+C init pour traceurs
+      call infotrac_init
+      ALLOCATE(q3d(iip1, jjp1, llm,nqtot))
 !      CALL phys_state_var_init()
       !
@@ -623,6 +639,4 @@
       phis(iip1,:) = phis(1,:)
 
-C init pour traceurs
-      call iniadvtrac(nq)
 C Ecriture
       CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
@@ -648,7 +662,7 @@
      *                phi,w, pbaru,pbarv,time+iday-dayref   )
        print*,'sortie caldyn0'     
-      CALL dynredem0("start.nc",dayref,phis,nqmx)
+      CALL dynredem0("start.nc",dayref,phis)
       print*,'sortie dynredem0'
-      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
+      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,masse ,
      .                            psol)
       print*,'sortie dynredem1' 
@@ -742,6 +756,6 @@
       visu_file='Etat0_visu.nc'
       CALL initdynav(visu_file,dayref,anneeref,time_step,
-     .              t_ops, t_wrt, nqmx, visuid)
-      CALL writedynav(visuid, nqmx, itau,vvent ,
+     .              t_ops, t_wrt, visuid)
+      CALL writedynav(visuid, itau,vvent ,
      .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
       else
@@ -750,4 +764,7 @@
       print*,'entree histclo'
       CALL histclo
+
+#endif 
+!#endif of #ifdef CPP_EARTH
       RETURN
       !
Index: LMDZ4/trunk/libf/dyn3dpar/filtreg_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/filtreg_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/filtreg_p.F	(revision 1146)
@@ -2,10 +2,13 @@
 
       SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv, 
-     .                       ifiltre, iaire, griscal ,iter)
-      USE Parallel, only : OMP_CHUNK 
+     &     ifiltre, iaire, griscal ,iter)
+      USE Parallel, only : OMP_CHUNK
       USE mod_filtre_fft
       USE timer_filtre
+      
+      USE filtreg_mod
+      
       IMPLICIT NONE
-
+      
 c=======================================================================
 c
@@ -50,5 +53,4 @@
 #include "dimensions.h"
 #include "paramet.h"
-#include "parafilt.h"
 #include "coefils.h"
 c
@@ -57,14 +59,6 @@
       INTEGER iim2,immjm
       INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
-
+      
       REAL  champ( iip1,nlat,nbniv)
-      REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs
-      COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus)
-     ,             , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs)
-     ,             ,  matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus)
-cym      REAL  eignq(iim), sdd1(iim),sdd2(iim)
-
-      REAL  eignq(iim)
-      REAL :: sdd1(iim),sdd2(iim)
       
       LOGICAL    griscal
@@ -74,23 +68,27 @@
       REAL :: champ_in(iip1,nlat,nbniv)
       
-      REAL,SAVE,TARGET :: sddu_loc(iim)
-      REAL,SAVE,TARGET :: sddv_loc(iim)
-      REAL,SAVE,TARGET :: unsddu_loc(iim)
-      REAL,SAVE,TARGET :: unsddv_loc(iim)
-c$OMP THREADPRIVATE(sddu_loc,sddv_loc,unsddu_loc,unsddv_loc)
       LOGICAL,SAVE     :: first=.TRUE.
 c$OMP THREADPRIVATE(first) 
 
+      REAL, DIMENSION(iip1,nlat,nbniv) :: champ_loc
+      INTEGER :: ll_nb, nbniv_loc
+      REAL, SAVE :: sdd12(iim,4)
+c$OMP THREADPRIVATE(sdd12) 
+
+      INTEGER, PARAMETER :: type_sddu=1
+      INTEGER, PARAMETER :: type_sddv=2
+      INTEGER, PARAMETER :: type_unsddu=3
+      INTEGER, PARAMETER :: type_unsddv=4
+
+      INTEGER :: sdd1_type, sdd2_type
+
       IF (first) THEN
-        sddu_loc(1:iim)=sddu(1:iim)
-        sddv_loc(1:iim)=sddv(1:iim)
-        unsddu_loc(1:iim)=unsddu(1:iim)
-        unsddv_loc(1:iim)=unsddv(1:iim)
-        CALL Init_timer
-        first=.FALSE.
-c	PRINT *,"----> sddu_loc=",sddu_loc
-c	PRINT *,"----> sddv_loc=",sddv_loc
-c	PRINT *,"----> unsddu_loc=",unsddu_loc
-c	PRINT *,"----> unsddv_loc=",unsddv_loc
+         sdd12(1:iim,type_sddu) = sddu(1:iim)
+         sdd12(1:iim,type_sddv) = sddv(1:iim)
+         sdd12(1:iim,type_unsddu) = unsddu(1:iim)
+         sdd12(1:iim,type_unsddv) = unsddv(1:iim)
+
+         CALL Init_timer
+         first=.FALSE.
       ENDIF
 
@@ -99,23 +97,25 @@
 c$OMP END MASTER
 
+c-------------------------------------------------------c
+
       IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) 
-     *    STOP'Pas de transformee simple dans cette version'
-
+     &     STOP'Pas de transformee simple dans cette version'
+      
       IF( iter.EQ. 2 )  THEN
-       PRINT *,' Pas d iteration du filtre dans cette version !'
-     * , ' Utiliser old_filtreg et repasser !'
-           STOP
+         PRINT *,' Pas d iteration du filtre dans cette version !'
+     &        , ' Utiliser old_filtreg et repasser !'
+         STOP
       ENDIF
 
       IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
-       PRINT *,' Cette routine ne calcule le filtre inverse que ',
-     * ' sur la grille des scalaires !'
-           STOP
+         PRINT *,' Cette routine ne calcule le filtre inverse que '
+     &        , ' sur la grille des scalaires !'
+         STOP
       ENDIF
 
       IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
-       PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
-     *,' corriger et repasser !'
-           STOP
+         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
+     &        , ' corriger et repasser !'
+         STOP
       ENDIF
 c
@@ -127,77 +127,53 @@
       IF( griscal )   THEN
          IF( nlat. NE. jjp1 )  THEN
-             PRINT  1111
-             STOP
+            PRINT  1111
+            STOP
          ELSE
-c
-             IF( iaire.EQ.1 )  THEN
-cym                CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 ) 
-cym                CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
-cym               sdd1=>sddv_loc
-cym               sdd2=>unsddv_loc
-               sdd1(1:iim)=sddv_loc(1:iim)
-               sdd2(1:iim)=unsddv_loc(1:iim)
-             ELSE
-cym                CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
-cym                CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
-               sdd1(1:iim)=unsddv_loc(1:iim)
-               sdd2(1:iim)=sddv_loc(1:iim)
-             END IF
-c
-             jdfil1 = 2
-             jffil1 = jfiltnu
-             jdfil2 = jfiltsu
-             jffil2 = jjm
-          END IF
+c     
+            IF( iaire.EQ.1 )  THEN
+               sdd1_type = type_sddv
+               sdd2_type = type_unsddv
+            ELSE
+               sdd1_type = type_unsddv
+               sdd2_type = type_sddv
+            ENDIF
+c
+            jdfil1 = 2
+            jffil1 = jfiltnu
+            jdfil2 = jfiltsu
+            jffil2 = jjm
+         ENDIF
       ELSE
-          IF( nlat.NE.jjm )  THEN
-             PRINT  2222
-             STOP
-          ELSE
-c
-             IF( iaire.EQ.1 )  THEN
-cym                CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 ) 
-cym                CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
-cym                sdd1=>sddu_loc
-cym                sdd2=>unsddu_loc
-                sdd1(1:iim)=sddu_loc(1:iim)
-                sdd2(1:iim)=unsddu_loc(1:iim)
-
-             ELSE
-cym                CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
-cym                CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
-cym               sdd1=>unsddu_loc
-cym               sdd2=>sddu_loc
-	       sdd1(1:iim)=unsddu_loc(1:iim)
-               sdd2(1:iim)=sddu_loc(1:iim)
-             END IF
-c
-             jdfil1 = 1
-             jffil1 = jfiltnv
-             jdfil2 = jfiltsv
-             jffil2 = jjm
-          END IF
-      END IF
-
-c      PRINT *,"APPEL a filtreg --> sdd1=",sdd1
-c      PRINT *,"APPEL a filtreg --> sdd2=",sdd2
-c      PRINT *,"----> sddu_loc=",sddu_loc
-c	PRINT *,"----> sddv_loc=",sddv_loc
-c	PRINT *,"----> unsddu_loc=",unsddu_loc
-c	PRINT *,"----> unsddv_loc=",unsddv_loc
- 
-c
-c
-      DO 100  hemisph = 1, 2
-c
-      IF ( hemisph.EQ.1 )  THEN
-c ym
-          jdfil = max(jdfil1,ibeg)
-          jffil = min(jffil1,iend)
-      ELSE
-c ym
-          jdfil = max(jdfil2,ibeg)
-          jffil = min(jffil2,iend)
-      END IF
+         IF( nlat.NE.jjm )  THEN
+            PRINT  2222
+            STOP
+         ELSE
+c
+            IF( iaire.EQ.1 )  THEN
+               sdd1_type = type_sddu
+               sdd2_type = type_unsddu
+            ELSE
+               sdd1_type = type_unsddu
+               sdd2_type = type_sddu
+            ENDIF
+c     
+            jdfil1 = 1
+            jffil1 = jfiltnv
+            jdfil2 = jfiltsv
+            jffil2 = jjm
+         ENDIF
+      ENDIF
+c      
+      DO hemisph = 1, 2
+c     
+         IF ( hemisph.EQ.1 )  THEN
+cym
+            jdfil = max(jdfil1,ibeg)
+            jffil = min(jffil1,iend)
+         ELSE
+cym
+            jdfil = max(jdfil2,ibeg)
+            jffil = min(jffil2,iend)
+         ENDIF
 
 
@@ -206,144 +182,199 @@
 cccccccccccccccccccccccccccccccccccccccccccc
 
-      IF (.NOT. use_filtre_fft) THEN
-      
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)    
-      DO 50  l = 1, nbniv
-        DO 30  j = jdfil,jffil
- 
- 
-          DO  5  i = 1, iim
-            champ(i,j,l) = champ(i,j,l) * sdd1(i)
-   5      CONTINUE
-c
-
-          IF( hemisph. EQ. 1 )      THEN
-
-            IF( ifiltre. EQ. -2 )   THEN
-
-
-              CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
-     .                     champ(1,j,l), 1, 0.0, eignq, 1)
-
-
-            ELSE IF ( griscal )     THEN
-
-              CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
-     .                    champ(1,j,l), 1, 0.0, eignq, 1)
-
-            ELSE 
-
-              CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
-     .                   champ(1,j,l), 1, 0.0, eignq, 1)
-            ENDIF
-
-          ELSE
-
-            IF( ifiltre. EQ. -2 )   THEN
-      
-              CALL SGEMV("N",iim,iim,1.0, matrinvs(1,1,j-jfiltsu+1),iim,
-     .                   champ(1,j,l), 1, 0.0, eignq, 1)
-      
-            ELSE IF ( griscal )     THEN
-      
-              CALL SGEMV("N",iim,iim,1.0,matriceus(1,1,j-jfiltsu+1),iim,
-     .                   champ(1,j,l), 1, 0.0, eignq, 1)
-            ELSE 
-          
-	      CALL SGEMV("N",iim,iim,1.0,matricevs(1,1,j-jfiltsv+1),iim,
-     .                    champ(1,j,l), 1, 0.0, eignq, 1)
-            ENDIF
-
-          ENDIF
-
-
-c
-          IF( ifiltre.EQ. 2 )  THEN
-         
-	    DO 15 i = 1, iim
-              champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
-  15        CONTINUE
-         
-	  ELSE
-        
-	    DO 16 i=1,iim
-               champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
-16          CONTINUE
-          
-	  ENDIF
-c
-          champ( iip1,j,l ) = champ( 1,j,l )
-c
-  30    CONTINUE
-c
-  50  CONTINUE
+         IF (.NOT. use_filtre_fft) THEN
+      
+c     !---------------------------------!
+c     ! Agregation des niveau verticaux !
+c     ! uniquement necessaire pour une  !
+c     ! execution OpenMP                !
+c     !---------------------------------!
+            ll_nb = 0
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+            DO l = 1, nbniv
+               ll_nb = ll_nb+1
+               DO j = jdfil,jffil
+                  DO i = 1, iim
+                     champ_loc(i,j,ll_nb) = 
+     &                    champ(i,j,l) * sdd12(i,sdd1_type)
+                  ENDDO
+               ENDDO
+            ENDDO
 c$OMP END DO NOWAIT
 
+            nbniv_loc = ll_nb
+
+            IF( hemisph.EQ.1 )      THEN
+               
+               IF( ifiltre.EQ.-2 )   THEN
+                  DO j = jdfil,jffil
+                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
+     &                    matrinvn(1,1,j), iim, 
+     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
+     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+                  ENDDO
+                  
+               ELSE IF ( griscal )     THEN
+                  DO j = jdfil,jffil
+                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
+     &                    matriceun(1,1,j), iim, 
+     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
+     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+                  ENDDO
+                  
+               ELSE 
+                  DO j = jdfil,jffil
+                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
+     &                    matricevn(1,1,j), iim, 
+     &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
+     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+                  ENDDO
+                  
+               ENDIF
+               
+            ELSE
+               
+               IF( ifiltre.EQ.-2 )   THEN
+                  DO j = jdfil,jffil
+                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
+     &                    matrinvs(1,1,j-jfiltsu+1), iim, 
+     &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
+     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+                  ENDDO
+                  
+               ELSE IF ( griscal )     THEN
+                  
+                  DO j = jdfil,jffil
+                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
+     &                    matriceus(1,1,j-jfiltsu+1), iim, 
+     &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
+     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+                  ENDDO
+                  
+               ELSE 
+                  
+                  DO j = jdfil,jffil
+                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
+     &                    matricevs(1,1,j-jfiltsv+1), iim, 
+     &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
+     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+                  ENDDO
+                  
+               ENDIF
+               
+            ENDIF
+!     c     
+            IF( ifiltre.EQ.2 )  THEN
+               
+c     !-------------------------------------!
+c     ! Dés-agregation des niveau verticaux !
+c     ! uniquement necessaire pour une      !
+c     ! execution OpenMP                    !
+c     !-------------------------------------!
+               ll_nb = 0
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+               DO l = 1, nbniv
+                  ll_nb = ll_nb + 1
+                  DO j = jdfil,jffil
+                     DO i = 1, iim
+                        champ( i,j,l ) = (champ_loc(i,j,ll_nb) 
+     &                       + champ_fft(i,j-jdfil+1,ll_nb))
+     &                       * sdd12(i,sdd2_type)
+                     ENDDO
+                  ENDDO
+               ENDDO
+c$OMP END DO NOWAIT
+               
+            ELSE
+               
+               ll_nb = 0
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+               DO l = 1, nbniv_loc
+                  ll_nb = ll_nb + 1
+                  DO j = jdfil,jffil
+                     DO i = 1, iim
+                        champ( i,j,l ) = (champ_loc(i,j,ll_nb) 
+     &                       - champ_fft(i,j-jdfil+1,ll_nb))
+     &                       * sdd12(i,sdd2_type)
+                     ENDDO
+                  ENDDO
+               ENDDO
+c$OMP END DO NOWAIT
+               
+            ENDIF
+            
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+            DO l = 1, nbniv
+               DO j = jdfil,jffil
+                  champ( iip1,j,l ) = champ( 1,j,l )
+               ENDDO
+            ENDDO
+c$OMP END DO NOWAIT
+            
 ccccccccccccccccccccccccccccccccccccccccccccc
 c Utilisation du filtre FFT
 ccccccccccccccccccccccccccccccccccccccccccccc
         
-       ELSE
+         ELSE
        
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-	  DO l=1,nbniv
-	    DO j=jdfil,jffil
-	      DO  i = 1, iim
-                champ( i,j,l)= champ(i,j,l)*sdd1(i)
-		champ_fft( i,j,l) = champ(i,j,l)
-              ENDDO
+            DO l=1,nbniv
+               DO j=jdfil,jffil
+                  DO  i = 1, iim
+                     champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
+                     champ_fft( i,j,l) = champ(i,j,l)
+                  ENDDO
+               ENDDO
             ENDDO
-	  ENDDO
 c$OMP END DO NOWAIT
 
-      IF (jdfil<=jffil) THEN
-        IF( ifiltre. EQ. -2 )   THEN
-          CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv) 
-        ELSE IF ( griscal )     THEN
-	  CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
-	ELSE
-	  CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
-	ENDIF
-      ENDIF
-
-
-        IF( ifiltre.EQ. 2 )  THEN
+            IF (jdfil<=jffil) THEN
+               IF( ifiltre. EQ. -2 )   THEN
+                  CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv) 
+               ELSE IF ( griscal )     THEN
+                  CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
+               ELSE
+                  CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
+               ENDIF
+            ENDIF
+
+
+            IF( ifiltre.EQ. 2 )  THEN
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)          
-	  DO l=1,nbniv
-	    DO j=jdfil,jffil
-	      DO  i = 1, iim
-                champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
-     .                             *sdd2(i)
-              ENDDO
-            ENDDO
-	  ENDDO
+               DO l=1,nbniv
+                  DO j=jdfil,jffil
+                     DO  i = 1, iim
+                        champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
+     &                       *sdd12(i,sdd2_type)
+                     ENDDO
+                  ENDDO
+               ENDDO
 c$OMP END DO NOWAIT	  
-	ELSE
+            ELSE
         
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)  
-	  DO l=1,nbniv
-	    DO j=jdfil,jffil
-	      DO  i = 1, iim
-                champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
-     .                            *sdd2(i)
-              ENDDO
+               DO l=1,nbniv
+                  DO j=jdfil,jffil
+                     DO  i = 1, iim
+                        champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
+     &                       *sdd12(i,sdd2_type)
+                     ENDDO
+                  ENDDO
+               ENDDO
+c$OMP END DO NOWAIT          
+            ENDIF
+c
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
+            DO l=1,nbniv
+               DO j=jdfil,jffil
+!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
+                  champ( iip1,j,l ) = champ( 1,j,l )
+               ENDDO
             ENDDO
-	  ENDDO
-c$OMP END DO NOWAIT          
-        ENDIF
-c
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
-	DO l=1,nbniv
-	  DO j=jdfil,jffil
-!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
-	    champ( iip1,j,l ) = champ( 1,j,l )
-          ENDDO
-	ENDDO
 c$OMP END DO NOWAIT          	
-      ENDIF 
+         ENDIF 
 c Fin de la zone de filtrage
 
 	
- 100  CONTINUE
+      ENDDO
 
 !      DO j=1,nlat
@@ -359,8 +390,8 @@
       
 c
-1111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a 
-     *filtrer, sur la grille des scalaires'/)
-2222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
-     *ltrer, sur la grille de V ou de Z'/)
+ 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a 
+     &     filtrer, sur la grille des scalaires'/)
+ 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
+     &     ltrer, sur la grille de V ou de Z'/)
 c$OMP MASTER      
       CALL stop_timer
Index: LMDZ4/trunk/libf/dyn3dpar/fluxstokenc_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/fluxstokenc_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/fluxstokenc_p.F	(revision 1146)
@@ -61,5 +61,5 @@
 	CALL initfluxsto_p( 'fluxstoke',
      .  time_step,istdyn* time_step,istdyn* time_step,
-     . nqmx, fluxid,fluxvid,fluxdid) 
+     .  fluxid,fluxvid,fluxdid) 
 	
         ijb=ij_begin
Index: LMDZ4/trunk/libf/dyn3dpar/gcm.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/gcm.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/gcm.F	(revision 1146)
@@ -9,14 +9,22 @@
       USE IOIPSL
 #endif
+
       USE mod_const_mpi, ONLY: init_const_mpi
       USE parallel
       USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
-      USE mod_grid_phy_lmdz
-      USE mod_phys_lmdz_omp_data, ONLY: klon_omp
-      USE dimphy
+      USE infotrac
       USE mod_interface_dyn_phys
-      USE comgeomphy
       USE mod_hallo
       USE Bands
+
+      USE filtreg_mod
+
+! Ehouarn: for now these only apply to Earth:
+#ifdef CPP_EARTH
+      USE mod_grid_phy_lmdz
+      USE mod_phys_lmdz_omp_data, ONLY: klon_omp 
+      USE dimphy
+      USE comgeomphy
+#endif
       IMPLICIT NONE
 
@@ -65,5 +73,4 @@
 #include "iniprint.h"
 #include "tracstoke.h"
-#include "advtrac.h"
 
       INTEGER         longcles
@@ -81,5 +88,5 @@
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
       REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
-      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: q ! champs advectes
       REAL ps(ip1jmp1)                       ! pression  au sol
 c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
@@ -136,5 +143,5 @@
 c    variables pour l'initialisation de la physique :
 c    ------------------------------------------------
-      INTEGER ngridmx,nq
+      INTEGER ngridmx
       PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
       REAL zcufi(ngridmx),zcvfi(ngridmx)
@@ -158,34 +165,4 @@
 
 
-c initialisation Anne
-      hadv_flg(:) = 0.
-      vadv_flg(:) = 0.
-      conv_flg(:) = 0.
-      pbl_flg(:)  = 0. 
-      tracnam(:)  = '        ' 
-      nprath = 1
-      nbtrac = 0
-      mmt_adj(:,:,:,:) = 1
-
-
-c--------------------------------------------------------------------------
-c   Iflag_phys controle l'appel a la physique :
-c   -------------------------------------------
-c      0 : pas de physique
-c      1 : Normale (appel a phylmd, phymars ...)
-c      2 : rappel Newtonien pour la temperature + friction au sol
-      iflag_phys=1
-
-c--------------------------------------------------------------------------
-c   Lecture de l'etat initial :
-c   ---------------------------
-c     T : on lit start.nc
-c     F : le modele s'autoinitialise avec un cas academique (iniacademic)
-#ifdef CPP_IOIPSL
-      read_start=.true.
-#else
-      read_start=.false.
-#endif
-
 c-----------------------------------------------------------------------
 c   Choix du calendrier
@@ -203,9 +180,10 @@
 c  ---------------------------------------
 c
-#ifdef CPP_IOIPSL
+! Ehouarn: dump possibility of using defrun
+!#ifdef CPP_IOIPSL
       CALL conf_gcm( 99, .TRUE. , clesphy0 )
-#else
-      CALL defrun( 99, .TRUE. , clesphy0 )
-#endif
+!#else
+!      CALL defrun( 99, .TRUE. , clesphy0 )
+!#endif
 c
 c
@@ -217,5 +195,10 @@
       call init_parallel
       call Read_Distrib
-      CALL Init_Phys_lmdz(iim,jjp1,llm,nqmx-2,mpi_size,distrib_phys)
+! Ehouarn : temporarily (?) keep this only for Earth
+      if (planet_type.eq."earth") then
+#ifdef CPP_EARTH
+        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
+#endif
+      endif ! of if (planet_type.eq."earth")
       CALL set_bands
       CALL Init_interface_dyn_phys
@@ -229,12 +212,17 @@
 c$OMP END PARALLEL
 
+! Ehouarn : temporarily (?) keep this only for Earth
+      if (planet_type.eq."earth") then
+#ifdef CPP_EARTH
 c$OMP PARALLEL
       call InitComgeomphy
 c$OMP END PARALLEL 
+#endif
+      endif ! of if (planet_type.eq."earth")
 
       IF (config_inca /= 'none') THEN
 #ifdef INCA
          call init_const_lmdz(
-     $        nbtrac,anneeref,dayref,
+     $        nbtr,anneeref,dayref,
      $        iphysiq,day_step,nday)
 
@@ -248,9 +236,11 @@
 c   Initialisation des traceurs
 c   ---------------------------
-c  Choix du schema pour l'advection
-c  dans fichier trac.def ou via INCA
-
-       call iniadvtrac(nq)
-c
+c  Choix du nombre de traceurs et du schema pour l'advection
+c  dans fichier traceur.def, par default ou via INCA
+      call infotrac_init
+
+c Allocation de la tableau q : champs advectes   
+      ALLOCATE(q(ip1jmp1,llm,nqtot))
+
 c-----------------------------------------------------------------------
 c   Lecture de l'etat initial :
@@ -259,7 +249,17 @@
 c  lecture du fichier start.nc
       if (read_start) then
-#ifdef CPP_IOIPSL
-         CALL dynetat0("start.nc",nqmx,vcov,ucov,
+      ! we still need to run iniacademic to initialize some
+      ! constants & fields, if we run the 'newtonian' case:
+        if (iflag_phys.eq.2) then
+          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
+        endif
+!#ifdef CPP_IOIPSL
+        if (planet_type.eq."earth") then
+#ifdef CPP_EARTH
+! Load an Earth-format start file
+         CALL dynetat0("start.nc",vcov,ucov,
      .              teta,q,masse,ps,phis, time_0)
+#endif
+        endif ! of if (planet_type.eq."earth")
 c       write(73,*) 'ucov',ucov
 c       write(74,*) 'vcov',vcov
@@ -268,14 +268,12 @@
 c       write(77,*) 'q',q
 
-#endif
-      endif
+      endif ! of if (read_start)
 
 c le cas echeant, creation d un etat initial
       IF (prt_level > 9) WRITE(lunout,*)
-     .                 'AVANT iniacademic AVANT AVANT AVANT AVANT'
+     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
       if (.not.read_start) then
-         CALL iniacademic(nqmx,vcov,ucov,teta,q,masse,ps,phis,time_0)
+         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
       endif
-
 
 c-----------------------------------------------------------------------
@@ -359,5 +357,4 @@
 c   Initialisation de la physique :
 c   -------------------------------
-#ifdef CPP_PHYS
       IF (call_iniphys.and.iflag_phys.eq.1) THEN
          latfi(1)=rlatu(1)
@@ -380,13 +377,14 @@
 
          WRITE(lunout,*)
-     .           'WARNING!!! vitesse verticale nulle dans la physique'
-
+     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
+! Earth:
+         if (planet_type.eq."earth") then
+#ifdef CPP_EARTH
          CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys ,
      ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
-
+#endif
+         endif ! of if (planet_type.eq."earth")
          call_iniphys=.false.
-
-      ENDIF
-#endif
+      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
 
 
@@ -404,4 +402,16 @@
 
 c-----------------------------------------------------------------------
+c   Initialisation des dimensions d'INCA :
+c   --------------------------------------
+      IF (config_inca /= 'none') THEN
+!$OMP PARALLEL
+#ifdef INCA
+         CALL init_inca_dim(klon_omp,llm,iim,jjm,
+     $        rlonu,rlatu,rlonv,rlatv)
+#endif
+!$OMP END PARALLEL
+      END IF
+
+c-----------------------------------------------------------------------
 c   Initialisation des I/O :
 c   ------------------------
@@ -410,10 +420,16 @@
       day_end = day_ini + nday
       WRITE(lunout,300)day_ini,day_end
+ 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
+
+!#ifdef CPP_IOIPSL
+      if (planet_type.eq."earth") then
+#ifdef CPP_EARTH
+      CALL dynredem0_p("restart.nc", day_end, phis)
+#endif
+      endif
+
+      ecripar = .TRUE.
 
 #ifdef CPP_IOIPSL
-      CALL dynredem0_p("restart.nc", day_end, phis, nqmx)
-
-      ecripar = .TRUE.
-
       if ( 1.eq.1) then
       time_step = zdtvr
@@ -421,11 +437,12 @@
       t_wrt = iecri * daysec
       CALL inithist_p(dynhist_file,day_ref,annee_ref,time_step,
-     .              t_ops, t_wrt, nqmx, histid, histvid)
-
-      t_ops = iperiod * time_step
-      t_wrt = periodav * daysec
-      CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
-     .              t_ops, t_wrt, nqmx, histaveid)
-
+     .              t_ops, t_wrt, histid, histvid)
+
+      IF (ok_dynzon) THEN 
+         t_ops = iperiod * time_step
+         t_wrt = periodav * daysec
+         CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
+     .        t_ops, t_wrt, histaveid)
+      END IF
       dtav = iperiod*dtvr/daysec
       endif
@@ -433,4 +450,5 @@
 
 #endif
+! #endif of #ifdef CPP_IOIPSL
 
 c  Choix des frequences de stokage pour le offline
@@ -453,11 +471,9 @@
 
 c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
-      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
+      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
      .              time_0)
 c$OMP END PARALLEL
 
 
- 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x,
-     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
       END
 
Index: LMDZ4/trunk/libf/dyn3dpar/groupeun_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/groupeun_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/groupeun_p.F	(revision 1146)
@@ -1,5 +1,5 @@
-      subroutine groupeun_p(jjmax,llmax,jjb,jje,q)
+      SUBROUTINE groupeun_p(jjmax,llmax,jjb,jje,q)
       USE parallel
-      implicit none
+      IMPLICIT NONE
 
 #include "dimensions.h"
@@ -8,70 +8,155 @@
 #include "comgeom2.h"
 
-      integer jjmax,llmax,jjb,jje
-      real q(iip1,jjmax,llmax)
+      INTEGER jjmax,llmax,jjb,jje
+      REAL q(iip1,jjmax,llmax)
 
-      integer ngroup
-      parameter (ngroup=3)
+      INTEGER ngroup
+      PARAMETER (ngroup=3)
 
-      real airen,airecn,qn
-      real aires,airecs,qs
+      REAL airecn,qn
+      REAL airecs,qs
 
-      integer i,j,l,ig,j1,j2,i0,jd
+      INTEGER i,j,l,ig,j1,j2,i0,jd
 
-Champs 3D
+c--------------------------------------------------------------------c 
+c Strategie d'optimisation                                           c
+c stocker les valeurs systematiquement recalculees                   c
+c et identiques d'un pas de temps sur l'autre. Il s'agit des         c
+c aires des cellules qui sont sommees. S'il n'y a pas de changement  c
+c de grille au cours de la simulation tout devrait bien se passer.   c
+c Autre optimisation : determination des bornes entre lesquelles "j" c
+c varie, au lieu de faire un test à chaque fois...
+c--------------------------------------------------------------------c 
+
+      INTEGER j_start, j_finish
+
+      REAL, SAVE :: airen_tab(iip1,jjp1,0:1)
+      REAL, SAVE :: aires_tab(iip1,jjp1,0:1)
+!$OMP THREADPRIVATE(airen_tab, aires_tab)
+
+      LOGICAL, SAVE :: first = .TRUE.
+!$OMP THREADPRIVATE(first)
+
+      IF (first) THEN
+         CALL INIT_GROUPEUN_P(airen_tab, aires_tab)
+         first = .FALSE.
+      ENDIF
+
+c Champs 3D
       jd=jjp1-jjmax
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      do l=1,llm
-      j1=1+jd
-      j2=2
-      do ig=1,ngroup
-         do j=j1-jd,j2-jd
-c           print*,'groupe ',ig,'  j= ',j,2**(ngroup-ig+1),'pts groupes'
-            if ( j >= jjb .AND. j <= jje) THEN
-             
-              do i0=1,iim,2**(ngroup-ig+1)
-                 
-                 airen=0.
-                 airecn=0.
-                 qn=0.
-                 
-                 do i=i0,i0+2**(ngroup-ig+1)-1
-                    airen=airen+aire(i,j)
-                    qn=qn+q(i,j,l)
-                 enddo
-                 airecn=0.
-                 do i=i0,i0+2**(ngroup-ig+1)-1
-                   q(i,j,l)=qn*aire(i,j)/airen
-                 enddo
-              enddo
-              q(iip1,j,l)=q(1,j,l)
-              
-            endif
+      DO l=1,llm
+         j1=1+jd
+         j2=2
+         DO ig=1,ngroup
+
+c     Concerne le pole nord
+            j_start  = MAX(jjb, j1-jd)
+            j_finish = MIN(jje, j2-jd)
+            DO j=j_start, j_finish
+               DO i0=1,iim,2**(ngroup-ig+1)
+                  qn=0.
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     qn=qn+q(i,j,l)
+                  ENDDO
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     q(i,j,l)=qn*airen_tab(i,j,jd)
+                  ENDDO
+               ENDDO
+               q(iip1,j,l)=q(1,j,l)
+            ENDDO
+        
+!c     Concerne le pole sud
+            j_start  = MAX(1+jjp1-jje-jd, j1-jd)
+            j_finish = MIN(1+jjp1-jjb-jd, j2-jd)
+            DO j=j_start, j_finish
+               DO i0=1,iim,2**(ngroup-ig+1)
+                  qs=0.
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     qs=qs+q(i,jjp1-j+1-jd,l)
+                  ENDDO
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     q(i,jjp1-j+1-jd,l)=qs*aires_tab(i,jjp1-j+1,jd)
+                  ENDDO
+               ENDDO
+               q(iip1,jjp1-j+1-jd,l)=q(1,jjp1-j+1-jd,l)
+            ENDDO
+        
+            j1=j2+1
+            j2=j2+2**ig
+         ENDDO
+      ENDDO
+!$OMP END DO NOWAIT
+
+      RETURN
+      END
+
+
+
+      SUBROUTINE INIT_GROUPEUN_P(airen_tab, aires_tab)
+
+      USE parallel
+      IMPLICIT NONE
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comgeom2.h"
+
+      INTEGER ngroup
+      PARAMETER (ngroup=3)
+
+      REAL airen,airecn
+      REAL aires,airecs
+
+      INTEGER i,j,l,ig,j1,j2,i0,jd
+
+      INTEGER j_start, j_finish
+
+      REAL :: airen_tab(iip1,jjp1,0:1)
+      REAL :: aires_tab(iip1,jjp1,0:1)
+
+      DO jd=0, 1
+         j1=1+jd
+         j2=2
+         DO ig=1,ngroup
             
-            if ( jjp1-j+1-jd >= jjb .AND. jjp1-j+1-jd <= jje) THEN
-             
-              do i0=1,iim,2**(ngroup-ig+1)
-                 aires=0.
-                 airecs=0.
-                 qs=0.
-                 do i=i0,i0+2**(ngroup-ig+1)-1
-                    aires=aires+aire(i,jjp1-j+1)
-                    qs=qs+q(i,jjp1-j+1-jd,l)
-                 enddo
-                 airecs=0.
-                 do i=i0,i0+2**(ngroup-ig+1)-1
-                   q(i,jjp1-j+1-jd,l)=qs*aire(i,jjp1-j+1)/aires
-                 enddo
-              enddo
-              q(iip1,jjp1-j+1-jd,l)=q(1,jjp1-j+1-jd,l)
-           
-            endif
-         enddo
-             
-           j1=j2+1
-           j2=j2+2**ig
-      enddo
-      enddo
-c$OMP END DO NOWAIT
-      return
-      end
+!     c     Concerne le pole nord
+            j_start = j1-jd
+            j_finish = j2-jd
+            DO j=j_start, j_finish
+               DO i0=1,iim,2**(ngroup-ig+1)
+                  airen=0.
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     airen = airen+aire(i,j)
+                  ENDDO
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     airen_tab(i,j,jd) = 
+     &                    aire(i,j) / airen
+                  ENDDO
+               ENDDO
+            ENDDO
+            
+!     c     Concerne le pole sud
+            j_start = j1-jd
+            j_finish = j2-jd
+            DO j=j_start, j_finish
+               DO i0=1,iim,2**(ngroup-ig+1)
+                  aires=0.
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     aires=aires+aire(i,jjp1-j+1)
+                  ENDDO
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     aires_tab(i,jjp1-j+1,jd) = 
+     &                    aire(i,jjp1-j+1) / aires
+                  ENDDO
+               ENDDO
+            ENDDO
+            
+            j1=j2+1
+            j2=j2+2**ig
+         ENDDO
+      ENDDO
+      
+      RETURN
+      END
Index: LMDZ4/trunk/libf/dyn3dpar/guide_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/guide_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/guide_p.F	(revision 1146)
@@ -4,4 +4,5 @@
       subroutine guide_pp(itau,ucov,vcov,teta,q,masse,ps)
       USE parallel
+      use netcdf
 
       IMPLICIT NONE
@@ -229,21 +230,26 @@
          IF (mpi_rank==0) THEN
           if (guide_modele) then
-            if (ncidpl.eq.-99) ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcod)
-          else	
-           if (guide_u) then
-            if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
-          endif
-c
-          if (guide_v) then
-            if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)
-          endif
-c
-          if (guide_T) then
-            if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)
-          endif
-c
-          if (guide_Q) then
-            if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)
-          endif
+             if (ncidpl.eq.-99) rcod=nf90_open('apbp.nc',Nf90_NOWRITe,
+     $            ncidpl)
+          else
+             if (guide_u) then
+                if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,
+     $               ncidpl)
+             endif
+c
+             if (guide_v) then
+                if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,
+     $               ncidpl)
+             endif
+c
+             if (guide_T) then
+                if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,
+     $               ncidpl)
+             endif
+c
+             if (guide_Q) then
+                if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite,
+     $               ncidpl)
+             endif
 c
           endif  !guide_modele
@@ -256,5 +262,5 @@
           status=NF_INQ_DIMLEN(ncidpl,rid,nlev)
          print *,'nlev guide', nlev 
-         call ncclos(ncidpl,rcod)
+         rcod = nf90_close(ncidpl)
 c   Lecture du premier etat des reanalyses.
          call Gather_Field(ps,ip1jmp1,1,0)
Index: LMDZ4/trunk/libf/dyn3dpar/infotrac.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/infotrac.F90	(revision 1146)
+++ LMDZ4/trunk/libf/dyn3dpar/infotrac.F90	(revision 1146)
@@ -0,0 +1,317 @@
+MODULE infotrac
+
+! nqtot : total number of tracers and higher order of moment, water vapor and liquid included
+  INTEGER, SAVE :: nqtot
+!!$OMP THREADPRIVATE(nqtot)   
+
+! nbtr : number of tracers not including higher order of moment or water vapor or liquid
+!        number of tracers used in the physics
+  INTEGER, SAVE :: nbtr
+!!$OMP THREADPRIVATE(nbtr)   
+
+! Name variables
+  CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics
+  CHARACTER(len=23), ALLOCATABLE, DIMENSION(:), SAVE :: ttext ! tracer long name for diagnostics
+!!$OMP THREADPRIVATE(tname,ttext)   
+
+! iadv  : index of trasport schema for each tracer
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: iadv
+!!$OMP THREADPRIVATE(iadv)   
+
+! niadv : vector keeping the coorspondance between all tracers(nqtot) treated in the 
+!         dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code. 
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
+!!$OMP THREADPRIVATE(niadv)   
+
+! Variables for INCA
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: conv_flg
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: pbl_flg
+!!$OMP THREADPRIVATE(conv_flg, pbl_flg)   
+
+CONTAINS
+
+  SUBROUTINE infotrac_init
+    IMPLICIT NONE
+!=======================================================================
+!
+!   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
+!   -------
+!   Modif special traceur F.Forget 05/94
+!   Modif M-A Filiberti 02/02 lecture de traceur.def
+!
+!   Objet:
+!   ------
+!   GCM LMD nouvelle grille
+!
+!=======================================================================
+!   ... modification de l'integration de q ( 26/04/94 ) ....
+!-----------------------------------------------------------------------
+! Declarations
+
+    INCLUDE "dimensions.h"
+    INCLUDE "control.h"
+    INCLUDE "iniprint.h"
+
+! Local variables
+    INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv  ! index of horizontal trasport schema
+    INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv  ! index of vertical trasport schema
+
+    CHARACTER(len=8), ALLOCATABLE, DIMENSION(:) :: tnom_0  ! tracer short name
+    CHARACTER(len=8), ALLOCATABLE, DIMENSION(:) :: tracnam ! name from INCA
+    CHARACTER(len=3), DIMENSION(30) :: descrq
+    CHARACTER(len=1), DIMENSION(3)  :: txts
+    CHARACTER(len=2), DIMENSION(9)  :: txtp
+    CHARACTER(len=13)               :: 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
+ 
+!-----------------------------------------------------------------------
+! Initialization :
+!
+    txts=(/'x','y','z'/)
+    txtp=(/'x ','y ','z ','xx','xy','xz','yy','yz','zz'/)
+
+    descrq(14)='VLH'
+    descrq(10)='VL1'
+    descrq(11)='VLP'
+    descrq(12)='FH1'
+    descrq(13)='FH2'
+    descrq(16)='PPM'
+    descrq(17)='PPS'
+    descrq(18)='PPP'
+    descrq(20)='SLP'
+    descrq(30)='PRA'
+
+!-----------------------------------------------------------------------
+!
+! 1) Get the true number of tracers + water vapor/liquid
+!    Here true tracers (nqtrue) means declared tracers (only first order)
+!
+!-----------------------------------------------------------------------
+    IF (config_inca == 'none') THEN
+       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
+       IF(ierr.EQ.0) THEN
+          WRITE(lunout,*) '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
+    ELSE
+       ! nbtr has been read from INCA by init_cont_lmdz() in gcm.F 
+       nqtrue=nbtr+2
+    END IF
+
+    IF (nqtrue < 2) THEN
+       WRITE(lunout,*) 'nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum'
+       CALL abort_gcm('infotrac_init','Not enough tracers',1)
+    END IF
+!
+! Allocate variables depending on nqtrue
+!
+    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue))
+
+    IF (config_inca /= 'none') THEN
+       ! Varaibles only needed in case of INCA
+       ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr))
+    END IF
+       
+!-----------------------------------------------------------------------
+! 2)     Choix  des schemas d'advection pour l'eau et les traceurs
+!
+!     iadv = 1    schema  transport type "humidite specifique LMD"
+!     iadv = 2    schema   amont
+!     iadv = 14   schema  Van-leer + humidite specifique 
+!                            Modif F.Codron
+!     iadv = 10   schema  Van-leer (retenu pour l'eau vapeur et liquide)
+!     iadv = 11   schema  Van-Leer pour hadv et version PPM (Monotone) pour vadv
+!     iadv = 12   schema  Frederic Hourdin I
+!     iadv = 13   schema  Frederic Hourdin II
+!     iadv = 16   schema  PPM Monotone(Collela & Woodward 1984)
+!     iadv = 17   schema  PPM Semi Monotone (overshoots autorisés)
+!     iadv = 18   schema  PPM Positif Defini (overshoots undershoots autorisés)
+!     iadv = 20   schema  Slopes
+!     iadv = 30   schema  Prather
+!
+!        Dans le tableau q(ij,l,iq) : iq = 1  pour l'eau vapeur
+!                                     iq = 2  pour l'eau liquide
+!       Et eventuellement             iq = 3,nqtot pour les autres traceurs
+!
+!        iadv(1): choix pour l'eau vap. et  iadv(2) : choix pour l'eau liq.
+!------------------------------------------------------------------------
+!
+!    Get choice of advection schema from file tracer.def or from INCA
+!---------------------------------------------------------------------
+    IF (config_inca == 'none') THEN
+       IF(ierr.EQ.0) THEN
+          ! Continue to read tracer.def
+          DO iq=1,nqtrue
+             READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
+          END DO
+          CLOSE(90)  
+       ELSE ! Without tracer.def
+          hadv(1) = 14
+          vadv(1) = 14
+          tnom_0(1) = 'H2Ov'
+          hadv(2) = 10
+          vadv(2) = 10
+          tnom_0(2) = 'H2Ol'
+          hadv(3) = 10
+          vadv(3) = 10
+          tnom_0(3) = 'RN'
+          hadv(4) = 10
+          vadv(4) = 10
+          tnom_0(4) = 'PB'
+       END IF
+       
+       WRITE(lunout,*) 'Valeur de traceur.def :'
+       WRITE(lunout,*) 'nombre de traceurs ',nqtrue
+       DO iq=1,nqtrue
+          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq)
+       END DO
+
+    ELSE  ! config_inca='aero' ou 'chem'
+! le module de chimie fournit les noms des traceurs
+! et les schemas d'advection associes.
+     
+#ifdef INCA
+       CALL init_transport( &
+            hadv, &
+            vadv, &
+            conv_flg, &
+            pbl_flg,  &
+            tracnam)
+#endif
+       tnom_0(1)='H2Ov'
+       tnom_0(2)='H2Ol'
+
+       DO iq =3,nqtrue
+          tnom_0(iq)=tracnam(iq-2)
+       END DO
+
+    END IF ! config_inca
+
+!-----------------------------------------------------------------------
+!
+! 3) Verify if advection schema 20 or 30 choosen
+!    Calculate total number of tracers needed: nqtot
+!    Allocate variables depending on total number of tracers
+!-----------------------------------------------------------------------
+    new_iq=0
+    DO iq=1,nqtrue
+       ! Add tracers for certain advection schema
+       IF (hadv(iq)<20 .AND. vadv(iq)<20 ) THEN
+          new_iq=new_iq+1  ! no tracers added
+       ELSE IF (hadv(iq)==20 .AND. vadv(iq)==20 ) THEN
+          new_iq=new_iq+4  ! 3 tracers added
+       ELSE IF (hadv(iq)==30 .AND. vadv(iq)==30 ) THEN
+          new_iq=new_iq+10 ! 9 tracers added
+       ELSE
+          WRITE(lunout,*) 'This choice of advection schema is not available'
+          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1)
+       END IF
+    END DO
+    
+    IF (new_iq /= nqtrue) THEN
+       ! The choice of advection schema imposes more tracers
+       ! Assigne total number of tracers
+       nqtot = new_iq
+
+       WRITE(lunout,*) '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'
+
+    ELSE
+       ! The true number of tracers is also the total number
+       nqtot = nqtrue
+    END IF
+
+!
+! Allocate variables with total number of tracers, nqtot
+!
+    ALLOCATE(tname(nqtot), ttext(nqtot))
+    ALLOCATE(iadv(nqtot), niadv(nqtot))
+
+!-----------------------------------------------------------------------
+!
+! 4) Determine iadv, long and short name
+!
+!-----------------------------------------------------------------------
+    new_iq=0
+    DO iq=1,nqtrue
+       new_iq=new_iq+1
+
+       ! Verify choice of advection schema
+       IF (hadv(iq)==vadv(iq)) THEN
+          iadv(new_iq)=hadv(iq)
+       ELSE IF (hadv(iq)==10 .AND. vadv(iq)==16) THEN
+          iadv(new_iq)=11
+       ELSE
+          WRITE(lunout,*)'This choice of advection schema is not available'
+          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1)
+       END IF
+      
+       str1=tnom_0(iq)
+       tname(new_iq)= tnom_0(iq)
+       IF (iadv(new_iq)==0) THEN
+          ttext(new_iq)=str1(1:lnblnk(str1))
+       ELSE
+          ttext(new_iq)=str1(1:lnblnk(str1))//descrq(iadv(new_iq))
+       END IF
+
+       ! schemas tenant compte des moments d'ordre superieur
+       str2=ttext(new_iq)
+       IF (iadv(new_iq)==20) THEN
+          DO jq=1,3
+             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)
+          END DO
+       ELSE IF (iadv(new_iq)==30) THEN
+          DO jq=1,9
+             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)
+          END DO
+       END IF
+    END DO
+
+!
+! Find vector keeping the correspodence between true and total tracers
+!
+    niadv(:)=0
+    iiq=0
+    DO iq=1,nqtot
+       IF(iadv(iq).GE.0) THEN
+          ! True tracer
+          iiq=iiq+1
+          niadv(iiq)=iq
+       ENDIF
+    END DO
+
+
+    WRITE(lunout,*) 'Information stored in dimtrac :'
+    WRITE(lunout,*) 'iadv  niadv tname  ttext :'
+    DO iq=1,nqtot
+       WRITE(lunout,*) iadv(iq),niadv(iq), tname(iq), ttext(iq)
+    END DO
+
+!-----------------------------------------------------------------------
+! Finalize :
+!
+    DEALLOCATE(tnom_0, hadv, vadv)
+    IF (config_inca /= 'none') DEALLOCATE(tracnam)
+
+999 FORMAT (i2,1x,i2,1x,a8)
+
+  END SUBROUTINE infotrac_init
+
+END MODULE infotrac
Index: LMDZ4/trunk/libf/dyn3dpar/iniacademic.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/iniacademic.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/iniacademic.F	(revision 1146)
@@ -4,5 +4,8 @@
 c
 c
-      SUBROUTINE iniacademic(nq,vcov,ucov,teta,q,masse,ps,phis,time_0)
+      SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
+
+      USE filtreg_mod
+      USE infotrac, ONLY : nqtot
 
 c%W%    %G%
@@ -42,9 +45,9 @@
 #include "temps.h"
 #include "control.h"
+#include "iniprint.h"
 
 c   Arguments:
 c   ----------
 
-      integer nq
       real time_0
 
@@ -52,21 +55,17 @@
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
       REAL teta(ip1jmp1,llm)                 ! temperature potentielle
-      REAL q(ip1jmp1,llm,nq)               ! champs advectes
+      REAL q(ip1jmp1,llm,nqtot)              ! champs advectes
       REAL ps(ip1jmp1)                       ! pression  au sol
       REAL masse(ip1jmp1,llm)                ! masse d'air
+      REAL phis(ip1jmp1)                     ! geopotentiel au sol
+
+c   Local:
+c   ------
+
       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
       REAL pks(ip1jmp1)                      ! exner au  sol
       REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
-      REAL phis(ip1jmp1)                     ! geopotentiel au sol
       REAL phi(ip1jmp1,llm)                  ! geopotentiel
-
-
-
-
-
-c   Local:
-c   ------
-
       REAL ddsin,tetarappelj,tetarappell,zsig
       real tetajl(jjp1,llm)
@@ -79,37 +78,38 @@
 
 c-----------------------------------------------------------------------
+! 1. Initializations for Earth-like case
+! --------------------------------------
+      if (planet_type=="earth") then
+c
+        time_0=0.
 
-c
-      time_0=0.
+        im         = iim
+        jm         = jjm
+        day_ini    = 0
+        omeg       = 4.*asin(1.)/86400.
+        rad    = 6371229.
+        g      = 9.8
+        daysec = 86400.
+        dtvr    = daysec/FLOAT(day_step)
+        zdtvr=dtvr
+        kappa  = 0.2857143
+        cpp    = 1004.70885
+        preff     = 101325.
+        pa        =  50000.
+        etot0      = 0.
+        ptot0      = 0.
+        ztot0      = 0.
+        stot0      = 0.
+        ang0       = 0.
 
-      im         = iim
-      jm         = jjm
-      day_ini    = 0
-      omeg       = 4.*asin(1.)/86400.
-      rad    = 6371229.
-      g      = 9.8
-      daysec = 86400.
-      dtvr    = daysec/FLOAT(day_step)
-      zdtvr=dtvr
-      kappa  = 0.2857143
-      cpp    = 1004.70885
-      preff     = 101325.
-      pa        =  50 000.
-      etot0      = 0.
-      ptot0      = 0.
-      ztot0      = 0.
-      stot0      = 0.
-      ang0       = 0.
-      pa         = 0.
+        CALL iniconst
+        CALL inigeom
+        CALL inifilr
 
-      CALL inicons0
-      CALL inigeom
-      CALL inifilr
-
-      ps=0.
-      phis=0.
+        ps=0.
+        phis=0.
 c---------------------------------------------------------------------
 
-      taurappel=10.*daysec
+        taurappel=10.*daysec
 
 c---------------------------------------------------------------------
@@ -117,63 +117,63 @@
 c   --------------------------------------
 
-      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
+        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
-      ENDDO
+           do j=1,jjp1
+             tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
+           enddo
+          endif ! of if (zsig.gt.0.3)
+        ENDDO ! of DO l=1,llm
 
-      do l=1,llm
-         do j=1,jjp1
-            do i=1,iip1
-               ij=(j-1)*iip1+i
-               tetarappel(ij,l)=tetajl(j,l)
-            enddo
-         enddo
-      enddo
+        do l=1,llm
+           do j=1,jjp1
+              do i=1,iip1
+                 ij=(j-1)*iip1+i
+                 tetarappel(ij,l)=tetajl(j,l)
+              enddo
+           enddo
+        enddo
 
-c     call dump2d(jjp1,llm,tetajl,'TEQ   ')
+c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
 
-      ps=1.e5
-      phis=0.
-      CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-      CALL massdair(p,masse)
+        ps=1.e5
+        phis=0.
+        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:nq)=0.
+        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 al\351atoire sur la temp\351rature
-      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
+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))
+           enddo
+        enddo
 
-      do l=1,llm
-         do ij=1,ip1jmp1,iip1
-            teta(ij+iim,l)=teta(ij,l)
-         enddo
-      enddo
+        do l=1,llm
+           do ij=1,ip1jmp1,iip1
+              teta(ij+iim,l)=teta(ij,l)
+           enddo
+        enddo
 
 
@@ -185,9 +185,14 @@
 
 c   initialisation d'un traceur sur une colonne
-      j=jjp1*3/4
-      i=iip1/2
-      ij=(j-1)*iip1+i
-      q(ij,:,3)=1.
-
+        j=jjp1*3/4
+        i=iip1/2
+        ij=(j-1)*iip1+i
+        q(ij,:,3)=1.
+      
+      else
+        write(lunout,*)"iniacademic: planet types other than earth",
+     &                 " not implemented (yet)."
+        stop
+      endif ! of if (planet_type=="earth")
       return
       END
Index: LMDZ4/trunk/libf/dyn3dpar/iniadvtrac.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/iniadvtrac.F	(revision 1122)
+++ 	(revision )
@@ -1,248 +1,0 @@
-!
-! $Header$
-!
-c
-c
-      subroutine iniadvtrac(nq)
-      USE ioipsl
-      IMPLICIT NONE
-c=======================================================================
-c
-c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
-c   -------
-c   Modif special traceur F.Forget 05/94
-c   Modif M-A Filiberti 02/02 lecture de traceur.def
-c
-c   Objet:
-c   ------
-c
-c   GCM LMD nouvelle grille
-c
-c=======================================================================
-c   ... modification de l'integration de q ( 26/04/94 ) ....
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
-C
-#include "dimensions.h"
-#include "advtrac.h"
-#include "control.h"
-
-c   local
-      character*3 descrq(30)
-      character*1 txts(3)
-      character*2 txtp(9)
-      character*13 str1,str2,str3
-
-      integer nq,iq,iiq,iiiq,ierr,ii
-      integer lnblnk
-      external lnblnk
-
-      data txts/'x','y','z'/
-      data txtp/'x','y','z','xx','xy','xz','yy','yz','zz'/
-
-c-----------------------------------------------------------------------
-c   Initialisations:
-c   ----------------
-      descrq(14)='VLH'
-      descrq(10)='VL1'
-      descrq(11)='VLP'
-      descrq(12)='FH1'
-      descrq(13)='FH2'
-      descrq(16)='PPM'
-      descrq(17)='PPS'
-      descrq(18)='PPP'
-      descrq(20)='SLP'
-      descrq(30)='PRA'
-
-      IF (config_inca /= 'none') THEN
-#ifdef INCA
-         CALL init_transport( 
-     $        hadv_flg,
-     $        vadv_flg,
-     $        conv_flg,
-     $        pbl_flg,
-     $        tracnam)
-#endif
-      END IF
-c-----------------------------------------------------------------------
-c        Choix  des schemas d'advection pour l'eau et les traceurs
-c
-c     iadv = 1    schema  transport type "humidite specifique LMD"
-c     iadv = 2    schema   amont
-c     iadv = 14    schema  Van-leer + humidite specifique 
-c                            Modif F.Codron
-c     iadv = 10   schema  Van-leer (retenu pour l'eau vapeur et liquide)
-c     iadv = 11   schema  Van-Leer pour hadv et version PPM (Monotone) pour vadv
-c     iadv = 12   schema  Frederic Hourdin I
-c     iadv = 13   schema  Frederic Hourdin II
-c     iadv = 16   schema  PPM Monotone(Collela & Woodward 1984)
-c     iadv = 17   schema  PPM Semi Monotone (overshoots autorisés)
-c     iadv = 18   schema  PPM Positif Defini (overshoots undershoots autorisés)
-c     iadv = 20   schema  Slopes
-c     iadv = 30   schema  Prather
-c
-c        Dans le tableau q(ij,l,iq) : iq = 1  pour l'eau vapeur
-c                                     iq = 2  pour l'eau liquide
-c        Et eventuellement            iq = 3,nqmx pour les autres traceurs
-c
-c        iadv(1): choix pour l'eau vap. et  iadv(2) : choix pour l'eau liq.
-C------------------------------------------------------------------------
-c     Choix du schema d'advection
-c------------------------------------------------------------------
-c choix par defaut = van leer pour tous les traceurs
-      do iq=1,nqmx
-       iadv(iq)=10
-       str1(1:1)='q'
-       if (nqmx.le.99) then
-       WRITE(str1(2:3),'(i2.2)') iq
-       else
-       WRITE(str1(2:4),'(i3.3)') iq
-       endif
-       tnom(iq)=str1
-       tname(iq)=tnom(iq) 
-       str2=tnom(iq) 
-       ttext(iq)=str2(1:lnblnk(str2))//descrq(iadv(iq))
-      end do
-      nq=nqmx
-c------------------------------------------------------------------
-c     Choix du schema pour l'advection
-c    dans fichier traceur.def
-c------------------------------------------------------------------
-
-      IF (config_inca /= 'none') THEN
-C le module de chimie fournit les noms des traceurs
-C et les schemas d'advection associes.
-         tnom(1)='H2Ov'
-         tnom(2)='H2Ol'
-         nq=nbtrac+2
-      
-       if (nq.gt.nqmx) then
-       print*,'nombre de traceurs incompatible INCA/LMDZT', nq, nbtrac
-       stop
-       endif
-      do iq =3,nq
-      tnom(iq)=tracnam(iq-2)
-      end do
-      do iq =1,nq
-      hadv(iq)= hadv_flg(iq) 
-      vadv(iq)= vadv_flg(iq) 
-      end do
-      ELSE  ! config_inca=none
-         print*,'ouverture de traceur.def'
-         open(90,file='traceur.def',form='formatted',status='old',
-     s        iostat=ierr)
-      if(ierr.eq.0) then
-        print*,'ouverture de traceur.def ok'
-        read(90,*) nq
-        print*,'nombre de traceurs ',nq
-        if (nq.gt.nqmx) then
-          print*,'nombre de traceurs trop important'
-          print*,'verifier traceur.def'
-          stop
-        endif
-C
-        do iq=1,nq
-          read(90,999) hadv(iq),vadv(iq),tnom(iq)
-        end do
-        close(90)  
-        PRINT*,'lecture de traceur.def :'   
-        do iq=1,nq
-          write(*,*) hadv(iq),vadv(iq),tnom(iq)
-        end do       
-      else
-        print*,'pb ouverture traceur.def'
-        print*,'ATTENTION on prend des valeurs par defaut'
-        nq = 4
-        hadv(1) = 14
-        vadv(1) = 14
-        tnom(1) = 'H2Ov'
-        hadv(2) = 10
-        vadv(2) = 10
-        tnom(2) = 'H2Ol'
-        hadv(3) = 10
-        vadv(3) = 10
-        tnom(3) = 'RN'
-        hadv(4) = 10
-        vadv(4) = 10
-        tnom(4) = 'PB'
-      ENDIF
-      PRINT*,'Valeur de traceur.def :'
-      do iq=1,nq
-        write(*,*) hadv(iq),vadv(iq),tnom(iq)
-      end do       
-
-      END IF ! config_inca
-
-c a partir du nom court du traceur et du schema d'advection au detemine le nom long.
-        iiq=0
-        ii=0
-        do iq=1,nq 
-         iiq=iiq+1
-         if (hadv(iq).ne.vadv(iq)) then
-           if (hadv(iq).eq.10.and.vadv(iq).eq.16) then
-             iadv(iiq)=11
-           else
-             print*,'le choix des schemas d''advection H et V'
-             print*, 'est non disponible actuellement'
-             stop 
-           endif
-         else
-          iadv(iiq)=hadv(iq)
-         endif
-c verification nombre de traceurs
-          if (iadv(iiq).lt.20) then
-             ii=ii+1
-          elseif (iadv(iiq).eq.20) then
-             ii=ii+4
-          elseif (iadv(iiq).eq.30) then
-             ii=ii+10
-          endif
- 
-         str1=tnom(iq)
-         tname(iiq)=tnom(iq)
-         IF (iadv(iiq).eq.0) THEN
-           ttext(iiq)=str1(1:lnblnk(str1))
-         ELSE
-           ttext(iiq)=str1(1:lnblnk(str1))//descrq(iadv(iiq))
-         ENDIF
-         str2=ttext(iiq)
-c   schemas tenant compte des moments d'ordre superieur.
-          if (iadv(iiq).eq.20) then
-             do iiiq=1,3
-               iiq=iiq+1
-               iadv(iiq)=-20
-               ttext(iiq)=str2(1:lnblnk(str2))//txts(iiiq)
-               tname(iiq)=str1(1:lnblnk(str1))//txts(iiiq)
-              enddo
-            elseif (iadv(iiq).eq.30) then
-              do iiiq=1,9
-               iiq=iiq+1
-               iadv(iiq)=-30
-               ttext(iiq)=str2(1:lnblnk(str2))//txtp(iiiq)
-               tname(iiq)=str1(1:lnblnk(str1))//txtp(iiiq)
-              enddo
-           endif
-        end do
-       if(ii.ne.nqmx) then
-       print*,'WARNING'
-       print*,'le nombre de traceurs et de moments eventuels'
-       print*,'est inferieur a nqmx '
-       endif
-       if (iiq.gt.nqmx) then
-       print*,'le choix des schemas est incompatible avec '
-       print*,'la dimension nqmx (nombre de traceurs)'
-       print*,'verifier traceur.def ou la namelist INCA'
-       print*,'ou recompiler avec plus de traceurs'
-       stop
-       endif
-      iiq=0
-      do iq=1,nqmx
-         if(iadv(iq).ge.0) then
-             iiq=iiq+1
-             niadv(iiq)=iq
-         endif
-      end do
-      return
-999   format (i2,1x,i2,1x,a8)
-      END
Index: LMDZ4/trunk/libf/dyn3dpar/inicons0.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/inicons0.F	(revision 1122)
+++ 	(revision )
@@ -1,55 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE inicons0
-      IMPLICIT NONE
-c
-c  P. Le Van
-c
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "temps.h"
-#include "control.h"
-#include "comvert.h"
-#include "comdissnew.h"
-
-
-c
-c
-c
-c-----------------------------------------------------------------------
-c   dimension des boucles:
-c   ----------------------
-
-      im      = iim
-      jm      = jjm
-      lllm    = llm
-      imp1    = iim 
-      jmp1    = jjm + 1
-      lllmm1  = llm - 1
-      lllmp1  = llm + 1
-
-c-----------------------------------------------------------------------
-
-      dtdiss  = idissip * dtvr
-      dtphys  = iphysiq * dtvr
-      unsim   = 1./iim
-      pi      = 2.*ASIN( 1. )
-
-c-----------------------------------------------------------------------
-      r       = cpp * kappa
-
-      PRINT*,' Cp R  kappa ',  cpp, r , kappa
-c
-c-----------------------------------------------------------------------
-
-       CALL disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
-c
-c
-       RETURN
-       END
Index: LMDZ4/trunk/libf/dyn3dpar/initdynav_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/initdynav_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/initdynav_p.F	(revision 1146)
@@ -4,6 +4,5 @@
 c
 c
-      subroutine initdynav_p(infile,day0,anne0,tstep,t_ops,t_wrt
-     .                     ,nq,fileid)
+      subroutine initdynav_p(infile,day0,anne0,tstep,t_ops,t_wrt,fileid)
 
        USE IOIPSL
@@ -11,4 +10,6 @@
        use Write_field
        use misc_mod
+       USE infotrac
+
       implicit none
 
@@ -30,5 +31,4 @@
 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:
@@ -50,5 +50,4 @@
 #include "description.h"
 #include "serre.h"
-#include "advtrac.h"
 
 C   Arguments
@@ -58,5 +57,4 @@
       real tstep, t_ops, t_wrt
       integer fileid
-      integer nq
       integer thoriid, zvertiid
 
@@ -82,5 +80,4 @@
       
       INTEGER :: dynave_domain_id
-      
       
       if (adjust) return
@@ -169,5 +166,5 @@
 C  Traceurs
 C
-        DO iq=1,nq
+        DO iq=1,nqtot
           call histdef(fileid, ttext(iq), ttext(iq), '-',
      .             iip1, jjn, thoriid, llm, 1, llm, zvertiid,
Index: LMDZ4/trunk/libf/dyn3dpar/initfluxsto_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/initfluxsto_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/initfluxsto_p.F	(revision 1146)
@@ -3,5 +3,5 @@
 !
       subroutine initfluxsto_p
-     .  (infile,tstep,t_ops,t_wrt,nq,
+     .  (infile,tstep,t_ops,t_wrt,
      .                    fileid,filevid,filedid)
 
@@ -30,5 +30,4 @@
 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:
@@ -58,5 +57,5 @@
       real tstep, t_ops, t_wrt
       integer fileid, filevid,filedid
-      integer nq,ndex(1)
+      integer ndex(1)
       real nivd(1)
 
@@ -87,5 +86,4 @@
       INTEGER :: dynu_domain_id
       INTEGER :: dynv_domain_id
-
 
 C
Index: LMDZ4/trunk/libf/dyn3dpar/inithist_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/inithist_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/inithist_p.F	(revision 1146)
@@ -2,5 +2,5 @@
 ! $Header$
 !
-      subroutine inithist_p(infile,day0,anne0,tstep,t_ops,t_wrt,nq,
+      subroutine inithist_p(infile,day0,anne0,tstep,t_ops,t_wrt,
      .                      fileid,filevid)
 
@@ -9,4 +9,5 @@
        use Write_field
        use misc_mod
+       USE infotrac
 
       implicit none
@@ -29,5 +30,4 @@
 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:
@@ -50,5 +50,4 @@
 #include "description.h"
 #include "serre.h"
-#include "advtrac.h"
 
 C   Arguments
@@ -58,5 +57,4 @@
       real tstep, t_ops, t_wrt
       integer fileid, filevid
-      integer nq
 
 C   Variables locales
@@ -83,4 +81,5 @@
       INTEGER :: dynu_domain_id
       INTEGER :: dynv_domain_id
+
 C
 C  Initialisations
@@ -217,5 +216,5 @@
 C  Traceurs
 C
-        DO iq=1,nq
+        DO iq=1,nqtot
           call histdef(fileid, ttext(iq),  ttext(iq), '-',
      .             iip1, jjn, thoriid, llm, 1, llm, zvertiid,
Index: LMDZ4/trunk/libf/dyn3dpar/integrd_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/integrd_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/integrd_p.F	(revision 1146)
@@ -32,5 +32,4 @@
 #include "temps.h"
 #include "serre.h"
-#include "advtrac.h"
 
 c   Arguments:
Index: LMDZ4/trunk/libf/dyn3dpar/ioipsl_getincom.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/ioipsl_getincom.F90	(revision 1146)
+++ LMDZ4/trunk/libf/dyn3dpar/ioipsl_getincom.F90	(revision 1146)
@@ -0,0 +1,2609 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Module and routines in this file are taken from IOIPSL
+! files getincom.f90
+! Module names has been changed to avoid problems
+! if compiling model with IOIPSL library
+! Ehouarn - March 2009
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+MODULE ioipsl_getincom
+!---------------------------------------------------------------------
+  USE ioipsl_stringop, &
+ &   ONLY : findpos,nocomma,cmpblank,strlowercase,gensig,find_sig
+!-
+  IMPLICIT NONE
+!-
+  PRIVATE
+  PUBLIC :: getin, getin_dump
+!-
+  INTERFACE getin
+    MODULE PROCEDURE getinrs, getinr1d, getinr2d, &
+ &                   getinis, getini1d, getini2d, &
+ &                   getincs, getinc1d, getinc2d, &
+ &                   getinls, getinl1d, getinl2d
+  END INTERFACE
+!-
+  INTEGER,PARAMETER :: max_files=100
+  CHARACTER(LEN=100),DIMENSION(max_files),SAVE :: filelist
+  INTEGER,SAVE      :: nbfiles
+!-
+  INTEGER,PARAMETER :: max_lines=4000
+  INTEGER,SAVE :: nb_lines
+  CHARACTER(LEN=100),DIMENSION(max_lines),SAVE :: fichier
+  INTEGER,DIMENSION(max_lines),SAVE :: targetsiglist,fromfile,compline
+  CHARACTER(LEN=30),DIMENSION(max_lines),SAVE  :: targetlist
+!-
+! The data base of parameters
+!-
+  INTEGER,PARAMETER :: memslabs=200
+  INTEGER,PARAMETER :: compress_lim = 20
+!-
+  INTEGER,SAVE :: nb_keys=0
+  INTEGER,SAVE :: keymemsize=0
+  INTEGER,SAVE,ALLOCATABLE :: keysig(:)
+  CHARACTER(LEN=30),SAVE,ALLOCATABLE :: keystr(:)
+!-
+! keystatus definition
+! keystatus = 1 : Value comes from run.def
+! keystatus = 2 : Default value is used
+! keystatus = 3 : Some vector elements were taken from default
+!-
+  INTEGER,SAVE,ALLOCATABLE :: keystatus(:)
+!-
+! keytype definition
+! keytype = 1 : Interger
+! keytype = 2 : Real
+! keytype = 3 : Character
+! keytype = 4 : Logical
+!-
+  INTEGER,SAVE,ALLOCATABLE :: keytype(:)
+!-
+! Allow compression for keys (only for integer and real)
+! keycompress < 0 : not compresses
+! keycompress > 0 : number of repeat of the value
+!-
+  INTEGER,SAVE,ALLOCATABLE :: keycompress(:)
+  INTEGER,SAVE,ALLOCATABLE :: keyfromfile(:)
+!-
+  INTEGER,SAVE,ALLOCATABLE :: keymemstart(:)
+  INTEGER,SAVE,ALLOCATABLE :: keymemlen(:)
+!-
+  INTEGER,SAVE,ALLOCATABLE :: intmem(:)
+  INTEGER,SAVE             :: intmemsize=0, intmempos=0
+  REAL,SAVE,ALLOCATABLE :: realmem(:)
+  INTEGER,SAVE          :: realmemsize=0, realmempos=0
+  CHARACTER(LEN=100),SAVE,ALLOCATABLE :: charmem(:)
+  INTEGER,SAVE             :: charmemsize=0, charmempos=0
+  LOGICAL,SAVE,ALLOCATABLE :: logicmem(:)
+  INTEGER,SAVE             :: logicmemsize=0, logicmempos=0
+!-
+CONTAINS
+!-
+!=== REAL INTERFACES
+!-
+SUBROUTINE getinrs (TARGET,ret_val)
+!---------------------------------------------------------------------
+!-  Get a real scalar. We first check if we find it
+!-  in the database and if not we get it from the run.def
+!-
+!-  getinr1d and getinr2d are written on the same pattern
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  REAL :: ret_val
+!-
+  REAL,DIMENSION(1) :: tmp_ret_val
+  INTEGER :: target_sig, pos, status=0, fileorig
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  tmp_ret_val(1) = ret_val
+!-
+  IF (pos < 0) THEN
+!-- Get the information out of the file
+    CALL getfilr (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwr (TARGET,target_sig,status,fileorig,1,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbrr (pos,1,TARGET,tmp_ret_val)
+  ENDIF
+  ret_val = tmp_ret_val(1)
+!---------------------
+END SUBROUTINE getinrs
+!-
+!===
+!-
+SUBROUTINE getinr1d (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- See getinrs for details. It is the same thing but for a vector
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  REAL,DIMENSION(:) :: ret_val
+!-
+  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: tmp_ret_val
+  INTEGER,SAVE :: tmp_ret_size = 0
+  INTEGER :: target_sig, pos, size_of_in, status=0, fileorig
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  size_of_in = SIZE(ret_val)
+  IF (.NOT.ALLOCATED(tmp_ret_val)) THEN
+    ALLOCATE (tmp_ret_val(size_of_in))
+  ELSE IF (size_of_in > tmp_ret_size) THEN
+    DEALLOCATE (tmp_ret_val)
+    ALLOCATE (tmp_ret_val(size_of_in))
+    tmp_ret_size = size_of_in
+  ENDIF
+  tmp_ret_val(1:size_of_in) = ret_val(1:size_of_in)
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfilr (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwr &
+ &   (TARGET,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbrr (pos,size_of_in,TARGET,tmp_ret_val)
+  ENDIF
+  ret_val(1:size_of_in) = tmp_ret_val(1:size_of_in)
+!----------------------
+END SUBROUTINE getinr1d
+!-
+!===
+!-
+SUBROUTINE getinr2d (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- See getinrs for details. It is the same thing but for a matrix
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  REAL,DIMENSION(:,:) :: ret_val
+!-
+  REAL,DIMENSION(:),ALLOCATABLE,SAVE :: tmp_ret_val
+  INTEGER,SAVE :: tmp_ret_size = 0
+  INTEGER :: target_sig,pos,size_of_in,size_1,size_2,status=0,fileorig
+  INTEGER :: jl, jj, ji
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  size_of_in = SIZE(ret_val)
+  size_1 = SIZE(ret_val,1)
+  size_2 = SIZE(ret_val,2)
+  IF (.NOT.ALLOCATED(tmp_ret_val)) THEN
+    ALLOCATE (tmp_ret_val(size_of_in))
+  ELSE IF (size_of_in > tmp_ret_size) THEN
+    DEALLOCATE (tmp_ret_val)
+    ALLOCATE (tmp_ret_val(size_of_in))
+    tmp_ret_size = size_of_in
+  ENDIF
+!-
+  jl=0
+  DO jj=1,size_2
+    DO ji=1,size_1
+      jl=jl+1
+      tmp_ret_val(jl) = ret_val(ji,jj)
+    ENDDO
+  ENDDO
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfilr (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwr &
+ &   (TARGET,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbrr (pos,size_of_in,TARGET,tmp_ret_val)
+  ENDIF
+!-
+  jl=0
+  DO jj=1,size_2
+    DO ji=1,size_1
+      jl=jl+1
+      ret_val(ji,jj) = tmp_ret_val(jl)
+    ENDDO
+  ENDDO
+!----------------------
+END SUBROUTINE getinr2d
+!-
+!===
+!-
+SUBROUTINE getfilr (TARGET,status,fileorig,ret_val)
+!---------------------------------------------------------------------
+!- Subroutine that will extract from the file the values
+!- attributed to the keyword target
+!-
+!- REALS
+!- -----
+!-
+!- target   : in  : CHARACTER(LEN=*)  target for which we will
+!-                                    look in the file
+!- status   : out : INTEGER tells us from where we obtained the data
+!- fileorig : out : The index of the file from which the key comes
+!- ret_val  : out : REAL(nb_to_ret) values read
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  INTEGER :: status, fileorig
+  REAL,DIMENSION(:) :: ret_val
+!-
+  INTEGER :: nb_to_ret
+  INTEGER :: it, pos, len_str, epos, ppos, int_tmp, status_cnt
+  CHARACTER(LEN=3)  :: cnt, tl, dl
+  CHARACTER(LEN=10) :: fmt
+  CHARACTER(LEN=30) :: full_target
+  CHARACTER(LEN=80) :: str_READ, str_READ_lower, str_tmp
+  INTEGER :: full_target_sig
+  REAL :: compvalue
+!-
+  INTEGER,SAVE :: max_len = 0
+  LOGICAL,SAVE,ALLOCATABLE :: found(:)
+  LOGICAL :: def_beha
+  LOGICAL :: compressed = .FALSE.
+!---------------------------------------------------------------------
+  nb_to_ret = SIZE(ret_val)
+  CALL getin_read
+!-
+! Get the variables and memory we need
+!-
+  IF (max_len == 0) THEN
+    ALLOCATE(found(nb_to_ret))
+    max_len = nb_to_ret
+  ENDIF
+  IF (max_len < nb_to_ret) THEN
+    DEALLOCATE(found)
+    ALLOCATE(found(nb_to_ret))
+    max_len = nb_to_ret
+  ENDIF
+  found(:) = .FALSE.
+!-
+! See what we find in the files read
+!-
+  DO it=1,nb_to_ret
+!---
+!-
+!-- First try the target as it is
+!---
+    full_target = TARGET(1:len_TRIM(target))
+    CALL gensig (full_target,full_target_sig)
+    CALL find_sig (nb_lines,targetlist,full_target, &
+ &                 targetsiglist,full_target_sig,pos)
+!---
+!-- Another try
+!---
+    IF (pos < 0) THEN
+      WRITE(cnt,'(I3.3)') it
+      full_target = TARGET(1:len_TRIM(target))//'__'//cnt
+      CALL gensig (full_target,full_target_sig)
+      CALL find_sig (nb_lines,targetlist,full_target, &
+ &                   targetsiglist,full_target_sig,pos)
+    ENDIF
+!---
+!-- A priori we dont know from which file the target could come.
+!-- Thus by default we attribute it to the first file :
+!---
+    fileorig = 1
+!--
+    IF (pos > 0) THEN
+!----
+      found(it) = .TRUE.
+      fileorig = fromfile(pos)
+!-----
+!---- DECODE
+!-----
+      str_READ = TRIM(ADJUSTL(fichier(pos)))
+      str_READ_lower = str_READ
+      CALL strlowercase (str_READ_lower)
+!----
+      IF (    (     (INDEX(str_READ_lower,'def') == 1)     &
+ &             .AND.(LEN_TRIM(str_READ_lower) == 3)   )    &
+ &        .OR.(     (INDEX(str_READ_lower,'default') == 1) &
+ &             .AND.(LEN_TRIM(str_READ_lower) == 7)   )   ) THEN
+        def_beha = .TRUE.
+      ELSE
+        def_beha = .FALSE.
+        len_str = LEN_TRIM(str_READ)
+        epos = INDEX(str_READ,'e')
+        ppos = INDEX(str_READ,'.')
+!------
+        IF (epos > 0) THEN
+          WRITE(tl,'(I3.3)') len_str
+          WRITE(dl,'(I3.3)') epos-ppos-1
+          fmt='(e'//tl//'.'//dl//')'
+          READ(str_READ,fmt) ret_val(it)
+        ELSE IF (ppos > 0) THEN
+          WRITE(tl,'(I3.3)') len_str
+          WRITE(dl,'(I3.3)') len_str-ppos
+          fmt='(f'//tl//'.'//dl//')'
+          READ(str_READ,fmt) ret_val(it)
+        ELSE
+          WRITE(tl,'(I3.3)') len_str
+          fmt = '(I'//tl//')'
+          READ(str_READ,fmt) int_tmp
+          ret_val(it) = REAL(int_tmp)
+        ENDIF
+      ENDIF
+!----
+      targetsiglist(pos) = -1
+!-----
+!---- Is this the value of a compressed field ?
+!-----
+      IF (compline(pos) > 0) THEN
+        IF (compline(pos) == nb_to_ret) THEN
+          compressed = .TRUE.
+          compvalue = ret_val(it)
+        ELSE
+          WRITE(*,*) 'WARNING from getfilr'
+          WRITE(*,*) 'For key ',TRIM(TARGET), &
+ & ' we have a compressed field but which does not have the right size.'
+          WRITE(*,*) 'We will try to fix that '
+          compressed = .TRUE.
+          compvalue = ret_val(it)
+        ENDIF
+      ENDIF
+    ELSE
+      found(it) = .FALSE.
+    ENDIF
+  ENDDO
+!--
+! If this is a compressed field then we will uncompress it
+!--
+  IF (compressed) THEN
+    DO it=1,nb_to_ret
+      IF (.NOT. found(it)) THEN
+        ret_val(it) = compvalue
+        found(it) = .TRUE.
+      ENDIF
+    ENDDO
+  ENDIF
+!-
+! Now we get the status for what we found
+!-
+  IF (def_beha) THEN
+    status = 2
+    WRITE(*,*) 'USING DEFAULT BEHAVIOUR FOR ',TRIM(TARGET)
+  ELSE
+    status_cnt = 0
+    DO it=1,nb_to_ret
+      IF (.NOT. found(it)) THEN
+        status_cnt = status_cnt+1
+        IF (nb_to_ret > 1) THEN
+          WRITE(str_tmp,'(a,"__",I3.3)') TRIM(TARGET),it
+        ELSE
+          str_tmp = TRIM(TARGET)
+        ENDIF
+        WRITE(*,*) 'USING DEFAULTS : ',TRIM(str_tmp),'=',ret_val(it)
+      ENDIF
+    ENDDO
+!---
+    IF (status_cnt == 0) THEN
+      status = 1
+    ELSE IF (status_cnt == nb_to_ret) THEN
+      status = 2
+    ELSE
+      status = 3
+    ENDIF
+  ENDIF
+!---------------------
+END SUBROUTINE getfilr
+!-
+!=== INTEGER INTERFACES
+!-
+SUBROUTINE getinis (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- Get a interer scalar. We first check if we find it
+!- in the database and if not we get it from the run.def
+!-
+!- getini1d and getini2d are written on the same pattern
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  INTEGER :: ret_val
+!-
+  INTEGER,DIMENSION(1) :: tmp_ret_val
+  INTEGER :: target_sig, pos, status=0, fileorig
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  tmp_ret_val(1) = ret_val
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfili (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwi (TARGET,target_sig,status,fileorig,1,tmp_ret_val)
+    ELSE
+!-- Get the value out of the database
+    CALL getdbri (pos,1,TARGET,tmp_ret_val)
+  ENDIF
+  ret_val = tmp_ret_val(1)
+!---------------------
+END SUBROUTINE getinis
+!-
+!===
+!-
+SUBROUTINE getini1d (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- See getinis for details. It is the same thing but for a vector
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  INTEGER,DIMENSION(:) :: ret_val
+!-
+  INTEGER,DIMENSION(:),ALLOCATABLE,SAVE :: tmp_ret_val
+  INTEGER,SAVE :: tmp_ret_size = 0
+  INTEGER :: target_sig, pos, size_of_in, status=0, fileorig
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  size_of_in = SIZE(ret_val)
+  IF (.NOT.ALLOCATED(tmp_ret_val)) THEN
+    ALLOCATE (tmp_ret_val(size_of_in))
+  ELSE IF (size_of_in > tmp_ret_size) THEN
+    DEALLOCATE (tmp_ret_val)
+    ALLOCATE (tmp_ret_val(size_of_in))
+    tmp_ret_size = size_of_in
+  ENDIF
+  tmp_ret_val(1:size_of_in) = ret_val(1:size_of_in)
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfili (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwi &
+ &   (TARGET,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbri (pos,size_of_in,TARGET,tmp_ret_val)
+  ENDIF
+  ret_val(1:size_of_in) = tmp_ret_val(1:size_of_in)
+!----------------------
+END SUBROUTINE getini1d
+!-
+!===
+!-
+SUBROUTINE getini2d (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- See getinis for details. It is the same thing but for a matrix
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  INTEGER,DIMENSION(:,:) :: ret_val
+!-
+  INTEGER,DIMENSION(:),ALLOCATABLE,SAVE :: tmp_ret_val
+  INTEGER,SAVE :: tmp_ret_size = 0
+  INTEGER :: target_sig,pos,size_of_in,size_1,size_2,status=0,fileorig
+  INTEGER :: jl, jj, ji
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  size_of_in = SIZE(ret_val)
+  size_1 = SIZE(ret_val,1)
+  size_2 = SIZE(ret_val,2)
+  IF (.NOT.ALLOCATED(tmp_ret_val)) THEN
+    ALLOCATE (tmp_ret_val(size_of_in))
+  ELSE IF (size_of_in > tmp_ret_size) THEN
+    DEALLOCATE (tmp_ret_val)
+    ALLOCATE (tmp_ret_val(size_of_in))
+    tmp_ret_size = size_of_in
+  ENDIF
+!-
+  jl=0
+  DO jj=1,size_2
+    DO ji=1,size_1
+      jl=jl+1
+      tmp_ret_val(jl) = ret_val(ji,jj)
+    ENDDO
+  ENDDO
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfili (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwi &
+ &   (TARGET,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbri (pos,size_of_in,TARGET,tmp_ret_val)
+  ENDIF
+!-
+  jl=0
+  DO jj=1,size_2
+    DO ji=1,size_1
+      jl=jl+1
+      ret_val(ji,jj) = tmp_ret_val(jl)
+    ENDDO
+  ENDDO
+!----------------------
+END SUBROUTINE getini2d
+!-
+!===
+!-
+SUBROUTINE getfili (TARGET,status,fileorig,ret_val)
+!---------------------------------------------------------------------
+!- Subroutine that will extract from the file the values
+!- attributed to the keyword target
+!-
+!- INTEGER
+!- -------
+!-
+!- target   : in  : CHARACTER(LEN=*)  target for which we will
+!-                                    look in the file
+!- status   : out : INTEGER tells us from where we obtained the data
+!- fileorig : out : The index of the file from which the key comes
+!- ret_val  : out : INTEGER(nb_to_ret) values read
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  INTEGER :: status, fileorig
+  INTEGER :: ret_val(:)
+!-
+  INTEGER :: nb_to_ret
+  INTEGER :: it, pos, len_str, status_cnt
+  CHARACTER(LEN=3)  :: cnt, chlen
+  CHARACTER(LEN=10) ::  fmt
+  CHARACTER(LEN=30) :: full_target
+  CHARACTER(LEN=80) :: str_READ, str_READ_lower, str_tmp
+  INTEGER :: full_target_sig
+  INTEGER :: compvalue
+!-
+  INTEGER,SAVE :: max_len = 0
+  LOGICAL,SAVE,ALLOCATABLE :: found(:)
+  LOGICAL :: def_beha
+  LOGICAL :: compressed = .FALSE.
+!---------------------------------------------------------------------
+  nb_to_ret = SIZE(ret_val)
+  CALL getin_read
+!-
+! Get the variables and memory we need
+!-
+  IF (max_len == 0) THEN
+    ALLOCATE(found(nb_to_ret))
+    max_len = nb_to_ret
+  ENDIF
+  IF (max_len < nb_to_ret) THEN
+    DEALLOCATE(found)
+    ALLOCATE(found(nb_to_ret))
+    max_len = nb_to_ret
+  ENDIF
+  found(:) = .FALSE.
+!-
+! See what we find in the files read
+!-
+  DO it=1,nb_to_ret
+!---
+!-- First try the target as it is
+!---
+    full_target = TARGET(1:len_TRIM(target))
+    CALL gensig (full_target,full_target_sig)
+    CALL find_sig (nb_lines,targetlist,full_target, &
+ &                 targetsiglist,full_target_sig,pos)
+!---
+!-- Another try
+!---
+    IF (pos < 0) THEN
+      WRITE(cnt,'(I3.3)') it
+      full_target = TARGET(1:len_TRIM(target))//'__'//cnt
+      CALL gensig (full_target,full_target_sig)
+      CALL find_sig (nb_lines,targetlist,full_target, &
+ &                   targetsiglist,full_target_sig,pos)
+    ENDIF
+!---
+!-- A priori we dont know from which file the target could come.
+!-- Thus by default we attribute it to the first file :
+!---
+    fileorig = 1
+!-
+    IF (pos > 0) THEN
+!-----
+      found(it) = .TRUE.
+      fileorig = fromfile(pos)
+!-----
+!---- DECODE
+!----
+      str_READ = TRIM(ADJUSTL(fichier(pos)))
+      str_READ_lower = str_READ
+      CALL strlowercase (str_READ_lower)
+!-----
+      IF (    (     (INDEX(str_READ_lower,'def') == 1)     &
+ &             .AND.(LEN_TRIM(str_READ_lower) == 3)   )    &
+ &        .OR.(     (INDEX(str_READ_lower,'default') == 1) &
+ &             .AND.(LEN_TRIM(str_READ_lower) == 7)   )   ) THEN
+        def_beha = .TRUE.
+      ELSE
+        def_beha = .FALSE.
+        len_str = LEN_TRIM(str_READ)
+        WRITE(chlen,'(I3.3)') len_str
+        fmt = '(I'//chlen//')'
+        READ(str_READ,fmt) ret_val(it)
+      ENDIF
+!-----
+      targetsiglist(pos) = -1
+!-----
+!---- Is this the value of a compressed field ?
+!-----
+      IF (compline(pos) > 0) THEN
+        IF (compline(pos) == nb_to_ret) THEN
+          compressed = .TRUE.
+          compvalue = ret_val(it)
+        ELSE
+          WRITE(*,*) 'WARNING from getfilr'
+          WRITE(*,*) 'For key ',TRIM(TARGET), &
+ & ' we have a compressed field but which does not have the right size.'
+          WRITE(*,*) 'We will try to fix that '
+          compressed = .TRUE.
+          compvalue = ret_val(it)
+        ENDIF
+      ENDIF
+    ELSE
+      found(it) = .FALSE.
+    ENDIF
+  ENDDO
+!-
+! If this is a compressed field then we will uncompress it
+!-
+  IF (compressed) THEN
+    DO it=1,nb_to_ret
+      IF (.NOT. found(it)) THEN
+        ret_val(it) = compvalue
+        found(it) = .TRUE.
+      ENDIF
+    ENDDO
+  ENDIF
+!-
+! Now we get the status for what we found
+!-
+  IF (def_beha) THEN
+    status = 2
+    WRITE(*,*) 'USING DEFAULT BEHAVIOUR FOR ',TRIM(TARGET)
+  ELSE
+    status_cnt = 0
+    DO it=1,nb_to_ret
+      IF (.NOT. found(it)) THEN
+        status_cnt = status_cnt+1
+        IF (nb_to_ret > 1) THEN
+          WRITE(str_tmp,'(a,"__",I3.3)') TRIM(TARGET),it
+        ELSE
+          str_tmp = TRIM(TARGET)
+        ENDIF
+        WRITE(*,*) 'USING DEFAULTS : ',TRIM(str_tmp),'=',ret_val(it)
+      ENDIF
+    ENDDO
+!---
+    IF (status_cnt == 0) THEN
+      status = 1
+    ELSE IF (status_cnt == nb_to_ret) THEN
+      status = 2
+    ELSE
+      status = 3
+    ENDIF
+  ENDIF
+!---------------------
+END SUBROUTINE getfili
+!-
+!=== CHARACTER INTERFACES
+!-
+SUBROUTINE getincs (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- Get a CHARACTER scalar. We first check if we find it
+!- in the database and if not we get it from the run.def
+!-
+!- getinc1d and getinc2d are written on the same pattern
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  CHARACTER(LEN=*) :: ret_val
+!-
+  CHARACTER(LEN=100),DIMENSION(1) :: tmp_ret_val
+  INTEGER :: target_sig, pos, status=0, fileorig
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  tmp_ret_val(1) = ret_val
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfilc (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwc (TARGET,target_sig,status,fileorig,1,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbrc (pos,1,TARGET,tmp_ret_val)
+  ENDIF
+  ret_val = tmp_ret_val(1)
+!---------------------
+END SUBROUTINE getincs
+!-
+!===
+!-
+SUBROUTINE getinc1d (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- See getincs for details. It is the same thing but for a vector
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  CHARACTER(LEN=*),DIMENSION(:) :: ret_val
+!-
+  CHARACTER(LEN=100),DIMENSION(:),ALLOCATABLE,SAVE :: tmp_ret_val
+  INTEGER,SAVE :: tmp_ret_size = 0
+  INTEGER :: target_sig, pos, size_of_in, status=0, fileorig
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  size_of_in = SIZE(ret_val)
+  IF (.NOT.ALLOCATED(tmp_ret_val)) THEN
+    ALLOCATE (tmp_ret_val(size_of_in))
+  ELSE IF (size_of_in > tmp_ret_size) THEN
+    DEALLOCATE (tmp_ret_val)
+    ALLOCATE (tmp_ret_val(size_of_in))
+    tmp_ret_size = size_of_in
+  ENDIF
+  tmp_ret_val(1:size_of_in) = ret_val(1:size_of_in)
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfilc (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwc &
+ &   (TARGET,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbrc (pos,size_of_in,TARGET,tmp_ret_val)
+  ENDIF
+  ret_val(1:size_of_in) = tmp_ret_val(1:size_of_in)
+!----------------------
+END SUBROUTINE getinc1d
+!-
+!===
+!-
+SUBROUTINE getinc2d (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- See getincs for details. It is the same thing but for a matrix
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  CHARACTER(LEN=*),DIMENSION(:,:) :: ret_val
+!-
+  CHARACTER(LEN=100),DIMENSION(:),ALLOCATABLE,SAVE :: tmp_ret_val
+  INTEGER,SAVE :: tmp_ret_size = 0
+  INTEGER :: target_sig,pos,size_of_in,size_1,size_2,status=0,fileorig
+  INTEGER :: jl,jj,ji
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  size_of_in = SIZE(ret_val)
+  size_1 = SIZE(ret_val,1)
+  size_2 = SIZE(ret_val,2)
+  IF (.NOT.ALLOCATED(tmp_ret_val)) THEN
+    ALLOCATE (tmp_ret_val(size_of_in))
+  ELSE IF (size_of_in > tmp_ret_size) THEN
+    DEALLOCATE (tmp_ret_val)
+    ALLOCATE (tmp_ret_val(size_of_in))
+    tmp_ret_size = size_of_in
+  ENDIF
+!-
+  jl=0
+  DO jj=1,size_2
+    DO ji=1,size_1
+      jl=jl+1
+      tmp_ret_val(jl) = ret_val(ji,jj)
+    ENDDO
+  ENDDO
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfilc (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwc &
+ &   (TARGET,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbrc (pos,size_of_in,TARGET,tmp_ret_val)
+  ENDIF
+!-
+  jl=0
+  DO jj=1,size_2
+    DO ji=1,size_1
+      jl=jl+1
+      ret_val(ji,jj) = tmp_ret_val(jl)
+    ENDDO
+  ENDDO
+!----------------------
+END SUBROUTINE getinc2d
+!-
+!===
+!-
+SUBROUTINE getfilc (TARGET,status,fileorig,ret_val)
+!---------------------------------------------------------------------
+!- Subroutine that will extract from the file the values
+!- attributed to the keyword target
+!-
+!- CHARACTER
+!- ---------
+!-
+!- target   : in  : CHARACTER(LEN=*)  target for which we will
+!-                                    look in the file
+!- status   : out : INTEGER tells us from where we obtained the data
+!- fileorig : out : The index of the file from which the key comes
+!- ret_val  : out : CHARACTER(nb_to_ret) values read
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+!-
+  CHARACTER(LEN=*) :: TARGET
+  INTEGER :: status, fileorig
+  CHARACTER(LEN=*),DIMENSION(:) :: ret_val
+!-
+  INTEGER :: nb_to_ret
+  INTEGER :: it, pos, len_str, status_cnt
+  CHARACTER(LEN=3)  :: cnt
+  CHARACTER(LEN=30) :: full_target
+  CHARACTER(LEN=80) :: str_READ, str_READ_lower, str_tmp
+  INTEGER :: full_target_sig
+!-
+  INTEGER,SAVE :: max_len = 0
+  LOGICAL,DIMENSION(:),SAVE,ALLOCATABLE :: found
+  LOGICAL :: def_beha
+!---------------------------------------------------------------------
+  nb_to_ret = SIZE(ret_val)
+  CALL getin_read
+!-
+! Get the variables and memory we need
+!-
+  IF (max_len == 0) THEN
+    ALLOCATE(found(nb_to_ret))
+    max_len = nb_to_ret
+  ENDIF
+  IF (max_len < nb_to_ret) THEN
+    DEALLOCATE(found)
+    ALLOCATE(found(nb_to_ret))
+    max_len = nb_to_ret
+  ENDIF
+  found(:) = .FALSE.
+!-
+! See what we find in the files read
+!-
+  DO it=1,nb_to_ret
+!---
+!-- First try the target as it is
+    full_target = TARGET(1:len_TRIM(target))
+    CALL gensig (full_target,full_target_sig)
+    CALL find_sig (nb_lines,targetlist,full_target, &
+ &                 targetsiglist,full_target_sig,pos)
+!---
+!-- Another try
+!---
+    IF (pos < 0) THEN
+      WRITE(cnt,'(I3.3)') it
+      full_target = TARGET(1:len_TRIM(target))//'__'//cnt
+      CALL gensig (full_target,full_target_sig)
+      CALL find_sig (nb_lines,targetlist,full_target, &
+ &                   targetsiglist,full_target_sig,pos)
+    ENDIF
+!---
+!-- A priori we dont know from which file the target could come.
+!-- Thus by default we attribute it to the first file :
+!---
+    fileorig = 1
+!---
+    IF (pos > 0) THEN
+!-----
+      found(it) = .TRUE.
+      fileorig = fromfile(pos)
+!-----
+!---- DECODE
+!-----
+      str_READ = TRIM(ADJUSTL(fichier(pos)))
+      str_READ_lower = str_READ
+      CALL strlowercase (str_READ_lower)
+!-----
+      IF (    (     (INDEX(str_READ_lower,'def') == 1)     &
+ &             .AND.(LEN_TRIM(str_READ_lower) == 3)   )    &
+ &        .OR.(     (INDEX(str_READ_lower,'default') == 1) &
+ &             .AND.(LEN_TRIM(str_READ_lower) == 7)   )   ) THEN
+        def_beha = .TRUE.
+      ELSE
+        def_beha = .FALSE.
+        len_str = LEN_TRIM(str_READ)
+        ret_val(it) = str_READ(1:len_str)
+      ENDIF
+!-----
+      targetsiglist(pos) = -1
+!-----
+    ELSE
+      found(it) = .FALSE.
+    ENDIF
+  ENDDO
+!-
+! Now we get the status for what we found
+!-
+  IF (def_beha) THEN
+    status = 2
+    WRITE(*,*) 'USING DEFAULT BEHAVIOUR FOR ',TRIM(TARGET)
+  ELSE
+    status_cnt = 0
+    DO it=1,nb_to_ret
+      IF (.NOT. found(it)) THEN
+        status_cnt = status_cnt+1
+        IF (nb_to_ret > 1) THEN
+          WRITE(str_tmp,'(a,"__",I3.3)') TRIM(TARGET),it
+        ELSE
+          str_tmp = TARGET(1:len_TRIM(target))
+        ENDIF
+        WRITE(*,*) 'USING DEFAULTS : ',TRIM(str_tmp),'=',ret_val(it)
+      ENDIF
+    ENDDO
+!-
+    IF (status_cnt == 0) THEN
+      status = 1
+    ELSE IF (status_cnt == nb_to_ret) THEN
+      status = 2
+    ELSE
+      status = 3
+    ENDIF
+  ENDIF
+!---------------------
+END SUBROUTINE getfilc
+!-
+!=== LOGICAL INTERFACES
+!-
+SUBROUTINE getinls (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- Get a logical scalar. We first check if we find it
+!- in the database and if not we get it from the run.def
+!-
+!- getinl1d and getinl2d are written on the same pattern
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  LOGICAL :: ret_val
+!-
+  LOGICAL,DIMENSION(1) :: tmp_ret_val
+  INTEGER :: target_sig, pos, status=0, fileorig
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  tmp_ret_val(1) = ret_val
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfill (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwl (TARGET,target_sig,status,fileorig,1,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbrl (pos,1,TARGET,tmp_ret_val)
+  ENDIF
+  ret_val = tmp_ret_val(1)
+!---------------------
+END SUBROUTINE getinls
+!-
+!===
+!-
+SUBROUTINE getinl1d (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- See getinls for details. It is the same thing but for a vector
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  LOGICAL,DIMENSION(:) :: ret_val
+!-
+  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: tmp_ret_val
+  INTEGER,SAVE :: tmp_ret_size = 0
+  INTEGER :: target_sig, pos, size_of_in, status=0, fileorig
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  size_of_in = SIZE(ret_val)
+  IF (.NOT.ALLOCATED(tmp_ret_val)) THEN
+    ALLOCATE (tmp_ret_val(size_of_in))
+  ELSE IF (size_of_in > tmp_ret_size) THEN
+    DEALLOCATE (tmp_ret_val)
+    ALLOCATE (tmp_ret_val(size_of_in))
+    tmp_ret_size = size_of_in
+  ENDIF
+  tmp_ret_val(1:size_of_in) = ret_val(1:size_of_in)
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfill (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwl &
+ &   (TARGET,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbrl (pos,size_of_in,TARGET,tmp_ret_val)
+  ENDIF
+  ret_val(1:size_of_in) = tmp_ret_val(1:size_of_in)
+!----------------------
+END SUBROUTINE getinl1d
+!-
+!===
+!-
+SUBROUTINE getinl2d (TARGET,ret_val)
+!---------------------------------------------------------------------
+!- See getinls for details. It is the same thing but for a matrix
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  LOGICAL,DIMENSION(:,:) :: ret_val
+!-
+  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: tmp_ret_val
+  INTEGER,SAVE :: tmp_ret_size = 0
+  INTEGER :: target_sig,pos,size_of_in,size_1,size_2,status=0,fileorig
+  INTEGER :: jl,jj,ji
+!---------------------------------------------------------------------
+!-
+! Compute the signature of the target
+!-
+  CALL gensig (TARGET,target_sig)
+!-
+! Do we have this target in our database ?
+!-
+  CALL find_sig (nb_keys,keystr,target,keysig,target_sig,pos)
+!-
+  size_of_in = SIZE(ret_val)
+  size_1 = SIZE(ret_val,1)
+  size_2 = SIZE(ret_val,2)
+  IF (.NOT.ALLOCATED(tmp_ret_val)) THEN
+    ALLOCATE (tmp_ret_val(size_of_in))
+  ELSE IF (size_of_in > tmp_ret_size) THEN
+    DEALLOCATE (tmp_ret_val)
+    ALLOCATE (tmp_ret_val(size_of_in))
+    tmp_ret_size = size_of_in
+  ENDIF
+!-
+  jl=0
+  DO jj=1,size_2
+    DO ji=1,size_1
+      jl=jl+1
+      tmp_ret_val(jl) = ret_val(ji,jj)
+    ENDDO
+  ENDDO
+!-
+  IF (pos < 0) THEN
+!-- Ge the information out of the file
+    CALL getfill (TARGET,status,fileorig,tmp_ret_val)
+!-- Put the data into the database
+    CALL getdbwl &
+ &   (TARGET,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+  ELSE
+!-- Get the value out of the database
+    CALL getdbrl (pos,size_of_in,TARGET,tmp_ret_val)
+  ENDIF
+!-
+  jl=0
+  DO jj=1,size_2
+    DO ji=1,size_1
+      jl=jl+1
+      ret_val(ji,jj) = tmp_ret_val(jl)
+    ENDDO
+  ENDDO
+!----------------------
+END SUBROUTINE getinl2d
+!-
+!===
+!-
+SUBROUTINE getfill (TARGET,status,fileorig,ret_val)
+!---------------------------------------------------------------------
+!- Subroutine that will extract from the file the values
+!- attributed to the keyword target
+!-
+!- LOGICAL
+!- -------
+!-
+!- target   : in  : CHARACTER(LEN=*)  target for which we will
+!-                                    look in the file
+!- status   : out : INTEGER tells us from where we obtained the data
+!- fileorig : out : The index of the file from which the key comes
+!- ret_val  : out : LOGICAL(nb_to_ret) values read
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: TARGET
+  INTEGER :: status, fileorig
+  LOGICAL,DIMENSION(:) :: ret_val
+!-
+  INTEGER :: nb_to_ret
+  INTEGER :: it, pos, len_str, ipos_tr, ipos_fl, status_cnt
+  CHARACTER(LEN=3)  :: cnt
+  CHARACTER(LEN=30) :: full_target
+  CHARACTER(LEN=80) :: str_READ, str_READ_lower, str_tmp
+  INTEGER :: full_target_sig
+!-
+  INTEGER,SAVE :: max_len = 0
+  LOGICAL,DIMENSION(:),SAVE,ALLOCATABLE :: found
+  LOGICAL :: def_beha
+!---------------------------------------------------------------------
+  nb_to_ret = SIZE(ret_val)
+  CALL getin_read
+!-
+! Get the variables and memory we need
+!-
+  IF (max_len == 0) THEN
+    ALLOCATE(found(nb_to_ret))
+    max_len = nb_to_ret
+  ENDIF
+  IF (max_len < nb_to_ret) THEN
+    DEALLOCATE(found)
+    ALLOCATE(found(nb_to_ret))
+    max_len = nb_to_ret
+  ENDIF
+  found(:) = .FALSE.
+!-
+! See what we find in the files read
+!-
+  DO it=1,nb_to_ret
+!---
+!-- First try the target as it is
+!---
+    full_target = TARGET(1:len_TRIM(target))
+    CALL gensig (full_target,full_target_sig)
+    CALL find_sig (nb_lines,targetlist,full_target, &
+ &                 targetsiglist,full_target_sig,pos)
+!---
+!-- Another try
+!---
+    IF (pos < 0) THEN
+      WRITE(cnt,'(I3.3)') it
+      full_target = TARGET(1:len_TRIM(target))//'__'//cnt
+      CALL gensig (full_target,full_target_sig)
+      CALL find_sig (nb_lines,targetlist,full_target, &
+ &                   targetsiglist,full_target_sig,pos)
+    ENDIF
+!---
+!-- A priori we dont know from which file the target could come.
+!-- Thus by default we attribute it to the first file :
+!---
+    fileorig = 1
+!---
+    IF (pos > 0) THEN
+!-----
+      found(it) = .TRUE.
+      fileorig = fromfile(pos)
+!-----
+!---- DECODE
+!-----
+      str_READ = TRIM(ADJUSTL(fichier(pos)))
+      str_READ_lower = str_READ
+      CALL strlowercase (str_READ_lower)
+!-----
+      IF (    (     (INDEX(str_READ_lower,'def') == 1)     &
+ &             .AND.(LEN_TRIM(str_READ_lower) == 3)   )    &
+ &        .OR.(     (INDEX(str_READ_lower,'default') == 1) &
+ &             .AND.(LEN_TRIM(str_READ_lower) == 7)   )   ) THEN
+        def_beha = .TRUE.
+      ELSE
+        def_beha = .FALSE.
+        len_str = LEN_TRIM(str_READ)
+        ipos_tr = -1
+        ipos_fl = -1
+!-------
+        ipos_tr = MAX(INDEX(str_READ,'tru'),INDEX(str_READ,'TRU'), &
+ &                    INDEX(str_READ,'y'),INDEX(str_READ,'Y'))
+        ipos_fl = MAX(INDEX(str_READ,'fal'),INDEX(str_READ,'FAL'), &
+ &                    INDEX(str_READ,'n'),INDEX(str_READ,'N'))
+!-------
+        IF (ipos_tr > 0) THEN
+          ret_val(it) = .TRUE.
+        ELSE IF (ipos_fl > 0) THEN
+          ret_val(it) = .FALSE.
+        ELSE
+          WRITE(*,*) "ERROR : getfill : TARGET ", &
+ &                   TRIM(TARGET)," is not of logical value"
+          STOP 'getinl'
+        ENDIF
+      ENDIF
+!-----
+      targetsiglist(pos) = -1
+!-----
+    ELSE
+!-
+      found(it) = .FALSE.
+!-
+    ENDIF
+!-
+  ENDDO
+!-
+! Now we get the status for what we found
+!-
+  IF (def_beha) THEN
+    status = 2
+    WRITE(*,*) 'USING DEFAULT BEHAVIOUR FOR ',TRIM(TARGET)
+  ELSE
+    status_cnt = 0
+    DO it=1,nb_to_ret
+      IF (.NOT. found(it)) THEN
+        status_cnt = status_cnt+1
+        IF (nb_to_ret > 1) THEN
+          WRITE(str_tmp,'(a,"__",I3.3)') TRIM(TARGET),it
+        ELSE
+          str_tmp = TRIM(TARGET)
+        ENDIF
+        WRITE(*,*) 'USING DEFAULTS : ',TRIM(str_tmp),'=',ret_val(it)
+      ENDIF
+    ENDDO
+!---
+    IF (status_cnt == 0) THEN
+      status = 1
+    ELSE IF (status_cnt == nb_to_ret) THEN
+      status = 2
+    ELSE
+      status = 3
+    ENDIF
+  ENDIF
+!---------------------
+END SUBROUTINE getfill
+!-
+!===
+!-
+SUBROUTINE getin_read
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER,SAVE :: allread=0
+  INTEGER,SAVE :: current,i
+!---------------------------------------------------------------------
+  IF (allread == 0) THEN
+!-- Allocate a first set of memory.
+    CALL getin_allockeys
+    CALL getin_allocmem (1,0)
+    CALL getin_allocmem (2,0)
+    CALL getin_allocmem (3,0)
+    CALL getin_allocmem (4,0)
+!-- Start with reading the files
+    nbfiles = 1
+    filelist(1) = 'run.def'
+    current = 1
+    nb_lines = 0
+!--
+    DO WHILE (current <= nbfiles)
+      CALL getin_readdef (current)
+      current = current+1
+    ENDDO
+    allread = 1
+    CALL getin_checkcohe ()
+  ENDIF
+!------------------------
+END SUBROUTINE getin_read
+!-
+!===
+!-
+  SUBROUTINE getin_readdef(current)
+!---------------------------------------------------------------------
+!- This subroutine will read the files and only keep the
+!- the relevant information. The information is kept as it
+!- found in the file. The data will be analysed later.
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: current
+!-
+  CHARACTER(LEN=100) :: READ_str, NEW_str, new_key, last_key, key_str
+  CHARACTER(LEN=3) :: cnt
+  INTEGER :: nb_lastkey
+!-
+  INTEGER :: eof, ptn, len_str, i, it, iund
+  LOGICAL :: check = .FALSE.
+!---------------------------------------------------------------------
+  eof = 0
+  ptn = 1
+  nb_lastkey = 0
+!-
+  IF (check) THEN
+    WRITE(*,*) 'getin_readdef : Open file ',TRIM(filelist(current))
+  ENDIF
+!-
+  OPEN (22,file=filelist(current),ERR=9997,STATUS="OLD")
+!-
+  DO WHILE (eof /= 1)
+!---
+    CALL getin_skipafew (22,READ_str,eof,nb_lastkey)
+    len_str = LEN_TRIM(READ_str)
+    ptn = INDEX(READ_str,'=')
+!---
+    IF (ptn > 0) THEN
+!---- Get the target
+      key_str = TRIM(ADJUSTL(READ_str(1:ptn-1)))
+!---- Make sure that if a vector keyword has the right length
+      iund =  INDEX(key_str,'__')
+      IF (iund > 0) THEN
+        SELECTCASE( len_trim(key_str)-iund )
+          CASE(2)
+            READ(key_str(iund+2:len_trim(key_str)),'(I1)') it
+          CASE(3)
+            READ(key_str(iund+2:len_trim(key_str)),'(I2)') it
+          CASE(4)
+            READ(key_str(iund+2:len_trim(key_str)),'(I3)') it
+          CASE DEFAULT
+            it = -1
+        END SELECT
+        IF (it > 0) THEN
+          WRITE(cnt,'(I3.3)') it
+          key_str = key_str(1:iund+1)//cnt
+        ELSE
+          WRITE(*,*) &
+ &          'getin_readdef : A very strange key has just been found'
+          WRITE(*,*) 'getin_readdef : ',key_str(1:len_TRIM(key_str))
+          STOP 'getin_readdef'
+        ENDIF
+      ENDIF
+!---- Prepare the content
+      NEW_str = TRIM(ADJUSTL(READ_str(ptn+1:len_str)))
+      CALL nocomma (NEW_str)
+      CALL cmpblank (NEW_str)
+      NEW_str  = TRIM(ADJUSTL(NEW_str))
+      IF (check) THEN
+        WRITE(*,*) &
+ &        '--> getin_readdef : ',TRIM(key_str),' :: ',TRIM(NEW_str)
+      ENDIF
+!---- Decypher the content of NEW_str
+!-
+!---- This has to be a new key word, thus :
+      nb_lastkey = 0
+!----
+      CALL getin_decrypt (current,key_str,NEW_str,last_key,nb_lastkey)
+!----
+    ELSE IF (len_str > 0) THEN
+!---- Prepare the key if we have an old one to which
+!---- we will add the line just read
+      IF (nb_lastkey > 0) THEN
+        iund =  INDEX(last_key,'__')
+        IF (iund > 0) THEN
+!-------- We only continue a keyword, thus it is easy
+          key_str = last_key(1:iund-1)
+        ELSE
+          IF (nb_lastkey /= 1) THEN
+            WRITE(*,*) &
+ &   'getin_readdef : An error has occured. We can not have a scalar'
+            WRITE(*,*) 'getin_readdef : keywod and a vector content'
+            STOP 'getin_readdef'
+          ENDIF
+!-------- The last keyword needs to be transformed into a vector.
+          targetlist(nb_lines) = &
+ &          last_key(1:MIN(len_trim(last_key),30))//'__001'
+          CALL gensig (targetlist(nb_lines),targetsiglist(nb_lines))
+          key_str = last_key(1:len_TRIM(last_key))
+        ENDIF
+      ENDIF
+!---- Prepare the content
+      NEW_str = TRIM(ADJUSTL(READ_str(1:len_str)))
+      CALL getin_decrypt (current,key_str,NEW_str,last_key,nb_lastkey)
+    ELSE
+!---- If we have an empty line the the keyword finishes
+      nb_lastkey = 0
+      IF (check) THEN
+        WRITE(*,*) 'getin_readdef : Have found an emtpy line '
+      ENDIF
+    ENDIF
+  ENDDO
+!-
+  CLOSE(22)
+!-
+  IF (check) THEN
+    OPEN (22,file='run.def.test')
+    DO i=1,nb_lines
+      WRITE(22,*) targetlist(i)," : ",fichier(i)
+    ENDDO
+    CLOSE(22)
+  ENDIF
+!-
+  RETURN
+!-
+9997 WRITE(*,*) "getin_readdef : Could not open file ", &
+          & TRIM(filelist(current))
+!---------------------------
+END SUBROUTINE getin_readdef
+!-
+!===
+!-
+SUBROUTINE getin_decrypt(current,key_str,NEW_str,last_key,nb_lastkey)
+!---------------------------------------------------------------------
+!- This subroutine is going to decypher the line.
+!- It essentialy checks how many items are included and
+!- it they can be attached to a key.
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+! ARGUMENTS
+!-
+  INTEGER :: current, nb_lastkey
+  CHARACTER(LEN=*) :: key_str, NEW_str, last_key
+!-
+! LOCAL
+!-
+  INTEGER :: len_str, blk, nbve, starpos
+  CHARACTER(LEN=100) :: tmp_str, new_key, mult
+  CHARACTER(LEN=3)   :: cnt, chlen
+  CHARACTER(LEN=10)  :: fmt
+!---------------------------------------------------------------------
+  len_str = LEN_TRIM(NEW_str)
+  blk = INDEX(NEW_str(1:len_str),' ')
+  tmp_str = NEW_str(1:len_str)
+!-
+! If the key is a new file then we take it up. Else
+! we save the line and go on.
+!-
+  IF (INDEX(key_str,'INCLUDEDEF') > 0) THEN
+    DO WHILE (blk > 0)
+      IF (nbfiles+1 > max_files) THEN
+        WRITE(*,*) 'FATAL ERROR : Too many files to include'
+        STOP 'getin_readdef'
+      ENDIF
+!-----
+      nbfiles = nbfiles+1
+      filelist(nbfiles) = tmp_str(1:blk)
+!-----
+      tmp_str = TRIM(ADJUSTL(tmp_str(blk+1:LEN_TRIM(tmp_str))))
+      blk = INDEX(tmp_str(1:LEN_TRIM(tmp_str)),' ')
+    ENDDO
+!---
+    IF (nbfiles+1 > max_files) THEN
+      WRITE(*,*) 'FATAL ERROR : Too many files to include'
+      STOP 'getin_readdef'
+    ENDIF
+!---
+    nbfiles =  nbfiles+1
+    filelist(nbfiles) = TRIM(ADJUSTL(tmp_str))
+!---
+    last_key = 'INCLUDEDEF'
+    nb_lastkey = 1
+  ELSE
+!-
+!-- We are working on a new line of input
+!-
+    nb_lines = nb_lines+1
+    IF (nb_lines > max_lines) THEN
+      WRITE(*,*) &
+ &      'Too many line in the run.def files. You need to increase'
+      WRITE(*,*) 'the parameter max_lines in the module getincom.'
+      STOP 'getin_decrypt'
+    ENDIF
+!-
+!-- First we solve the issue of conpressed information. Once
+!-- this is done all line can be handled in the same way.
+!-
+    starpos = INDEX(NEW_str(1:len_str),'*')
+    IF ( (starpos > 0).AND.(tmp_str(1:1) /= '"') &
+ &                    .AND.(tmp_str(1:1) /= "'") ) THEN
+!-----
+      IF (INDEX(key_str(1:len_TRIM(key_str)),'__') > 0) THEN
+        WRITE(*,*) 'ERROR : getin_decrypt'
+        WRITE(*,*) &
+&         'We can not have a compressed field of values for in a'
+        WRITE(*,*) &
+&         'vector notation. If a target is of the type TARGET__1'
+        WRITE(*,*) 'then only a scalar value is allowed'
+        WRITE(*,*) 'The key at fault : ',key_str(1:len_TRIM(key_str))
+        STOP 'getin_decrypt'
+      ENDIF
+!-
+!---- Read the multiplied
+!-
+      mult = TRIM(ADJUSTL(NEW_str(1:starpos-1)))
+!---- Construct the new string and its parameters
+      NEW_str = TRIM(ADJUSTL(NEW_str(starpos+1:len_str)))
+      len_str = LEN_TRIM(NEW_str)
+      blk = INDEX(NEW_str(1:len_str),' ')
+      IF (blk > 1) THEN
+        WRITE(*,*) &
+ &       'This is a strange behavior of getin_decrypt you could report'
+      ENDIF
+      WRITE(chlen,'(I3.3)') LEN_TRIM(mult)
+      fmt = '(I'//chlen//')'
+      READ(mult,fmt) compline(nb_lines)
+!---
+    ELSE
+      compline(nb_lines) = -1
+    ENDIF
+!-
+!-- If there is no space wthin the line then the target is a scalar
+!-- or the element of a properly written vector.
+!-- (ie of the type TARGET__1)
+!-
+    IF (    (blk <= 1) &
+ &      .OR.(tmp_str(1:1) == '"') &
+ &      .OR.(tmp_str(1:1) == "'") ) THEN
+!-
+      IF (nb_lastkey == 0) THEN
+!------ Save info of current keyword as a scalar
+!------ if it is not a continuation
+        targetlist(nb_lines) = key_str(1:MIN(len_trim(key_str),30))
+        last_key = key_str(1:MIN(len_trim(key_str),30))
+        nb_lastkey = 1
+      ELSE
+!------ We are continuing a vector so the keyword needs
+!------ to get the underscores
+        WRITE(cnt,'(I3.3)') nb_lastkey+1
+        targetlist(nb_lines) = &
+ &        key_str(1:MIN(len_trim(key_str),25))//'__'//cnt
+        last_key = key_str(1:MIN(len_trim(key_str),25))//'__'//cnt
+        nb_lastkey = nb_lastkey+1
+      ENDIF
+!-----
+      CALL gensig (targetlist(nb_lines),targetsiglist(nb_lines))
+      fichier(nb_lines) = NEW_str(1:len_str)
+      fromfile(nb_lines) = current
+    ELSE
+!-
+!---- If there are blanks whithin the line then we are dealing
+!---- with a vector and we need to split it in many entries
+!---- with the TRAGET__1 notation.
+!----
+!---- Test if the targer is not already a vector target !
+!-
+      IF (INDEX(TRIM(key_str),'__') > 0) THEN
+        WRITE(*,*) 'ERROR : getin_decrypt'
+        WRITE(*,*) 'We have found a mixed vector notation'
+        WRITE(*,*) 'If a target is of the type TARGET__1'
+        WRITE(*,*) 'then only a scalar value is allowed'
+        WRITE(*,*) 'The key at fault : ',key_str(1:len_TRIM(key_str))
+        STOP 'getin_decrypt'
+      ENDIF
+!-
+      nbve = nb_lastkey
+      nbve = nbve+1
+      WRITE(cnt,'(I3.3)') nbve
+!-
+      DO WHILE (blk > 0)
+!-
+!------ Save the content of target__nbve
+!-
+        fichier(nb_lines) = tmp_str(1:blk)
+        new_key =  key_str(1:MIN(len_trim(key_str),25))//'__'//cnt
+        targetlist(nb_lines) = new_key(1:MIN(len_trim(new_key),30))
+        CALL gensig (targetlist(nb_lines),targetsiglist(nb_lines))
+        fromfile(nb_lines) = current
+!-
+        tmp_str = TRIM(ADJUSTL(tmp_str(blk+1:LEN_TRIM(tmp_str))))
+        blk = INDEX(TRIM(tmp_str),' ')
+!-
+        nb_lines = nb_lines+1
+        IF (nb_lines > max_lines) THEN
+          WRITE(*,*) &
+ &          'Too many line in the run.def files. You need to increase'
+          WRITE(*,*) 'the parameter max_lines in the module getincom.'
+          STOP 'getin_decrypt'
+        ENDIF
+        nbve = nbve+1
+        WRITE(cnt,'(I3.3)') nbve
+!-
+      ENDDO
+!-
+!---- Save the content of the last target
+!-
+      fichier(nb_lines) = tmp_str(1:LEN_TRIM(tmp_str))
+      new_key = key_str(1:MIN(len_trim(key_str),25))//'__'//cnt
+      targetlist(nb_lines) = new_key(1:MIN(len_trim(new_key),30))
+      CALL gensig (targetlist(nb_lines),targetsiglist(nb_lines))
+      fromfile(nb_lines) = current
+!-
+      last_key = key_str(1:MIN(len_trim(key_str),25))//'__'//cnt
+      nb_lastkey = nbve
+!-
+    ENDIF
+!-
+  ENDIF
+!---------------------------
+END SUBROUTINE getin_decrypt
+!-
+!===
+!-
+SUBROUTINE getin_checkcohe ()
+!---------------------------------------------------------------------
+!- This subroutine checks for redundancies.
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+! Arguments
+!-
+!-
+! LOCAL
+!-
+  INTEGER :: line,i,sig
+  INTEGER :: found
+  CHARACTER(LEN=30) :: str
+!---------------------------------------------------------------------
+  DO line=1,nb_lines-1
+!-
+    CALL find_sig &
+ &    (nb_lines-line,targetlist(line+1:nb_lines),targetlist(line), &
+ &     targetsiglist(line+1:nb_lines),targetsiglist(line),found)
+!---
+!-- IF we have found it we have a problem to solve.
+!---
+    IF (found > 0) THEN
+      WRITE(*,*) 'COUNT : ', &
+ &  COUNT(ABS(targetsiglist(line+1:nb_lines)-targetsiglist(line)) < 1)
+!-----
+      WRITE(*,*) &
+ & 'getin_checkcohe : Found a problem on key ',targetlist(line)
+      WRITE(*,*) &
+ & 'getin_checkcohe : The following values were encoutered :'
+      WRITE(*,*) &
+ & '                ',TRIM(targetlist(line)), &
+ &               targetsiglist(line),' == ',fichier(line)
+      WRITE(*,*) &
+ & '                ',TRIM(targetlist(line+found)), &
+ &               targetsiglist(line+found),' == ',fichier(line+found)
+      WRITE(*,*) &
+ & 'getin_checkcohe : We will keep only the last value'
+!-----
+      targetsiglist(line) = 1
+    ENDIF
+  ENDDO
+!-
+END SUBROUTINE getin_checkcohe
+!-
+!===
+!-
+SUBROUTINE getin_skipafew (unit,out_string,eof,nb_lastkey)
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: unit, eof, nb_lastkey
+  CHARACTER(LEN=100) :: dummy
+  CHARACTER(LEN=100) :: out_string
+  CHARACTER(LEN=1) :: first
+!---------------------------------------------------------------------
+  first="#"
+  eof = 0
+  out_string = "    "
+!-
+  DO WHILE (first == "#")
+    READ (unit,'(a100)',ERR=9998,END=7778) dummy
+    dummy = TRIM(ADJUSTL(dummy))
+    first=dummy(1:1)
+    IF (first == "#") THEN
+      nb_lastkey = 0
+    ENDIF
+  ENDDO
+  out_string=dummy
+!-
+  RETURN
+!-
+9998 WRITE(*,*) " GETIN_SKIPAFEW : Error while reading file "
+  STOP 'getin_skipafew'
+!-
+7778 eof = 1
+!----------------------------
+END SUBROUTINE getin_skipafew
+!-
+!=== INTEGER database INTERFACE
+!-
+SUBROUTINE getdbwi &
+ &  (target,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+!---------------------------------------------------------------------
+!- Write the INTEGER data into the data base
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: target
+  INTEGER :: target_sig, status, fileorig, size_of_in
+  INTEGER,DIMENSION(:) :: tmp_ret_val
+!---------------------------------------------------------------------
+!-
+! First check if we have sufficiant space for the new key
+!-
+  IF (nb_keys+1 > keymemsize) THEN
+    CALL getin_allockeys ()
+  ENDIF
+!-
+! Fill out the items of the data base
+!-
+  nb_keys = nb_keys+1
+  keysig(nb_keys) = target_sig
+  keystr(nb_keys) = target(1:MIN(len_trim(target),30))
+  keystatus(nb_keys) = status
+  keytype(nb_keys) = 1
+  keyfromfile(nb_keys) = fileorig
+!-
+! Can we compress the data base entry ?
+!-
+  IF (     (MINVAL(tmp_ret_val) == MAXVAL(tmp_ret_val)) &
+ &    .AND.(size_of_in > compress_lim)) THEN
+    keymemstart(nb_keys) = intmempos+1
+    keycompress(nb_keys) = size_of_in
+    keymemlen(nb_keys) = 1
+  ELSE
+    keymemstart(nb_keys) = intmempos+1
+    keycompress(nb_keys) = -1
+    keymemlen(nb_keys) = size_of_in
+  ENDIF
+!-
+! Before writing the actual size lets see if we have the space
+!-
+  IF (keymemstart(nb_keys)+keymemlen(nb_keys) > intmemsize) THEN
+    CALL getin_allocmem (1,keymemlen(nb_keys))
+  ENDIF
+!-
+  intmem(keymemstart(nb_keys): &
+ &       keymemstart(nb_keys)+keymemlen(nb_keys)-1) = &
+ &  tmp_ret_val(1:keymemlen(nb_keys))
+  intmempos = keymemstart(nb_keys)+keymemlen(nb_keys)-1
+!---------------------
+END SUBROUTINE getdbwi
+!-
+!===
+!-
+SUBROUTINE getdbri (pos,size_of_in,target,tmp_ret_val)
+!---------------------------------------------------------------------
+!- Read the required variables in the database for INTEGERS
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: pos, size_of_in
+  CHARACTER(LEN=*) :: target
+  INTEGER,DIMENSION(:) :: tmp_ret_val
+!---------------------------------------------------------------------
+  IF (keytype(pos) /= 1) THEN
+    WRITE(*,*) 'FATAL ERROR : Wrong data type for keyword ',target
+    STOP 'getdbri'
+  ENDIF
+!-
+  IF (keycompress(pos) > 0) THEN
+    IF ( keycompress(pos) /= size_of_in .OR. keymemlen(pos) /= 1 ) THEN
+      WRITE(*,*) &
+ &      'FATAL ERROR : Wrong compression length for keyword ',target
+      STOP 'getdbri'
+    ELSE
+      tmp_ret_val(1:size_of_in) = intmem(keymemstart(pos))
+    ENDIF
+  ELSE
+    IF (keymemlen(pos) /= size_of_in) THEN
+      WRITE(*,*) 'FATAL ERROR : Wrong array length for keyword ',target
+      STOP 'getdbri'
+    ELSE
+      tmp_ret_val(1:size_of_in) = &
+ &      intmem(keymemstart(pos):keymemstart(pos)+keymemlen(pos)-1)
+    ENDIF
+  ENDIF
+!---------------------
+END SUBROUTINE getdbri
+!-
+!=== REAL database INTERFACE
+!-
+SUBROUTINE getdbwr &
+ &  (target,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+!---------------------------------------------------------------------
+!- Write the REAL data into the data base
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: target
+  INTEGER :: target_sig, status, fileorig, size_of_in
+  REAL,DIMENSION(:) :: tmp_ret_val
+!---------------------------------------------------------------------
+!-
+! First check if we have sufficiant space for the new key
+!-
+  IF (nb_keys+1 > keymemsize) THEN
+    CALL getin_allockeys ()
+  ENDIF
+!-
+! Fill out the items of the data base
+!-
+  nb_keys = nb_keys+1
+  keysig(nb_keys) = target_sig
+  keystr(nb_keys) = target(1:MIN(len_trim(target),30))
+  keystatus(nb_keys) = status
+  keytype(nb_keys) = 2
+  keyfromfile(nb_keys) = fileorig
+!-
+! Can we compress the data base entry ?
+!-
+  IF (     (MINVAL(tmp_ret_val) == MAXVAL(tmp_ret_val)) &
+&     .AND.(size_of_in > compress_lim)) THEN
+    keymemstart(nb_keys) = realmempos+1
+    keycompress(nb_keys) = size_of_in
+    keymemlen(nb_keys) = 1
+  ELSE
+    keymemstart(nb_keys) = realmempos+1
+    keycompress(nb_keys) = -1
+    keymemlen(nb_keys) = size_of_in
+  ENDIF
+!-
+! Before writing the actual size lets see if we have the space
+!-
+  IF (keymemstart(nb_keys)+keymemlen(nb_keys) > realmemsize) THEN
+    CALL getin_allocmem (2,keymemlen(nb_keys))
+  ENDIF
+!-
+  realmem(keymemstart(nb_keys): &
+ &        keymemstart(nb_keys)+keymemlen(nb_keys)-1) = &
+ &  tmp_ret_val(1:keymemlen(nb_keys))
+  realmempos = keymemstart(nb_keys)+keymemlen(nb_keys)-1
+!---------------------
+END SUBROUTINE getdbwr
+!-
+!===
+!-
+SUBROUTINE getdbrr (pos,size_of_in,target,tmp_ret_val)
+!---------------------------------------------------------------------
+!- Read the required variables in the database for REALS
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: pos, size_of_in
+  CHARACTER(LEN=*) :: target
+  REAL,DIMENSION(:) :: tmp_ret_val
+!---------------------------------------------------------------------
+  IF (keytype(pos) /= 2) THEN
+    WRITE(*,*) 'FATAL ERROR : Wrong data type for keyword ',target
+    STOP 'getdbrr'
+  ENDIF
+!-
+  IF (keycompress(pos) > 0) THEN
+    IF (    (keycompress(pos) /= size_of_in) &
+ &      .OR.(keymemlen(pos) /= 1) ) THEN
+      WRITE(*,*) &
+ &      'FATAL ERROR : Wrong compression length for keyword ',target
+      STOP 'getdbrr'
+    ELSE
+      tmp_ret_val(1:size_of_in) = realmem(keymemstart(pos))
+    ENDIF
+  ELSE
+    IF (keymemlen(pos) /= size_of_in) THEN
+      WRITE(*,*) 'FATAL ERROR : Wrong array length for keyword ',target
+      STOP 'getdbrr'
+    ELSE
+      tmp_ret_val(1:size_of_in) = &
+ &      realmem(keymemstart(pos):keymemstart(pos)+keymemlen(pos)-1)
+    ENDIF
+  ENDIF
+!---------------------
+END SUBROUTINE getdbrr
+!-
+!=== CHARACTER database INTERFACE
+!-
+SUBROUTINE getdbwc &
+ &  (target,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+!---------------------------------------------------------------------
+!- Write the CHARACTER data into the data base
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: target
+  INTEGER :: target_sig,status,fileorig,size_of_in
+  CHARACTER(LEN=*),DIMENSION(:) :: tmp_ret_val
+!---------------------------------------------------------------------
+!-
+! First check if we have sufficiant space for the new key
+!-
+  IF (nb_keys+1 > keymemsize) THEN
+    CALL getin_allockeys ()
+  ENDIF
+!-
+! Fill out the items of the data base
+!-
+  nb_keys = nb_keys+1
+  keysig(nb_keys) = target_sig
+  keystr(nb_keys) = target(1:MIN(len_trim(target),30))
+  keystatus(nb_keys) = status
+  keytype(nb_keys) = 3
+  keyfromfile(nb_keys) = fileorig
+  keymemstart(nb_keys) = charmempos+1
+  keymemlen(nb_keys) = size_of_in
+!-
+! Before writing the actual size lets see if we have the space
+!-
+  IF (keymemstart(nb_keys)+keymemlen(nb_keys) > realmemsize) THEN
+    CALL getin_allocmem (3,keymemlen(nb_keys))
+  ENDIF
+!-
+  charmem(keymemstart(nb_keys): &
+ &        keymemstart(nb_keys)+keymemlen(nb_keys)-1) = &
+ &  tmp_ret_val(1:keymemlen(nb_keys))
+  charmempos = keymemstart(nb_keys)+keymemlen(nb_keys)-1
+!---------------------
+END SUBROUTINE getdbwc
+!-
+!===
+!-
+SUBROUTINE getdbrc(pos,size_of_in,target,tmp_ret_val)
+!---------------------------------------------------------------------
+!- Read the required variables in the database for CHARACTER
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: pos, size_of_in
+  CHARACTER(LEN=*) :: target
+  CHARACTER(LEN=*),DIMENSION(:) :: tmp_ret_val
+!---------------------------------------------------------------------
+  IF (keytype(pos) /= 3) THEN
+    WRITE(*,*) 'FATAL ERROR : Wrong data type for keyword ',target
+    STOP 'getdbrc'
+  ENDIF
+!-
+  IF (keymemlen(pos) /= size_of_in) THEN
+    WRITE(*,*) 'FATAL ERROR : Wrong array length for keyword ',target
+    STOP 'getdbrc'
+  ELSE
+    tmp_ret_val(1:size_of_in) = &
+ &    charmem(keymemstart(pos):keymemstart(pos)+keymemlen(pos)-1)
+  ENDIF
+!---------------------
+END SUBROUTINE getdbrc
+!-
+!=== LOGICAL database INTERFACE
+!-
+SUBROUTINE getdbwl &
+ &  (target,target_sig,status,fileorig,size_of_in,tmp_ret_val)
+!---------------------------------------------------------------------
+!- Write the LOGICAL data into the data base
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(LEN=*) :: target
+  INTEGER :: target_sig, status, fileorig, size_of_in
+  LOGICAL,DIMENSION(:) :: tmp_ret_val
+!---------------------------------------------------------------------
+!-
+! First check if we have sufficiant space for the new key
+!-
+  IF (nb_keys+1 > keymemsize) THEN
+    CALL getin_allockeys ()
+  ENDIF
+!-
+! Fill out the items of the data base
+!-
+  nb_keys = nb_keys+1
+  keysig(nb_keys) = target_sig
+  keystr(nb_keys) = target(1:MIN(len_trim(target),30))
+  keystatus(nb_keys) = status
+  keytype(nb_keys) = 4
+  keyfromfile(nb_keys) = fileorig
+  keymemstart(nb_keys) = logicmempos+1
+  keymemlen(nb_keys) = size_of_in
+!-
+! Before writing the actual size lets see if we have the space
+!-
+  IF (keymemstart(nb_keys)+keymemlen(nb_keys) > logicmemsize) THEN
+    CALL getin_allocmem (4,keymemlen(nb_keys))
+  ENDIF
+!-
+  logicmem(keymemstart(nb_keys): &
+ &         keymemstart(nb_keys)+keymemlen(nb_keys)-1) = &
+ &  tmp_ret_val(1:keymemlen(nb_keys))
+  logicmempos = keymemstart(nb_keys)+keymemlen(nb_keys)-1
+!---------------------
+END SUBROUTINE getdbwl
+!-
+!===
+!-
+SUBROUTINE getdbrl(pos,size_of_in,target,tmp_ret_val)
+!---------------------------------------------------------------------
+!- Read the required variables in the database for LOGICALS
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: pos, size_of_in
+  CHARACTER(LEN=*) :: target
+  LOGICAL,DIMENSION(:) :: tmp_ret_val
+!---------------------------------------------------------------------
+  IF (keytype(pos) /= 4) THEN
+    WRITE(*,*) 'FATAL ERROR : Wrong data type for keyword ',target
+    STOP 'getdbrl'
+  ENDIF
+!-
+  IF (keymemlen(pos) /= size_of_in) THEN
+    WRITE(*,*) 'FATAL ERROR : Wrong array length for keyword ',target
+    STOP 'getdbrl'
+  ELSE
+    tmp_ret_val(1:size_of_in) = &
+ &    logicmem(keymemstart(pos):keymemstart(pos)+keymemlen(pos)-1)
+  ENDIF
+!---------------------
+END SUBROUTINE getdbrl
+!-
+!===
+!-
+SUBROUTINE getin_allockeys ()
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER,ALLOCATABLE :: tmp_int(:)
+  CHARACTER(LEN=100),ALLOCATABLE :: tmp_str(:)
+!-
+  INTEGER :: ier
+!---------------------------------------------------------------------
+!-
+! Either nothing exists in these arrays and it is easy to do
+!-
+  IF (keymemsize == 0) THEN
+!-
+    ALLOCATE(keysig(memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keysig to size ', &
+ &      memslabs
+      STOP
+    ENDIF
+!-
+    ALLOCATE(keystr(memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keystr to size ', &
+ &      memslabs
+      STOP
+    ENDIF
+!-
+    ALLOCATE(keystatus(memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keystatus to size ', &
+ &      memslabs
+      STOP
+    ENDIF
+!-
+    ALLOCATE(keytype(memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keytype to size ', &
+ &      memslabs
+      STOP
+    ENDIF
+!-
+    ALLOCATE(keycompress(memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keycompress to size ', &
+ &      memslabs
+      STOP
+    ENDIF
+!-
+    ALLOCATE(keyfromfile(memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keyfromfile to size ', &
+ &      memslabs
+      STOP
+    ENDIF
+!-
+    ALLOCATE(keymemstart(memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keymemstart to size ', &
+ &      memslabs
+      STOP
+    ENDIF
+!-
+    ALLOCATE(keymemlen(memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keymemlen to size ', &
+ &      memslabs
+      STOP
+    ENDIF
+!-
+    nb_keys = 0
+    keymemsize = memslabs
+    keycompress(:) = -1
+!-
+  ELSE
+!-
+!-- There is something already in the memory,
+!-- we need to transfer and reallocate.
+!-
+    ALLOCATE(tmp_str(keymemsize),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate tmp_str to size ', &
+ &      keymemsize
+      STOP
+    ENDIF
+!-
+    ALLOCATE(tmp_int(keymemsize),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate tmp_int to size ', &
+ &      keymemsize
+      STOP
+    ENDIF
+!-
+    tmp_int(1:keymemsize) = keysig(1:keymemsize)
+    DEALLOCATE(keysig)
+    ALLOCATE(keysig(keymemsize+memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keysig to size ', &
+ &      keymemsize+memslabs
+      STOP
+    ENDIF
+    keysig(1:keymemsize) = tmp_int(1:keymemsize)
+!-
+    tmp_str(1:keymemsize) = keystr(1:keymemsize)
+    DEALLOCATE(keystr)
+    ALLOCATE(keystr(keymemsize+memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keystr to size ', &
+ &      keymemsize+memslabs
+      STOP
+    ENDIF
+    keystr(1:keymemsize) = tmp_str(1:keymemsize)
+!-
+    tmp_int(1:keymemsize) = keystatus(1:keymemsize)
+    DEALLOCATE(keystatus)
+    ALLOCATE(keystatus(keymemsize+memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keystatus to size ', &
+ &      keymemsize+memslabs
+      STOP
+    ENDIF
+    keystatus(1:keymemsize) = tmp_int(1:keymemsize)
+!-
+    tmp_int(1:keymemsize) = keytype(1:keymemsize)
+    DEALLOCATE(keytype)
+    ALLOCATE(keytype(keymemsize+memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keytype to size ', &
+ &      keymemsize+memslabs
+      STOP
+    ENDIF
+    keytype(1:keymemsize) = tmp_int(1:keymemsize)
+!-
+    tmp_int(1:keymemsize) = keycompress(1:keymemsize)
+    DEALLOCATE(keycompress)
+    ALLOCATE(keycompress(keymemsize+memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keycompress to size ', &
+ &      keymemsize+memslabs
+      STOP
+    ENDIF
+    keycompress(:) = -1
+    keycompress(1:keymemsize) = tmp_int(1:keymemsize)
+!-
+    tmp_int(1:keymemsize) = keyfromfile(1:keymemsize)
+    DEALLOCATE(keyfromfile)
+    ALLOCATE(keyfromfile(keymemsize+memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keyfromfile to size ', &
+ &      keymemsize+memslabs
+      STOP
+    ENDIF
+    keyfromfile(1:keymemsize) = tmp_int(1:keymemsize)
+!-
+    tmp_int(1:keymemsize) = keymemstart(1:keymemsize)
+    DEALLOCATE(keymemstart)
+    ALLOCATE(keymemstart(keymemsize+memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keymemstart to size ', &
+ &      keymemsize+memslabs
+      STOP
+    ENDIF
+    keymemstart(1:keymemsize) = tmp_int(1:keymemsize)
+!-
+    tmp_int(1:keymemsize) = keymemlen(1:keymemsize)
+    DEALLOCATE(keymemlen)
+    ALLOCATE(keymemlen(keymemsize+memslabs),stat=ier)
+    IF (ier /= 0) THEN
+      WRITE(*,*) &
+ &      'getin_allockeys : Can not allocate keymemlen to size ', &
+ &      keymemsize+memslabs
+      STOP
+    ENDIF
+    keymemlen(1:keymemsize) = tmp_int(1:keymemsize)
+!-
+    keymemsize = keymemsize+memslabs
+!-
+    DEALLOCATE(tmp_int)
+    DEALLOCATE(tmp_str)
+  ENDIF
+!-----------------------------
+END SUBROUTINE getin_allockeys
+!-
+!===
+!-
+SUBROUTINE getin_allocmem (type,len_wanted)
+!---------------------------------------------------------------------
+!- Allocate the memory of the data base for all 4 types of memory
+!-
+!- 1 = INTEGER
+!- 2 = REAL
+!- 3 = CHAR
+!- 4 = LOGICAL
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: type, len_wanted
+!-
+  INTEGER,ALLOCATABLE :: tmp_int(:)
+  CHARACTER(LEN=100),ALLOCATABLE :: tmp_char(:)
+  REAL,ALLOCATABLE :: tmp_real(:)
+  LOGICAL,ALLOCATABLE :: tmp_logic(:)
+  INTEGER :: ier
+!---------------------------------------------------------------------
+  SELECT CASE (type)
+  CASE(1)
+    IF (intmemsize == 0) THEN
+      ALLOCATE(intmem(memslabs),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to allocate db-memory intmem to ', &
+ &    memslabs
+        STOP
+      ENDIF
+      intmemsize=memslabs
+    ELSE
+      ALLOCATE(tmp_int(intmemsize),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to allocate tmp_int to ', &
+ &    intmemsize
+        STOP
+      ENDIF
+      tmp_int(1:intmemsize) = intmem(1:intmemsize)
+      DEALLOCATE(intmem)
+      ALLOCATE(intmem(intmemsize+MAX(memslabs,len_wanted)),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to re-allocate db-memory intmem to ', &
+ &    intmemsize+MAX(memslabs,len_wanted)
+        STOP
+      ENDIF
+      intmem(1:intmemsize) = tmp_int(1:intmemsize)
+      intmemsize = intmemsize+MAX(memslabs,len_wanted)
+      DEALLOCATE(tmp_int)
+    ENDIF
+  CASE(2)
+    IF (realmemsize == 0) THEN
+      ALLOCATE(realmem(memslabs),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to allocate db-memory realmem to ', &
+ &    memslabs
+        STOP
+      ENDIF
+      realmemsize =  memslabs
+    ELSE
+      ALLOCATE(tmp_real(realmemsize),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to allocate tmp_real to ', &
+ &    realmemsize
+        STOP
+      ENDIF
+      tmp_real(1:realmemsize) = realmem(1:realmemsize)
+      DEALLOCATE(realmem)
+      ALLOCATE(realmem(realmemsize+MAX(memslabs,len_wanted)),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to re-allocate db-memory realmem to ', &
+ &    realmemsize+MAX(memslabs,len_wanted)
+        STOP
+      ENDIF
+      realmem(1:realmemsize) = tmp_real(1:realmemsize)
+      realmemsize = realmemsize+MAX(memslabs,len_wanted)
+      DEALLOCATE(tmp_real)
+    ENDIF
+  CASE(3)
+    IF (charmemsize == 0) THEN
+      ALLOCATE(charmem(memslabs),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to allocate db-memory charmem to ', &
+ &    memslabs
+        STOP
+      ENDIF
+      charmemsize = memslabs
+    ELSE
+      ALLOCATE(tmp_char(charmemsize),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to allocate tmp_char to ', &
+ &    charmemsize
+        STOP
+      ENDIF
+      tmp_char(1:charmemsize) = charmem(1:charmemsize)
+      DEALLOCATE(charmem)
+      ALLOCATE(charmem(charmemsize+MAX(memslabs,len_wanted)),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to re-allocate db-memory charmem to ', &
+ &    charmemsize+MAX(memslabs,len_wanted)
+        STOP
+      ENDIF
+      charmem(1:charmemsize) = tmp_char(1:charmemsize)
+      charmemsize = charmemsize+MAX(memslabs,len_wanted)
+      DEALLOCATE(tmp_char)
+    ENDIF
+  CASE(4)
+    IF (logicmemsize == 0) THEN
+      ALLOCATE(logicmem(memslabs),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to allocate db-memory logicmem to ', &
+ &    memslabs
+        STOP
+      ENDIF
+      logicmemsize = memslabs
+    ELSE
+      ALLOCATE(tmp_logic(logicmemsize),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to allocate tmp_logic to ', &
+ &    logicmemsize
+        STOP
+      ENDIF
+      tmp_logic(1:logicmemsize) = logicmem(1:logicmemsize)
+      DEALLOCATE(logicmem)
+      ALLOCATE(logicmem(logicmemsize+MAX(memslabs,len_wanted)),stat=ier)
+      IF (ier /= 0) THEN
+        WRITE(*,*) &
+ &    'getin_allocmem : Unable to re-allocate db-memory logicmem to ', &
+ &    logicmemsize+MAX(memslabs,len_wanted)
+        STOP
+      ENDIF
+      logicmem(1:logicmemsize) = tmp_logic(1:logicmemsize)
+      logicmemsize = logicmemsize+MAX(memslabs,len_wanted)
+      DEALLOCATE(tmp_logic)
+    ENDIF
+  CASE DEFAULT
+    WRITE(*,*) 'getin_allocmem : Unknown type : ',type
+    STOP
+  END SELECT
+!----------------------------
+END SUBROUTINE getin_allocmem
+!-
+!===
+!-
+SUBROUTINE getin_dump (fileprefix)
+!---------------------------------------------------------------------
+!- This subroutine will dump the content of the database into  file
+!- which has the same format as the run.def. The idea is that the user
+!- can see which parameters were used and re-use the file for another
+!- run.
+!-
+!- The argument file allows the user to change the name of the file
+!- in which the data will be archived
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  CHARACTER(*),OPTIONAL :: fileprefix
+!-
+  CHARACTER(LEN=80) :: usedfileprefix = "used"
+  INTEGER :: ikey,if,iff,iv
+  CHARACTER(LEN=3) :: tmp3
+  CHARACTER(LEN=100) :: tmp_str, used_filename
+  LOGICAL :: check = .FALSE.
+!---------------------------------------------------------------------
+  IF (PRESENT(fileprefix)) THEN
+    usedfileprefix = fileprefix(1:MIN(len_trim(fileprefix),80))
+  ENDIF
+!-
+  DO if=1,nbfiles
+!---
+    used_filename = TRIM(usedfileprefix)//'_'//TRIM(filelist(if))
+    IF (check) THEN
+      WRITE(*,*) &
+ &      'GETIN_DUMP : opens file : ',TRIM(used_filename),' if = ',if
+      WRITE(*,*) 'GETIN_DUMP : NUMBER OF KEYS : ',nb_keys
+    ENDIF
+    OPEN(unit=76,file=used_filename)
+!-
+!-- If this is the first file we need to add the list
+!-- of file which belong to it
+!-
+    IF ( (if == 1) .AND. (nbfiles > 1) ) THEN
+      WRITE(76,*) '# '
+      WRITE(76,*) '# This file is linked to the following files :'
+      WRITE(76,*) '# '
+      DO iff=2,nbfiles
+        WRITE(76,*) 'INCLUDEDEF = ',TRIM(filelist(iff))
+      ENDDO
+      WRITE(76,*) '# '
+    ENDIF
+!---
+    DO ikey=1,nb_keys
+!-
+!---- Is this key form this file ?
+!-
+      IF (keyfromfile(ikey) == if) THEN
+!-
+!---- Write some comments
+!-
+        WRITE(76,*) '#'
+        SELECT CASE (keystatus(ikey))
+        CASE(1)
+          WRITE(76,*) '# Values of ', &
+ &          TRIM(keystr(ikey)),' comes from the run.def.'
+        CASE(2)
+          WRITE(76,*) '# Values of ', &
+ &          TRIM(keystr(ikey)),' are all defaults.'
+        CASE(3)
+          WRITE(76,*) '# Values of ', &
+ &          TRIM(keystr(ikey)),' are a mix of run.def and defaults.'
+        CASE DEFAULT
+          WRITE(76,*) '# Dont know from where the value of ', &
+ &          TRIM(keystr(ikey)),' comes.'
+        END SELECT
+        WRITE(76,*) '#'
+!-
+!---- Write the values
+!-
+        SELECT CASE (keytype(ikey))
+!-
+        CASE(1)
+          IF (keymemlen(ikey) == 1) THEN
+            IF (keycompress(ikey) < 0) THEN
+              WRITE(76,*) &
+ &              TRIM(keystr(ikey)),' = ',intmem(keymemstart(ikey))
+            ELSE
+              WRITE(76,*) &
+ &              TRIM(keystr(ikey)),' = ',keycompress(ikey), &
+ &              ' * ',intmem(keymemstart(ikey))
+            ENDIF
+          ELSE
+            DO iv=0,keymemlen(ikey)-1
+              WRITE(tmp3,'(I3.3)') iv+1
+              WRITE(76,*) &
+ &              TRIM(keystr(ikey)),'__',tmp3, &
+ &              ' = ',intmem(keymemstart(ikey)+iv)
+            ENDDO
+          ENDIF
+!-
+        CASE(2)
+          IF (keymemlen(ikey) == 1) THEN
+            IF (keycompress(ikey) < 0) THEN
+              WRITE(76,*) &
+ &              TRIM(keystr(ikey)),' = ',realmem(keymemstart(ikey))
+            ELSE
+              WRITE(76,*) &
+ &              TRIM(keystr(ikey)),' = ',keycompress(ikey),&
+                   & ' * ',realmem(keymemstart(ikey))
+            ENDIF
+          ELSE
+            DO iv=0,keymemlen(ikey)-1
+              WRITE(tmp3,'(I3.3)') iv+1
+              WRITE(76,*) &
+ &              TRIM(keystr(ikey)),'__',tmp3, &
+ &              ' = ',realmem(keymemstart(ikey)+iv)
+            ENDDO
+          ENDIF
+        CASE(3)
+          IF (keymemlen(ikey) == 1) THEN
+            tmp_str = charmem(keymemstart(ikey))
+            WRITE(76,*) TRIM(keystr(ikey)),' = ',TRIM(tmp_str)
+          ELSE
+            DO iv=0,keymemlen(ikey)-1
+              WRITE(tmp3,'(I3.3)') iv+1
+              tmp_str = charmem(keymemstart(ikey)+iv)
+              WRITE(76,*) &
+ &              TRIM(keystr(ikey)),'__',tmp3,' = ',TRIM(tmp_str)
+            ENDDO
+          ENDIF
+        CASE(4)
+          IF (keymemlen(ikey) == 1) THEN
+            IF (logicmem(keymemstart(ikey))) THEN
+              WRITE(76,*) TRIM(keystr(ikey)),' = TRUE '
+            ELSE
+              WRITE(76,*) TRIM(keystr(ikey)),' = FALSE '
+            ENDIF
+          ELSE
+            DO iv=0,keymemlen(ikey)-1
+              WRITE(tmp3,'(I3.3)') iv+1
+              IF (logicmem(keymemstart(ikey)+iv)) THEN
+                WRITE(76,*) TRIM(keystr(ikey)),'__',tmp3,' = TRUE '
+              ELSE
+                WRITE(76,*) TRIM(keystr(ikey)),'__',tmp3,' = FALSE '
+              ENDIF
+            ENDDO
+          ENDIF
+!-
+        CASE DEFAULT
+          WRITE(*,*) &
+ &          'FATAL ERROR : Unknown type for variable ', &
+ &          TRIM(keystr(ikey))
+          STOP 'getin_dump'
+        END SELECT
+      ENDIF
+    ENDDO
+!-
+    CLOSE(unit=76)
+!-
+  ENDDO
+!------------------------
+END SUBROUTINE getin_dump
+!-
+!===
+!-
+END MODULE ioipsl_getincom
Index: LMDZ4/trunk/libf/dyn3dpar/ioipsl_stringop.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/ioipsl_stringop.F90	(revision 1146)
+++ LMDZ4/trunk/libf/dyn3dpar/ioipsl_stringop.F90	(revision 1146)
@@ -0,0 +1,256 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Module and routines in this file are taken from IOIPSL
+! files stringop.f90
+! Module names has been changed to avoid problems
+! if compiling model with IOIPSL library
+! Ehouarn - March 2009
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+MODULE ioipsl_stringop
+!---------------------------------------------------------------------
+!-
+  INTEGER, DIMENSION(30) :: &
+       & prime=(/1,2,3,5,7,11,13,17,19,23,29,31,37,41,43, &
+       & 47,53,59,61,67,71,73,79,83,89,97,101,103,107,109/)
+!-
+!---------------------------------------------------------------------
+CONTAINS
+!=
+   SUBROUTINE cmpblank (str)
+!---------------------------------------------------------------------
+!-
+!---------------------------------------------------------------------
+   CHARACTER(LEN=*),INTENT(inout) :: str
+!-
+   INTEGER :: lcc,ipb
+!---------------------------------------------------------------------
+   lcc = LEN_TRIM(str)
+   ipb = 1
+   DO
+     IF (ipb >= lcc)   EXIT
+     IF (str(ipb:ipb+1) == '  ') THEN
+       str(ipb+1:) = str(ipb+2:lcc)
+       lcc = lcc-1
+     ELSE
+       ipb = ipb+1
+     ENDIF
+   ENDDO
+!-------------------------
+   END SUBROUTINE cmpblank
+!=
+   INTEGER FUNCTION cntpos (c_c,l_c,c_r,l_r)
+!---------------------------------------------------------------------
+!- Finds number of occurences of c_r in c_c
+!---------------------------------------------------------------------
+   IMPLICIT NONE
+!-
+   CHARACTER(LEN=*),INTENT(in) :: c_c
+   INTEGER,INTENT(IN) :: l_c
+   CHARACTER(LEN=*),INTENT(in) :: c_r
+   INTEGER,INTENT(IN) :: l_r
+!-
+   INTEGER :: ipos,indx,ires
+!---------------------------------------------------------------------
+   cntpos = 0
+   ipos   = 1
+   DO
+     indx = INDEX(c_c(ipos:l_c),c_r(1:l_r))
+     IF (indx > 0) THEN
+       cntpos = cntpos+1
+       ipos   = ipos+indx+l_r-1
+     ELSE
+       EXIT
+     ENDIF
+   ENDDO
+!---------------------
+   END FUNCTION cntpos
+!=
+   INTEGER FUNCTION findpos (c_c,l_c,c_r,l_r)
+!---------------------------------------------------------------------
+!- Finds position of c_r in c_c
+!---------------------------------------------------------------------
+   IMPLICIT NONE
+!-
+   CHARACTER(LEN=*),INTENT(in) :: c_c
+   INTEGER,INTENT(IN) :: l_c
+   CHARACTER(LEN=*),INTENT(in) :: c_r
+   INTEGER,INTENT(IN) :: l_r
+!---------------------------------------------------------------------
+    findpos = INDEX(c_c(1:l_c),c_r(1:l_r))
+    IF (findpos == 0)   findpos=-1
+!----------------------
+   END FUNCTION findpos
+!=
+   SUBROUTINE find_str (nb_str,str_tab,str_len_tab,str,pos)
+!---------------------------------------------------------------------
+!- This subroutine looks for a string in a table
+!---------------------------------------------------------------------
+!- INPUT
+!-   nb_str      : length of table
+!-   str_tab     : Table  of strings
+!-   str_len_tab : Table  of string-length
+!-   str         : Target we are looking for
+!- OUTPUT
+!-   pos         : -1 if str not found, else value in the table
+!---------------------------------------------------------------------
+   IMPLICIT NONE
+!-
+   INTEGER :: nb_str
+   CHARACTER(LEN=*),DIMENSION(nb_str) :: str_tab
+   INTEGER,DIMENSION(nb_str) :: str_len_tab
+   CHARACTER(LEN=*) :: str
+   INTEGER :: pos
+!-
+   INTEGER :: i,il
+!---------------------------------------------------------------------
+   pos = -1
+   il = LEN_TRIM(str)
+   IF ( nb_str > 0 ) THEN
+      DO i=1,nb_str
+         IF (     (INDEX(str_tab(i),str(1:il)) > 0) &
+              .AND.(str_len_tab(i) == il) ) THEN
+            pos = i
+            EXIT
+         ENDIF
+      ENDDO
+   ENDIF
+!-------------------------
+   END SUBROUTINE find_str
+!=
+   SUBROUTINE nocomma (str)
+!---------------------------------------------------------------------
+!-
+!---------------------------------------------------------------------
+   IMPLICIT NONE
+!-
+   CHARACTER(LEN=*) :: str
+!-
+   INTEGER :: i
+!---------------------------------------------------------------------
+   DO i=1,LEN_TRIM(str)
+     IF (str(i:i) == ',')   str(i:i) = ' '
+   ENDDO
+!------------------------
+   END SUBROUTINE nocomma
+!=
+   SUBROUTINE strlowercase (str)
+!---------------------------------------------------------------------
+!- Converts a string into lowercase
+!---------------------------------------------------------------------
+   IMPLICIT NONE
+!-
+   CHARACTER(LEN=*) :: str
+!-
+   INTEGER :: i,ic
+!---------------------------------------------------------------------
+   DO i=1,LEN_TRIM(str)
+     ic = IACHAR(str(i:i))
+     IF ( (ic >= 65) .AND. (ic <= 90) )   str(i:i) = ACHAR(ic+32)
+   ENDDO
+!-----------------------------
+   END SUBROUTINE strlowercase
+!=
+   SUBROUTINE struppercase (str)
+!---------------------------------------------------------------------
+!- Converts a string into uppercase
+!---------------------------------------------------------------------
+   IMPLICIT NONE
+!-
+   CHARACTER(LEN=*) :: str
+!-
+   INTEGER :: i,ic
+!---------------------------------------------------------------------
+   DO i=1,LEN_TRIM(str)
+     ic = IACHAR(str(i:i))
+     IF ( (ic >= 97) .AND. (ic <= 122) )   str(i:i) = ACHAR(ic-32)
+   ENDDO
+!-----------------------------
+   END SUBROUTINE struppercase
+!=
+!------------------
+   SUBROUTINE gensig (str, sig)
+!---------------------------------------------------------------------
+!- Generate a signature from the first 30 characters of the string
+!- This signature is not unique and thus when one looks for the
+!- one needs to also verify the string.
+!---------------------------------------------------------------------
+   IMPLICIT NONE
+!-
+   CHARACTER(LEN=*) :: str
+   INTEGER          :: sig
+!-
+   INTEGER :: i
+!---------------------------------------------------------------------
+   sig = 0
+   DO i=1,MIN(len_trim(str),30)
+      sig = sig  + prime(i)*IACHAR(str(i:i))
+   ENDDO
+!-----------------------------
+ END SUBROUTINE gensig
+!=
+!------------------
+   SUBROUTINE find_sig (nb_sig, str_tab, str, sig_tab, sig, pos)
+!---------------------------------------------------------------------
+!- Find the string signature in a list of signatures
+!---------------------------------------------------------------------
+!- INPUT
+!-   nb_sig      : length of table of signatures
+!-   str_tab     : Table of strings
+!-   str         : Target string we are looking for 
+!-   sig_tab     : Table of signatures
+!-   sig         : Target signature we are looking for
+!- OUTPUT
+!-   pos         : -1 if str not found, else value in the table
+!---------------------------------------------------------------------
+   IMPLICIT NONE
+!-
+   INTEGER :: nb_sig
+   CHARACTER(LEN=*),DIMENSION(nb_sig) :: str_tab
+   CHARACTER(LEN=*) :: str
+   INTEGER, DIMENSION(nb_sig) :: sig_tab
+   INTEGER :: sig
+!-
+   INTEGER :: pos
+   INTEGER, DIMENSION(nb_sig) :: loczeros
+!-
+   INTEGER :: il, len
+   INTEGER, DIMENSION(1) :: minpos
+!---------------------------------------------------------------------
+!-
+   pos = -1
+   il = LEN_TRIM(str)
+!-
+   IF ( nb_sig > 0 ) THEN
+      !
+      loczeros = ABS(sig_tab(1:nb_sig)-sig)
+      !
+      IF ( COUNT(loczeros < 1) == 1 ) THEN
+         !
+         minpos = MINLOC(loczeros)
+         len = LEN_TRIM(str_tab(minpos(1)))
+         IF ( (INDEX(str_tab(minpos(1)),str(1:il)) > 0) &
+                 .AND.(len == il) ) THEN
+            pos = minpos(1)
+         ENDIF
+         !
+      ELSE IF ( COUNT(loczeros < 1) > 1 ) THEN
+         !
+         DO WHILE (COUNT(loczeros < 1) >= 1 .AND. pos < 0 )
+            minpos = MINLOC(loczeros)
+            len = LEN_TRIM(str_tab(minpos(1)))
+            IF ( (INDEX(str_tab(minpos(1)),str(1:il)) > 0) &
+                 .AND.(len == il) ) THEN
+               pos = minpos(1)
+            ELSE
+               loczeros(minpos(1)) = 99999
+            ENDIF
+         ENDDO
+         !
+      ENDIF
+      !
+   ENDIF
+!-
+ END SUBROUTINE find_sig
+!=
+!------------------
+END MODULE ioipsl_stringop
Index: LMDZ4/trunk/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 1146)
@@ -4,10 +4,6 @@
 c
 c
-#define IO_DEBUG
-
-#undef CPP_IOIPSL
-#define CPP_IOIPSL
-
-      SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
+
+      SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
      &                    time_0)
 
@@ -21,4 +17,5 @@
        USE vampir
        USE timer_filtre, ONLY : print_filtre_timer
+       USE infotrac
 
       IMPLICIT NONE
@@ -69,13 +66,6 @@
 #include "com_io_dyn.h"
 #include "iniprint.h"
-
-c#include "tracstoke.h"
-
 #include "academic.h"
-!#include "clesphys.h"
-#include "advtrac.h"
       
-      integer nq
-
       INTEGER         longcles
       PARAMETER     ( longcles = 20 )
@@ -88,5 +78,5 @@
       REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
       REAL :: teta(ip1jmp1,llm)                 ! temperature potentielle 
-      REAL :: q(ip1jmp1,llm,nqmx)               ! champs advectes
+      REAL :: q(ip1jmp1,llm,nqtot)              ! champs advectes
       REAL :: ps(ip1jmp1)                       ! pression  au sol
       REAL,SAVE :: p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
@@ -109,5 +99,6 @@
 c   tendances dynamiques
       REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
-      REAL,SAVE :: dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
+      REAL,SAVE :: dteta(ip1jmp1,llm),dp(ip1jmp1)
+      REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq
 
 c   tendances de la dissipation
@@ -118,5 +109,6 @@
       REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
       REAL,SAVE :: dtetafi(ip1jmp1,llm)
-      REAL,SAVE :: dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
+      REAL,SAVE :: dpfi(ip1jmp1)
+      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
 
 c   variables pour le fichier histoire
@@ -186,5 +178,6 @@
       type(Request) :: Request_physic
       REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
-      REAL,SAVE :: dtetafi_tmp(iip1,llm),dqfi_tmp(iip1,llm,nqmx)
+      REAL,SAVE :: dtetafi_tmp(iip1,llm)
+      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi_tmp
       REAL,SAVE :: dpfi_tmp(iip1)
 
@@ -195,4 +188,5 @@
       INTEGER :: var_time
       LOGICAL :: ok_start_timer=.FALSE.
+      LOGICAL, SAVE :: firstcall=.TRUE.
 
 c$OMP MASTER
@@ -208,5 +202,5 @@
       itaufin   = nday*day_step
       itaufinp1 = itaufin +1
-
+      modname="leapfrog_p"
 
       itau = 0
@@ -217,4 +211,15 @@
           iday = iday+1
          ENDIF
+
+c Allocate variables depending on dynamic variable nqtot
+c$OMP MASTER
+         IF (firstcall) THEN
+            firstcall=.FALSE.
+            ALLOCATE(dq(ip1jmp1,llm,nqtot))
+            ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
+            ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
+         END IF
+c$OMP END MASTER      
+c$OMP BARRIER
 
 c-----------------------------------------------------------------------
@@ -276,5 +281,5 @@
 c$OMP BARRIER
        else
-          
+! Save fields obtained at previous time step as '...m1'
          ijb=ij_begin
          ije=ij_end
@@ -303,5 +308,5 @@
      .                    llm, -2,2, .TRUE., 1 )
 
-       endif
+       endif ! of if (FirstCaldyn)
        
       forward = .TRUE.
@@ -347,9 +352,9 @@
          IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
          IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 
-     s          .and. iflag_phys.NE.0                 ) apphys = .TRUE.
+     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
       ELSE
          IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
          IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
-         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.NE.0) apphys=.TRUE.
+         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
       END IF
 
@@ -455,5 +460,5 @@
      &                                jj_Nb_caldyn,0,0,TestRequest)
  
-        do j=1,nqmx
+        do j=1,nqtot
          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
      &                                jj_nb_caldyn,0,0,TestRequest)
@@ -490,5 +495,5 @@
        call Register_Hallo(p,ip1jmp1,llmp1,1,1,1,1,TestRequest)
        
-c       do j=1,nqmx
+c       do j=1,nqtot
 c         call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
 c     *                       TestRequest)
@@ -516,5 +521,5 @@
         call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
         call WriteField_p('phis',reshape(phis,(/iip1,jmp1/)))
-        do j=1,nqmx
+        do j=1,nqtot
           call WriteField_p('q'//trim(int2str(j)),
      .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
@@ -528,5 +533,7 @@
 
 c$OMP MASTER
-      print*,"Iteration No",True_itau
+      IF (prt_level>9) THEN
+        WRITE(lunout,*)"leapfrog_p: Iteration No",True_itau
+      ENDIF
 
 
@@ -585,7 +592,7 @@
 
 
-         ENDIF
-c
-      ENDIF
+         ENDIF ! of IF (offline)
+c
+      ENDIF ! of IF( forward. OR . leapf )
 cc$OMP END PARALLEL
 
@@ -608,4 +615,5 @@
 c$OMP BARRIER
 !       CALL FTRACE_REGION_BEGIN("integrd")
+
        CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
      $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
@@ -625,5 +633,5 @@
 c
 c      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
-c      do j=1,nqmx
+c      do j=1,nqtot
 c        call WriteField_p('q'//trim(int2str(j)),
 c     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
@@ -663,5 +671,7 @@
 c$OMP MASTER
          call suspend_timer(timer_caldyn)
-         print*,'Entree dans la physique : Iteration No ',true_itau
+
+         write(lunout,*)
+     &   'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
 c$OMP END MASTER
 
@@ -669,5 +679,4 @@
 
 c$OMP BARRIER
-
          CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
 c$OMP BARRIER
@@ -683,5 +692,4 @@
 c   -----------------------------------------------------
 
-#ifdef CPP_PHYS
 c+jld
 
@@ -689,6 +697,9 @@
       IF (ip_ebil_dyn.ge.1 ) THEN 
           ztit='bil dyn'
-          CALL diagedyn(ztit,2,1,1,dtphys
-     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
+! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
+           IF (planet_type.eq."earth") THEN
+            CALL diagedyn(ztit,2,1,1,dtphys
+     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
+           ENDIF
       ENDIF 
 c-jld
@@ -725,5 +736,5 @@
         
 c        call SetDistrib(jj_nb_vanleer)
-        do j=1,nqmx
+        do j=1,nqtot
  
           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
@@ -756,5 +767,5 @@
 cc$OMP BARRIER
 !        CALL FTRACE_REGION_BEGIN("calfis")
-        CALL calfis_p( nq, lafin ,rdayvrai,time  ,
+        CALL calfis_p(lafin ,rdayvrai,time  ,
      $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
      $               du,dv,dteta,dq,
@@ -777,5 +788,5 @@
           dpfi_tmp(1:iip1)     = dpfi(ijb:ijb+iim)  
 c$OMP END MASTER
-        endif
+        endif ! of if ( .not. pole_nord)
 
 c$OMP BARRIER
@@ -799,5 +810,5 @@
      *                      1,0,0,1,Request_physic)
 
-        do j=1,nqmx
+        do j=1,nqtot
           call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
      *                        1,0,0,1,Request_physic)
@@ -833,5 +844,5 @@
 c$OMP END MASTER
           
-        endif
+        endif ! of if (.not. pole_nord)
 c$OMP BARRIER
 cc$OMP MASTER        
@@ -842,5 +853,5 @@
 cc$OMP END MASTER
 c      
-c      do j=1,nqmx
+c      do j=1,nqtot
 c        call WriteField_p('dqfi'//trim(int2str(j)),
 c     .                reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
@@ -853,5 +864,5 @@
          ENDIF
        
-          CALL addfi_p( nqmx, dtphys, leapf, forward   ,
+          CALL addfi_p( dtphys, leapf, forward   ,
      $                  ucov, vcov, teta , q   ,ps ,
      $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
@@ -889,5 +900,5 @@
      *                               jj_Nb_caldyn,Request_physic)
 
-        do j=1,nqmx
+        do j=1,nqtot
         
           call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
@@ -922,6 +933,16 @@
 cc$OMP END MASTER
 
-#else
-
+
+c-jld
+c$OMP MASTER
+         call resume_timer(timer_caldyn)
+         if (FirstPhysic) then
+           ok_start_timer=.TRUE.
+           FirstPhysic=.false.
+         endif
+c$OMP END MASTER
+       ENDIF ! of IF( apphys )
+
+      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
 c   Calcul academique de la physique = Rappel Newtonien + fritcion 
 c   --------------------------------------------------------------
@@ -939,21 +960,9 @@
 
        call friction_p(ucov,vcov,iphysiq*dtvr)
-
-#endif
-
-c-jld
-c$OMP MASTER
-         call resume_timer(timer_caldyn)
-         if (FirstPhysic) then
-           ok_start_timer=.TRUE.
-           FirstPhysic=.false.
-         endif
-c$OMP END MASTER
-       ENDIF
+      ENDIF ! of IF(iflag_phys.EQ.2)
+
 
         CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
 c$OMP BARRIER
-
-
         CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
 c$OMP BARRIER
@@ -1285,4 +1294,5 @@
                ENDIF
 #ifdef CPP_IOIPSL
+             IF (ok_dynzon) THEN 
              call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
              call SendRequest(TestRequest)
@@ -1291,8 +1301,12 @@
 c$OMP BARRIER
 c$OMP MASTER
-              CALL writedynav_p(histaveid, nqmx, itau,vcov ,
+              CALL writedynav_p(histaveid, itau,vcov ,
      ,                          ucov,teta,pk,phi,q,masse,ps,phis)
-c$OMP END MASTER
-
+
+c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
+              CALL bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav, 
+     ,             ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q) 
+c$OMP END MASTER
+              ENDIF !ok_dynzon
 #endif
             ENDIF
@@ -1304,45 +1318,59 @@
 c      IF( MOD(itau,iecri         ).EQ.0) THEN
 
-           IF( MOD(itau,iecri*day_step).EQ.0) THEN
-c$OMP BARRIER
-c$OMP MASTER
-               nbetat = nbetatdem
-        CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi )
+            IF( MOD(itau,iecri*day_step).EQ.0) THEN
+c$OMP BARRIER
+c$OMP MASTER
+              nbetat = nbetatdem
+              CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
        
 cym        unat=0.
         
-        ijb=ij_begin
-        ije=ij_end
-        
-        if (pole_nord) then
-          ijb=ij_begin+iip1
-          unat(1:iip1,:)=0.
-        endif
-        
-        if (pole_sud) then 
-          ije=ij_end-iip1
-          unat(ij_end-iip1+1:ij_end,:)=0.
-        endif
+              ijb=ij_begin
+              ije=ij_end
+        
+              if (pole_nord) then
+                ijb=ij_begin+iip1
+                unat(1:iip1,:)=0.
+              endif
+        
+              if (pole_sud) then 
+                ije=ij_end-iip1
+                unat(ij_end-iip1+1:ij_end,:)=0.
+              endif
             
-        do l=1,llm
-           unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
-        enddo
-
-        ijb=ij_begin
-        ije=ij_end
-        if (pole_sud) ije=ij_end-iip1
-        
-        do l=1,llm
-           vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
-        enddo
+              do l=1,llm
+                unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
+              enddo
+
+              ijb=ij_begin
+              ije=ij_end
+              if (pole_sud) ije=ij_end-iip1
+        
+              do l=1,llm
+                vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
+              enddo
         
 #ifdef CPP_IOIPSL
  
-        CALL writehist_p(histid,histvid, nqmx,itau,vcov, 
-     s                       ucov,teta,phi,q,masse,ps,phis)
+              CALL writehist_p(histid,histvid, itau,vcov, 
+     &                         ucov,teta,phi,q,masse,ps,phis)
 
 #endif
-c$OMP END MASTER
-           ENDIF
+! For some Grads outputs of fields
+              if (output_grads_dyn) then
+! Ehouarn: hope this works the way I think it does:
+                  call Gather_Field(unat,ip1jmp1,llm,0)
+                  call Gather_Field(vnat,ip1jm,llm,0)
+                  call Gather_Field(teta,ip1jmp1,llm,0)
+                  call Gather_Field(ps,ip1jmp1,1,0)
+                  do iq=1,nqtot
+                    call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
+                  enddo
+                  if (mpi_rank==0) then
+#include "write_grads_dyn.h"
+                  endif
+              endif ! of if (output_grads_dyn)
+c$OMP END MASTER
+            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
 
             IF(itau.EQ.itaufin) THEN
@@ -1351,13 +1379,16 @@
 c$OMP MASTER
 
-c#ifdef CPP_IOIPSL
-
-       CALL dynredem1_p("restart.nc",0.0,
-     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
-c#endif
+              if (planet_type.eq."earth") then
+#ifdef CPP_EARTH
+! Write an Earth-format restart file
+                CALL dynredem1_p("restart.nc",0.0,
+     &                           vcov,ucov,teta,q,masse,ps)
+
+#endif
+              endif ! of if (planet_type.eq."earth")
 
               CLOSE(99)
 c$OMP END MASTER
-            ENDIF
+            ENDIF ! of IF (itau.EQ.itaufin)
 
 c-----------------------------------------------------------------------
@@ -1390,7 +1421,9 @@
                  dt  = 2.*dtvr
                  GO TO 2
-            END IF
-
-      ELSE
+            END IF ! of IF (MOD(itau,iperiod).EQ.0)
+                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
+
+
+      ELSE ! of IF (.not.purmats)
 
 c       ........................................................
@@ -1419,7 +1452,7 @@
                GO TO 2
 
-            ELSE
-
-            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
+            ELSE ! of IF(forward)
+
+              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
                IF(itau.EQ.itaufin) THEN
                   iav=1
@@ -1428,20 +1461,23 @@
                ENDIF
 #ifdef CPP_IOIPSL
-c$OMP BARRIER
-
-              call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
-              call SendRequest(TestRequest)
-c$OMP BARRIER
-              call WaitRequest(TestRequest)
-
-c$OMP BARRIER
-c$OMP MASTER
-              CALL writedynav_p(histaveid, nqmx, itau,vcov ,
+               IF (ok_dynzon) THEN
+c$OMP BARRIER
+
+               call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
+               call SendRequest(TestRequest)
+c$OMP BARRIER
+               call WaitRequest(TestRequest)
+
+c$OMP BARRIER
+c$OMP MASTER
+               CALL writedynav_p(histaveid, itau,vcov ,
      ,                          ucov,teta,pk,phi,q,masse,ps,phis)
-               call bilan_dyn_p (2,dtvr*iperiod,dtvr*day_step*periodav,
+               CALL bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
      ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
 c$OMP END MASTER
+               END IF !ok_dynzon
 #endif
-            ENDIF
+              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
+
 
 c               IF(MOD(itau,iecri         ).EQ.0) THEN
@@ -1449,69 +1485,76 @@
 c$OMP BARRIER
 c$OMP MASTER
-                  nbetat = nbetatdem
-       CALL geopot_p( ip1jmp1, teta  , pk , pks,  phis  , phi   )
+                nbetat = nbetatdem
+                CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
 
 cym        unat=0.
-        ijb=ij_begin
-        ije=ij_end
-        
-        if (pole_nord) then
-          ijb=ij_begin+iip1
-          unat(1:iip1,:)=0.
-        endif
-        
-        if (pole_sud) then 
-          ije=ij_end-iip1
-          unat(ij_end-iip1+1:ij_end,:)=0.
-        endif
+                ijb=ij_begin
+                ije=ij_end
+        
+                if (pole_nord) then
+                  ijb=ij_begin+iip1
+                  unat(1:iip1,:)=0.
+                endif
+        
+                if (pole_sud) then 
+                  ije=ij_end-iip1
+                  unat(ij_end-iip1+1:ij_end,:)=0.
+                endif
             
-        do l=1,llm
-           unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
-        enddo
-
-        ijb=ij_begin
-        ije=ij_end
-        if (pole_sud) ije=ij_end-iip1
-        
-        do l=1,llm
-           vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
-        enddo
+                do l=1,llm
+                  unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
+                enddo
+
+                ijb=ij_begin
+                ije=ij_end
+                if (pole_sud) ije=ij_end-iip1
+        
+                do l=1,llm
+                  vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
+                enddo
 
 #ifdef CPP_IOIPSL
 
-       CALL writehist_p( histid, histvid, nqmx, itau,vcov , 
-     ,                           ucov,teta,phi,q,masse,ps,phis)
-c#else
-c      call Gather_Field(unat,ip1jmp1,llm,0)
-c      call Gather_Field(vnat,ip1jm,llm,0)
-c      call Gather_Field(teta,ip1jmp1,llm,0)
-c      call Gather_Field(ps,ip1jmp1,1,0)
-c      do iq=1,nqmx
-c        call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
-c      enddo
+                CALL writehist_p(histid, histvid, itau,vcov , 
+     &                           ucov,teta,phi,q,masse,ps,phis)
+#endif
+! For some Grads output (but does it work?)
+                if (output_grads_dyn) then
+                  call Gather_Field(unat,ip1jmp1,llm,0)
+                  call Gather_Field(vnat,ip1jm,llm,0)
+                  call Gather_Field(teta,ip1jmp1,llm,0)
+                  call Gather_Field(ps,ip1jmp1,1,0)
+                  do iq=1,nqtot
+                    call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
+                  enddo
 c      
-c      if (mpi_rank==0) then
-c#include "write_grads_dyn.h"
-c      endif
+                  if (mpi_rank==0) then
+#include "write_grads_dyn.h"
+                  endif
+                endif ! of if (output_grads_dyn)
+
+c$OMP END MASTER
+              ENDIF ! of IF(MOD(itau,iecri*day_step).EQ.0)
+
+              IF(itau.EQ.itaufin) THEN
+                if (planet_type.eq."earth") then
+#ifdef CPP_EARTH
+c$OMP MASTER
+                   CALL dynredem1_p("restart.nc",0.0,
+     .                               vcov,ucov,teta,q,masse,ps)
+c$OMP END MASTER
 #endif
-
-c$OMP END MASTER
-               ENDIF
-
-                 IF(itau.EQ.itaufin) THEN
-c$OMP MASTER
-                   CALL dynredem1_p("restart.nc",0.0,
-     .                               vcov,ucov,teta,q,nqmx,masse,ps)
-c$OMP END MASTER
-                 ENDIF
-                 forward = .TRUE.
-                 GO TO  1
-
-            ENDIF
-
-      END IF
-c$OMP MASTER
-        call finalize_parallel
-c$OMP END MASTER
-        RETURN
+                endif ! of if (planet_type.eq."earth")
+              ENDIF ! of IF(itau.EQ.itaufin)
+
+              forward = .TRUE.
+              GO TO  1
+
+            ENDIF ! of IF (forward)
+
+      END IF ! of IF(.not.purmats)
+c$OMP MASTER
+      call finalize_parallel
+c$OMP END MASTER
+      RETURN
       END
Index: LMDZ4/trunk/libf/dyn3dpar/parallel.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/parallel.F90	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/parallel.F90	(revision 1146)
@@ -51,4 +51,11 @@
 #else
        using_mpi=.FALSE.
+#endif
+      
+
+#ifdef CPP_OMP
+       using_OMP=.TRUE.
+#else
+       using_OMP=.FALSE.
 #endif
       
Index: LMDZ4/trunk/libf/dyn3dpar/qminimum_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/qminimum_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/qminimum_p.F	(revision 1146)
@@ -50,7 +50,8 @@
       DO 1000 k = 1, llm
       DO 1040 i = ijb, ije
-            zx_defau      = AMAX1( seuil_liq - q(i,k,iq_liq), 0.0 )
-            q(i,k,iq_vap) = q(i,k,iq_vap) - zx_defau
-            q(i,k,iq_liq) = q(i,k,iq_liq) + zx_defau
+            if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
+               q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq
+               q(i,k,iq_liq) = seuil_liq
+            endif
  1040 CONTINUE
  1000 CONTINUE
@@ -69,8 +70,9 @@
 c$OMP DO SCHEDULE(STATIC)
       DO i = ijb, ije
-         zx_abc = deltap(i,k)/deltap(i,k-1)
-         zx_defau    = AMAX1( seuil_vap - q(i,k,iq), 0.0 )
-         q(i,k-1,iq) =  q(i,k-1,iq) - zx_defau * zx_abc
-         q(i,k,iq)   =  q(i,k,iq)   + zx_defau  
+         if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then
+            q(i,k-1,iq) =  q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) *
+     &           deltap(i,k) / deltap(i,k-1)
+            q(i,k,iq)   =  seuil_vap  
+         endif
       ENDDO
 c$OMP END DO NOWAIT
Index: LMDZ4/trunk/libf/dyn3dpar/read_reanalyse.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/read_reanalyse.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/read_reanalyse.F	(revision 1146)
@@ -11,4 +11,5 @@
 
        USE parallel
+       use netcdf
 c -----------------------------------------------------------------
 c   Declarations
@@ -72,7 +73,7 @@
             print *,'Vous êtes entrain de lire des données sur 
      .               niveaux modèle'
-            ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcode)
-            varidap=NCVID(ncidpl,'AP',rcode)
-            varidbp=NCVID(ncidpl,'BP',rcode)
+            rcode=nf90_open('apbp.nc',nf90_nowrite,ncidpl)
+            rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
+            rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
             print*,'ncidpl,varidap',ncidpl,varidap
             endif
@@ -80,6 +81,6 @@
 c Vent zonal
             if (guide_u) then
-               ncidu=NCOPN('u.nc',NCNOWRIT,rcode)
-               varidu=NCVID(ncidu,'UWND',rcode)
+               rcode=nf90_open('u.nc',nf90_nowrite,ncidu)
+               rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
                print*,'ncidu,varidu',ncidu,varidu
                if (ncidpl.eq.-99) ncidpl=ncidu
@@ -88,6 +89,6 @@
 c Vent meridien
             if (guide_v) then
-               ncidv=NCOPN('v.nc',NCNOWRIT,rcode)
-               varidv=NCVID(ncidv,'VWND',rcode)
+               rcode=nf90_open('v.nc',nf90_nowrite,ncidv)
+               rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
                print*,'ncidv,varidv',ncidv,varidv
                if (ncidpl.eq.-99) ncidpl=ncidv
@@ -96,6 +97,6 @@
 c Temperature
             if (guide_T) then
-               ncidt=NCOPN('T.nc',NCNOWRIT,rcode)
-               varidt=NCVID(ncidt,'AIR',rcode)
+               rcode=nf90_open('T.nc',nf90_nowrite,ncidt)
+               rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
                print*,'ncidt,varidt',ncidt,varidt
                if (ncidpl.eq.-99) ncidpl=ncidt
@@ -104,6 +105,6 @@
 c Humidite
             if (guide_Q) then
-               ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)
-               varidQ=NCVID(ncidQ,'RH',rcode)
+               rcode=nf90_open('hur.nc',nf90_nowrite,ncidQ)
+               rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
                print*,'ncidQ,varidQ',ncidQ,varidQ
                if (ncidpl.eq.-99) ncidpl=ncidQ
@@ -112,6 +113,6 @@
 c Pression de surface
             if ((guide_P).OR.(guide_modele)) then
-               ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
-               varidps=NCVID(ncidps,'SP',rcode)
+               rcode=nf90_open('ps.nc',nf90_nowrite,ncidps)
+               rcode = nf90_inq_varid(ncidps, 'SP', varidps)
                print*,'ncidps,varidps',ncidps,varidps
             endif
@@ -119,12 +120,12 @@
 c Coordonnee verticale
             if (.not.guide_modele) then
-              if (ncep) then
-               print*,'Vous etes entrain de lire des donnees NCEP'
-               varidpl=NCVID(ncidpl,'LEVEL',rcode)
-            else
-               print*,'Vous etes entrain de lire des donnees ECMWF'
-               varidpl=NCVID(ncidpl,'PRESSURE',rcode)
-              endif
-              print*,'ncidpl,varidpl',ncidpl,varidpl
+               if (ncep) then
+                  print*,'Vous etes entrain de lire des donnees NCEP'
+                  rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
+               else
+                  print*,'Vous etes entrain de lire des donnees ECMWF'
+                  rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
+               endif
+               print*,'ncidpl,varidpl',ncidpl,varidpl
             endif
             print*,'ncidu,varidpl',ncidu,varidpl
Index: LMDZ4/trunk/libf/dyn3dpar/serre.h
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/serre.h	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/serre.h	(revision 1146)
@@ -2,10 +2,10 @@
 ! $Header$
 !
-c
-c
-c..include serre.h
-c
-       REAL clon,clat,transx,transy,alphax,alphay,pxo,pyo,
-     ,  grossismx, grossismy, dzoomx, dzoomy,taux,tauy
-       COMMON/serre/clon,clat,transx,transy,alphax,alphay,pxo,pyo , 
-     ,  grossismx, grossismy, dzoomx, dzoomy,taux,tauy
+!c
+!c
+!c..include serre.h
+!c
+       REAL clon,clat,transx,transy,alphax,alphay,pxo,pyo,              &
+     &  grossismx, grossismy, dzoomx, dzoomy,taux,tauy
+       COMMON/serre/clon,clat,transx,transy,alphax,alphay,pxo,pyo ,     &
+     &  grossismx, grossismy, dzoomx, dzoomy,taux,tauy
Index: LMDZ4/trunk/libf/dyn3dpar/test_period.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/test_period.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/test_period.F	(revision 1146)
@@ -3,4 +3,5 @@
 !
       SUBROUTINE test_period ( ucov, vcov, teta, q, p, phis )
+      USE infotrac, ONLY : nqtot
 c
 c     Auteur : P. Le Van  
@@ -17,5 +18,5 @@
 c
       REAL ucov(ip1jmp1,llm), vcov(ip1jm,llm), teta(ip1jmp1,llm) ,
-     ,      q(ip1jmp1,llm,nqmx), p(ip1jmp1,llmp1), phis(ip1jmp1)
+     ,      q(ip1jmp1,llm,nqtot), p(ip1jmp1,llmp1), phis(ip1jmp1)
 c
 c   .....  Variables  locales  .....
@@ -68,5 +69,5 @@
       
 c
-      DO nq =1, nqmx
+      DO nq =1, nqtot
         DO l =1, llm
           DO ij = 1, ip1jmp1, iip1
Index: LMDZ4/trunk/libf/dyn3dpar/times.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/times.F90	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/times.F90	(revision 1146)
@@ -207,5 +207,5 @@
       endif
       
-    ENDIF  ! using_mpî
+    ENDIF  ! using_mp�
   end subroutine allgather_timer_average
   
@@ -222,5 +222,5 @@
   end subroutine InitTime
   
-  function DiffTime
+  function DiffTime()
   implicit none
     double precision :: DiffTime
@@ -236,5 +236,5 @@
   end function DiffTime
   
-  function DiffCpuTime
+  function DiffCpuTime()
   implicit none
     real :: DiffCpuTime
Index: LMDZ4/trunk/libf/dyn3dpar/vlspltgen_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/vlspltgen_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/vlspltgen_p.F	(revision 1146)
@@ -26,4 +26,5 @@
       USE Write_Field_p
       USE VAMPIR
+      USE infotrac, ONLY : nqtot
       IMPLICIT NONE
 
@@ -38,8 +39,8 @@
 c   Arguments:
 c   ----------
-      INTEGER iadv(nqmx)
+      INTEGER iadv(nqtot)
       REAL masse(ip1jmp1,llm),pente_max
       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
-      REAL q(ip1jmp1,llm,nqmx)
+      REAL q(ip1jmp1,llm,nqtot)
       REAL w(ip1jmp1,llm),pdt
       REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
@@ -51,9 +52,9 @@
 c
       REAL,SAVE :: qsat(ip1jmp1,llm)
-      REAL,SAVE :: zm(ip1jmp1,llm,nqmx)
+      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zm
       REAL,SAVE :: mu(ip1jmp1,llm)
       REAL,SAVE :: mv(ip1jm,llm)
       REAL,SAVE :: mw(ip1jmp1,llm+1)
-      REAL,SAVE :: zq(ip1jmp1,llm,nqmx)
+      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zq
       REAL zzpbar, zzw
 
@@ -67,4 +68,6 @@
       REAL tempe(ip1jmp1)
       INTEGER ijb,ije,iq
+      LOGICAL, SAVE :: firstcall=.TRUE.
+!$OMP THREADPRIVATE(firstcall)
       type(request) :: MyRequest1
       type(request) :: MyRequest2
@@ -84,4 +87,14 @@
         rtt  = 273.16
 
+c Allocate variables depending on dynamic variable nqtot
+
+         IF (firstcall) THEN
+            firstcall=.FALSE.
+!$OMP MASTER
+            ALLOCATE(zm(ip1jmp1,llm,nqtot))
+            ALLOCATE(zq(ip1jmp1,llm,nqtot))
+!$OMP END MASTER
+!$OMP BARRIER
+         END IF
 c-- Calcul de Qsat en chaque point
 c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
@@ -164,5 +177,5 @@
        ije=ij_end
 
-      DO iq=1,nqmx
+      DO iq=1,nqtot
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
         DO l=1,llm
@@ -175,5 +188,5 @@
 
 c$OMP BARRIER           
-      DO iq=1,nqmx
+      DO iq=1,nqtot
 
         if(iadv(iq) == 0) then
@@ -245,5 +258,5 @@
 c$OMP END MASTER       
 c$OMP BARRIER
-      do iq=1,nqmx
+      do iq=1,nqtot
 
         if(iadv(iq) == 0) then
@@ -285,5 +298,5 @@
 c$OMP BARRIER
  
-      do iq=1,nqmx
+      do iq=1,nqtot
 
         if(iadv(iq) == 0) then
@@ -308,5 +321,5 @@
 
 
-      do iq=1,nqmx
+      do iq=1,nqtot
 
         if(iadv(iq) == 0) then 
@@ -359,5 +372,5 @@
 
 c$OMP BARRIER
-      do iq=1,nqmx
+      do iq=1,nqtot
 
         if(iadv(iq) == 0) then
@@ -398,5 +411,5 @@
 
 
-      do iq=1,nqmx
+      do iq=1,nqtot
 
         if(iadv(iq) == 0) then
@@ -420,5 +433,5 @@
        enddo
 
-      do iq=1,nqmx
+      do iq=1,nqtot
 
         if(iadv(iq) == 0) then 
@@ -450,5 +463,5 @@
 
 
-      DO iq=1,nqmx
+      DO iq=1,nqtot
 
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)          
Index: LMDZ4/trunk/libf/dyn3dpar/write_grads_dyn.h
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/write_grads_dyn.h	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/write_grads_dyn.h	(revision 1146)
@@ -24,5 +24,5 @@
       string10='teta'
       CALL wrgrads(1,llm,teta,string10,string10)
-      do iq=1,nqmx
+      do iq=1,nqtot
          string10='q'
          write(string10(2:2),'(i1)') iq
Index: LMDZ4/trunk/libf/dyn3dpar/writedynav_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/writedynav_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/writedynav_p.F	(revision 1146)
@@ -2,5 +2,5 @@
 ! $Header$
 !
-      subroutine writedynav_p( histid, nq, time, vcov, 
+      subroutine writedynav_p( histid, time, vcov, 
      ,                          ucov,teta,ppk,phi,q,masse,ps,phis)
 
@@ -8,4 +8,5 @@
       USE parallel
       USE misc_mod
+      USE infotrac
       implicit none
 
@@ -17,5 +18,4 @@
 C   Entree:
 C      histid: ID du fichier histoire
-C      nqmx: nombre maxi de traceurs
 C      time: temps de l'ecriture
 C      vcov: vents v covariants
@@ -47,5 +47,4 @@
 #include "description.h"
 #include "serre.h"
-#include "advtrac.h"
 
 C
@@ -53,10 +52,10 @@
 C
 
-      INTEGER histid, nq
+      INTEGER histid
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) 
       REAL teta(ip1jmp1,llm),phi(ip1jmp1,llm),ppk(ip1jmp1,llm)                  
       REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
       REAL phis(ip1jmp1)                  
-      REAL q(ip1jmp1,llm,nq)
+      REAL q(ip1jmp1,llm,nqtot)
       integer time
 
@@ -105,6 +104,4 @@
 C  Vents V scalaire
 C
-      if (pole_sud) ije=ij_end-iip1
-      if (pole_sud) jjn=jj_nb-1
       
       call gr_v_scal_p(llm, vnat, vs)
@@ -114,8 +111,5 @@
 C  Temperature potentielle moyennee
 C
-      ijb=ij_begin
-      ije=ij_end
-      jjn=jj_nb
-      
+     
       call histwrite(histid, 'theta', itau_w, teta(ijb:ije,:), 
      .                iip1*jjn*llm, ndex3d)
@@ -139,5 +133,5 @@
 C  Traceurs
 C
-        DO iq=1,nq
+        DO iq=1,nqtot
           call histwrite(histid, ttext(iq), itau_w, q(ijb:ije,:,iq), 
      .                   iip1*jjn*llm, ndex3d)
Index: LMDZ4/trunk/libf/dyn3dpar/writehist_p.F
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/writehist_p.F	(revision 1122)
+++ LMDZ4/trunk/libf/dyn3dpar/writehist_p.F	(revision 1146)
@@ -2,5 +2,5 @@
 ! $Header$
 !
-      subroutine writehist_p( histid, histvid, nq, time, vcov, 
+      subroutine writehist_p( histid, histvid, time, vcov, 
      ,                          ucov,teta,phi,q,masse,ps,phis)
 
@@ -8,4 +8,5 @@
       USE parallel
       USE misc_mod
+      USE infotrac
       implicit none
 
@@ -18,5 +19,4 @@
 C      histid: ID du fichier histoire
 C      histvid:ID du fichier histoire pour les vents V (appele a disparaitre)
-C      nqmx: nombre maxi de traceurs
 C      time: temps de l'ecriture
 C      vcov: vents v covariants
@@ -48,5 +48,4 @@
 #include "description.h"
 #include "serre.h"
-#include "advtrac.h"
 
 C
@@ -54,10 +53,10 @@
 C
 
-      INTEGER histid, nq, histvid
+      INTEGER histid, histvid
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) 
       REAL teta(ip1jmp1,llm),phi(ip1jmp1,llm)                   
       REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
       REAL phis(ip1jmp1)                  
-      REAL q(ip1jmp1,llm,nq)
+      REAL q(ip1jmp1,llm,nqtot)
       integer time
 
@@ -119,5 +118,5 @@
 C  Traceurs
 C
-        DO iq=1,nq
+        DO iq=1,nqtot
           call histwrite(histid, ttext(iq), itau_w, q(ijb:ije,:,iq), 
      .                   iip1*jjn*llm, ndexu)
