Index: LMDZ5/branches/testing/libf/dyn3dmem/advtrac_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/advtrac_loc.F	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/advtrac_loc.F	(revision 2298)
@@ -24,5 +24,5 @@
       USE Vampir
       USE times
-      USE infotrac, ONLY: nqtot, iadv
+      USE infotrac, ONLY: nqtot, iadv, ok_iso_verif
       USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
       USE advtrac_mod, ONLY: finmasse
@@ -82,5 +82,5 @@
 !$OMP THREADPRIVATE(testRequest)
 
-c  test sur l'eventuelle creation de valeurs negatives de la masse
+c  test sur l''eventuelle creation de valeurs negatives de la masse
          ijb=ij_begin
          ije=ij_end
@@ -155,7 +155,13 @@
 c$OMP BARRIER
                  
+          !write(*,*) 'advtrac 157: appel de vlspltgen_loc'
           call vlspltgen_loc( q,iadv, 2., massem, wg ,
      *                        pbarug,pbarvg,dtvr,p,
      *                        pk,teta )
+
+          !write(*,*) 'advtrac 162: apres appel vlspltgen_loc'
+      if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'advtrac 162')
+      endif !if (ok_iso_verif) then
 
 #ifdef DEBUG_IO     
@@ -356,5 +362,6 @@
 c$OMP END DO
 
-       CALL qminimum_loc( q, 2, finmasse )
+        ! CRisi: on passe nqtot et non nq
+       CALL qminimum_loc( q, nqtot, finmasse )
 
       endif ! of if (planet_type=="earth")
Index: LMDZ5/branches/testing/libf/dyn3dmem/caladvtrac_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/caladvtrac_loc.F	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/caladvtrac_loc.F	(revision 2298)
@@ -56,5 +56,5 @@
 !$OMP THREADPRIVATE(Request_vanleer)
 
-           
+      !write(*,*) 'caladvtrac 58: entree'     
       ijbu=ij_begin
       ijeu=ij_end
@@ -109,4 +109,5 @@
 
       IF ( iadvtr.EQ.iapp_tracvl ) THEN
+      !write(*,*) 'caladvtrac 133'
 c$OMP MASTER
         call suspend_timer(timer_caldyn)
@@ -183,7 +184,9 @@
          CALL WriteField_u('wg1',wg_adv)
 #endif        
+      !write(*,*) 'caladvtrac 185'  
       CALL advtrac_loc( pbarug_adv,pbarvg_adv,wg_adv, 
      *             p_adv,  massem_adv,q_adv, teta_adv,
-     .             pk_adv)
+     .             pk_adv)     
+      !write(*,*) 'caladvtrac 189'
 
 
Index: LMDZ5/branches/testing/libf/dyn3dmem/calfis_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/calfis_loc.F	(revision 2258)
+++ 	(revision )
@@ -1,1189 +1,0 @@
-!
-! $Id$
-!
-C
-C
-      SUBROUTINE calfis_loc(lafin,
-     $                  jD_cur, jH_cur,
-     $                  pucov,
-     $                  pvcov,
-     $                  pteta,
-     $                  pq,
-     $                  pmasse,
-     $                  pps,
-     $                  pp,
-     $                  ppk,
-     $                  pphis,
-     $                  pphi,
-     $                  pducov,
-     $                  pdvcov,
-     $                  pdteta,
-     $                  pdq,
-     $                  flxw,
-     $                  pdufi,
-     $                  pdvfi,
-     $                  pdhfi,
-     $                  pdqfi,
-     $                  pdpsfi)
-#ifdef CPP_PHYS
-! If using physics
-c
-c    Auteur :  P. Le Van, F. Hourdin 
-c   .........
-      USE dimphy
-      USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root 
-      USE mod_interface_dyn_phys
-      USE IOPHY
-#endif
-      USE parallel_lmdz,ONLY:omp_chunk,using_mpi,jjb_u,jje_u,jjb_v,jje_v
-      USE Write_Field
-      Use Write_field_p
-      USE Times
-      USE infotrac, ONLY: nqtot, niadv, tname
-      USE control_mod, ONLY: planet_type, nsplit_phys
-
-      IMPLICIT NONE
-c=======================================================================
-c
-c   1. rearrangement des tableaux et transformation
-c      variables dynamiques  >  variables physiques
-c   2. calcul des termes physiques
-c   3. retransformation des tendances physiques en tendances dynamiques
-c
-c   remarques:
-c   ----------
-c
-c    - les vents sont donnes dans la physique par leurs composantes 
-c      naturelles.
-c    - la variable thermodynamique de la physique est une variable
-c      intensive :   T 
-c      pour la dynamique on prend    T * ( preff / p(l) ) **kappa
-c    - les deux seules variables dependant de la geometrie necessaires
-c      pour la physique sont la latitude pour le rayonnement et 
-c      l'aire de la maille quand on veut integrer une grandeur 
-c      horizontalement.
-c    - les points de la physique sont les points scalaires de la 
-c      la dynamique; numerotation:
-c          1 pour le pole nord
-c          (jjm-1)*iim pour l'interieur du domaine
-c          ngridmx pour le pole sud
-c      ---> ngridmx=2+(jjm-1)*iim
-c
-c     Input :
-c     -------
-c       ecritphy        frequence d'ecriture (en jours)de histphy
-c       pucov           covariant zonal velocity
-c       pvcov           covariant meridional velocity 
-c       pteta           potential temperature
-c       pps             surface pressure
-c       pmasse          masse d'air dans chaque maille
-c       pts             surface temperature  (K)
-c       callrad         clef d'appel au rayonnement
-c
-c    Output :
-c    --------
-c        pdufi          tendency for the natural zonal velocity (ms-1)
-c        pdvfi          tendency for the natural meridional velocity 
-c        pdhfi          tendency for the potential temperature
-c        pdtsfi         tendency for the surface temperature
-c
-c        pdtrad         radiative tendencies  \  both input
-c        pfluxrad       radiative fluxes      /  and output
-c
-c=======================================================================
-c
-c-----------------------------------------------------------------------
-c
-c    0.  Declarations :
-c    ------------------
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "temps.h"
-
-      INTEGER ngridmx
-      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
-
-#include "comconst.h"
-#include "comvert.h"
-#include "comgeom2.h"
-#include "iniprint.h"
-#ifdef CPP_MPI
-      include 'mpif.h'
-#endif
-c    Arguments :
-c    -----------
-      LOGICAL,INTENT(IN) ::  lafin ! .true. for the very last call to physics
-      REAL,INTENT(IN):: jD_cur, jH_cur
-      REAL,INTENT(IN):: pvcov(iip1,jjb_v:jje_v,llm) ! covariant meridional velocity
-      REAL,INTENT(IN):: pucov(iip1,jjb_u:jje_u,llm) ! covariant zonal velocity
-      REAL,INTENT(IN):: pteta(iip1,jjb_u:jje_u,llm) ! potential temperature
-      REAL,INTENT(IN):: pmasse(iip1,jjb_u:jje_u,llm) ! mass in each cell ! not used
-      REAL,INTENT(IN):: pq(iip1,jjb_u:jje_u,llm,nqtot) ! tracers
-      REAL,INTENT(IN):: pphis(iip1,jjb_u:jje_u) ! surface geopotential
-      REAL,INTENT(IN):: pphi(iip1,jjb_u:jje_u,llm) ! geopotential
-
-      REAL,INTENT(IN) :: pdvcov(iip1,jjb_v:jje_v,llm) ! dynamical tendency on vcov ! not used
-      REAL,INTENT(IN) :: pducov(iip1,jjb_u:jje_u,llm) ! dynamical tendency on ucov
-      REAL,INTENT(IN) :: pdteta(iip1,jjb_u:jje_u,llm) ! dynamical tendency on teta ! not used
-      REAL,INTENT(IN) :: pdq(iip1,jjb_u:jje_u,llm,nqtot) ! dynamical tendency on tracers ! not used
-
-      REAL,INTENT(IN) :: pps(iip1,jjb_u:jje_u) ! surface pressure (Pa)
-      REAL,INTENT(IN) :: pp(iip1,jjb_u:jje_u,llmp1) ! pressure at mesh interfaces (Pa)
-      REAL,INTENT(IN) :: ppk(iip1,jjb_u:jje_u,llm) ! Exner at mid-layer
-      REAL,INTENT(IN) :: flxw(iip1,jjb_u:jje_u,llm)  ! Vertical mass flux on dynamics grid
-
-      ! tendencies (in */s) from the physics
-      REAL,INTENT(OUT) :: pdvfi(iip1,jjb_v:jje_v,llm) ! tendency on covariant meridional wind
-      REAL,INTENT(OUT) :: pdufi(iip1,jjb_u:jje_u,llm) ! tendency on covariant zonal wind
-      REAL,INTENT(OUT) :: pdhfi(iip1,jjb_u:jje_u,llm) ! tendency on potential temperature (K/s)
-      REAL,INTENT(OUT) :: pdqfi(iip1,jjb_u:jje_u,llm,nqtot) ! tendency on tracers
-      REAL,INTENT(OUT) :: pdpsfi(iip1,jjb_u:jje_u) ! tendency on surface pressure (Pa/s)
-
-#ifdef CPP_PHYS
-! Ehouarn: for now calfis_p needs some informations from physics to compile
-c    Local variables :
-c    -----------------
-
-      INTEGER i,j,l,ig0,ig,iq,iiq
-      REAL,ALLOCATABLE,SAVE :: zpsrf(:)
-      REAL,ALLOCATABLE,SAVE :: zplev(:,:),zplay(:,:)
-      REAL,ALLOCATABLE,SAVE :: zphi(:,:),zphis(:)
-c
-      REAL,ALLOCATABLE,SAVE :: zufi(:,:), zvfi(:,:)
-      REAL,ALLOCATABLE,SAVE :: ztfi(:,:),zqfi(:,:,:)
-c
-      REAL,ALLOCATABLE,SAVE :: pcvgu(:,:), pcvgv(:,:)
-      REAL,ALLOCATABLE,SAVE :: pcvgt(:,:), pcvgq(:,:,:)
-c
-      REAL,ALLOCATABLE,SAVE :: zdufi(:,:),zdvfi(:,:)
-      REAL,ALLOCATABLE,SAVE :: zdtfi(:,:),zdqfi(:,:,:)
-      REAL,ALLOCATABLE,SAVE :: zdpsrf(:)
-      REAL,SAVE,ALLOCATABLE ::  flxwfi(:,:)     ! Flux de masse verticale sur la grille physiq
-
-c
-      REAL,ALLOCATABLE,SAVE :: zplev_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zplay_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zphi_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zphis_omp(:)
-      REAL,ALLOCATABLE,SAVE :: presnivs_omp(:)
-      REAL,ALLOCATABLE,SAVE :: zufi_omp(:,:) 
-      REAL,ALLOCATABLE,SAVE :: zvfi_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: ztfi_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zqfi_omp(:,:,:)
-      REAL,ALLOCATABLE,SAVE :: zdufi_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zdvfi_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zdtfi_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zdqfi_omp(:,:,:)
-      REAL,ALLOCATABLE,SAVE :: zdpsrf_omp(:)
-      REAL,SAVE,ALLOCATABLE ::  flxwfi_omp(:,:)     ! Flux de masse verticale sur la grille physiq
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Introduction du splitting (FH)
-! Question pour Yann :
-! J'ai été surpris au début que les tableaux zufi_omp, zdufi_omp n'co soitent
-! en SAVE. Je crois comprendre que c'est parce que tu voulais qu'il
-! soit allocatable (plutot par exemple que de passer une dimension
-! dépendant du process en argument des routines) et que, du coup,
-! le SAVE évite d'avoir à refaire l'allocation à chaque appel.
-! Tu confirmes ?
-! J'ai suivi le même principe pour les zdufic_omp
-! Mais c'est surement bien que tu controles.
-! 
-
-      REAL,ALLOCATABLE,SAVE :: zdufic_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zdvfic_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zdtfic_omp(:,:)
-      REAL,ALLOCATABLE,SAVE :: zdqfic_omp(:,:,:)
-      REAL jH_cur_split,zdt_split
-      LOGICAL debut_split,lafin_split
-      INTEGER isplit
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-c$OMP THREADPRIVATE(zplev_omp,zplay_omp,zphi_omp,zphis_omp,
-c$OMP+                 presnivs_omp,zufi_omp,zvfi_omp,ztfi_omp,
-c$OMP+                 zqfi_omp,zdufi_omp,zdvfi_omp,
-c$OMP+                 zdtfi_omp,zdqfi_omp,zdpsrf_omp,flxwfi_omp,
-c$OMP+                 zdufic_omp,zdvfic_omp,zdtfic_omp,zdqfic_omp)       
-
-      LOGICAL,SAVE :: first_omp=.true.
-c$OMP THREADPRIVATE(first_omp)
-      
-      REAL zsin(iim),zcos(iim),z1(iim)
-      REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
-      REAL unskap, pksurcp
-c
-      REAL SSUM
-
-      LOGICAL,SAVE :: firstcal=.true., debut=.true.
-c$OMP THREADPRIVATE(firstcal,debut)
-      
-      REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
-      INTEGER :: ierr
-#ifdef CPP_MPI
-      INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
-#else
-      INTEGER,dimension(1,4) :: Status
-#endif
-      INTEGER, dimension(4) :: Req
-      REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
-      integer :: k,kstart,kend
-      INTEGER :: offset
-
-      LOGICAL tracerdyn  
-c
-c-----------------------------------------------------------------------
-c
-c    1. Initialisations :
-c    --------------------
-c
-
-      klon=klon_mpi
-      
-c
-      IF ( firstcal )  THEN
-        debut = .TRUE.
-        IF (ngridmx.NE.2+(jjm-1)*iim) THEN
-          write(lunout,*) 'STOP dans calfis' 
-          write(lunout,*)  
-     &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim' 
-          write(lunout,*) '  ngridmx  jjm   iim   ' 
-          write(lunout,*) ngridmx,jjm,iim
-          STOP
-        ENDIF
-c$OMP MASTER
-      ALLOCATE(zpsrf(klon))
-      ALLOCATE(zplev(klon,llm+1),zplay(klon,llm))
-      ALLOCATE(zphi(klon,llm),zphis(klon))
-      ALLOCATE(zufi(klon,llm), zvfi(klon,llm))
-      ALLOCATE(ztfi(klon,llm),zqfi(klon,llm,nqtot))
-      ALLOCATE(pcvgu(klon,llm), pcvgv(klon,llm))
-      ALLOCATE(pcvgt(klon,llm), pcvgq(klon,llm,2))
-      ALLOCATE(zdufi(klon,llm),zdvfi(klon,llm))
-      ALLOCATE(zdtfi(klon,llm),zdqfi(klon,llm,nqtot))
-      ALLOCATE(zdpsrf(klon))
-      ALLOCATE(zdufi2(klon+iim,llm),zdvfi2(klon+iim,llm))
-      ALLOCATE(flxwfi(klon,llm))
-c$OMP END MASTER
-c$OMP BARRIER	  
-      ELSE
-          debut = .FALSE.
-      ENDIF
-
-c
-c
-c-----------------------------------------------------------------------
-c   40. transformation des variables dynamiques en variables physiques:
-c   ---------------------------------------------------------------
-
-c   41. pressions au sol (en Pascals)
-c   ----------------------------------
-
-c$OMP MASTER
-      call start_timer(timer_physic)
-c$OMP END MASTER
-
-c$OMP MASTER             
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-      do ig0=1,klon
-        i=index_i(ig0)
-        j=index_j(ig0)
-        zpsrf(ig0)=pps(i,j)
-      enddo
-c$OMP END MASTER
-
-
-c   42. pression intercouches :
-c
-c   -----------------------------------------------------------------
-c     .... zplev  definis aux (llm +1) interfaces des couches  ....
-c     .... zplay  definis aux (  llm )    milieux des couches  .... 
-c   -----------------------------------------------------------------
-
-c    ...    Exner = cp * ( p(l) / preff ) ** kappa     ....
-c
-       unskap   = 1./ kappa
-c
-c      print *,omp_rank,'klon--->',klon
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l = 1, llmp1
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-        do ig0=1,klon
-          i=index_i(ig0)
-          j=index_j(ig0)
-          zplev( ig0,l ) = pp(i,j,l)
-        enddo
-      ENDDO
-c$OMP END DO NOWAIT
-c
-c
-
-c   43. temperature naturelle (en K) et pressions milieux couches .
-c   ---------------------------------------------------------------
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-        do ig0=1,klon
-          i=index_i(ig0)
-          j=index_j(ig0)
-          pksurcp        = ppk(i,j,l) / cpp
-          zplay(ig0,l)   = preff * pksurcp ** unskap
-          ztfi(ig0,l)    = pteta(i,j,l)  * pksurcp
-        enddo
-
-      ENDDO
-c$OMP END DO NOWAIT
-
-c   43.bis traceurs
-c   ---------------
-c
-
-      DO iq=1,nqtot
-         iiq=niadv(iq)
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-         DO l=1,llm
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-           do ig0=1,klon
-             i=index_i(ig0)
-             j=index_j(ig0)
-             zqfi(ig0,l,iq)  = pq(i,j,l,iiq)
-           enddo
-         ENDDO
-c$OMP END DO NOWAIT	 
-      ENDDO
-
-
-c   Geopotentiel calcule par rapport a la surface locale:
-c   -----------------------------------------------------
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-         DO l=1,llm
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-           do ig0=1,klon
-             i=index_i(ig0)
-             j=index_j(ig0)
-             zphi(ig0,l)  = pphi(i,j,l)
-           enddo
-         ENDDO
-c$OMP END DO NOWAIT	 
-
-c      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,pphi,zphi)
-
-c$OMP MASTER
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-           do ig0=1,klon
-             i=index_i(ig0)
-             j=index_j(ig0)
-             zphis(ig0)  = pphis(i,j)
-           enddo
-c$OMP END MASTER
-
-
-c      CALL gr_dyn_fi_p(1,iip1,jjp1,klon,pphis,zphis)
-
-c$OMP BARRIER
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-	 DO ig=1,klon
-	   zphi(ig,l)=zphi(ig,l)-zphis(ig)
-	 ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
-      
-
-c
-c   45. champ u:
-c   ------------
-
-      kstart=1
-      kend=klon
-      
-      if (is_north_pole) kstart=2
-      if (is_south_pole) kend=klon-1
-      
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-      DO l=1,llm
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-!CDIR SPARSE
-        do ig0=kstart,kend
-          i=index_i(ig0)
-          j=index_j(ig0)
-          if (i==1) then
-            zufi(ig0,l)= 0.5 *(  pucov(iim,j,l)/cu(iim,j)
-     $                         + pucov(1,j,l)/cu(1,j) )
-          else
-            zufi(ig0,l)= 0.5*(  pucov(i-1,j,l)/cu(i-1,j) 
-     $                       + pucov(i,j,l)/cu(i,j) )
-          endif
-        enddo
-      ENDDO
-c$OMP END DO NOWAIT
-
-c   46.champ v:
-c   -----------
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-        DO ig0=kstart,kend
-          i=index_i(ig0)
-          j=index_j(ig0)
-          zvfi(ig0,l)= 0.5 *(  pvcov(i,j-1,l)/cv(i,j-1) 
-     $                       + pvcov(i,j,l)/cv(i,j) )
-    
-         ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
-
-c   47. champs de vents aux pole nord   
-c   ------------------------------
-c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
-c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
-
-      if (is_north_pole) then
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-        DO l=1,llm
-
-           z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1)
-           DO i=2,iim
-              z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
-           ENDDO
-  
-           DO i=1,iim
-              zcos(i)   = COS(rlonv(i))*z1(i)
-              zsin(i)   = SIN(rlonv(i))*z1(i)
-           ENDDO
-  
-           zufi(1,l)  = SSUM(iim,zcos,1)/pi
-           zvfi(1,l)  = SSUM(iim,zsin,1)/pi
-  
-        ENDDO
-c$OMP END DO NOWAIT      
-      endif
-
-
-c   48. champs de vents aux pole sud:
-c   ---------------------------------
-c        U = 1 / pi  *  integrale [ v * cos(long) * d long ]
-c        V = 1 / pi  *  integrale [ v * sin(long) * d long ]
-
-      if (is_south_pole) then
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-        DO l=1,llm
-  
-         z1(1)   =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm)
-           DO i=2,iim
-             z1(i)   =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
-	   ENDDO
-  
-           DO i=1,iim
-              zcos(i)    = COS(rlonv(i))*z1(i)
-              zsin(i)    = SIN(rlonv(i))*z1(i)
-	   ENDDO
-  
-           zufi(klon,l)  = SSUM(iim,zcos,1)/pi
-           zvfi(klon,l)  = SSUM(iim,zsin,1)/pi
-        ENDDO
-c$OMP END DO NOWAIT       
-      endif
-
-c On change de grille, dynamique vers physiq, pour le flux de masse verticale
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-         DO l=1,llm
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-           do ig0=1,klon
-             i=index_i(ig0)
-             j=index_j(ig0)
-             flxwfi(ig0,l)  = flxw(i,j,l)
-           enddo
-         ENDDO
-c$OMP END DO NOWAIT
-
-c      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
-
-c-----------------------------------------------------------------------
-c   Appel de la physique:
-c   ---------------------
-
-
-c$OMP BARRIER
-      if (first_omp) then
-        klon=klon_omp
-
-        allocate(zplev_omp(klon,llm+1))
-        allocate(zplay_omp(klon,llm))
-        allocate(zphi_omp(klon,llm))
-        allocate(zphis_omp(klon))
-        allocate(presnivs_omp(llm))
-        allocate(zufi_omp(klon,llm))
-        allocate(zvfi_omp(klon,llm))
-        allocate(ztfi_omp(klon,llm))
-        allocate(zqfi_omp(klon,llm,nqtot))
-        allocate(zdufi_omp(klon,llm))
-        allocate(zdvfi_omp(klon,llm))
-        allocate(zdtfi_omp(klon,llm))
-        allocate(zdqfi_omp(klon,llm,nqtot))
-        allocate(zdufic_omp(klon,llm))
-        allocate(zdvfic_omp(klon,llm))
-        allocate(zdtfic_omp(klon,llm))
-        allocate(zdqfic_omp(klon,llm,nqtot))
-        allocate(zdpsrf_omp(klon))
-        allocate(flxwfi_omp(klon,llm))
-	first_omp=.false.
-      endif
-       
-	   
-      klon=klon_omp
-      offset=klon_omp_begin-1
-      
-      do l=1,llm+1
-        do i=1,klon
-          zplev_omp(i,l)=zplev(offset+i,l)
-	enddo 
-      enddo
-	  
-       do l=1,llm
-        do i=1,klon  
-	  zplay_omp(i,l)=zplay(offset+i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  zphi_omp(i,l)=zphi(offset+i,l)
-	enddo 
-      enddo
-	
-      do i=1,klon
-	zphis_omp(i)=zphis(offset+i)
-      enddo 
-     
-	
-      do l=1,llm
-        presnivs_omp(l)=presnivs(l)
-      enddo 
-	
-      do l=1,llm
-        do i=1,klon
-	  zufi_omp(i,l)=zufi(offset+i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  zvfi_omp(i,l)=zvfi(offset+i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  ztfi_omp(i,l)=ztfi(offset+i,l)
-	enddo 
-      enddo
-	
-      do iq=1,nqtot
-        do l=1,llm
-          do i=1,klon
-            zqfi_omp(i,l,iq)=zqfi(offset+i,l,iq)
-	  enddo
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  zdufi_omp(i,l)=zdufi(offset+i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  zdvfi_omp(i,l)=zdvfi(offset+i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-          zdtfi_omp(i,l)=zdtfi(offset+i,l)
-	enddo 
-      enddo
-	
-      do iq=1,nqtot
-        do l=1,llm
-          do i=1,klon
-	    zdqfi_omp(i,l,iq)=zdqfi(offset+i,l,iq)
-	  enddo 
-        enddo
-      enddo
-      	
-      do i=1,klon
-	zdpsrf_omp(i)=zdpsrf(offset+i)
-      enddo 
-
-      do l=1,llm
-        do i=1,klon
-          flxwfi_omp(i,l)=flxwfi(offset+i,l)
-	enddo 
-      enddo
-      
-c$OMP BARRIER
-      
-
-!$OMP MASTER 
-!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys 
-!$OMP END MASTER
-      zdt_split=dtphys/nsplit_phys
-      zdufic_omp(:,:)=0.
-      zdvfic_omp(:,:)=0.
-      zdtfic_omp(:,:)=0.
-      zdqfic_omp(:,:,:)=0.
-
-#ifdef CPP_PHYS
-      do isplit=1,nsplit_phys
-
-         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
-         debut_split=debut.and.isplit==1
-         lafin_split=lafin.and.isplit==nsplit_phys
-
-      if (planet_type=="earth") then
-
-      CALL physiq (klon,
-     .             llm,
-     .             debut_split,
-     .             lafin_split,
-     .             jD_cur,
-     .             jH_cur_split,
-     .             zdt_split,
-     .             zplev_omp,
-     .             zplay_omp,
-     .             zphi_omp,
-     .             zphis_omp,
-     .             presnivs_omp,
-     .             zufi_omp,
-     .             zvfi_omp,
-     .             ztfi_omp,
-     .             zqfi_omp,
-     .             flxwfi_omp,
-     .             zdufi_omp,
-     .             zdvfi_omp,
-     .             zdtfi_omp,
-     .             zdqfi_omp,
-     .             zdpsrf_omp,
-     .             pducov)
-
-      else if ( planet_type=="generic" ) then
-
-      CALL physiq (klon,     !! ngrid
-     .             llm,            !! nlayer
-     .             nqtot,          !! nq
-     .             tname,          !! tracer names from dynamical core (given in infotrac)
-     .             debut_split,    !! firstcall 
-     .             lafin_split,    !! lastcall
-     .             jD_cur,         !! pday. see leapfrog_p 
-     .             jH_cur_split,   !! ptime "fraction of day"
-     .             zdt_split,      !! ptimestep
-     .             zplev_omp,  !! pplev
-     .             zplay_omp,  !! pplay
-     .             zphi_omp,   !! pphi
-     .             zufi_omp,   !! pu
-     .             zvfi_omp,   !! pv
-     .             ztfi_omp,   !! pt
-     .             zqfi_omp,   !! pq
-     .             flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
-     .             zdufi_omp,  !! pdu
-     .             zdvfi_omp,  !! pdv
-     .             zdtfi_omp,  !! pdt
-     .             zdqfi_omp,  !! pdq
-     .             zdpsrf_omp, !! pdpsrf
-     .             tracerdyn)      !! tracerdyn <-- utilite ???
-
-      endif ! of if (planet_type=="earth")
-
-
-         zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split
-         zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split
-         ztfi_omp(:,:)=ztfi_omp(:,:)+zdtfi_omp(:,:)*zdt_split
-         zqfi_omp(:,:,:)=zqfi_omp(:,:,:)+zdqfi_omp(:,:,:)*zdt_split
-
-         zdufic_omp(:,:)=zdufic_omp(:,:)+zdufi_omp(:,:)
-         zdvfic_omp(:,:)=zdvfic_omp(:,:)+zdvfi_omp(:,:)
-         zdtfic_omp(:,:)=zdtfic_omp(:,:)+zdtfi_omp(:,:)
-         zdqfic_omp(:,:,:)=zdqfic_omp(:,:,:)+zdqfi_omp(:,:,:)
-
-      enddo
-
-#endif
-! of #ifdef CPP_PHYS
-
-
-      zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
-      zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
-      zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
-      zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
-
-c$OMP BARRIER
-
-      do l=1,llm+1
-        do i=1,klon
-          zplev(offset+i,l)=zplev_omp(i,l)
-	enddo 
-      enddo
-	  
-       do l=1,llm
-        do i=1,klon  
-	  zplay(offset+i,l)=zplay_omp(i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  zphi(offset+i,l)=zphi_omp(i,l)
-	enddo 
-      enddo
-	
-
-      do i=1,klon
-	zphis(offset+i)=zphis_omp(i)
-      enddo 
-     
-	
-      do l=1,llm
-        presnivs(l)=presnivs_omp(l)
-      enddo 
-	
-      do l=1,llm
-        do i=1,klon
-	  zufi(offset+i,l)=zufi_omp(i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  zvfi(offset+i,l)=zvfi_omp(i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  ztfi(offset+i,l)=ztfi_omp(i,l)
-	enddo 
-      enddo
-	
-      do iq=1,nqtot
-        do l=1,llm
-          do i=1,klon
-            zqfi(offset+i,l,iq)=zqfi_omp(i,l,iq)
-	  enddo
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  zdufi(offset+i,l)=zdufi_omp(i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-	  zdvfi(offset+i,l)=zdvfi_omp(i,l)
-	enddo 
-      enddo
-	
-      do l=1,llm
-        do i=1,klon
-          zdtfi(offset+i,l)=zdtfi_omp(i,l)
-	enddo 
-      enddo
-	
-      do iq=1,nqtot
-        do l=1,llm
-          do i=1,klon
-	    zdqfi(offset+i,l,iq)=zdqfi_omp(i,l,iq)
-	  enddo 
-        enddo
-      enddo
-      	
-      do i=1,klon
-	zdpsrf(offset+i)=zdpsrf_omp(i)
-      enddo 
-      
-
-      klon=klon_mpi
-500   CONTINUE
-c$OMP BARRIER
-
-c$OMP MASTER
-      call stop_timer(timer_physic)
-c$OMP END MASTER
-
-      IF (using_mpi) THEN
-            
-      if (MPI_rank>0) then
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
-       DO l=1,llm      
-        du_send(1:iim,l)=zdufi(1:iim,l)
-        dv_send(1:iim,l)=zdvfi(1:iim,l)
-       ENDDO
-c$OMP END DO NOWAIT       
-
-c$OMP BARRIER
-#ifdef CPP_MPI 
-c$OMP MASTER
-!$OMP CRITICAL (MPI)
-        call MPI_ISSEND(du_send,iim*llm,MPI_REAL8,MPI_Rank-1,401,
-     &                   COMM_LMDZ,Req(1),ierr)
-        call MPI_ISSEND(dv_send,iim*llm,MPI_REAL8,MPI_Rank-1,402,
-     &                  COMM_LMDZ,Req(2),ierr)
-!$OMP END CRITICAL (MPI)
-c$OMP END MASTER
-#endif
-c$OMP BARRIER
-     
-      endif
-   
-      if (MPI_rank<MPI_Size-1) then
-c$OMP BARRIER
-#ifdef CPP_MPI 
-c$OMP MASTER      
-!$OMP CRITICAL (MPI)
-        call MPI_IRECV(du_recv,iim*llm,MPI_REAL8,MPI_Rank+1,401,
-     &                 COMM_LMDZ,Req(3),ierr)
-        call MPI_IRECV(dv_recv,iim*llm,MPI_REAL8,MPI_Rank+1,402,
-     &                 COMM_LMDZ,Req(4),ierr)
-!$OMP END CRITICAL (MPI)
-c$OMP END MASTER
-#endif
-      endif
-
-c$OMP BARRIER
-
-
-#ifdef CPP_MPI 
-c$OMP MASTER    
-!$OMP CRITICAL (MPI)
-      if (MPI_rank>0 .and. MPI_rank< MPI_Size-1) then
-        call MPI_WAITALL(4,Req(1),Status,ierr)
-      else if (MPI_rank>0) then
-        call MPI_WAITALL(2,Req(1),Status,ierr)
-      else if (MPI_rank <MPI_Size-1) then
-        call MPI_WAITALL(2,Req(3),Status,ierr)
-      endif
-!$OMP END CRITICAL (MPI)
-c$OMP END MASTER
-#endif
-
-c$OMP BARRIER     
-
-      ENDIF ! using_mpi
-      
-      
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-            
-        zdufi2(1:klon,l)=zdufi(1:klon,l)
-        zdufi2(klon+1:klon+iim,l)=du_recv(1:iim,l)
-            
-        zdvfi2(1:klon,l)=zdvfi(1:klon,l)
-        zdvfi2(klon+1:klon+iim,l)=dv_recv(1:iim,l) 
-
-        pdhfi(:,jj_begin,l)=0
-        pdqfi(:,jj_begin,l,:)=0
-        pdufi(:,jj_begin,l)=0
-        pdvfi(:,jj_begin,l)=0
-                
-        if (.not. is_south_pole) then
-          pdhfi(:,jj_end:jj_end+1,l)=0
-          pdqfi(:,jj_end:jj_end+1,l,:)=0
-          pdufi(:,jj_end:jj_end+1,l)=0
-          pdvfi(:,jj_end:jj_end+1,l)=0
-        endif
-      
-       ENDDO 
-c$OMP END DO NOWAIT
-
-c$OMP MASTER
-        pdpsfi(:,jj_begin)=0    
-       
-       if (.not. is_south_pole) then
-	 pdpsfi(:,jj_end:jj_end+1)=0
-       endif
-c$OMP END MASTER
-c-----------------------------------------------------------------------
-c   transformation des tendances physiques en tendances dynamiques:
-c   ---------------------------------------------------------------
-
-c  tendance sur la pression :
-c  -----------------------------------
-c      CALL gr_fi_dyn_p(1,klon,iip1,jjp1,zdpsrf,pdpsfi)
-
-c$OMP MASTER
-      kstart=1
-      kend=klon
-
-      if (is_north_pole) kstart=2
-      if (is_south_pole)  kend=klon-1
-
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-!cdir NODEP
-        do ig0=kstart,kend
-          i=index_i(ig0)
-          j=index_j(ig0)
-          pdpsfi(i,j) = zdpsrf(ig0)
-          if (i==1) pdpsfi(iip1,j) =  zdpsrf(ig0)
-         enddo          
-
-        if (is_north_pole) then
-            DO i=1,iip1
-              pdpsfi(i,1)    = zdpsrf(1)
-            enddo
-        endif
-        
-        if (is_south_pole) then
-            DO i=1,iip1
-              pdpsfi(i,jjp1) = zdpsrf(klon)
-            ENDDO
-        endif
-c$OMP END MASTER
-cc$OMP BARRIER
-
-c
-c   62. enthalpie potentielle
-c   ---------------------
-      
-      kstart=1
-      kend=klon
-
-      if (is_north_pole) kstart=2
-      if (is_south_pole)  kend=klon-1
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-!cdir NODEP
-        do ig0=kstart,kend
-          i=index_i(ig0)
-          j=index_j(ig0)
-          pdhfi(i,j,l) = cpp * zdtfi(ig0,l) / ppk(i,j,l)
-          if (i==1) pdhfi(iip1,j,l) =  cpp * zdtfi(ig0,l) / ppk(i,j,l)
-         enddo          
-
-        if (is_north_pole) then
-            DO i=1,iip1
-              pdhfi(i,1,l)    = cpp *  zdtfi(1,l)      / ppk(i, 1  ,l)
-            enddo
-        endif
-        
-        if (is_south_pole) then
-            DO i=1,iip1
-              pdhfi(i,jjp1,l) = cpp *  zdtfi(klon,l)/ ppk(i,jjp1,l)
-            ENDDO
-        endif
-      ENDDO
-c$OMP END DO NOWAIT
-      
-c   62. humidite specifique
-c   ---------------------
-! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
-!      DO iq=1,nqtot
-!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-!         DO l=1,llm
-!!!cdir NODEP 
-!           do ig0=kstart,kend
-!             i=index_i(ig0)
-!             j=index_j(ig0)
-!             pdqfi(i,j,l,iq) = zdqfi(ig0,l,iq) 
-!             if (i==1) pdqfi(iip1,j,l,iq) = zdqfi(ig0,l,iq) 
-!           enddo
-!           
-!           if (is_north_pole) then
-!             do i=1,iip1
-!               pdqfi(i,1,l,iq)    = zdqfi(1,l,iq)             
-!             enddo
-!           endif
-!           
-!           if (is_south_pole) then
-!             do i=1,iip1
-!               pdqfi(i,jjp1,l,iq) = zdqfi(klon,l,iq) 
-!             enddo
-!           endif
-!         ENDDO
-!c$OMP END DO NOWAIT
-!      ENDDO
-
-c   63. traceurs
-c   ------------
-C     initialisation des tendances
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-        pdqfi(:,jj_begin:jj_end,l,:)=0.
-      ENDDO
-c$OMP END DO NOWAIT	 
-
-C
-!cdir NODEP
-      DO iq=1,nqtot
-         iiq=niadv(iq)
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-         DO l=1,llm
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-!cdir NODEP           
-	     DO ig0=kstart,kend
-              i=index_i(ig0)
-              j=index_j(ig0)
-              pdqfi(i,j,l,iiq) = zdqfi(ig0,l,iq)
-              if (i==1) pdqfi(iip1,j,l,iiq) = zdqfi(ig0,l,iq)
-            ENDDO
-	    
-	    IF (is_north_pole) then
-	      DO i=1,iip1
-                pdqfi(i,1,l,iiq)    = zdqfi(1,l,iq)
-	      ENDDO
-	    ENDIF
-	    
-	    IF (is_south_pole) then
-	      DO i=1,iip1
-                pdqfi(i,jjp1,l,iiq) = zdqfi(klon,l,iq)
-	      ENDDO
-	    ENDIF
-	    
-         ENDDO
-c$OMP END DO NOWAIT	 
-      ENDDO
-      
-c   65. champ u:
-c   ------------
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,llm
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-!cdir NODEP
-         do ig0=kstart,kend
-           i=index_i(ig0)
-           j=index_j(ig0)
-           
-           if (i/=iim) then
-             pdufi(i,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
-           endif
-           
-           if (i==1) then
-              pdufi(iim,j,l)=0.5*(  zdufi2(ig0,l)
-     $                            + zdufi2(ig0+iim-1,l))*cu(iim,j)
-             pdufi(iip1,j,l)=0.5*(zdufi2(ig0,l)+zdufi2(ig0+1,l))*cu(i,j)
-           endif
-         
-         enddo
-         
-         if (is_north_pole) then
-           DO i=1,iip1
-            pdufi(i,1,l)    = 0.
-           ENDDO
-         endif
-         
-         if (is_south_pole) then
-           DO i=1,iip1
-            pdufi(i,jjp1,l) = 0.
-           ENDDO
-         endif
-         
-      ENDDO
-c$OMP END DO NOWAIT
-
-c   67. champ v:
-c   ------------
-
-      kstart=1
-      kend=klon
-
-      if (is_north_pole) kstart=2
-      if (is_south_pole)  kend=klon-1-iim
-      
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-      DO l=1,llm
-!CDIR ON_ADB(index_i)
-!CDIR ON_ADB(index_j) 
-!cdir NODEP
-        do ig0=kstart,kend
-           i=index_i(ig0)
-           j=index_j(ig0)
-           pdvfi(i,j,l)=0.5*(zdvfi2(ig0,l)+zdvfi2(ig0+iim,l))*cv(i,j)
-           if (i==1) pdvfi(iip1,j,l) = 0.5*(zdvfi2(ig0,l)+
-     $	                                    zdvfi2(ig0+iim,l))
-     $				          *cv(i,j)
-        enddo
-         
-      ENDDO
-c$OMP END DO NOWAIT
-
-
-c   68. champ v pres des poles:
-c   ---------------------------
-c      v = U * cos(long) + V * SIN(long)
-
-      if (is_north_pole) then
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-        DO l=1,llm
-
-          DO i=1,iim
-            pdvfi(i,1,l)=
-     $      zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i))
-       
-            pdvfi(i,1,l)=
-     $      0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1)
-          ENDDO
-
-          pdvfi(iip1,1,l)  = pdvfi(1,1,l)
-
-        ENDDO
-c$OMP END DO NOWAIT
-
-      endif    
-      
-      if (is_south_pole) then
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
-         DO l=1,llm
-  
-           DO i=1,iim
-              pdvfi(i,jjm,l)=zdufi(klon,l)*COS(rlonv(i))
-     $        +zdvfi(klon,l)*SIN(rlonv(i))
-
-              pdvfi(i,jjm,l)=
-     $        0.5*(pdvfi(i,jjm,l)+zdvfi(klon-iip1+i,l))*cv(i,jjm)
-           ENDDO
-
-           pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
-
-        ENDDO
-c$OMP END DO NOWAIT
-     
-      endif
-c-----------------------------------------------------------------------
-
-700   CONTINUE
- 
-      firstcal = .FALSE.
-
-#else 
-      write(lunout,*) 
-     & "calfis_p: for now can only work with parallel physics" 
-      stop 
-#endif 
-! of #ifdef CPP_PHYS
-      RETURN
-      END
Index: LMDZ5/branches/testing/libf/dyn3dmem/call_calfis_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/call_calfis_mod.F90	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/call_calfis_mod.F90	(revision 2298)
@@ -227,9 +227,10 @@
   !$OMP BARRIER
 
+#ifdef CPP_PHYS
     CALL calfis_loc(lafin ,jD_cur, jH_cur,                       &
                      ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,   &
                      du,dv,dteta,dq,                             &
                      flxw, dufi,dvfi,dtetafi,dqfi,dpfi  )
-
+#endif
     ijb=ij_begin
     ije=ij_end  
Index: LMDZ5/branches/testing/libf/dyn3dmem/check_isotopes_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/check_isotopes_loc.F	(revision 2298)
+++ LMDZ5/branches/testing/libf/dyn3dmem/check_isotopes_loc.F	(revision 2298)
@@ -0,0 +1,215 @@
+        subroutine check_isotopes(q,ijb,ije,err_msg)
+        USE infotrac
+        USE parallel_lmdz
+        implicit none
+
+#include "dimensions.h"
+
+        ! inputs
+        integer ijb,ije ! peut être local et différent de ijb_u,ije_u, ex: dans qminimum
+        real q(ijb_u:ije_u,llm,nqtot)
+        character*(*) err_msg ! message d''erreur à afficher
+
+        ! locals
+        integer ixt,phase,k,i,iq,iiso,izone,ieau,iqeau
+        real xtractot,xiiso
+        real borne
+        real qmin
+        real errmax ! erreur maximale en absolu.
+        real errmaxrel ! erreur maximale en relatif autorisée
+        real deltaDmax,deltaDmin
+        real ridicule
+        parameter (borne=1e19)
+        parameter (errmax=1e-8)
+        parameter (errmaxrel=1e-3)
+        parameter (qmin=1e-11)
+        parameter (deltaDmax=200.0,deltaDmin=-999.9)
+        parameter (ridicule=1e-12)
+        real deltaD
+
+        if (ok_isotopes) then
+
+        !write(*,*) 'check_isotopes 31: err_msg=',err_msg
+        ! verifier que rien n'est NaN
+        do ixt=1,ntraciso
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+            do k=1,llm
+              DO i = ijb,ije
+                if ((q(i,k,iq).gt.-borne).and.
+     :            (q(i,k,iq).lt.borne)) then
+                else !if ((x(ixt,i,j).gt.-borne).and.
+                  write(*,*) 'erreur detectee par iso_verif_noNaN:'
+                  write(*,*) err_msg
+                  write(*,*) 'q,i,k,iq=',q(i,k,iq),i,k,iq
+                  write(*,*) 'borne=',borne
+                  stop
+                endif  !if ((x(ixt,i,j).gt.-borne).and.
+              enddo !DO i = ijb,ije
+            enddo !do k=1,llm
+c$OMP END DO NOWAIT
+          enddo !do phase=1,nqo
+        enddo !do ixt=1,ntraciso
+
+        !write(*,*) 'check_isotopes 52'
+        ! verifier que l'eau normale est OK
+        if (use_iso(1)) then
+          ixt=indnum_fn_num(1)
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+            do k=1,llm
+            DO i = ijb,ije  
+              if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
+     :          (abs((q(i,k,phase)-q(i,k,iq))/
+     :           max(max(abs(q(i,k,phase)),abs(q(i,k,iq))),1e-18))
+     :           .gt.errmaxrel)) then
+                  write(*,*) 'erreur detectee par iso_verif_egalite:'
+                  write(*,*) err_msg
+                  write(*,*) 'ixt,phase,ijb=',ixt,phase,ijb
+                  write(*,*) 'q,iq,i,k=',q(i,k,iq),iq,i,k
+                  write(*,*) 'q(i,k,phase)=',q(i,k,phase)
+                  stop
+              endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
+              ! bidouille pour éviter divergence:
+              q(i,k,iq)= q(i,k,phase) 
+            enddo ! DO i = ijb,ije
+            enddo !do k=1,llm
+c$OMP END DO NOWAIT
+          enddo ! do phase=1,nqo 
+        endif !if (use_iso(1)) then
+        
+        !write(*,*) 'check_isotopes 78'
+        ! verifier que HDO est raisonable
+        if (use_iso(2)) then
+          ixt=indnum_fn_num(2)
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+            do k=1,llm
+            DO i = ijb,ije
+            if (q(i,k,iq).gt.qmin) then
+             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(2)-1)*1000
+             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur detectee par iso_verif_aberrant:'
+                  write(*,*) err_msg
+                  write(*,*) 'ixt,phase=',ixt,phase
+                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
+                  write(*,*) 'q=',q(i,k,:)
+                  write(*,*) 'deltaD=',deltaD
+                  stop
+             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+            endif !if (q(i,k,iq).gt.qmin) then
+            enddo !DO i = ijb,ije
+            enddo !do k=1,llm
+c$OMP END DO NOWAIT
+          enddo ! do phase=1,nqo 
+        endif !if (use_iso(2)) then
+
+        !write(*,*) 'check_isotopes 103'
+        ! verifier que O18 est raisonable
+        if (use_iso(3)) then
+          ixt=indnum_fn_num(3)
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+            do k=1,llm
+            DO i = ijb,ije
+            if (q(i,k,iq).gt.qmin) then
+             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(3)-1)*1000
+             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur detectee iso_verif_aberrant O18:'
+                  write(*,*) err_msg
+                  write(*,*) 'ixt,phase=',ixt,phase
+                  write(*,*) 'q,iq,i,k,=',q(i,k,phase),iq,i,k
+                  write(*,*) 'xt=',q(i,k,:)
+                  write(*,*) 'deltaO18=',deltaD
+                  stop
+             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+            endif !if (q(i,k,iq).gt.qmin) then
+            enddo !DO i = ijb,ije
+            enddo !do k=1,llm
+c$OMP END DO NOWAIT
+          enddo ! do phase=1,nqo 
+        endif !if (use_iso(2)) then
+
+
+        !write(*,*) 'check_isotopes 129'
+        if (ok_isotrac) then
+
+          if (use_iso(2).and.use_iso(1)) then
+            do izone=1,ntraceurs_zone
+             ixt=index_trac(izone,indnum_fn_num(2))
+             ieau=index_trac(izone,indnum_fn_num(1))
+             do phase=1,nqo
+               iq=iqiso(ixt,phase)
+               iqeau=iqiso(ieau,phase)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+               do k=1,llm
+                DO i = ijb,ije
+                if (q(i,k,iq).gt.qmin) then
+                 deltaD=(q(i,k,iq)/q(i,k,iqeau)/tnat(2)-1)*1000
+                 if ((deltaD.gt.deltaDmax).or.
+     &                   (deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur dans iso_verif_aberrant trac:'
+                  write(*,*) err_msg
+                  write(*,*) 'izone,phase=',izone,phase
+                  write(*,*) 'ixt,ieau=',ixt,ieau
+                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
+                  write(*,*) 'deltaD=',deltaD
+                  stop
+                 endif !if ((deltaD.gt.deltaDmax).or.
+                endif !if (q(i,k,iq).gt.qmin) then
+                enddo !DO i = ijb,ije
+                enddo  ! do k=1,llm
+c$OMP END DO NOWAIT
+              enddo ! do phase=1,nqo    
+            enddo !do izone=1,ntraceurs_zone
+          endif !if (use_iso(2).and.use_iso(1)) then
+
+          do iiso=1,niso
+           do phase=1,nqo
+              iq=iqiso(iiso,phase)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+              do k=1,llm
+                DO i = ijb,ije
+                   xtractot=0.0
+                   xiiso=q(i,k,iq)
+                   do izone=1,ntraceurs_zone
+                      iq=iqiso(index_trac(izone,iiso),phase)
+                      xtractot=xtractot+ q(i,k,iq)
+                   enddo !do izone=1,ntraceurs_zone
+                   if ((abs(xtractot-xiiso).gt.errmax).and.
+     :                  (abs(xtractot-xiiso)/
+     :                  max(max(abs(xtractot),abs(xiiso)),1e-18)
+     :                  .gt.errmaxrel)) then
+                  write(*,*) 'erreur detectee par iso_verif_traceurs:'
+                  write(*,*) err_msg
+                  write(*,*) 'iiso,phase=',iiso,phase
+                  write(*,*) 'i,k,=',i,k
+                  write(*,*) 'q(i,k,:)=',q(i,k,:)
+                  stop
+                 endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
+                  
+                 ! bidouille pour éviter divergence:
+                 if (abs(xtractot).gt.ridicule) then
+                   do izone=1,ntraceurs_zone
+                     ixt=index_trac(izone,iiso) 
+                     q(i,k,iq)=q(i,k,iq)/xtractot*xiiso
+                   enddo !do izone=1,ntraceurs_zone                
+                  endif !if ((abs(xtractot).gt.ridicule) then
+                enddo !DO i = ijb,ije
+              enddo !do k=1,llm
+c$OMP END DO NOWAIT
+           enddo !do phase=1,nqo
+          enddo !do iiso=1,niso
+
+        endif !if (ok_isotrac) then
+
+        endif ! if (ok_isotopes)
+        !write(*,*) 'check_isotopes 198'
+        
+        end
+
+
Index: LMDZ5/branches/testing/libf/dyn3dmem/dynetat0_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/dynetat0_loc.F	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/dynetat0_loc.F	(revision 2298)
@@ -366,4 +366,19 @@
            write(lunout,*)"Il est donc initialise a zero"
            q(:,:,iq)=0.
+
+           ! CRisi: pour les isotopes, on peut faire init théorique
+           ! distill de Rayleigh très simplifiée
+           if (ok_isotopes) then
+              if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.0)) then
+                q(:,:,iq)=q(:,:,iqpere(iq))                             &
+     &                   *tnat(iso_num(iq))                             &
+     &                   *(q(:,:,iqpere(iq))/30.e-3)                    &
+     &                   **(alpha_ideal(iso_num(iq))-1)
+              endif
+              if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.1)) then
+                  q(:,:,iq)=q(:,:,iqiso(iso_indnum(iq),phase_num(iq)))
+              endif  
+           endif !if (ok_isotopes) then        
+
         ELSE
 #ifdef NC_DOUBLE
@@ -380,5 +395,9 @@
 
         ENDIF
-      ENDDO
+      ENDDO !DO iq=1,nqtot
+
+      if (ok_iso_verif) then
+         call check_isotopes(q,ijb_u,ije_u,'dynetat0_loc')
+      endif !if (ok_iso_verif) then
 
       DEALLOCATE(q_glo)
Index: LMDZ5/branches/testing/libf/dyn3dmem/gcm.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/gcm.F	(revision 2258)
+++ 	(revision )
@@ -1,504 +1,0 @@
-!
-! $Id: $
-!
-c
-c
-      PROGRAM gcm
-
-#ifdef CPP_IOIPSL
-      USE IOIPSL
-#endif
-
-      USE mod_const_mpi, ONLY: init_const_mpi
-      USE parallel_lmdz
-      USE infotrac
-      USE mod_interface_dyn_phys
-      USE mod_hallo
-      USE Bands
-      USE getparam
-      USE filtreg_mod
-      USE control_mod
-
-#ifdef INCA
-! Only INCA needs these informations (from the Earth's physics)
-      USE indice_sol_mod
-      USE mod_phys_lmdz_omp_data, ONLY: klon_omp 
-#endif
-
-#ifdef CPP_PHYS
-!      USE mod_grid_phy_lmdz
-!      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
-!      USE dimphy
-!      USE comgeomphy
-#endif
-      IMPLICIT NONE
-
-c      ......   Version  du 10/01/98    ..........
-
-c             avec  coordonnees  verticales hybrides 
-c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
-
-c=======================================================================
-c
-c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
-c   -------
-c
-c   Objet:
-c   ------
-c
-c   GCM LMD nouvelle grille
-c
-c=======================================================================
-c
-c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
-c      et possibilite d'appeler une fonction f(y)  a derivee tangente
-c      hyperbolique a la  place de la fonction a derivee sinusoidale.
-c  ... Possibilite de choisir le schema pour l'advection de
-c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
-c
-c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
-c      Pour Van-Leer iadv=10
-c
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comdissnew.h"
-#include "comvert.h"
-#include "comgeom.h"
-#include "logic.h"
-#include "temps.h"
-#include "ener.h"
-#include "description.h"
-#include "serre.h"
-!#include "com_io_dyn.h"
-#include "iniprint.h"
-#include "tracstoke.h"
-
-#ifdef INCA
-! Only INCA needs these informations (from the Earth's physics)
-!#include "indicesol.h"
-#endif
-
-      REAL zdtvr
-
-c   variables dynamiques
-      REAL,ALLOCATABLE,SAVE  :: vcov(:,:),ucov(:,:) ! vents covariants
-      REAL,ALLOCATABLE,SAVE  :: teta(:,:)     ! temperature potentielle 
-      REAL, ALLOCATABLE,SAVE :: q(:,:,:)      ! champs advectes
-      REAL,ALLOCATABLE,SAVE  :: ps(:)         ! pression  au sol
-c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-      REAL,ALLOCATABLE,SAVE  :: masse(:,:)    ! masse d'air
-      REAL,ALLOCATABLE,SAVE  :: phis(:)       ! geopotentiel au sol
-c      REAL phi(ip1jmp1,llm)                  ! geopotentiel
-c      REAL w(ip1jmp1,llm)                    ! vitesse verticale
-
-c variables dynamiques intermediaire pour le transport
-
-c   variables pour le fichier histoire
-      REAL dtav      ! intervalle de temps elementaire
-
-      REAL time_0
-
-      LOGICAL lafin
-
-      real time_step, t_wrt, t_ops
-
-c+jld variables test conservation energie
-c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
-C     Tendance de la temp. potentiel d (theta)/ d t due a la 
-C     tansformation d'energie cinetique en energie thermique
-C     cree par la dissipation
-c      REAL dhecdt(ip1jmp1,llm)
-c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
-c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
-c      CHARACTER (len=15) :: ztit
-c-jld 
-
-
-      character (len=80) :: dynhist_file, dynhistave_file
-      character (len=20) :: modname
-      character (len=80) :: abort_message
-! locales pour gestion du temps
-      INTEGER :: an, mois, jour
-      REAL :: heure
-
-
-c-----------------------------------------------------------------------
-c   Initialisations:
-c   ----------------
-
-      abort_message = 'last timestep reached'
-      modname = 'gcm'
-      descript = 'Run GCM LMDZ'
-      lafin    = .FALSE.
-      dynhist_file = 'dyn_hist'
-      dynhistave_file = 'dyn_hist_ave'
-
-
-
-c----------------------------------------------------------------------
-c  lecture des fichiers gcm.def ou run.def
-c  ---------------------------------------
-c
-      CALL conf_gcm( 99, .TRUE. )
-      if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm",
-     s "iphysiq must be a multiple of iperiod", 1)
-c
-c
-c------------------------------------
-c   Initialisation partie parallele
-c------------------------------------
-      CALL init_const_mpi
-
-      call init_parallel
-      call ini_getparam("out.def")
-      call Read_Distrib
-
-#ifdef CPP_PHYS
-        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
-!#endif
-!      CALL set_bands
-!#ifdef CPP_PHYS
-      CALL Init_interface_dyn_phys
-#endif
-      CALL barrier
-
-      CALL set_bands
-      if (mpi_rank==0) call WriteBands
-      call Set_Distrib(distrib_caldyn)
-
-c$OMP PARALLEL
-      call Init_Mod_hallo
-c$OMP END PARALLEL
-
-!#ifdef CPP_PHYS
-!c$OMP PARALLEL
-!      call InitComgeomphy ! now done in iniphysiq
-!c$OMP END PARALLEL 
-!#endif
-
-c-----------------------------------------------------------------------
-c   Choix du calendrier
-c   -------------------
-
-c      calend = 'earth_365d'
-
-#ifdef CPP_IOIPSL
-      if (calend == 'earth_360d') then
-        call ioconf_calendar('360d')
-        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
-      else if (calend == 'earth_365d') then
-        call ioconf_calendar('noleap')
-        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
-      else if (calend == 'gregorian') then
-        call ioconf_calendar('gregorian')
-        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
-      else
-        abort_message = 'Mauvais choix de calendrier'
-        call abort_gcm(modname,abort_message,1)
-      endif
-#endif
-
-      IF (type_trac == 'inca') THEN
-#ifdef INCA
-         call init_const_lmdz(
-     $        nbtr,anneeref,dayref,
-     $        iphysiq,day_step,nday, 
-     $        nbsrf, is_oce,is_sic,
-     $        is_ter,is_lic, calend)
-
-         call init_inca_para(
-     $        iim,jjm+1,llm,klon_glo,mpi_size,
-     $        distrib_phys,COMM_LMDZ)
-#endif
-      END IF
-
-c-----------------------------------------------------------------------
-c   Initialisation des traceurs
-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(ucov(ijb_u:ije_u,llm))
-      ALLOCATE(vcov(ijb_v:ije_v,llm))
-      ALLOCATE(teta(ijb_u:ije_u,llm))
-      ALLOCATE(masse(ijb_u:ije_u,llm))
-      ALLOCATE(ps(ijb_u:ije_u))
-      ALLOCATE(phis(ijb_u:ije_u))
-      ALLOCATE(q(ijb_u:ije_u,llm,nqtot))
-
-c-----------------------------------------------------------------------
-c   Lecture de l'etat initial :
-c   ---------------------------
-
-c  lecture du fichier start.nc
-      if (read_start) then
-      ! we still need to run iniacademic to initialize some
-      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
-        if (iflag_phys.ne.1) then
-          CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
-        endif
-
-!        if (planet_type.eq."earth") then
-! Load an Earth-format start file
-         CALL dynetat0_loc("start.nc",vcov,ucov,
-     &              teta,q,masse,ps,phis, time_0)
-!        endif ! of if (planet_type.eq."earth")
-
-c       write(73,*) 'ucov',ucov
-c       write(74,*) 'vcov',vcov
-c       write(75,*) 'teta',teta
-c       write(76,*) 'ps',ps
-c       write(77,*) 'q',q
-
-      endif ! of if (read_start)
-
-c le cas echeant, creation d un etat initial
-      IF (prt_level > 9) WRITE(lunout,*)
-     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
-      if (.not.read_start) then
-         CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
-      endif
-
-c-----------------------------------------------------------------------
-c   Lecture des parametres de controle pour la simulation :
-c   -------------------------------------------------------
-c  on recalcule eventuellement le pas de temps
-
-      IF(MOD(day_step,iperiod).NE.0) THEN
-        abort_message = 
-     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
-        call abort_gcm(modname,abort_message,1)
-      ENDIF
-
-      IF(MOD(day_step,iphysiq).NE.0) THEN
-        abort_message = 
-     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
-        call abort_gcm(modname,abort_message,1)
-      ENDIF
-
-      zdtvr    = daysec/REAL(day_step)
-        IF(dtvr.NE.zdtvr) THEN
-         WRITE(lunout,*)
-     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
-        ENDIF
-
-C
-C on remet le calendrier \`a zero si demande
-c
-      IF (start_time /= starttime) then
-        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
-     &,' fichier restart ne correspond pas a celle lue dans le run.def'
-        IF (raz_date == 1) then
-          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
-          start_time = starttime
-        ELSE
-          WRITE(lunout,*)'Je m''arrete'
-          CALL abort
-        ENDIF
-      ENDIF
-      IF (raz_date == 1) THEN
-        annee_ref = anneeref
-        day_ref = dayref
-        day_ini = dayref
-        itau_dyn = 0
-        itau_phy = 0
-        time_0 = 0.
-        write(lunout,*)
-     .   'GCM: On reinitialise a la date lue dans gcm.def'
-      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
-        write(lunout,*)
-     .  'GCM: Attention les dates initiales lues dans le fichier'
-        write(lunout,*)
-     .  ' restart ne correspondent pas a celles lues dans '
-        write(lunout,*)' gcm.def'
-        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
-        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
-        write(lunout,*)' Pas de remise a zero'
-      ENDIF
-c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
-c        write(lunout,*)
-c     .  'GCM: Attention les dates initiales lues dans le fichier'
-c        write(lunout,*)
-c     .  ' restart ne correspondent pas a celles lues dans '
-c        write(lunout,*)' gcm.def'
-c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
-c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
-c        if (raz_date .ne. 1) then
-c          write(lunout,*)
-c     .    'GCM: On garde les dates du fichier restart'
-c        else
-c          annee_ref = anneeref
-c          day_ref = dayref
-c          day_ini = dayref
-c          itau_dyn = 0
-c          itau_phy = 0
-c          time_0 = 0.
-c          write(lunout,*)
-c     .   'GCM: On reinitialise a la date lue dans gcm.def'
-c        endif
-c      ELSE
-c        raz_date = 0
-c      endif
-
-#ifdef CPP_IOIPSL
-      mois = 1
-      heure = 0.
-      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
-      jH_ref = jD_ref - int(jD_ref)
-      jD_ref = int(jD_ref)
-
-      call ioconf_startdate(INT(jD_ref), jH_ref)
-
-      write(lunout,*)'DEBUG'
-      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
-      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
-      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
-      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
-      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
-#else
-! Ehouarn: we still need to define JD_ref and JH_ref
-! and since we don't know how many days there are in a year
-! we set JD_ref to 0 (this should be improved ...)
-      jD_ref=0
-      jH_ref=0
-#endif
-
-      if (iflag_phys.eq.1) then
-      ! these initialisations have already been done (via iniacademic)
-      ! if running in SW or Newtonian mode
-c-----------------------------------------------------------------------
-c   Initialisation des constantes dynamiques :
-c   ------------------------------------------
-        dtvr = zdtvr
-        CALL iniconst
-
-c-----------------------------------------------------------------------
-c   Initialisation de la geometrie :
-c   --------------------------------
-        CALL inigeom
-
-c-----------------------------------------------------------------------
-c   Initialisation du filtre :
-c   --------------------------
-        CALL inifilr
-      endif ! of if (iflag_phys.eq.1)
-c
-c-----------------------------------------------------------------------
-c   Initialisation de la dissipation :
-c   ----------------------------------
-
-      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
-     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
-
-c-----------------------------------------------------------------------
-c   Initialisation de la physique :
-c   -------------------------------
-      IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
-! Physics:
-#ifdef CPP_PHYS
-         CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
-     &                rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,
-     &                iflag_phys)
-#endif
-      ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
-
-
-c-----------------------------------------------------------------------
-c   Initialisation des dimensions d'INCA :
-c   --------------------------------------
-      IF (type_trac == 'inca') 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   ------------------------
-
-
-      if (nday>=0) then
-         day_end = day_ini + nday
-      else
-         day_end = day_ini - nday/day_step
-      endif
-  
-      WRITE(lunout,300)day_ini,day_end
- 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
-
-#ifdef CPP_IOIPSL
-      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
-      write (lunout,301)jour, mois, an
-      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
-      write (lunout,302)jour, mois, an
- 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
- 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
-#endif
-
-!      if (planet_type.eq."earth") then
-! Write an Earth-format restart file
-        CALL dynredem0_loc("restart.nc", day_end, phis)
-!      endif
-
-      ecripar = .TRUE.
-
-#ifdef CPP_IOIPSL
-      time_step = zdtvr
-      IF (mpi_rank==0) then
-        if (ok_dyn_ins) then
-          ! initialize output file for instantaneous outputs
-          ! t_ops = iecri * daysec ! do operations every t_ops
-          t_ops =((1.0*iecri)/day_step) * daysec  
-          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
-          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
-          CALL inithist(day_ref,annee_ref,time_step,
-     &                  t_ops,t_wrt)
-        endif
-
-      IF (ok_dyn_ave) THEN 
-         ! initialize output file for averaged outputs
-         t_ops = iperiod * time_step ! do operations every t_ops
-         t_wrt = periodav * daysec   ! write output every t_wrt
-         CALL initdynav_loc(day_ref,annee_ref,time_step,t_ops,t_wrt)
-        END IF
-      ENDIF
-      dtav = iperiod*dtvr/daysec
-#endif
-! #endif of #ifdef CPP_IOIPSL
-
-c  Choix des frequences de stokage pour le offline
-c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
-c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
-      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
-      istphy=istdyn/iphysiq     
-
-
-c
-c-----------------------------------------------------------------------
-c   Integration temporelle du modele :
-c   ----------------------------------
-
-c       write(78,*) 'ucov',ucov
-c       write(78,*) 'vcov',vcov
-c       write(78,*) 'teta',teta
-c       write(78,*) 'ps',ps
-c       write(78,*) 'q',q
-
-c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
-      CALL leapfrog_loc(ucov,vcov,teta,ps,masse,phis,q,time_0)
-c$OMP END PARALLEL
-
-!      OPEN(unit=5487,file='ok_lmdz',status='replace')
-!      WRITE(5487,*) 'ok_lmdz'
-!      CLOSE(5487)
-      END
-
Index: LMDZ5/branches/testing/libf/dyn3dmem/gcm.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/gcm.F90	(revision 2298)
+++ LMDZ5/branches/testing/libf/dyn3dmem/gcm.F90	(revision 2298)
@@ -0,0 +1,500 @@
+! $Id: $
+
+PROGRAM gcm
+
+#ifdef CPP_IOIPSL
+  USE IOIPSL
+#endif
+
+  USE mod_const_mpi, ONLY: init_const_mpi
+  USE parallel_lmdz
+  USE infotrac
+#ifdef CPP_PHYS
+  USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys
+#endif
+  USE mod_hallo
+  USE Bands
+  USE filtreg_mod
+  USE control_mod
+
+#ifdef INCA
+  ! Only INCA needs these informations (from the Earth's physics)
+  USE indice_sol_mod
+  USE mod_phys_lmdz_omp_data, ONLY: klon_omp 
+#endif
+
+#ifdef CPP_PHYS
+  !      USE mod_grid_phy_lmdz
+  !      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
+  !      USE dimphy
+  !      USE comgeomphy
+#endif
+  IMPLICIT NONE
+
+  !      ......   Version  du 10/01/98    ..........
+
+  !             avec  coordonnees  verticales hybrides 
+  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
+
+  !=======================================================================
+  !
+  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
+  !   -------
+  !
+  !   Objet:
+  !   ------
+  !
+  !   GCM LMD nouvelle grille
+  !
+  !=======================================================================
+  !
+  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
+  !      et possibilite d'appeler une fonction f(y)  a derivee tangente
+  !      hyperbolique a la  place de la fonction a derivee sinusoidale.
+  !  ... Possibilite de choisir le schema pour l'advection de
+  !        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
+  !
+  !      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
+  !      Pour Van-Leer iadv=10
+  !
+  !-----------------------------------------------------------------------
+  !   Declarations:
+  !   -------------
+  include "dimensions.h"
+  include "paramet.h"
+  include "comconst.h"
+  include "comdissnew.h"
+  include "comvert.h"
+  include "comgeom.h"
+  include "logic.h"
+  include "temps.h"
+  include "ener.h"
+  include "description.h"
+  include "serre.h"
+  !include "com_io_dyn.h"
+  include "iniprint.h"
+  include "tracstoke.h"
+
+#ifdef INCA
+  ! Only INCA needs these informations (from the Earth's physics)
+  !include "indicesol.h"
+#endif
+
+  REAL zdtvr
+
+  !   variables dynamiques
+  REAL,ALLOCATABLE,SAVE  :: vcov(:,:),ucov(:,:) ! vents covariants
+  REAL,ALLOCATABLE,SAVE  :: teta(:,:)     ! temperature potentielle 
+  REAL, ALLOCATABLE,SAVE :: q(:,:,:)      ! champs advectes
+  REAL,ALLOCATABLE,SAVE  :: ps(:)         ! pression  au sol
+  !      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
+  REAL,ALLOCATABLE,SAVE  :: masse(:,:)    ! masse d'air
+  REAL,ALLOCATABLE,SAVE  :: phis(:)       ! geopotentiel au sol
+  !      REAL phi(ip1jmp1,llm)                  ! geopotentiel
+  !      REAL w(ip1jmp1,llm)                    ! vitesse verticale
+
+  ! variables dynamiques intermediaire pour le transport
+
+  !   variables pour le fichier histoire
+  REAL dtav      ! intervalle de temps elementaire
+
+  REAL time_0
+
+  LOGICAL lafin
+
+  real time_step, t_wrt, t_ops
+
+  !+jld variables test conservation energie
+  !      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
+  !     Tendance de la temp. potentiel d (theta)/ d t due a la 
+  !     tansformation d'energie cinetique en energie thermique
+  !     cree par la dissipation
+  !      REAL dhecdt(ip1jmp1,llm)
+  !      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
+  !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
+  !      CHARACTER (len=15) :: ztit
+  !-jld 
+
+
+  character (len=80) :: dynhist_file, dynhistave_file
+  character (len=20) :: modname
+  character (len=80) :: abort_message
+  ! locales pour gestion du temps
+  INTEGER :: an, mois, jour
+  REAL :: heure
+
+
+  !-----------------------------------------------------------------------
+  !   Initialisations:
+  !   ----------------
+
+  abort_message = 'last timestep reached'
+  modname = 'gcm'
+  descript = 'Run GCM LMDZ'
+  lafin    = .FALSE.
+  dynhist_file = 'dyn_hist'
+  dynhistave_file = 'dyn_hist_ave'
+
+
+
+  !----------------------------------------------------------------------
+  !  lecture des fichiers gcm.def ou run.def
+  !  ---------------------------------------
+  !
+  CALL conf_gcm( 99, .TRUE. )
+  if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm", &
+       "iphysiq must be a multiple of iperiod", 1)
+  !
+  !
+  !------------------------------------
+  !   Initialisation partie parallele
+  !------------------------------------
+  CALL init_const_mpi
+
+  call init_parallel
+  call Read_Distrib
+
+#ifdef CPP_PHYS
+  CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
+  !#endif
+  !      CALL set_bands
+  !#ifdef CPP_PHYS
+  CALL Init_interface_dyn_phys
+#endif
+  CALL barrier
+
+  CALL set_bands
+  if (mpi_rank==0) call WriteBands
+  call Set_Distrib(distrib_caldyn)
+
+  !$OMP PARALLEL
+  call Init_Mod_hallo
+  !$OMP END PARALLEL
+
+  !#ifdef CPP_PHYS
+  !c$OMP PARALLEL
+  !      call InitComgeomphy ! now done in iniphysiq
+  !c$OMP END PARALLEL 
+  !#endif
+
+  !-----------------------------------------------------------------------
+  !   Choix du calendrier
+  !   -------------------
+
+  !      calend = 'earth_365d'
+
+#ifdef CPP_IOIPSL
+  if (calend == 'earth_360d') then
+     call ioconf_calendar('360d')
+     write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
+  else if (calend == 'earth_365d') then
+     call ioconf_calendar('noleap')
+     write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
+  else if (calend == 'gregorian') then
+     call ioconf_calendar('gregorian')
+     write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
+  else
+     abort_message = 'Mauvais choix de calendrier'
+     call abort_gcm(modname,abort_message,1)
+  endif
+#endif
+
+  IF (type_trac == 'inca') THEN
+#ifdef INCA
+     call init_const_lmdz( &
+          nbtr,anneeref,dayref, &
+          iphysiq,day_step,nday,  &
+          nbsrf, is_oce,is_sic, &
+          is_ter,is_lic, calend)
+
+     call init_inca_para( &
+          iim,jjm+1,llm,klon_glo,mpi_size, &
+          distrib_phys,COMM_LMDZ)
+#endif
+  END IF
+
+  !-----------------------------------------------------------------------
+  !   Initialisation des traceurs
+  !   ---------------------------
+  !  Choix du nombre de traceurs et du schema pour l'advection
+  !  dans fichier traceur.def, par default ou via INCA
+  call infotrac_init
+
+  ! Allocation de la tableau q : champs advectes   
+  ALLOCATE(ucov(ijb_u:ije_u,llm))
+  ALLOCATE(vcov(ijb_v:ije_v,llm))
+  ALLOCATE(teta(ijb_u:ije_u,llm))
+  ALLOCATE(masse(ijb_u:ije_u,llm))
+  ALLOCATE(ps(ijb_u:ije_u))
+  ALLOCATE(phis(ijb_u:ije_u))
+  ALLOCATE(q(ijb_u:ije_u,llm,nqtot))
+
+  !-----------------------------------------------------------------------
+  !   Lecture de l'etat initial :
+  !   ---------------------------
+
+  !  lecture du fichier start.nc
+  if (read_start) then
+     ! we still need to run iniacademic to initialize some
+     ! constants & fields, if we run the 'newtonian' or 'SW' cases:
+     if (iflag_phys.ne.1) then
+        CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
+     endif
+
+     !        if (planet_type.eq."earth") then
+     ! Load an Earth-format start file
+     CALL dynetat0_loc("start.nc",vcov,ucov, &
+          teta,q,masse,ps,phis, time_0)
+     !        endif ! of if (planet_type.eq."earth")
+
+     !       write(73,*) 'ucov',ucov
+     !       write(74,*) 'vcov',vcov
+     !       write(75,*) 'teta',teta
+     !       write(76,*) 'ps',ps
+     !       write(77,*) 'q',q
+
+  endif ! of if (read_start)
+
+  ! le cas echeant, creation d un etat initial
+  IF (prt_level > 9) WRITE(lunout,*) &
+       'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
+  if (.not.read_start) then
+     CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
+  endif
+
+  !-----------------------------------------------------------------------
+  !   Lecture des parametres de controle pour la simulation :
+  !   -------------------------------------------------------
+  !  on recalcule eventuellement le pas de temps
+
+  IF(MOD(day_step,iperiod).NE.0) THEN
+     abort_message =  &
+          'Il faut choisir un nb de pas par jour multiple de iperiod'
+     call abort_gcm(modname,abort_message,1)
+  ENDIF
+
+  IF(MOD(day_step,iphysiq).NE.0) THEN
+     abort_message =  &
+          'Il faut choisir un nb de pas par jour multiple de iphysiq'
+     call abort_gcm(modname,abort_message,1)
+  ENDIF
+
+  zdtvr    = daysec/REAL(day_step)
+  IF(dtvr.NE.zdtvr) THEN
+     WRITE(lunout,*) &
+          'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
+  ENDIF
+
+  !
+  ! on remet le calendrier \`a zero si demande
+  !
+  IF (start_time /= starttime) then
+     WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' &
+          ,' fichier restart ne correspond pas a celle lue dans le run.def'
+     IF (raz_date == 1) then
+        WRITE(lunout,*)'Je prends l''heure lue dans run.def'
+        start_time = starttime
+     ELSE
+        WRITE(lunout,*)'Je m''arrete'
+        CALL abort
+     ENDIF
+  ENDIF
+  IF (raz_date == 1) THEN
+     annee_ref = anneeref
+     day_ref = dayref
+     day_ini = dayref
+     itau_dyn = 0
+     itau_phy = 0
+     time_0 = 0.
+     write(lunout,*) &
+          'GCM: On reinitialise a la date lue dans gcm.def'
+  ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
+     write(lunout,*) &
+          'GCM: Attention les dates initiales lues dans le fichier'
+     write(lunout,*) &
+          ' restart ne correspondent pas a celles lues dans '
+     write(lunout,*)' gcm.def'
+     write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
+     write(lunout,*)' day_ref=',day_ref," dayref=",dayref
+     write(lunout,*)' Pas de remise a zero'
+  ENDIF
+  !      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
+  !        write(lunout,*)
+  !     .  'GCM: Attention les dates initiales lues dans le fichier'
+  !        write(lunout,*)
+  !     .  ' restart ne correspondent pas a celles lues dans '
+  !        write(lunout,*)' gcm.def'
+  !        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
+  !        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
+  !        if (raz_date .ne. 1) then
+  !          write(lunout,*)
+  !     .    'GCM: On garde les dates du fichier restart'
+  !        else
+  !          annee_ref = anneeref
+  !          day_ref = dayref
+  !          day_ini = dayref
+  !          itau_dyn = 0
+  !          itau_phy = 0
+  !          time_0 = 0.
+  !          write(lunout,*)
+  !     .   'GCM: On reinitialise a la date lue dans gcm.def'
+  !        endif
+  !      ELSE
+  !        raz_date = 0
+  !      endif
+
+#ifdef CPP_IOIPSL
+  mois = 1
+  heure = 0.
+  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
+  jH_ref = jD_ref - int(jD_ref)
+  jD_ref = int(jD_ref)
+
+  call ioconf_startdate(INT(jD_ref), jH_ref)
+
+  write(lunout,*)'DEBUG'
+  write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
+  write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
+  call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
+  write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
+  write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
+#else
+  ! Ehouarn: we still need to define JD_ref and JH_ref
+  ! and since we don't know how many days there are in a year
+  ! we set JD_ref to 0 (this should be improved ...)
+  jD_ref=0
+  jH_ref=0
+#endif
+
+  if (iflag_phys.eq.1) then
+     ! these initialisations have already been done (via iniacademic)
+     ! if running in SW or Newtonian mode
+     !-----------------------------------------------------------------------
+     !   Initialisation des constantes dynamiques :
+     !   ------------------------------------------
+     dtvr = zdtvr
+     CALL iniconst
+
+     !-----------------------------------------------------------------------
+     !   Initialisation de la geometrie :
+     !   --------------------------------
+     CALL inigeom
+
+     !-----------------------------------------------------------------------
+     !   Initialisation du filtre :
+     !   --------------------------
+     CALL inifilr
+  endif ! of if (iflag_phys.eq.1)
+  !
+  !-----------------------------------------------------------------------
+  !   Initialisation de la dissipation :
+  !   ----------------------------------
+
+  CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   , &
+       tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
+
+  !-----------------------------------------------------------------------
+  !   Initialisation de la physique :
+  !   -------------------------------
+  IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
+     ! Physics:
+#ifdef CPP_PHYS
+     CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, &
+          rlatu,rlonv,aire,cu,cv,rad,g,r,cpp, &
+          iflag_phys)
+#endif
+  ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
+
+
+  !-----------------------------------------------------------------------
+  !   Initialisation des dimensions d'INCA :
+  !   --------------------------------------
+  IF (type_trac == 'inca') THEN
+     !$OMP PARALLEL
+#ifdef INCA
+     CALL init_inca_dim(klon_omp,llm,iim,jjm, &
+          rlonu,rlatu,rlonv,rlatv)
+#endif
+     !$OMP END PARALLEL
+  END IF
+
+  !-----------------------------------------------------------------------
+  !   Initialisation des I/O :
+  !   ------------------------
+
+
+  if (nday>=0) then
+     day_end = day_ini + nday
+  else
+     day_end = day_ini - nday/day_step
+  endif
+
+  WRITE(lunout,300)day_ini,day_end
+300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
+
+#ifdef CPP_IOIPSL
+  call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
+  write (lunout,301)jour, mois, an
+  call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
+  write (lunout,302)jour, mois, an
+301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
+302 FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
+#endif
+
+  !      if (planet_type.eq."earth") then
+  ! Write an Earth-format restart file
+  CALL dynredem0_loc("restart.nc", day_end, phis)
+  !      endif
+
+  ecripar = .TRUE.
+
+#ifdef CPP_IOIPSL
+  time_step = zdtvr
+  IF (mpi_rank==0) then
+     if (ok_dyn_ins) then
+        ! initialize output file for instantaneous outputs
+        ! t_ops = iecri * daysec ! do operations every t_ops
+        t_ops =((1.0*iecri)/day_step) * daysec  
+        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
+        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
+        CALL inithist(day_ref,annee_ref,time_step, &
+             t_ops,t_wrt)
+     endif
+
+     IF (ok_dyn_ave) THEN 
+        ! initialize output file for averaged outputs
+        t_ops = iperiod * time_step ! do operations every t_ops
+        t_wrt = periodav * daysec   ! write output every t_wrt
+        CALL initdynav_loc(day_ref,annee_ref,time_step,t_ops,t_wrt)
+     END IF
+  ENDIF
+  dtav = iperiod*dtvr/daysec
+#endif
+  ! #endif of #ifdef CPP_IOIPSL
+
+  !  Choix des frequences de stokage pour le offline
+  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
+  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
+  istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
+  istphy=istdyn/iphysiq     
+
+
+  !
+  !-----------------------------------------------------------------------
+  !   Integration temporelle du modele :
+  !   ----------------------------------
+
+  !       write(78,*) 'ucov',ucov
+  !       write(78,*) 'vcov',vcov
+  !       write(78,*) 'teta',teta
+  !       write(78,*) 'ps',ps
+  !       write(78,*) 'q',q
+
+  !$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
+  CALL leapfrog_loc(ucov,vcov,teta,ps,masse,phis,q,time_0)
+  !$OMP END PARALLEL
+
+  !      OPEN(unit=5487,file='ok_lmdz',status='replace')
+  !      WRITE(5487,*) 'ok_lmdz'
+  !      CLOSE(5487)
+END PROGRAM gcm
Index: LMDZ5/branches/testing/libf/dyn3dmem/gr_dyn_fi_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/gr_dyn_fi_p.F	(revision 2258)
+++ 	(revision )
@@ -1,44 +1,0 @@
-!
-! $Id$
-!
-      SUBROUTINE gr_dyn_fi_p(nfield,im,jm,ngrid,pdyn,pfi)
-#ifdef CPP_PHYS
-! Interface with parallel physics,
-      USE mod_interface_dyn_phys
-      USE dimphy
-      USE parallel_lmdz
-      IMPLICIT NONE
-c=======================================================================
-c   passage d'un champ de la grille scalaire a la grille physique
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-      INTEGER im,jm,ngrid,nfield
-      REAL pdyn(im,jm,nfield)
-      REAL pfi(ngrid,nfield)
-
-      INTEGER i,j,ig,l
-
-c-----------------------------------------------------------------------
-c   calcul:
-c   -------
-
-c      IF(ngrid.NE.2+(jm-2)*(im-1)) STOP 'probleme de dim'
-c   traitement des poles
-c   traitement des point normaux
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l=1,nfield    
-       DO ig=1,klon
-         i=index_i(ig)
-         j=index_j(ig)
-         pfi(ig,l)=pdyn(i,j,l)
-       ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
-#endif
-! of #ifdef CPP_PHYS
-      RETURN
-      END
Index: LMDZ5/branches/testing/libf/dyn3dmem/gr_fi_dyn_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/gr_fi_dyn_p.F	(revision 2258)
+++ 	(revision )
@@ -1,56 +1,0 @@
-!
-! $Id$
-!
-      SUBROUTINE gr_fi_dyn_p(nfield,ngrid,im,jm,pfi,pdyn)
-#ifdef CPP_PHYS
-! Interface with parallel physics,
-      USE mod_interface_dyn_phys
-      USE dimphy
-      USE parallel_lmdz
-      IMPLICIT NONE
-c=======================================================================
-c   passage d'un champ de la grille scalaire a la grille physique
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-      INTEGER im,jm,ngrid,nfield
-      REAL pdyn(im,jm,nfield)
-      REAL pfi(ngrid,nfield)
-
-      INTEGER i,j,ifield,ig
-
-c-----------------------------------------------------------------------
-c   calcul:
-c   -------
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO ifield=1,nfield
-
-        do ig=1,klon
-          i=index_i(ig)
-          j=index_j(ig)
-          pdyn(i,j,ifield)=pfi(ig,ifield)
-          if (i==1) pdyn(im,j,ifield)=pdyn(i,j,ifield)
-	enddo
-
-c   traitement des poles
-      if (pole_nord) then
-        do i=1,im
-	  pdyn(i,1,ifield)=pdyn(1,1,ifield)
-	enddo
-      endif
-       
-      if (pole_sud) then
-        do i=1,im
-	  pdyn(i,jm,ifield)=pdyn(1,jm,ifield)
-	enddo
-      endif
-      
-      ENDDO
-c$OMP END DO NOWAIT
-#endif
-! of #ifdef CPP_PHYS
-      RETURN
-      END
Index: LMDZ5/branches/testing/libf/dyn3dmem/guide_loc_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/guide_loc_mod.F90	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/guide_loc_mod.F90	(revision 2298)
@@ -87,4 +87,5 @@
 ! Lecture des parametres:  
 ! ---------------------------------------------
+    call ini_getparam("nudging_parameters_out.txt")
 ! Variables guidees
     CALL getpar('guide_u',.true.,guide_u,'guidage de u')
@@ -159,4 +160,6 @@
     CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
 
+    call fin_getparam
+    
 ! ---------------------------------------------
 ! Determination du nombre de niveaux verticaux
Index: LMDZ5/branches/testing/libf/dyn3dmem/iniacademic_loc.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/iniacademic_loc.F90	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/iniacademic_loc.F90	(revision 2298)
@@ -7,5 +7,6 @@
   use exner_hyb_m, only: exner_hyb
   use exner_milieu_m, only: exner_milieu
-  USE infotrac, ONLY : nqtot
+  USE infotrac, ONLY: nqtot,niso_possibles,ok_isotopes,iqpere,ok_iso_verif,tnat,alpha_ideal, &
+        & iqiso,phase_num,iso_indnum,iso_num,zone_num
   USE control_mod, ONLY: day_step,planet_type
   USE parallel_lmdz, ONLY: ijb_u, ije_u, ijb_v, ije_v
@@ -110,5 +111,5 @@
   ztot0      = 0.
   stot0      = 0.
-  ang0       = 0.
+  ang0       = 0.      
 
   if (llm == 1) then
@@ -269,12 +270,32 @@
         if (planet_type=="earth") then
            ! Earth: first two tracers will be water
+
            do i=1,nqtot
               if (i == 1) q(ijb_u:ije_u,:,i)=1.e-10
               if (i == 2) q(ijb_u:ije_u,:,i)=1.e-15
               if (i.gt.2) q(ijb_u:ije_u,:,i)=0.
+
+              ! CRisi: init des isotopes
+              ! distill de Rayleigh très simplifiée
+              if (ok_isotopes) then
+                if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then          
+                   q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqpere(i))       &
+      &                  *tnat(iso_num(i))                             &
+      &                  *(q(ijb_u:ije_u,:,iqpere(i))/30.e-3)                              &
+     &                   **(alpha_ideal(iso_num(i))-1)
+                endif                
+                if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then
+                  q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqiso(iso_indnum(i),phase_num(i)))
+                endif
+              endif !if (ok_isotopes) then
+
            enddo
         else
            q(ijb_u:ije_u,:,:)=0
         endif ! of if (planet_type=="earth")
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'iniacademic_loc')
+        endif !if (ok_iso_verif) then
 
         ! add random perturbation to temperature
Index: LMDZ5/branches/testing/libf/dyn3dmem/integrd_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/integrd_loc.F	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/integrd_loc.F	(revision 2298)
@@ -11,4 +11,5 @@
       USE write_field
       USE integrd_mod
+      USE infotrac, ONLY: ok_iso_verif ! ajout CRisi
       IMPLICIT NONE
 
@@ -86,5 +87,7 @@
       INTEGER :: ierr
 
+      !write(*,*) 'integrd 88: entree, nq=',nq
 c-----------------------------------------------------------------------
+
 c$OMP BARRIER     
       if (pole_nord) THEN
@@ -125,6 +128,8 @@
       DO 2 ij = ijb,ije
        pscr (ij)    = ps0(ij)
-       ps (ij)      = psm1(ij) + dt * dp(ij)
+       ps (ij)      = psm1(ij) + dt * dp(ij)     
+
    2  CONTINUE
+
 c$OMP END DO  
 c$OMP BARRIER
@@ -159,4 +164,5 @@
 c$OMP END MASTER
 c$OMP BARRIER
+        !write(*,*) 'integrd 170'
       IF (.NOT. Checksum_all) THEN
         call WriteField_v('int_vcov',vcov)
@@ -188,4 +194,5 @@
        
 c
+        !write(*,*) 'integrd 200'
 C$OMP MASTER
       if (pole_nord) THEN
@@ -214,4 +221,5 @@
 c$OMP END MASTER
 c$OMP BARRIER
+      !write(*,*) 'integrd 217'  
 c
 c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
@@ -219,4 +227,5 @@
 
       CALL pression_loc ( ip1jmp1, ap, bp, ps, p )
+
 c$OMP BARRIER
       CALL massdair_loc (     p  , masse         )
@@ -276,4 +285,5 @@
 c
 c
+        !write(*,*) 'integrd 291'
       IF (pole_nord) THEN
        
@@ -334,8 +344,19 @@
            ENDDO
           ENDDO
+          
 c$OMP END DO NOWAIT
 c$OMP BARRIER
 
-          CALL qminimum_loc( q, nq, deltap )
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'integrd 342')
+        endif !if (ok_iso_verif) then
+
+        !write(*,*) 'integrd 341'
+        CALL qminimum_loc( q, nq, deltap )
+        !write(*,*) 'integrd 343'
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'integrd 346')
+        endif !if (ok_iso_verif) then
 c
 c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
@@ -387,4 +408,8 @@
       
       ENDIF
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'integrd 409')
+        endif !if (ok_iso_verif) then
       
 ! Ehouarn: forget about finvmaold
@@ -404,4 +429,5 @@
 
 15    continue
+          !write(*,*) 'integrd 410'
 
 c$OMP DO SCHEDULE(STATIC)
Index: LMDZ5/branches/testing/libf/dyn3dmem/leapfrog_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/leapfrog_loc.F	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/leapfrog_loc.F	(revision 2298)
@@ -200,4 +200,8 @@
       LOGICAL, SAVE :: firstcall=.TRUE.
       TYPE(distrib),SAVE :: new_dist
+
+      if (ok_iso_verif) then
+         call check_isotopes(q0,ijb_u,ije_u,'leapfrog204: debut')
+      endif !if (ok_iso_verif) then
       
 c$OMP MASTER
@@ -219,4 +223,8 @@
       itaufinp1 = itaufin +1
 
+      if (ok_iso_verif) then
+        call check_isotopes(q0,ijb_u,ije_u,'leapfrog 226')
+      endif !if (ok_iso_verif) then
+
       itau = 0
       physic=.true.
@@ -231,4 +239,8 @@
       phis=phis0
       q=q0
+
+      if (ok_iso_verif) then
+        call check_isotopes(q,ijb_u,ije_u,'leapfrog 239')
+      endif !if (ok_iso_verif) then
       
 !      iday = day_ini+itau/day_step
@@ -296,5 +308,5 @@
 
    1  CONTINUE ! Matsuno Forward step begins here
-
+      !write(*,*) 'leapfrog 298: itau=',itau
       jD_cur = jD_ref + day_ini - day_ref +                             &
      &          itau/day_step
@@ -306,4 +318,7 @@
       endif
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 321')
+        endif !if (ok_iso_verif) then
 
 #ifdef CPP_IOIPSL
@@ -384,5 +399,15 @@
 cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
 
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 400')
+        endif !if (ok_iso_verif) then
+
    2  CONTINUE ! Matsuno backward or leapfrog step begins here
+
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 402')
+        endif !if (ok_iso_verif) then
 
 c$OMP MASTER
@@ -455,4 +480,8 @@
 c$OMP END MASTER      
 
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 471')
+        endif !if (ok_iso_verif) then
 
 !ym  PAS D'AJUSTEMENT POUR LE MOMENT     
@@ -574,4 +603,7 @@
      
       
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 589')
+        endif !if (ok_iso_verif) then
       
 c-----------------------------------------------------------------------
@@ -635,5 +667,8 @@
       ! compute geopotential phi()
       CALL geopot_loc  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
-
+        
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 651')
+        endif !if (ok_iso_verif) then
       
       call VTb(VTcaldyn)
@@ -644,4 +679,5 @@
 !      CALL FTRACE_REGION_BEGIN("caldyn")
       time = jD_cur + jH_cur 
+
       CALL caldyn_loc 
      $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
@@ -670,10 +706,21 @@
 c   -------------------------------------------------------------
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,
+     &           'leapfrog 686: avant caladvtrac')
+        endif !if (ok_iso_verif) then
       
       IF( forward. OR . leapf )  THEN
 ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
+        !write(*,*) 'leapfrog 679: avant CALL caladvtrac_loc'
          CALL caladvtrac_loc(q,pbaru,pbarv,
      *        p, masse, dq,  teta,
      .        flxw,pk, iapptrac)
+
+         !write(*,*) 'leapfrog 719'
+         if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,
+     &           'leapfrog 698: apres caladvtrac')
+         endif !if (ok_iso_verif) then
 
 !      do j=1,nqtot
@@ -708,8 +755,19 @@
 !       CALL FTRACE_REGION_BEGIN("integrd")
 
-       CALL integrd_loc ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
+       !write(*,*) 'leapfrog 720' 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 756')
+        endif !if (ok_iso_verif) then
+
+       ! CRisi: pourquoi aller jusqu'à 2 et non pas jusqu'à nqtot?? 
+       CALL integrd_loc ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
      $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis)
 !     $              finvmaold                                    )
 
+       !write(*,*) 'leapfrog 724'       
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 762')
+        endif !if (ok_iso_verif) then 
+ 
 !       CALL FTRACE_REGION_END("integrd")
 c$OMP BARRIER
@@ -724,4 +782,9 @@
       call WriteField_u('ps_int',ps)
 #endif    
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 775')
+        endif !if (ok_iso_verif) then
+
 c      do j=1,nqtot
 c        call WriteField_p('q'//trim(int2str(j)),
@@ -1082,4 +1145,9 @@
        ENDIF ! of IF( apphys )
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1132')
+        endif !if (ok_iso_verif) then
+        !write(*,*) 'leapfrog 1134: iflag_phys=',iflag_phys
+
       IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
 c$OMP MASTER
@@ -1146,9 +1214,12 @@
 
 cc$OMP END PARALLEL
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1196')
+        endif !if (ok_iso_verif) then
 
 c-----------------------------------------------------------------------
 c   dissipation horizontale et verticale  des petites echelles:
 c   ----------------------------------------------------------
-
+      !write(*,*) 'leapfrog 1163: apdiss=',apdiss
       IF(apdiss) THEN
       
@@ -1379,5 +1450,9 @@
 c                call abort_gcm(modname,abort_message,0)
 c              ENDIF
-        
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1430')
+        endif !if (ok_iso_verif) then     
+ 
 c   ********************************************************************
 c   ********************************************************************
@@ -1455,4 +1530,8 @@
       ENDIF
       
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1509')
+        endif !if (ok_iso_verif) then
+
       IF ( .NOT.purmats ) THEN
 c       ........................................................
@@ -1526,4 +1605,8 @@
             ENDIF
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1584')
+        endif !if (ok_iso_verif) then
+
 c-----------------------------------------------------------------------
 c   ecriture de la bande histoire:
@@ -1562,4 +1645,8 @@
             ENDIF ! of IF (itau.EQ.itaufin)
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1624')
+        endif !if (ok_iso_verif) then
+
 c-----------------------------------------------------------------------
 c   gestion de l'integration temporelle:
@@ -1596,4 +1683,9 @@
 
       ELSE ! of IF (.not.purmats)
+
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1664')
+        endif !if (ok_iso_verif) then
 
 c       ........................................................
@@ -1631,4 +1723,9 @@
 
             ELSE ! of IF(forward) i.e. backward step
+
+              
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1698')
+        endif !if (ok_iso_verif) then  
 
               IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
@@ -1683,4 +1780,9 @@
             ENDIF ! of IF (forward)
 
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1750')
+        endif !if (ok_iso_verif) then
+
       END IF ! of IF(.not.purmats)
 c$OMP MASTER
Index: LMDZ5/branches/testing/libf/dyn3dmem/mod_interface_dyn_phys.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/mod_interface_dyn_phys.F90	(revision 2258)
+++ 	(revision )
@@ -1,58 +1,0 @@
-! 
-! $Id$
-!
-MODULE mod_interface_dyn_phys
-  INTEGER,SAVE,dimension(:),allocatable :: index_i
-  INTEGER,SAVE,dimension(:),allocatable :: index_j
-  
-  
-#ifdef CPP_PHYS
-! Interface with parallel physics,
-CONTAINS
-  
-  SUBROUTINE Init_interface_dyn_phys
-    USE mod_phys_lmdz_mpi_data
-    IMPLICIT NONE
-    include 'dimensions.h'    
-    
-    INTEGER :: i,j,k
-    
-    ALLOCATE(index_i(klon_mpi))
-    ALLOCATE(index_j(klon_mpi))
-    
-    k=1
-    IF (is_north_pole) THEN
-      index_i(k)=1
-      index_j(k)=1
-      k=2
-    ELSE
-      DO i=ii_begin,iim
-	index_i(k)=i
-	index_j(k)=jj_begin
-	k=k+1
-       ENDDO
-    ENDIF
-    
-    DO j=jj_begin+1,jj_end-1
-      DO i=1,iim
-	index_i(k)=i
-	index_j(k)=j
-	k=k+1
-      ENDDO
-    ENDDO
-    
-    IF (is_south_pole) THEN
-      index_i(k)=1
-      index_j(k)=jj_end
-    ELSE
-      DO i=1,ii_end
-	index_i(k)=i
-	index_j(k)=jj_end
-	k=k+1
-       ENDDO
-    ENDIF
-  
-  END SUBROUTINE Init_interface_dyn_phys 
-#endif
-! of #ifdef CPP_PHYS
-END MODULE mod_interface_dyn_phys
Index: LMDZ5/branches/testing/libf/dyn3dmem/qminimum_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/qminimum_loc.F	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/qminimum_loc.F	(revision 2298)
@@ -1,4 +1,5 @@
-      SUBROUTINE qminimum_loc( q,nq,deltap )
+      SUBROUTINE qminimum_loc( q,nqtot,deltap )
       USE parallel_lmdz
+      USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif
       IMPLICIT none
 c
@@ -10,6 +11,6 @@
 #include "comvert.h"
 c
-      INTEGER nq
-      REAL q(ijb_u:ije_u,llm,nq), deltap(ijb_u:ije_u,llm)
+      INTEGER nqtot ! CRisi: on remplace nq par nqtot
+      REAL q(ijb_u:ije_u,llm,nqtot), deltap(ijb_u:ije_u,llm)
 c
       INTEGER iq_vap, iq_liq
@@ -27,4 +28,7 @@
       INTEGER i, k, iq
       REAL zx_defau, zx_abc, zx_pump(ijb_u:ije_u), pompe
+
+      real zx_defau_diag(ijb_u:ije_u,llm,2) 
+      real q_follow(ijb_u:ije_u,llm,2) 
 c
       REAL SSUM
@@ -38,4 +42,6 @@
       INTEGER Index_pump(ij_end-ij_begin+1)
       INTEGER nb_pump
+      INTEGER ixt
+      INTEGER iso_verif_noNaN_nostop
 c
 c Quand l'eau liquide est trop petite (ou negative), on prend
@@ -44,11 +50,26 @@
 c
 
+        !write(*,*) 'qminimum 52: entree'
+        if (ok_iso_verif) then
+           call check_isotopes(q,ij_begin,ij_end,'qminimum 52')   
+        endif !if (ok_iso_verif) then     
+
       ijb=ij_begin
       ije=ij_end
 
+      zx_defau_diag(ijb:ije,:,:)=0.0
+      q_follow(ijb:ije,:,1:2)=q(ijb:ije,:,1:2)  
+
+      !write(*,*) 'qminimum 57'
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
       DO 1000 k = 1, llm
       DO 1040 i = ijb, ije
             if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
+
+              if (ok_isotopes) then
+                 zx_defau_diag(i,k,iq_liq)=AMAX1
+     :               ( seuil_liq - q(i,k,iq_liq), 0.0 )
+              endif !if (ok_isotopes) 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
@@ -60,8 +81,10 @@
 c --->  SYNCHRO OPENMP ICI
 
+
 c
 c Quand l'eau vapeur est trop faible (ou negative), on complete
 c le defaut en prennant de l'eau vapeur de la couche au-dessous.
 c
+      !write(*,*) 'qminimum 81'
       iq = iq_vap
 c
@@ -70,8 +93,15 @@
 c$OMP DO SCHEDULE(STATIC)
       DO i = ijb, ije
+
          if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then
+
+            if (ok_isotopes) then
+              zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 )
+            endif !if (ok_isotopes) 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
@@ -79,8 +109,10 @@
       ENDDO
 c$OMP BARRIER
+
 c
 c Quand il s'agit de la premiere couche au-dessus du sol, on
 c doit imprimer un message d'avertissement (saturation possible).
 c
+      !write(*,*) 'qminimum 106'
       nb_pump=0
 c$OMP DO SCHEDULE(STATIC)
@@ -103,4 +135,114 @@
          ENDDO
       ENDIF
+
+      !write(*,*) 'qminimum 128'
+      if (ok_isotopes) then
+      ! CRisi: traiter de même les traceurs d'eau
+      ! Mais il faut les prendre à l'envers pour essayer de conserver la
+      ! masse.
+      ! 1) pompage dans le sol  
+      ! On suppose que ce pompage se fait sans isotopes -> on ne modifie
+      ! rien ici et on croise les doigts pour que ça ne soit pas trop
+      ! génant
+      DO i = ijb, ije
+        if (zx_pump(i).gt.0.0) then
+          q_follow(i,1,iq_vap)=q_follow(i,1,iq_vap)+zx_pump(i)
+        endif !if (zx_pump(i).gt.0.0) then
+      enddo !DO i = ijb, ije  
+
+      ! 2) transfert de vap vers les couches plus hautes
+      !write(*,*) 'qminimum 139'
+      do k=2,llm
+        DO i = ijb, ije
+          if (zx_defau_diag(i,k,iq_vap).gt.0.0) then             
+              ! on ajoute la vapeur en k              
+              do ixt=1,ntraciso
+               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
+     :              +zx_defau_diag(i,k,iq_vap)
+     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
+                
+              if (ok_iso_verif) then
+                if (iso_verif_noNaN_nostop(q(i,k,iqiso(ixt,iq_vap)),
+     :                   'qminimum 155').eq.1) then
+                   write(*,*) 'i,k,ixt=',i,k,ixt
+                   write(*,*) 'q_follow(i,k-1,iq_vap)=',
+     :                   q_follow(i,k-1,iq_vap)
+                   write(*,*) 'q(i,k,iqiso(ixt,iq_vap))=',
+     :                   q(i,k,iqiso(ixt,iq_vap))
+                   write(*,*) 'zx_defau_diag(i,k,iq_vap)=',
+     :                   zx_defau_diag(i,k,iq_vap)
+                   write(*,*) 'q(i,k-1,iqiso(ixt,iq_vap))=',
+     :                   q(i,k-1,iqiso(ixt,iq_vap))
+                   stop
+                endif
+              endif
+
+              ! et on la retranche en k-1
+               q(i,k-1,iqiso(ixt,iq_vap))=q(i,k-1,iqiso(ixt,iq_vap))
+     :              -zx_defau_diag(i,k,iq_vap)
+     :              *deltap(i,k)/deltap(i,k-1)
+     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
+
+               if (ok_iso_verif) then
+                if (iso_verif_noNaN_nostop(q(i,k-1,iqiso(ixt,iq_vap)),
+     :                   'qminimum 175').eq.1) then
+                   write(*,*) 'k,i,ixt=',k,i,ixt
+                   write(*,*) 'q_follow(i,k-1,iq_vap)=',
+     :                   q_follow(i,k-1,iq_vap)
+                   write(*,*) 'q(i,k,iqiso(ixt,iq_vap))=',
+     :                   q(i,k,iqiso(ixt,iq_vap))
+                   write(*,*) 'zx_defau_diag(i,k,iq_vap)=',
+     :                   zx_defau_diag(i,k,iq_vap)
+                   write(*,*) 'q(i,k-1,iqiso(ixt,iq_vap))=',
+     :                   q(i,k-1,iqiso(ixt,iq_vap))
+                   stop
+                endif
+              endif 
+
+              enddo !do ixt=1,niso
+              q_follow(i,k,iq_vap)=   q_follow(i,k,iq_vap)
+     :               +zx_defau_diag(i,k,iq_vap)
+              q_follow(i,k-1,iq_vap)=   q_follow(i,k-1,iq_vap)
+     :               -zx_defau_diag(i,k,iq_vap)
+     :              *deltap(i,k)/deltap(i,k-1)
+          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
+        enddo !DO i = 1, ip1jmp1        
+       enddo !do k=2,llm
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'qminimum 168')
+        endif !if (ok_iso_verif) then
+        
+      
+        ! 3) transfert d'eau de la vapeur au liquide
+        !write(*,*) 'qminimum 164'
+        do k=1,llm
+        DO i = ijb, ije
+          if (zx_defau_diag(i,k,iq_liq).gt.0.0) then
+
+              ! on ajoute eau liquide en k en k              
+              do ixt=1,ntraciso
+               q(i,k,iqiso(ixt,iq_liq))=q(i,k,iqiso(ixt,iq_liq))
+     :              +zx_defau_diag(i,k,iq_liq)
+     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)
+              ! et on la retranche à la vapeur en k
+               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
+     :              -zx_defau_diag(i,k,iq_liq)
+     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)   
+              enddo !do ixt=1,niso
+              q_follow(i,k,iq_liq)=   q_follow(i,k,iq_liq)
+     :               +zx_defau_diag(i,k,iq_liq)
+              q_follow(i,k,iq_vap)=   q_follow(i,k,iq_vap)
+     :               -zx_defau_diag(i,k,iq_liq)
+          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
+        enddo !DO i = 1, ip1jmp1
+       enddo !do k=2,llm  
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'qminimum 197')
+        endif !if (ok_iso_verif) then
+
+      endif !if (ok_isotopes) then
+      !write(*,*) 'qminimum 188'
 c
       RETURN
Index: LMDZ5/branches/testing/libf/dyn3dmem/vlsplt_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/vlsplt_loc.F	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/vlsplt_loc.F	(revision 2298)
@@ -2,5 +2,5 @@
 ! $Id$
 !
-      SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x)
+      RECURSIVE SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x,iq)
 
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -14,4 +14,5 @@
 c   --------------------------------------------------------------------
       USE parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
       IMPLICIT NONE
 c
@@ -25,8 +26,9 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
-      REAL u_m( ijb_u:ije_u,llm ),pbarv( iip1,jjb_v:jje_v,llm)
-      REAL q(ijb_u:ije_u,llm)
-      REAL w(ijb_u:ije_u,llm)
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
+      REAL u_m( ijb_u:ije_u,llm),pbarv( iip1,jjb_v:jje_v,llm)
+      REAL q(ijb_u:ije_u,llm,nqtot) ! CRisi: ajout dimension nqtot
+      REAL w(ijb_u:ije_u,llm) 
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -42,4 +44,7 @@
       REAL u_mq(ijb_u:ije_u,llm)
 
+      REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
+
       Logical extremum
 
@@ -51,4 +56,6 @@
       INTEGER ijb,ije,ijb_x,ije_x
       
+      !write(*,*) 'vlsplt 58: entree dans vlx_loc, iq,ijb_x=',
+!     &   iq,ijb_x
 c   calcul de la pente a droite et a gauche de la maille
 
@@ -64,5 +71,5 @@
 c   calcul des pentes avec limitation, Van Leer scheme I:
 c   -----------------------------------------------------
-
+      ! on a besoin de q entre ijb et ije
 c   calcul de la pente aux points u
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
@@ -70,7 +77,7 @@
             
             DO ij=ijb,ije-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
 c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
-c              sigu(ij)=u_m(ij,l)/masse(ij,l)
+c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
             ENDDO
             DO ij=ijb+iip1-1,ije,iip1
@@ -126,5 +133,5 @@
          DO l = 1, llm
             DO ij=ijb,ije-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
             ENDDO
             DO ij=ijb+iip1-1,ije,iip1
@@ -147,4 +154,6 @@
       ENDIF ! (pente_max.lt.-1.e-5)
 
+      !write(*,*) 'vlx 156: iq,ijb_x=',iq,ijb_x
+
 c   bouclage de la pente en iip1:
 c   -----------------------------
@@ -168,11 +177,11 @@
       DO l=1,llm
        DO ij=ijb,ije-1
-          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
-     ,                     1.+u_m(ij,l)/masse(ij+1,l),
-     ,                     u_m(ij,l))
+          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
+     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
+     ,                     u_m(ij,l,iq))
           zdum(ij,l)=0.5*zdum(ij,l)
           u_mq(ij,l)=cvmgp(
-     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
-     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
+     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
+     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
      ,                u_m(ij,l))
           u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
@@ -185,13 +194,16 @@
 c	print*,'Cumule ....'
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        ! on a besoin de masse entre ijb et ije
       DO l=1,llm
        DO ij=ijb,ije-1
-c	print*,'masse(',ij,')=',masse(ij,l)
+c	print*,'masse(',ij,')=',masse(ij,l,iq)
           IF (u_m(ij,l).gt.0.) THEN
-             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
-             u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
+             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
+             u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)
+     :           +0.5*zdum(ij,l)*dxq(ij,l))
           ELSE
-             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
-             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
+             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
+             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)
+     :           -0.5*zdum(ij,l)*dxq(ij+1,l))
           ENDIF
        ENDDO
@@ -215,4 +227,5 @@
 c$OMP END DO NOWAIT
 c	print*,'Ok test 1'
+
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
@@ -223,5 +236,5 @@
 c$OMP END DO NOWAIT
 c	 print*,'Ok test 2'
-
+        
 
 c   traitement special pour le cas ou on advecte en longitude plus que le
@@ -247,4 +260,6 @@
 c     &       ,'contenu de la maille : ',n0
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+
+
          DO l=1,llm
             IF(nl(l).gt.0) THEN
@@ -258,5 +273,5 @@
                ENDDO
                niju=iju
-c              PRINT*,'niju,nl',niju,nl(l)
+               !PRINT*,'vlx 278, niju,nl',niju,nl(l)
 
 c  traitement des mailles
@@ -270,7 +285,8 @@
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m-masse(ijq,l)
+                     do while(zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)
+     &                          +q(ijq,l,iq)*masse(ijq,l,iq)
+                        zu_m=zu_m-masse(ijq,l,iq)
                         i=mod(i-2+iim,iim)+1
                         ijq=(j-1)*iip1+i
@@ -278,18 +294,20 @@
 c   ajout de la maille non completement advectee
                      u_mq(ij,l)=u_mq(ij,l)+zu_m*
-     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
+     &               (q(ijq,l,iq)+0.5*
+     &               (1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ELSE
                      ijq=ij+1
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(-zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m+masse(ijq,l)
+                     do while(-zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
+     &                           *masse(ijq,l,iq)
+                        zu_m=zu_m+masse(ijq,l,iq)
                         i=mod(i,iim)+1
                         ijq=(j-1)*iip1+i
                      ENDDO
 c   ajout de la maille non completement advectee
-                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
-     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
+                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
+     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ENDIF
                ENDDO
@@ -299,5 +317,4 @@
 cym      ENDIF  ! n0.gt.0 
 9999    continue
-
 
 c   bouclage en latitude
@@ -311,21 +328,71 @@
 c$OMP END DO NOWAIT
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+
+      !write(*,*) 'vlsplt 326: iq,ijb_x,nqfils(iq)=',iq,ijb_x,nqfils(iq) 
+
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO l=1,llm
+          DO ij=ijb,ije
+           ! On a besoin de q et masse seulement entre ijb et ije. On ne
+           ! les calcule donc que de ijb à ije
+           masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
+
+      !write(*,*) 'vlsplt 360: iq,ijb_x=',iq,ijb_x
+
 c   calcul des tENDances
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
          DO ij=ijb+1,ije
-            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+
-     &      u_mq(ij-1,l)-u_mq(ij,l))
-     &      /new_m
-            masse(ij,l)=new_m
+            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
+     &        u_mq(ij-1,l)-u_mq(ij,l))
+     &        /new_m
+            masse(ij,l,iq)=new_m
          ENDDO
 c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
          DO ij=ijb+iip1-1,ije,iip1
-            q(ij-iim,l)=q(ij,l)
-            masse(ij-iim,l)=masse(ij,l)
-         ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
+            q(ij-iim,l,iq)=q(ij,l,iq)
+            masse(ij-iim,l,iq)=masse(ij,l,iq)
+         ENDDO
+      ENDDO
+c$OMP END DO NOWAIT
+      !write(*,*) 'vlsplt 380: iq,ijb_x=',iq,ijb_x
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio 
+      ! puis on boucle en longitude
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb+1,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+          DO ij=ijb+iip1-1,ije,iip1
+             q(ij-iim,l,iq2)=q(ij,l,iq2)
+          enddo ! DO ij=ijb+iip1-1,ije,iip1
+         enddo !DO l=1,llm
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+      !write(*,*) 'vlsplt 399: iq,ijb_x=',iq,ijb_x
 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
 c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
@@ -336,5 +403,5 @@
 
 
-      SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v)
+      RECURSIVE SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -349,4 +416,5 @@
 c   --------------------------------------------------------------------
       USE parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
       IMPLICIT NONE
 c
@@ -361,7 +429,8 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
       REAL masse_adv_v( ijb_v:ije_v,llm)
-      REAL q(ijb_u:ije_u,llm), dq( ijb_u:ije_u,llm)
+      REAL q(ijb_u:ije_u,llm,nqtot), dq( ijb_u:ije_u,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -392,4 +461,7 @@
       SAVE airej2,airejjm
 c$OMP THREADPRIVATE(airej2,airejjm)
+
+      REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
 c
 c
@@ -401,6 +473,11 @@
       INTEGER ijb,ije
 
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1  
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end
+
       IF(first) THEN
-c         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
+         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
          first=.false.
          do i=2,iip1
@@ -434,5 +511,5 @@
       if (pole_nord) then
         DO i = 1, iim
-          airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
+          airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
         ENDDO
         qpns   = SSUM( iim,  airescb ,1 ) / airej2
@@ -441,11 +518,9 @@
       if (pole_sud) then
         DO i = 1, iim
-          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
+          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
         ENDDO
         qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
       endif
       
-      
-
 c   calcul des pentes aux points v
 
@@ -455,8 +530,12 @@
       if (pole_sud)  ije=ij_end-iip1
       
+      ! on a besoin de q entre ij_begin-2*iip1 et ij_end+2*iip1
+      ! Si pole sud, entre ij_begin-2*iip1 et ij_end
+      ! Si pole Nord, entre ij_begin et ij_end+2*iip1
       DO ij=ijb,ije
-         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
+         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
          adyqv(ij)=abs(dyqv(ij))
       ENDDO
+ 
 
 c   calcul des pentes aux points scalaires
@@ -475,5 +554,5 @@
       IF (pole_nord) THEN
         DO ij=1,iip1
-           dyq(ij,l)=qpns-q(ij+iip1,l)
+           dyq(ij,l)=qpns-q(ij+iip1,l,iq)
         ENDDO
         
@@ -497,5 +576,5 @@
 
         DO ij=1,iip1
-           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
+           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
         ENDDO
 
@@ -633,9 +712,10 @@
        DO ij=ijb,ije
           IF(masse_adv_v(ij,l).gt.0) THEN
-              qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
-     ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
+              qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)*
+     ,                   0.5*(1.-masse_adv_v(ij,l)
+     ,                   /masse(ij+iip1,l,iq))
           ELSE
-              qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
-     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
+              qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
+     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq))
           ENDIF
           qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
@@ -643,4 +723,33 @@
       ENDDO
 c$OMP END DO NOWAIT
+
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      !write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq) 
+
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end
+   
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
+         DO l=1,llm
+         DO ij=ijb,ije
+           masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vly_loc(Ratio,pente_max,masse,qbyv,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
       
       ijb=ij_begin
@@ -652,11 +761,15 @@
       DO l=1,llm
          DO ij=ijb,ije
-            newmasse=masse(ij,l)
-     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
-     
-            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
-     &         /newmasse
-            masse(ij,l)=newmasse
-         ENDDO
+            newmasse=masse(ij,l,iq)
+     &         +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
+
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
+     &         -qbyv(ij-iip1,l))/newmasse
+
+            masse(ij,l,iq)=newmasse
+
+         ENDDO
+
+
 c.-. ancienne version
 c        convpn=SSUM(iim,qbyv(1,l),1)/apoln
@@ -665,12 +778,12 @@
            convpn=SSUM(iim,qbyv(1,l),1)
            convmpn=ssum(iim,masse_adv_v(1,l),1)
-           massepn=ssum(iim,masse(1,l),1)
+           massepn=ssum(iim,masse(1,l,iq),1)
            qpn=0.
            do ij=1,iim
-              qpn=qpn+masse(ij,l)*q(ij,l)
+              qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
            enddo
            qpn=(qpn+convpn)/(massepn+convmpn)
            do ij=1,iip1
-              q(ij,l)=qpn
+              q(ij,l,iq)=qpn
            enddo
          endif
@@ -683,12 +796,12 @@
            convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
            convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
-           masseps=ssum(iim, masse(ip1jm+1,l),1)
+           masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
            qps=0.
            do ij = ip1jm+1,ip1jmp1-1
-              qps=qps+masse(ij,l)*q(ij,l)
+              qps=qps+masse(ij,l,iq)*q(ij,l,iq)
            enddo
            qps=(qps+convps)/(masseps+convmps)
            do ij=ip1jm+1,ip1jmp1
-              q(ij,l)=qps
+              q(ij,l,iq)=qps
            enddo
          endif
@@ -704,5 +817,5 @@
 c        DO ij = 1,iip1
 c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
@@ -714,5 +827,5 @@
 c        DO ij = ip1jm+1,ip1jmp1
 c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c._. fin nouvelle version
@@ -720,4 +833,24 @@
 c$OMP END DO NOWAIT
 
+! retablir les fils en rapport de melange par rapport a l'air:
+      ijb=ij_begin
+      ije=ij_end
+!      if (pole_nord) ijb=ij_begin
+!      if (pole_sud)  ije=ij_end
+
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+
       RETURN
       END
@@ -725,5 +858,5 @@
       
       
-      SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x)
+      RECURSIVE SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -739,4 +872,5 @@
       USE parallel_lmdz
       USE vlz_mod
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi  
       IMPLICIT NONE
 c
@@ -750,7 +884,8 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
-      REAL q(ijb_u:ije_u,llm)
-      REAL w(ijb_u:ije_u,llm+1)
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
+      REAL q(ijb_u:ije_u,llm,nqtot)
+      REAL w(ijb_u:ije_u,llm+1,nqtot)
+      INTEGER iq
 c
 c      Local 
@@ -779,5 +914,9 @@
       LOGICAL,SAVE :: first=.TRUE.
 !$OMP THREADPRIVATE(first)
-      
+
+      !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      ! Ces varibles doivent être déclarées en pointer et en save dans
+      ! vlz_loc si on veut qu'elles soient vues par tous les threads.  
+      INTEGER ifils,iq2 ! CRisi
 
       IF (first) THEN
@@ -787,4 +926,5 @@
 c    sens de W
 
+      !write(*,*) 'vlsplt 926: entree dans vlz_loc, iq=',iq
 #ifdef BIDON
       IF(testcpu) THEN
@@ -799,5 +939,5 @@
       DO l=2,llm
          DO ij=ijb,ije
-            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
+            dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
             adzqw(ij,l)=abs(dzqw(ij,l))
          ENDDO
@@ -842,37 +982,86 @@
 c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
 
+       !write(*,*) 'vlz 982,ijb,ije=',ijb,ije
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1,llm-1
          do  ij = ijb,ije
-          IF(w(ij,l+1).gt.0.) THEN
-             sigw=w(ij,l+1)/masse(ij,l+1)
-             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
+          IF(w(ij,l+1,iq).gt.0.) THEN
+             sigw=w(ij,l+1,iq)/masse(ij,l+1,iq)
+             wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l+1,iq)
+     :           +0.5*(1.-sigw)*dzq(ij,l+1))
           ELSE
-             sigw=w(ij,l+1)/masse(ij,l)
-             wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
+             sigw=w(ij,l+1,iq)/masse(ij,l,iq)
+             wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l,iq)
+     :           -0.5*(1.+sigw)*dzq(ij,l))
           ENDIF
          ENDDO
        ENDDO
-c$OMP END DO NOWAIT
+c$OMP END DO NOWAIT    
+       !write(*,*) 'vlz 1001'    
 
 c$OMP MASTER
        DO ij=ijb,ije
-          wq(ij,llm+1)=0.
-          wq(ij,1)=0.
+          wq(ij,llm+1,iq)=0.
+          wq(ij,1,iq)=0.
        ENDDO
 c$OMP END MASTER
 c$OMP BARRIER
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO l=1,llm
+          DO ij=ijb,ije
+           masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+           !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015
+           w(ij,l,iq2)=wq(ij,l,iq)
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+c$OMP BARRIER
+
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vlz_loc(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi  
+
+! CRisi: On rajoute ici une barrière car on veut être sur que tous les
+! wq soient synchronisés
+
+c$OMP BARRIER
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
          DO ij=ijb,ije
-            newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
+            newmasse=masse(ij,l,iq)+w(ij,l+1,iq)-w(ij,l,iq)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)
+     &         +wq(ij,l+1,iq)-wq(ij,l,iq))
      &         /newmasse
-            masse(ij,l)=newmasse
-         ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
-
+            masse(ij,l,iq)=newmasse
+         ENDDO
+      ENDDO
+c$OMP END DO NOWAIT
+
+      
+! retablir les fils en rapport de melange par rapport a l'air:
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
 
       RETURN
Index: LMDZ5/branches/testing/libf/dyn3dmem/vlspltgen_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/vlspltgen_loc.F	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/vlspltgen_loc.F	(revision 2298)
@@ -27,5 +27,7 @@
       USE Write_Field_loc
       USE VAMPIR
-      USE infotrac, ONLY : nqtot
+      ! CRisi: on rajoute variables utiles d'infotrac  
+      USE infotrac, ONLY : nqtot,nqperes,nqdesc,nqfils,iqfils, 
+     &    ok_iso_verif
       USE vlspltgen_mod
       IMPLICIT NONE
@@ -64,5 +66,5 @@
       REAL ptarg,pdelarg,foeew,zdelta
       REAL tempe(ijb_u:ije_u)
-      INTEGER ijb,ije,iq
+      INTEGER ijb,ije,iq,iq2,ifils
       LOGICAL, SAVE :: firstcall=.TRUE.
 !$OMP THREADPRIVATE(firstcall)
@@ -150,17 +152,21 @@
       ije=ij_end
 
+      DO iq=1,nqtot
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
       DO l=1,llm
          DO ij=ijb,ije
-            mw(ij,l)=w(ij,l) * zzw
+            mw(ij,l,iq)=w(ij,l) * zzw
          ENDDO
       ENDDO
 c$OMP END DO NOWAIT
-
+      ENDDO
+
+      DO iq=1,nqtot  
 c$OMP MASTER
       DO ij=ijb,ije
-         mw(ij,llm+1)=0.
+         mw(ij,llm+1,iq)=0.
       ENDDO
 c$OMP END MASTER
+      ENDDO
 
 c      CALL SCOPY(ijp1llm,q,1,zq,1)
@@ -170,5 +176,5 @@
        ije=ij_end
 
-      DO iq=1,nqtot
+      DO iq=1,nqtot       
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
         DO l=1,llm
@@ -179,5 +185,5 @@
       ENDDO
 
-#ifdef DEBUG_IO    
+#ifdef DEBUG_IO     
        CALL WriteField_u('mu',mu)
        CALL WriteField_v('mv',mv)
@@ -186,11 +192,169 @@
 #endif
 
+      ! verif temporaire
+      ijb=ij_begin
+      ije=ij_end  
+      if (ok_iso_verif) then
+        call check_isotopes(zq,ijb,ije,'vlspltgen_loc 191')
+      endif !if (ok_iso_verif) then   
+
 c$OMP BARRIER           
-      DO iq=1,nqtot
-
+!      DO iq=1,nqtot
+      DO iq=1,nqperes ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
+       !write(*,*) 'vlspltgen 192: iq,iadv=',iq,iadv(iq)
+#ifdef DEBUG_IO    
+       CALL WriteField_u('zq',zq(:,:,iq))
+       CALL WriteField_u('zm',zm(:,:,iq)) 
+#endif
+        if(iadv(iq) == 0) then
+	
+	  cycle 
+	
+	else if (iadv(iq)==10) then
+
+#ifdef _ADV_HALO        
+! CRisi: on ajoute les nombres de fils et tableaux des fils
+! On suppose qu'on ne peut advecter les fils que par le schéma 10.  
+	  call vlx_loc(zq,pente_max,zm,mu,
+     &	             ij_begin,ij_begin+2*iip1-1,iq)
+          call vlx_loc(zq,pente_max,zm,mu,
+     &               ij_end-2*iip1+1,ij_end,iq)
+#else
+	  call vlx_loc(zq,pente_max,zm,mu,
+     &	             ij_begin,ij_end,iq)
+#endif
+
+c$OMP MASTER
+          call VTb(VTHallo)
+c$OMP END MASTER
+          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
+          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
+! CRisi
+          do ifils=1,nqdesc(iq)
+            iq2=iqfils(ifils,iq)
+            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1)
+            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1)
+          enddo
+
+c$OMP MASTER
+          call VTe(VTHallo)
+c$OMP END MASTER
+	else if (iadv(iq)==14) then
+
+#ifdef _ADV_HALO           
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                   qsat,ij_begin,ij_begin+2*iip1-1,iq)
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                   qsat,ij_end-2*iip1+1,ij_end,iq)
+#else
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                   qsat,ij_begin,ij_end,iq)
+#endif
+
+c$OMP MASTER
+          call VTb(VTHallo)
+c$OMP END MASTER
+
+          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
+          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
+          do ifils=1,nqdesc(iq)
+            iq2=iqfils(ifils,iq)
+            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1)
+            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1)
+          enddo
+
+c$OMP MASTER
+          call VTe(VTHallo)
+c$OMP END MASTER 
+        else
+	
+	  stop 'vlspltgen_p : schema non parallelise'
+      
+        endif
+      
+      enddo !DO iq=1,nqperes
+      
+      
+c$OMP BARRIER      
+c$OMP MASTER      
+      call VTb(VTHallo)
+c$OMP END MASTER
+
+      call SendRequest(MyRequest1)
+
+c$OMP MASTER
+      call VTe(VTHallo)
+c$OMP END MASTER       
+c$OMP BARRIER
+
+      ! verif temporaire
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1  
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end  
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 280')
+      endif !if (ok_iso_verif) then
+
+      do iq=1,nqperes
+        !write(*,*) 'vlspltgen 279: iq=',iq
+
+        if(iadv(iq) == 0) then
+	
+	  cycle 
+	
+	else if (iadv(iq)==10) then
+
+#ifdef _ADV_HALLO
+          call vlx_loc(zq,pente_max,zm,mu,
+     &                 ij_begin+2*iip1,ij_end-2*iip1,iq)
+#endif        
+	else if (iadv(iq)==14) then
+#ifdef _ADV_HALLO
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1,iq)
+#endif    
+        else
+	
+	  stop 'vlspltgen_p : schema non parallelise'
+      
+        endif
+      
+      enddo
+c$OMP BARRIER      
+c$OMP MASTER
+      call VTb(VTHallo)
+c$OMP END MASTER
+
+!      call WaitRecvRequest(MyRequest1)
+!      call WaitSendRequest(MyRequest1)
+c$OMP BARRIER
+       call WaitRequest(MyRequest1)
+
+
+c$OMP MASTER
+      call VTe(VTHallo)
+c$OMP END MASTER
+c$OMP BARRIER
+
+      
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326')
+      endif !if (ok_iso_verif) then        
+      if (ok_iso_verif) then
+           ijb=ij_begin-2*iip1
+           ije=ij_end+2*iip1
+           if (pole_nord) ijb=ij_begin
+           if (pole_sud)  ije=ij_end
+           call check_isotopes(zq,ijb,ije,'vlspltgen_loc 336')
+      endif !if (ok_iso_verif) then  
+
+      do iq=1,nqperes
+       !write(*,*) 'vlspltgen 321: iq=',iq
 #ifdef DEBUG_IO    
        CALL WriteField_u('zq',zq(:,:,iq))
        CALL WriteField_u('zm',zm(:,:,iq))
 #endif
+
         if(iadv(iq) == 0) then
 	
@@ -198,47 +362,12 @@
 	
 	else if (iadv(iq)==10) then
-
-#ifdef _ADV_HALO        
-	  call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &	             ij_begin,ij_begin+2*iip1-1)
-          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &               ij_end-2*iip1+1,ij_end)
-#else
-	  call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &	             ij_begin,ij_end)
-#endif
-
-c$OMP MASTER
-          call VTb(VTHallo)
-c$OMP END MASTER
-          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
-          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
-
-c$OMP MASTER
-          call VTe(VTHallo)
-c$OMP END MASTER
+        
+          call vly_loc(zq,pente_max,zm,mv,iq)
+  
 	else if (iadv(iq)==14) then
-
-#ifdef _ADV_HALO           
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                   qsat,ij_begin,ij_begin+2*iip1-1)
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                   qsat,ij_end-2*iip1+1,ij_end)
-#else
-
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                   qsat,ij_begin,ij_end)
-#endif
-
-c$OMP MASTER
-          call VTb(VTHallo)
-c$OMP END MASTER
-
-          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
-          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
-
-c$OMP MASTER
-          call VTe(VTHallo)
-c$OMP END MASTER 
+      
+          call vlyqs_loc(zq,pente_max,zm,mv,
+     &                   qsat,iq)
+ 
         else
 	
@@ -246,78 +375,92 @@
       
         endif
-      
-      enddo
-      
-      
-c$OMP BARRIER      
-c$OMP MASTER      
-      call VTb(VTHallo)
-c$OMP END MASTER
-
-      call SendRequest(MyRequest1)
-
-c$OMP MASTER
-      call VTe(VTHallo)
-c$OMP END MASTER       
-c$OMP BARRIER
-      do iq=1,nqtot
-
-        if(iadv(iq) == 0) then
-	
-	  cycle 
-	
-	else if (iadv(iq)==10) then
-
-#ifdef _ADV_HALLO
-          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                 ij_begin+2*iip1,ij_end-2*iip1)
-#endif        
-	else if (iadv(iq)==14) then
-#ifdef _ADV_HALLO
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1)
-#endif    
-        else
-	
-	  stop 'vlspltgen_p : schema non parallelise'
-      
-        endif
-      
-      enddo
-c$OMP BARRIER      
-c$OMP MASTER
-      call VTb(VTHallo)
-c$OMP END MASTER
-
-!      call WaitRecvRequest(MyRequest1)
-!      call WaitSendRequest(MyRequest1)
-c$OMP BARRIER
-       call WaitRequest(MyRequest1)
-
-
-c$OMP MASTER
-      call VTe(VTHallo)
-c$OMP END MASTER
-c$OMP BARRIER
-
-
-      do iq=1,nqtot
+       
+       enddo
+
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357')
+      endif !if (ok_iso_verif) then
+
+      do iq=1,nqperes
+      !write(*,*) 'vlspltgen 349: iq=',iq
 #ifdef DEBUG_IO    
        CALL WriteField_u('zq',zq(:,:,iq))
        CALL WriteField_u('zm',zm(:,:,iq))
 #endif
+        if(iadv(iq) == 0) then 
+	  
+	  cycle 
+	
+	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
+
+c$OMP BARRIER        
+#ifdef _ADV_HALLO
+          call vlz_loc(zq,pente_max,zm,mw,
+     &               ij_begin,ij_begin+2*iip1-1,iq)
+          call vlz_loc(zq,pente_max,zm,mw,
+     &               ij_end-2*iip1+1,ij_end,iq)
+#else
+          call vlz_loc(zq,pente_max,zm,mw,
+     &               ij_begin,ij_end,iq)
+#endif
+c$OMP BARRIER
+
+c$OMP MASTER
+          call VTb(VTHallo)
+c$OMP END MASTER
+
+          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)
+          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)
+          ! CRisi
+          do ifils=1,nqdesc(iq)
+            iq2=iqfils(ifils,iq)
+            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest2)
+            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest2)
+          enddo      
+c$OMP MASTER
+          call VTe(VTHallo)
+c$OMP END MASTER	
+c$OMP BARRIER
+        else
+	
+	  stop 'vlspltgen_p : schema non parallelise'
+      
+        endif
+      
+      enddo
+c$OMP BARRIER      
+
+c$OMP MASTER        
+      call VTb(VTHallo)
+c$OMP END MASTER
+
+      call SendRequest(MyRequest2)
+
+c$OMP MASTER
+      call VTe(VTHallo)
+c$OMP END MASTER	
+
+
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429')
+      endif !if (ok_iso_verif) then 
+
+c$OMP BARRIER
+      do iq=1,nqperes
+      !write(*,*) 'vlspltgen 409: iq=',iq
+
         if(iadv(iq) == 0) then
-	
+	  
 	  cycle 
 	
-	else if (iadv(iq)==10) then
-        
-          call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv)
-  
-	else if (iadv(iq)==14) then
-      
-          call vlyqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv,
-     &                   qsat)
- 
+	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
+c$OMP BARRIER        
+
+#ifdef _ADV_HALLO
+          call vlz_loc(zq,pente_max,zm,mw,
+     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
+#endif
+
+c$OMP BARRIER        
         else
 	
@@ -325,162 +468,102 @@
       
         endif
-       
-       enddo
-
-
-      do iq=1,nqtot
+      
+      enddo
+      !write(*,*) 'vlspltgen_loc 476'
+
+c$OMP BARRIER
+      !write(*,*) 'vlspltgen_loc 477'
+c$OMP MASTER
+      call VTb(VTHallo)
+c$OMP END MASTER
+
+!      call WaitRecvRequest(MyRequest2)
+!      call WaitSendRequest(MyRequest2)
+c$OMP BARRIER
+       CALL WaitRequest(MyRequest2)
+
+c$OMP MASTER
+      call VTe(VTHallo)
+c$OMP END MASTER
+c$OMP BARRIER
+
+
+      !write(*,*) 'vlspltgen_loc 494'
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461')
+      endif !if (ok_iso_verif) then
+
+      do iq=1,nqperes
+      !write(*,*) 'vlspltgen 449: iq=',iq
 #ifdef DEBUG_IO    
        CALL WriteField_u('zq',zq(:,:,iq))
        CALL WriteField_u('zm',zm(:,:,iq))
 #endif
+        if(iadv(iq) == 0) then
+	
+	  cycle 
+	
+	else if (iadv(iq)==10) then
+        
+          call vly_loc(zq,pente_max,zm,mv,iq)
+  
+	else if (iadv(iq)==14) then
+      
+          call vlyqs_loc(zq,pente_max,zm,mv,
+     &                   qsat,iq)
+ 
+        else
+	
+	  stop 'vlspltgen_p : schema non parallelise'
+      
+        endif
+       
+       enddo !do iq=1,nqperes
+
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493')
+      endif !if (ok_iso_verif) then
+
+      do iq=1,nqperes
+      !write(*,*) 'vlspltgen 477: iq=',iq
+#ifdef DEBUG_IO    
+       CALL WriteField_u('zq',zq(:,:,iq))
+       CALL WriteField_u('zm',zm(:,:,iq))
+#endif
         if(iadv(iq) == 0) then 
 	  
 	  cycle 
 	
-	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
-
-c$OMP BARRIER        
-#ifdef _ADV_HALLO
-          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
-     &               ij_begin,ij_begin+2*iip1-1)
-          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
-     &               ij_end-2*iip1+1,ij_end)
-#else
-          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
-     &               ij_begin,ij_end)
-#endif
-c$OMP BARRIER
-
-c$OMP MASTER
-          call VTb(VTHallo)
-c$OMP END MASTER
-
-          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)
-          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)
-
-c$OMP MASTER
-          call VTe(VTHallo)
-c$OMP END MASTER	
-c$OMP BARRIER
-        else
-	
-	  stop 'vlspltgen_p : schema non parallelise'
-      
-        endif
-      
-      enddo
-c$OMP BARRIER      
-
-c$OMP MASTER        
-      call VTb(VTHallo)
-c$OMP END MASTER
-
-      call SendRequest(MyRequest2)
-
-c$OMP MASTER
-      call VTe(VTHallo)
-c$OMP END MASTER	
-
-c$OMP BARRIER
-      do iq=1,nqtot
-
-        if(iadv(iq) == 0) then
-	  
-	  cycle 
-	
-	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
-c$OMP BARRIER        
-
-#ifdef _ADV_HALLO
-          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
-     &               ij_begin+2*iip1,ij_end-2*iip1)
-#endif
-
-c$OMP BARRIER        
-        else
-	
-	  stop 'vlspltgen_p : schema non parallelise'
-      
-        endif
-      
-      enddo
-
-c$OMP BARRIER
-c$OMP MASTER
-      call VTb(VTHallo)
-c$OMP END MASTER
-
-!      call WaitRecvRequest(MyRequest2)
-!      call WaitSendRequest(MyRequest2)
-c$OMP BARRIER
-       CALL WaitRequest(MyRequest2)
-
-c$OMP MASTER
-      call VTe(VTHallo)
-c$OMP END MASTER
-c$OMP BARRIER
-
-
-      do iq=1,nqtot
-#ifdef DEBUG_IO    
-       CALL WriteField_u('zq',zq(:,:,iq))
-       CALL WriteField_u('zm',zm(:,:,iq))
-#endif
-        if(iadv(iq) == 0) then
-	
-	  cycle 
-	
 	else if (iadv(iq)==10) then
         
-          call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv)
+          call vlx_loc(zq,pente_max,zm,mu,
+     &               ij_begin,ij_end,iq)
   
 	else if (iadv(iq)==14) then
       
-          call vlyqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv,
-     &                   qsat)
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                 qsat, ij_begin,ij_end,iq)
  
         else
 	
-	  stop 'vlspltgen_p : schema non parallelise'
+          stop 'vlspltgen_p : schema non parallelise'
       
         endif
        
-       enddo
-
-
-      do iq=1,nqtot
-#ifdef DEBUG_IO    
-       CALL WriteField_u('zq',zq(:,:,iq))
-       CALL WriteField_u('zm',zm(:,:,iq))
-#endif
-        if(iadv(iq) == 0) then 
-	  
-	  cycle 
-	
-	else if (iadv(iq)==10) then
-        
-          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &               ij_begin,ij_end)
-  
-	else if (iadv(iq)==14) then
-      
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                 qsat, ij_begin,ij_end)
- 
-        else
-	
-          stop 'vlspltgen_p : schema non parallelise'
-      
-        endif
-       
-       enddo
-
+       enddo !do iq=1,nqperes
+
+      !write(*,*) 'vlspltgen 550: apres derniere serie de call vlx'
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521')
+      endif !if (ok_iso_verif) then
      
       ijb=ij_begin
       ije=ij_end
+      !write(*,*) 'vlspltgen_loc 557'
 c$OMP BARRIER      
 
-
+      !write(*,*) 'vlspltgen_loc 559'  
       DO iq=1,nqtot
+       !write(*,*) 'vlspltgen_loc 561, iq=',iq  
 #ifdef DEBUG_IO    
        CALL WriteField_u('zq',zq(:,:,iq))
@@ -495,5 +578,6 @@
            ENDDO
         ENDDO
-c$OMP END DO NOWAIT          
+c$OMP END DO NOWAIT    
+      !write(*,*) 'vlspltgen_loc 575'      
 
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
@@ -504,8 +588,11 @@
         ENDDO
 c$OMP END DO NOWAIT  
-
-      ENDDO
+      !write(*,*) 'vlspltgen_loc 583'  
+      ENDDO !DO iq=1,nqtot
         
-      
+      if (ok_iso_verif) then
+           call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557')
+      endif !if (ok_iso_verif) then
+
 c$OMP BARRIER
 
@@ -516,4 +603,5 @@
 cc$OMP BARRIER
 
+      !write(*,*) 'vlspltgen 597: sortie'  
       RETURN
       END
Index: LMDZ5/branches/testing/libf/dyn3dmem/vlspltgen_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/vlspltgen_mod.F90	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/vlspltgen_mod.F90	(revision 2298)
@@ -2,7 +2,7 @@
 
   REAL,POINTER,SAVE :: qsat(:,:)
-  REAL,POINTER,SAVE :: mu(:,:)
+  REAL,POINTER,SAVE :: mu(:,:) ! CRisi: on ajoute une dimension
   REAL,POINTER,SAVE :: mv(:,:)
-  REAL,POINTER,SAVE :: mw(:,:)
+  REAL,POINTER,SAVE :: mw(:,:,:)
   REAL,POINTER,SAVE :: zm(:,:,:)
   REAL,POINTER,SAVE :: zq(:,:,:)
@@ -25,5 +25,5 @@
     CALL allocate_u(mu,llm,d)
     CALL allocate_v(mv,llm,d)
-    CALL allocate_u(mw,llm+1,d)
+    CALL allocate_u(mw,llm+1,nqtot,d)
     CALL allocate_u(zm,llm,nqtot,d)
     CALL allocate_u(zq,llm,nqtot,d)
Index: LMDZ5/branches/testing/libf/dyn3dmem/vlspltqs_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/vlspltqs_loc.F	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/vlspltqs_loc.F	(revision 2298)
@@ -1,12 +1,13 @@
-      SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
+      SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
 c
 c    ********************************************************************
-c     Shema  d'advection " pseudo amont " .
+c     Shema  d''advection " pseudo amont " .
 c    ********************************************************************
 c
 c   --------------------------------------------------------------------
       USE parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi  
       IMPLICIT NONE
 c
@@ -20,8 +21,9 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
       REAL u_m( ijb_u:ije_u,llm )
-      REAL q(ijb_u:ije_u,llm)
+      REAL q(ijb_u:ije_u,llm,nqtot)
       REAL qsat(ijb_u:ije_u,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -36,4 +38,7 @@
       REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
       REAL u_mq(ijb_u:ije_u,llm)
+      REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
+
 
       REAL      SSUM
@@ -42,4 +47,6 @@
       INTEGER ijb,ije,ijb_x,ije_x
       
+      !write(*,*) 'vlspltqs 58: entree vlxqs_loc, iq,ijb_x=',
+!     &   iq,ijb_x
 
 c   calcul de la pente a droite et a gauche de la maille
@@ -65,7 +72,7 @@
          DO l = 1, llm
             DO ij=ijb,ije-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
 c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
-c              sigu(ij)=u_m(ij,l)/masse(ij,l)
+c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
             ENDDO
             DO ij=ijb+iip1-1,ije,iip1
@@ -120,5 +127,5 @@
          DO l = 1, llm
             DO ij=ijb,ije-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
             ENDDO
             DO ij=ijb+iip1-1,ije,iip1
@@ -179,11 +186,11 @@
       DO l=1,llm
        DO ij=ijb,ije-1
-          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
-     ,                     1.+u_m(ij,l)/masse(ij+1,l),
+          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
+     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
      ,                     u_m(ij,l))
           zdum(ij,l)=0.5*zdum(ij,l)
           u_mq(ij,l)=cvmgp(
-     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
-     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
+     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
+     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
      ,                u_m(ij,l))
           u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
@@ -195,16 +202,16 @@
 c   on cumule le flux correspondant a toutes les mailles dont la masse
 c   au travers de la paroi pENDant le pas de temps.
-c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
+c   le rapport de melange de l''air advecte est min(q_vanleer, Qsat_downwind)
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
        DO ij=ijb,ije-1
           IF (u_m(ij,l).gt.0.) THEN
-             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
+             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
              u_mq(ij,l)=u_m(ij,l)*
-     $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
+     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
           ELSE
-             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
+             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
              u_mq(ij,l)=u_m(ij,l)*
-     $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
+     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
           ENDIF
        ENDDO
@@ -273,5 +280,5 @@
                ENDDO
                niju=iju
-c              PRINT*,'niju,nl',niju,nl(l)
+               !PRINT*,'vlxqs 280: niju,nl',niju,nl(l)
 
 c  traitement des mailles
@@ -285,26 +292,28 @@
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m-masse(ijq,l)
+                     do while(zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
+     &                     *masse(ijq,l,iq) 
+                        zu_m=zu_m-masse(ijq,l,iq)
                         i=mod(i-2+iim,iim)+1
                         ijq=(j-1)*iip1+i
                      ENDDO
 c   ajout de la maille non completement advectee
-                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
-     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
+                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)
+     &                 +0.5*(1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ELSE
                      ijq=ij+1
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(-zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m+masse(ijq,l)
+                     do while(-zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
+     &                   *masse(ijq,l,iq)
+                        zu_m=zu_m+masse(ijq,l,iq)
                         i=mod(i,iim)+1
                         ijq=(j-1)*iip1+i
                      ENDDO
 c   ajout de la maille non completement advectee
-                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
-     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
+                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
+     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ENDIF
                ENDDO
@@ -325,28 +334,77 @@
 c$OMP END DO NOWAIT
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      !write(*,*) 'vlspltqs 336: iq,ijb_x,nqfils(iq)=',
+!     &     iq,ijb_x,nqfils(iq)  
+
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+           masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)      
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqfils(iq)
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         !write(*,*) 'vlxqs 349: on appelle vlx pour iq2=',iq2
+         call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
+
+      !write(*,*) 'vlspltqs 360: iq,ijb_x=',iq,ijb_x   
+
 c   calcul des tendances
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
          DO ij=ijb+1,ije
-            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+
+            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
      &      u_mq(ij-1,l)-u_mq(ij,l))
      &      /new_m
-            masse(ij,l)=new_m
-         ENDDO
-c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
+            masse(ij,l,iq)=new_m
+         ENDDO
+c   Modif Fred 22 03 96 correction d''un bug (les scopy ci-dessous)
          DO ij=ijb+iip1-1,ije,iip1
-            q(ij-iim,l)=q(ij,l)
-            masse(ij-iim,l)=masse(ij,l)
-         ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
+            q(ij-iim,l,iq)=q(ij,l,iq)
+            masse(ij-iim,l,iq)=masse(ij,l,iq)
+         ENDDO
+      ENDDO
+c$OMP END DO NOWAIT
+
+      !write(*,*) 'vlspltqs 380: iq,ijb_x=',iq,ijb_x
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb+1,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+          DO ij=ijb+iip1-1,ije,iip1
+             q(ij-iim,l,iq2)=q(ij,l,iq2)
+          enddo ! DO ij=ijb+iip1-1,ije,iip1
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+      !write(*,*) 'vlspltqs 399: iq,ijb_x=',iq,ijb_x
+
 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
-c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
+c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1,iq),iip1,masse(iip2,1,iq),iip1)
 
 
       RETURN
       END
-      SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat)
+      SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -361,4 +419,5 @@
 c   --------------------------------------------------------------------
       USE parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi  
       IMPLICIT NONE
 c
@@ -373,8 +432,9 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
       REAL masse_adv_v( ijb_v:ije_v,llm)
-      REAL q(ijb_u:ije_u,llm)
+      REAL q(ijb_u:ije_u,llm,nqtot)
       REAL qsat(ijb_u:ije_u,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -386,5 +446,5 @@
       REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v)
       REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
-      REAL qbyv(ijb_v:ije_v,llm)
+      REAL qbyv(ijb_v:ije_v,llm,nqtot)
 
       REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
@@ -402,4 +462,7 @@
 c
 c
+      REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
+
       REAL      SSUM
 
@@ -407,6 +470,18 @@
       INTEGER ijb,ije
 
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1  
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end
+      ij=3525
+      l=3
+      if ((ij.ge.ijb).and.(ij.le.ije)) then
+        !write(*,*) 'vlyqs 480: ij,l,iq,ijb,q(ij,l,:)=',
+!     &             ij,l,iq,ijb,q(ij,l,:)
+      endif  
+
       IF(first) THEN
          PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
+         PRINT*,'vlyqs_loc, iq=',iq
          first=.false.
          do i=2,iip1
@@ -439,5 +514,5 @@
       if (pole_nord) then
         DO i = 1, iim
-          airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
+          airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
         ENDDO
         qpns   = SSUM( iim,  airescb ,1 ) / airej2
@@ -446,5 +521,5 @@
       if (pole_sud) then
         DO i = 1, iim
-          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
+          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
         ENDDO
         qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
@@ -460,5 +535,5 @@
       
       DO ij=ijb,ije
-         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
+         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
          adyqv(ij)=abs(dyqv(ij))
       ENDDO
@@ -482,5 +557,5 @@
 c   calcul des pentes aux poles
         DO ij=1,iip1
-           dyq(ij,l)=qpns-q(ij+iip1,l)
+           dyq(ij,l)=qpns-q(ij+iip1,l,iq)
         ENDDO
 
@@ -513,5 +588,5 @@
 
         DO ij=1,iip1
-           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
+           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
         ENDDO
 
@@ -636,14 +711,45 @@
        DO ij=ijb,ije
          IF( masse_adv_v(ij,l).GT.0. ) THEN
-           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
-     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
+           qbyv(ij,l,iq)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
+     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
+     ,      /masse(ij+iip1,l,iq)))
          ELSE
-              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
-     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
+              qbyv(ij,l,iq)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
+     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
          ENDIF
-          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
+          qbyv(ij,l,iq) = masse_adv_v(ij,l)*qbyv(ij,l,iq)
        ENDDO
       ENDDO
 c$OMP END DO NOWAIT
+
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      !write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
+      
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end  
+
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+           masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)   
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vly_loc(Ratio,pente_max,masse,qbyv,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+        
+! end CRisi
 
       ijb=ij_begin
@@ -655,9 +761,9 @@
       DO l=1,llm
          DO ij=ijb,ije
-            newmasse=masse(ij,l)
+            newmasse=masse(ij,l,iq)
      &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
-     &         /newmasse
-            masse(ij,l)=newmasse
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l,iq)
+     &         -qbyv(ij-iip1,l,iq))/newmasse
+            masse(ij,l,iq)=newmasse
          ENDDO
 c.-. ancienne version
@@ -665,11 +771,11 @@
          IF (pole_nord) THEN
 
-           convpn=SSUM(iim,qbyv(1,l),1)/apoln
+           convpn=SSUM(iim,qbyv(1,l,iq),1)/apoln
            convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
            DO ij = 1,iip1
-              newmasse=masse(ij,l)+convmpn*aire(ij)
-              q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
+              newmasse=masse(ij,l,iq)+convmpn*aire(ij)
+              q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
      &                 newmasse
-              masse(ij,l)=newmasse
+              masse(ij,l,iq)=newmasse
            ENDDO
 	 
@@ -678,11 +784,11 @@
 	 IF (pole_sud) THEN
 	 
-	   convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
+	   convps  = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols
            convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
            DO ij = ip1jm+1,ip1jmp1
-              newmasse=masse(ij,l)+convmps*aire(ij)
-              q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
+              newmasse=masse(ij,l,iq)+convmps*aire(ij)
+              q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
      &                 newmasse
-              masse(ij,l)=newmasse
+              masse(ij,l,iq)=newmasse
            ENDDO
 	 
@@ -691,27 +797,48 @@
 
 c._. nouvelle version
-c        convpn=SSUM(iim,qbyv(1,l),1)
+c        convpn=SSUM(iim,qbyv(1,l,iq),1)
 c        convmpn=ssum(iim,masse_adv_v(1,l),1)
-c        oldmasse=ssum(iim,masse(1,l),1)
+c        oldmasse=ssum(iim,masse(1,l,iq),1)
 c        newmasse=oldmasse+convmpn
-c        newq=(q(1,l)*oldmasse+convpn)/newmasse
+c        newq=(q(1,l,iq)*oldmasse+convpn)/newmasse
 c        newmasse=newmasse/apoln
 c        DO ij = 1,iip1
-c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           q(ij,l,iq)=newq
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
-c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
+c        convps=-SSUM(iim,qbyv(ip1jm-iim,l,iq),1)
 c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
-c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
+c        oldmasse=ssum(iim,masse(ip1jm-iim,l,iq),1)
 c        newmasse=oldmasse+convmps
-c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
+c        newq=(q(ip1jmp1,l,iq)*oldmasse+convps)/newmasse
 c        newmasse=newmasse/apols
 c        DO ij = ip1jm+1,ip1jmp1
-c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           q(ij,l,iq)=newq
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c._. fin nouvelle version
       ENDDO
 c$OMP END DO NOWAIT
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      ijb=ij_begin
+      ije=ij_end
+!      if (pole_nord) ijb=ij_begin+iip1
+!      if (pole_sud)  ije=ij_end-iip1
+  
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+
       RETURN
       END
Index: LMDZ5/branches/testing/libf/dyn3dmem/vlz_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/vlz_mod.F90	(revision 2258)
+++ LMDZ5/branches/testing/libf/dyn3dmem/vlz_mod.F90	(revision 2298)
@@ -1,8 +1,11 @@
 MODULE vlz_mod
 
-  REAL,POINTER,SAVE :: wq(:,:)
+  REAL,POINTER,SAVE :: wq(:,:,:)
   REAL,POINTER,SAVE :: dzq(:,:)
   REAL,POINTER,SAVE :: dzqw(:,:)
   REAL,POINTER,SAVE :: adzqw(:,:)
+  ! CRisi: pour les traceurs:  
+  !REAL,POINTER,SAVE :: masseq(:,:,:)
+  REAL,POINTER,SAVE :: Ratio(:,:,:)
   
 CONTAINS
@@ -18,8 +21,12 @@
     
     d=>distrib_vanleer
-    CALL allocate_u(wq,llm+1,d)
+    CALL allocate_u(wq,llm+1,nqtot,d)
     CALL allocate_u(dzq,llm,d)
     CALL allocate_u(dzqw,llm,d)
     CALL allocate_u(adzqw,llm,d)
+    if (nqdesc_tot.gt.0) then
+    !CALL allocate_u(masseq,llm,nqtot,d)
+    CALL allocate_u(Ratio,llm,nqtot,d)
+    endif !if (nqdesc_tot.gt.0) then
 
   END SUBROUTINE vlz_allocate
@@ -29,4 +36,5 @@
   USE bands
   USE parallel_lmdz
+  USE infotrac
   IMPLICIT NONE
     TYPE(distrib),INTENT(IN) :: dist
@@ -36,4 +44,9 @@
     CALL switch_u(dzqw,distrib_vanleer,dist)
     CALL switch_u(adzqw,distrib_vanleer,dist)
+    ! CRisi:
+    if (nqdesc_tot.gt.0) then    
+    !CALL switch_u(masseq,distrib_vanleer,dist)
+    CALL switch_u(Ratio,distrib_vanleer,dist)
+    endif !if (nqdesc_tot.gt.0) then     
 
   END SUBROUTINE vlz_switch_vanleer  
