Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/advtrac.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/advtrac.F90	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/advtrac.F90	(revision 4368)
@@ -1,408 +1,354 @@
 ! $Id$
 
-SUBROUTINE advtrac(pbaru,pbarv , p,  masse,q,iapptrac,teta, flxw, pk)
-  !     Auteur :  F. Hourdin
-  !
-  !     Modif. P. Le Van     (20/12/97)
-  !            F. Codron     (10/99)
-  !            D. Le Croller (07/2001)
-  !            M.A Filiberti (04/2002)
-  !
-  USE infotrac, ONLY: nqtot, iadv,nqperes,ok_iso_verif
-  USE control_mod, ONLY: iapp_tracvl, day_step
-  USE comconst_mod, ONLY: dtvr
-
-  IMPLICIT NONE
-  !
-  include "dimensions.h"
-  include "paramet.h"
-  include "comdissip.h"
-  include "comgeom2.h"
-  include "description.h"
-  include "iniprint.h"
-
-  !-------------------------------------------------------------------
-  !     Arguments
-  !-------------------------------------------------------------------
-  INTEGER,INTENT(OUT) :: iapptrac
-  REAL,INTENT(IN) :: pbaru(ip1jmp1,llm)
-  REAL,INTENT(IN) :: pbarv(ip1jm,llm)
-  REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
-  REAL,INTENT(IN) :: masse(ip1jmp1,llm)
-  REAL,INTENT(IN) :: p( ip1jmp1,llmp1 )
-  REAL,INTENT(IN) :: teta(ip1jmp1,llm)
-  REAL,INTENT(IN) :: pk(ip1jmp1,llm)
-  REAL,INTENT(OUT) :: flxw(ip1jmp1,llm)
-  !-------------------------------------------------------------------
-  !     Ajout PPM
-  !--------------------------------------------------------
-  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
-  !-------------------------------------------------------------
-  !     Variables locales
-  !-------------------------------------------------------------
-
-  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
-  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
-  REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 
-  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
-  INTEGER iadvtr
-  INTEGER ij,l,iq,iiq
-  REAL zdpmin, zdpmax
-  EXTERNAL  minmax
-  SAVE iadvtr, massem, pbaruc, pbarvc
-  DATA iadvtr/0/
-  !----------------------------------------------------------
-  !     Rajouts pour PPM
-  !----------------------------------------------------------
-  INTEGER indice,n
-  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
-  REAL CFLmaxz,aaa,bbb ! CFL maximum
-  REAL psppm(iim,jjp1) ! pression  au sol
-  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
-  REAL qppm(iim*jjp1,llm,nqtot)
-  REAL fluxwppm(iim,jjp1,llm)
-  REAL apppm(llmp1), bpppm(llmp1)
-  LOGICAL dum,fill
-  DATA fill/.true./
-  DATA dum/.true./
-
-  integer,save :: countcfl=0
-  real cflx(ip1jmp1,llm)
-  real cfly(ip1jm,llm)
-  real cflz(ip1jmp1,llm)
-  real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
-
-  IF(iadvtr.EQ.0) THEN
-     pbaruc(:,:)=0
-     pbarvc(:,:)=0
-  ENDIF
-
-  !   accumulation des flux de masse horizontaux
-  DO l=1,llm
-     DO ij = 1,ip1jmp1
-        pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
-     ENDDO
-     DO ij = 1,ip1jm
-        pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
-     ENDDO
-  ENDDO
-
-  !   selection de la masse instantannee des mailles avant le transport.
-  IF(iadvtr.EQ.0) THEN
-
+#define DEBUG_IO
+#undef DEBUG_IO
+SUBROUTINE advtrac(pbaru, pbarv, p, masse,q,iapptrac,teta, flxw, pk)
+   !     Auteur :  F. Hourdin
+   !
+   !     Modif. P. Le Van     (20/12/97)
+   !            F. Codron     (10/99)
+   !            D. Le Croller (07/2001)
+   !            M.A Filiberti (04/2002)
+   !
+   USE infotrac,     ONLY: nqtot, tracers, isoCheck
+   USE control_mod,  ONLY: iapp_tracvl, day_step
+   USE comconst_mod, ONLY: dtvr
+
+   IMPLICIT NONE
+   !
+   include "dimensions.h"
+   include "paramet.h"
+   include "comdissip.h"
+   include "comgeom2.h"
+   include "description.h"
+   include "iniprint.h"
+
+   !---------------------------------------------------------------------------
+   !     Arguments
+   !---------------------------------------------------------------------------
+   INTEGER, INTENT(OUT) :: iapptrac
+   REAL, INTENT(IN) :: pbaru(ip1jmp1,llm)
+   REAL, INTENT(IN) :: pbarv(ip1jm,  llm)
+   REAL, INTENT(INOUT) ::  q(ip1jmp1,llm,nqtot)
+   REAL, INTENT(IN) :: masse(ip1jmp1,llm)
+   REAL, INTENT(IN) ::     p(ip1jmp1,llmp1 )
+   REAL, INTENT(IN) ::  teta(ip1jmp1,llm)
+   REAL, INTENT(IN) ::    pk(ip1jmp1,llm)
+   REAL, INTENT(OUT) :: flxw(ip1jmp1,llm)
+   !---------------------------------------------------------------------------
+   !     Ajout PPM
+   !---------------------------------------------------------------------------
+   REAL :: massebx(ip1jmp1,llm), masseby(ip1jm,llm)
+   !---------------------------------------------------------------------------
+   !     Variables locales
+   !---------------------------------------------------------------------------
+   INTEGER :: ij, l, iq, iadv
+!   REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu
+   REAL :: zdp(ip1jmp1), zdpmin, zdpmax
+   INTEGER, SAVE :: iadvtr=0
+   REAL, DIMENSION(ip1jmp1,llm) :: pbaruc, pbarug, massem, wg
+   REAL, DIMENSION(ip1jm,  llm) :: pbarvc, pbarvg
+   EXTERNAL  minmax
+   SAVE massem, pbaruc, pbarvc
+   !---------------------------------------------------------------------------
+   !     Rajouts pour PPM
+   !---------------------------------------------------------------------------
+   INTEGER indice, n
+   REAL :: dtbon                       ! Pas de temps adaptatif pour que CFL<1
+   REAL :: CFLmaxz, aaa, bbb           ! CFL maximum
+   REAL, DIMENSION(iim,jjp1,llm) :: unatppm, vnatppm, fluxwppm
+   REAL ::    qppm(iim*jjp1,llm,nqtot)
+   REAL ::   psppm(iim,jjp1)           ! pression  au sol
+   REAL, DIMENSION(llmp1) :: apppm, bpppm
+   LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE.
+
+   INTEGER, SAVE :: countcfl=0
+   REAL, DIMENSION(ip1jmp1,llm) :: cflx, cflz
+   REAL, DIMENSION(ip1jm  ,llm) :: cfly
+   REAL, DIMENSION(llm), SAVE :: cflxmax, cflymax, cflzmax
+
+   IF(iadvtr == 0) THEN
+      pbaruc(:,:)=0
+      pbarvc(:,:)=0
+   END IF
+
+   !--- Accumulation des flux de masse horizontaux
+   DO l=1,llm
+      DO ij = 1,ip1jmp1
+         pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
+      END DO
+      DO ij = 1,ip1jm
+         pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
+      END DO
+   END DO
+
+   !--- Selection de la masse instantannee des mailles avant le transport.
+   IF(iadvtr == 0) THEN
      CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
-     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
-     !
-  ENDIF
-
-  iadvtr   = iadvtr+1
-  iapptrac = iadvtr
-
-
-  !   Test pour savoir si on advecte a ce pas de temps
-  IF ( iadvtr.EQ.iapp_tracvl ) THEN
-
-     !c   ..  Modif P.Le Van  ( 20/12/97 )  ....
-     !c
-
-     !   traitement des flux de masse avant advection.
-     !     1. calcul de w
-     !     2. groupement des mailles pres du pole.
-
-     CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
-
-     ! ... Flux de masse diaganostiques traceurs
-     flxw = wg / REAL(iapp_tracvl)
-
-     !  test sur l'eventuelle creation de valeurs negatives de la masse
-     DO l=1,llm-1
-        DO ij = iip2+1,ip1jm
-           zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l) &
-                - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
-                +       wg(ij,l+1)  - wg(ij,l)
-        ENDDO
-        CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
-        DO ij = iip2,ip1jm
-           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 
-        ENDDO
-
-
-        CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
-
-        IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
-           PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin, &
-                '   MAX:', zdpmax
-        ENDIF
-
-     ENDDO
-
-
-     !-------------------------------------------------------------------
-     ! Calcul des criteres CFL en X, Y et Z
-     !-------------------------------------------------------------------
-
-     if (countcfl == 0. ) then
-        cflxmax(:)=0.
-        cflymax(:)=0.
-        cflzmax(:)=0.
-     endif
-
-     countcfl=countcfl+iapp_tracvl
-     cflx(:,:)=0.
-     cfly(:,:)=0.
-     cflz(:,:)=0.
-     do l=1,llm
-        do ij=iip2,ip1jm-1
-           if (pbarug(ij,l)>=0.) then
-              cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
-           else
-              cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
-           endif
-        enddo
-     enddo
-     do l=1,llm
-        do ij=iip2,ip1jm-1,iip1
-           cflx(ij+iip1,l)=cflx(ij,l)
-        enddo
-     enddo
-
-     do l=1,llm
-        do ij=1,ip1jm
-           if (pbarvg(ij,l)>=0.) then
-              cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
-           else
-              cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
-           endif
-        enddo
-     enddo
-
-     do l=2,llm
-        do ij=1,ip1jm
-           if (wg(ij,l)>=0.) then
-              cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
-           else
-              cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
-           endif
-        enddo
-     enddo
-
-     do l=1,llm
-        cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
-        cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
-        cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
-     enddo
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     ! Par defaut, on sort le diagnostic des CFL tous les jours.
-     ! Si on veut le sortir a chaque pas d'advection en cas de plantage 
-     !     if (countcfl==iapp_tracvl) then
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     if (countcfl==day_step) then
-        do l=1,llm
-           write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), &
-                cflzmax(l)
-        enddo
-        countcfl=0
-     endif
-
-     !-------------------------------------------------------------------
-     !   Advection proprement dite (Modification Le Croller (07/2001)
-     !-------------------------------------------------------------------
-
-     !----------------------------------------------------
-     !        Calcul des moyennes basées sur la masse
-     !----------------------------------------------------
-     call massbar(massem,massebx,masseby)          
-
-     !-----------------------------------------------------------
-     !     Appel des sous programmes d'advection
-     !-----------------------------------------------------------
-
-     if (ok_iso_verif) then
-           write(*,*) 'advtrac 227'
-           call check_isotopes_seq(q,ip1jmp1,'advtrac 162')
-     endif !if (ok_iso_verif) then
-
-     do iq=1,nqperes
-        !        call clock(t_initial)
-        if(iadv(iq) == 0) cycle 
-        !   ----------------------------------------------------------------
-        !   Schema de Van Leer I MUSCL
-        !   ----------------------------------------------------------------
-        if(iadv(iq).eq.10) THEN
-           ! CRisi: on fait passer tout q pour avoir acces aux fils
-           
-           !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
-           call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
-           
-           !   ----------------------------------------------------------------
-           !   Schema "pseudo amont" + test sur humidite specifique
-           !    pour la vapeur d'eau. F. Codron
-           !   ----------------------------------------------------------------
-        else if(iadv(iq).eq.14) then
-           !
-           !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
-           CALL vlspltqs( q, 2., massem, wg , &
-                pbarug,pbarvg,dtvr,p,pk,teta,iq)
-           
-           !   ----------------------------------------------------------------
-           !   Schema de Frederic Hourdin
-           !   ----------------------------------------------------------------
-        else if(iadv(iq).eq.12) then
-           !            Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
-                   dtvr,'n=',n
-           endif
-           do indice=1,n
-              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
-           end do
-        else if(iadv(iq).eq.13) then
-           !            Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
-                   dtvr,'n=',n
-           endif
-           do indice=1,n
-              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
-           end do
-           !   ----------------------------------------------------------------
-           !   Schema de pente SLOPES
-           !   ----------------------------------------------------------------
-        else if (iadv(iq).eq.20) then
-           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
-
-           !   ----------------------------------------------------------------
-           !   Schema de Prather
-           !   ----------------------------------------------------------------
-        else if (iadv(iq).eq.30) then
-           !            Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
-                   dtvr,'n=',n
-           endif
-           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
-                n,dtbon)
-
-           !   ----------------------------------------------------------------
-           !   Schemas PPM Lin et Rood
-           !   ----------------------------------------------------------------
-        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
-             iadv(iq).LE.18)) then
-
-           !        Test sur le flux horizontal
-           !        Pas de temps adaptatif
-           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
-           if (n.GT.1) then
-              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
-                   dtvr,'n=',n
-           endif
-           !        Test sur le flux vertical
-           CFLmaxz=0.
-           do l=2,llm
-              do ij=iip2,ip1jm
-                 aaa=wg(ij,l)*dtvr/massem(ij,l)
-                 CFLmaxz=max(CFLmaxz,aaa)
-                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
-                 CFLmaxz=max(CFLmaxz,bbb)
-              enddo
-           enddo
-           if (CFLmaxz.GE.1) then
-              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
-           endif
-
-           !-----------------------------------------------------------
-           !        Ss-prg interface LMDZ.4->PPM3d
-           !-----------------------------------------------------------
-
-           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
-                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
-                unatppm,vnatppm,psppm)
-
-           do indice=1,n
-              !----------------------------------------------------------------
-              !                         VL (version PPM) horiz. et PPM vert.
-              !----------------------------------------------------------------
-              if (iadv(iq).eq.11) then
-                 !                  Ss-prg PPM3d de Lin
-                 call ppm3d(1,qppm(1,1,iq), &
-                      psppm,psppm, &
-                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
-                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
-                      fill,dum,220.)
-
-                 !-------------------------------------------------------------
-                 !                           Monotonic PPM
-                 !-------------------------------------------------------------
-              else if (iadv(iq).eq.16) then
-                 !                  Ss-prg PPM3d de Lin
-                 call ppm3d(1,qppm(1,1,iq), &
-                      psppm,psppm, &
-                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
-                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
-                      fill,dum,220.)
-                 !-------------------------------------------------------------
-
-                 !-------------------------------------------------------------
-                 !                           Semi Monotonic PPM
-                 !-------------------------------------------------------------
-              else if (iadv(iq).eq.17) then
-                 !                  Ss-prg PPM3d de Lin
-                 call ppm3d(1,qppm(1,1,iq), &
-                      psppm,psppm, &
-                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
-                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
-                      fill,dum,220.)
-                 !-------------------------------------------------------------
-
-                 !-------------------------------------------------------------
-                 !                         Positive Definite PPM
-                 !-------------------------------------------------------------
-              else if (iadv(iq).eq.18) then
-                 !                  Ss-prg PPM3d de Lin
-                 call ppm3d(1,qppm(1,1,iq), &
-                      psppm,psppm, &
-                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
-                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
-                      fill,dum,220.)
-                 !-------------------------------------------------------------
-              endif
-           enddo
-           !-----------------------------------------------------------------
-           !               Ss-prg interface PPM3d-LMDZ.4
-           !-----------------------------------------------------------------
-           call interpost(q(1,1,iq),qppm(1,1,iq))
-        endif
-        !----------------------------------------------------------------------
-
-        !-----------------------------------------------------------------
-        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
-        ! et Nord j=1
-        !-----------------------------------------------------------------
-
-        !                  call traceurpole(q(1,1,iq),massem)
-
-        ! calcul du temps cpu pour un schema donne
-
-        !                  call clock(t_final)
-        !ym                  tps_cpu=t_final-t_initial
-        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
-
-     end DO
-
-     if (ok_iso_verif) then
-           write(*,*) 'advtrac 402'
-           call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
-     endif !if (ok_iso_verif) then
-
-     !------------------------------------------------------------------
-     !   on reinitialise a zero les flux de masse cumules
-     !---------------------------------------------------
-     iadvtr=0
-
-  ENDIF ! if iadvtr.EQ.iapp_tracvl
+   ! CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
+   END IF
+
+   iadvtr   = iadvtr+1
+   iapptrac = iadvtr
+
+   !--- Test pour savoir si on advecte a ce pas de temps
+   IF(iadvtr /= iapp_tracvl) RETURN
+
+   !   ..  Modif P.Le Van  ( 20/12/97 )  ....
+   !
+   !   traitement des flux de masse avant advection.
+   !       1. calcul de w
+   !       2. groupement des mailles pres du pole.
+
+   CALL groupe(massem, pbaruc, pbarvc, pbarug, pbarvg, wg)
+
+   !--- Flux de masse diaganostiques traceurs
+   flxw = wg / REAL(iapp_tracvl)
+
+   !--- Test sur l'eventuelle creation de valeurs negatives de la masse
+   DO l=1,llm-1
+      DO ij = iip2+1,ip1jm
+         zdp(ij) = pbarug(ij-1,l)    - pbarug(ij,l) &
+                 - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
+                 +     wg(ij,l+1)    -     wg(ij,l)
+      END DO
+! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
+      CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
+      DO ij = iip2,ip1jm
+         zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
+      END DO
+
+      CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
+
+      IF(MAX(ABS(zdpmin),ABS(zdpmax)) > 0.5) &
+         WRITE(*,*)'WARNING DP/P l=',l,'  MIN:',zdpmin,' MAX:', zdpmax
+
+   END DO
+
+   !-------------------------------------------------------------------------
+   ! Calcul des criteres CFL en X, Y et Z
+   !-------------------------------------------------------------------------
+   IF(countcfl == 0. ) then
+      cflxmax(:)=0.
+      cflymax(:)=0.
+      cflzmax(:)=0.
+   END IF
+
+   countcfl=countcfl+iapp_tracvl
+   cflx(:,:)=0.
+   cfly(:,:)=0.
+   cflz(:,:)=0.
+   DO l=1,llm
+      DO ij=iip2,ip1jm-1
+         IF(pbarug(ij,l)>=0.) then
+            cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
+         ELSE
+            cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
+         END IF
+      END DO
+   END DO
+
+   DO l=1,llm
+      DO ij=iip2,ip1jm-1,iip1
+         cflx(ij+iip1,l)=cflx(ij,l)
+      END DO
+   END DO
+
+   DO l=1,llm
+      DO ij=1,ip1jm
+         IF(pbarvg(ij,l)>=0.) then
+            cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
+         ELSE
+            cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
+         END IF
+      END DO
+   END DO
+
+   DO l=2,llm
+      DO ij=1,ip1jm
+         IF(wg(ij,l) >= 0.) THEN
+            cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
+         ELSE
+            cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
+         END IF
+      END DO
+   END DO
+
+   DO l=1,llm
+      cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
+      cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
+      cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
+   END DO
+
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Par defaut, on sort le diagnostic des CFL tous les jours.
+   ! Si on veut le sortir a chaque pas d'advection en cas de plantage 
+   !       IF(countcfl==iapp_tracvl) then
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   IF(countcfl==day_step) then
+      DO l=1,llm
+         WRITE(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), cflzmax(l)
+      END DO
+      countcfl=0
+   END IF
+
+   !---------------------------------------------------------------------------
+   !   Advection proprement dite (Modification Le Croller (07/2001)
+   !---------------------------------------------------------------------------
+
+   !---------------------------------------------------------------------------
+   !   Calcul des moyennes basees sur la masse
+   !---------------------------------------------------------------------------
+   CALL massbar(massem,massebx,masseby)
+
+#ifdef DEBUG_IO    
+   CALL WriteField_u('massem',massem)
+   CALL WriteField_u('wg',wg)
+   CALL WriteField_u('pbarug',pbarug)
+   CALL WriteField_v('pbarvg',pbarvg)
+   CALL WriteField_u('p_tmp',p)
+   CALL WriteField_u('pk_tmp',pk)
+   CALL WriteField_u('teta_tmp',teta)
+   DO iq=1,nqtot
+      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
+   END DO
+#endif
+
+   IF(isoCheck) WRITE(*,*) 'advtrac 227'
+   CALL check_isotopes_seq(q,ip1jmp1,'advtrac 162')
+
+   !-------------------------------------------------------------------------
+   !       Appel des sous programmes d'advection
+   !-------------------------------------------------------------------------
+   DO iq = 1, nqtot
+!     CALL clock(t_initial)
+      IF(tracers(iq)%parent /= 'air') CYCLE
+      iadv = tracers(iq)%iadv
+      !-----------------------------------------------------------------------
+      SELECT CASE(iadv)
+      !-----------------------------------------------------------------------
+         CASE(0); CYCLE
+         !--------------------------------------------------------------------
+         CASE(10)  !--- Schema de Van Leer I MUSCL
+         !--------------------------------------------------------------------
+!           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
+            CALL vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
+
+         !--------------------------------------------------------------------
+         CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
+                   !--- pour la vapeur d'eau. F. Codron
+         !--------------------------------------------------------------------
+!           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
+            CALL vlspltqs(q,2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta,iq)
+
+         !--------------------------------------------------------------------
+         CASE(12)  !--- Schema de Frederic Hourdin
+         !--------------------------------------------------------------------
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            DO indice=1,n
+              CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
+            END DO
+
+         !--------------------------------------------------------------------
+         CASE(13)  !--- Pas de temps adaptatif
+         !--------------------------------------------------------------------
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            DO indice=1,n
+               CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
+            END DO
+
+         !--------------------------------------------------------------------
+         CASE(20)  !--- Schema de pente SLOPES
+         !--------------------------------------------------------------------
+            CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
+
+         !--------------------------------------------------------------------
+         CASE(30)  !--- Schema de Prather
+         !--------------------------------------------------------------------
+            ! Pas de temps adaptatif
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon)
+
+         !--------------------------------------------------------------------
+         CASE(11,16,17,18)   !--- Schemas PPM Lin et Rood
+         !--------------------------------------------------------------------
+            ! Test sur le flux horizontal
+            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
+            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
+            ! Test sur le flux vertical
+            CFLmaxz=0.
+            DO l=2,llm
+               DO ij=iip2,ip1jm
+                  aaa=wg(ij,l)*dtvr/massem(ij,l)
+                  CFLmaxz=max(CFLmaxz,aaa)
+                  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
+                  CFLmaxz=max(CFLmaxz,bbb)
+               END DO
+            END DO
+            IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
+            !----------------------------------------------------------------
+            !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
+            !----------------------------------------------------------------
+            CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
+                 apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
+                 unatppm,vnatppm,psppm)
+
+            !----------------------------------------------------------------
+            DO indice=1,n     !--- VL (version PPM) horiz. et PPM vert.
+            !----------------------------------------------------------------
+               SELECT CASE(iadv)
+                  !----------------------------------------------------------
+                  CASE(11)
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+                                2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
+                  !----------------------------------------------------------
+                  CASE(16) !--- Monotonic PPM
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+                                3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
+                  !----------------------------------------------------------
+                  CASE(17) !--- Semi monotonic PPM
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+                                4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.)
+                  !----------------------------------------------------------
+                  CASE(18) !--- Positive Definite PPM
+                  !----------------------------------------------------------
+                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
+                                5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
+               END SELECT
+            !----------------------------------------------------------------
+            END DO
+            !----------------------------------------------------------------
+            !     Ss-prg interface PPM3d-LMDZ.4
+            !----------------------------------------------------------------
+            CALL interpost(q(1,1,iq),qppm(1,1,iq))
+      !----------------------------------------------------------------------
+      END SELECT
+      !----------------------------------------------------------------------
+
+      !----------------------------------------------------------------------
+      ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
+      !----------------------------------------------------------------------
+      !  CALL traceurpole(q(1,1,iq),massem)
+
+      !--- Calcul du temps cpu pour un schema donne
+      !  CALL clock(t_final)
+      !ym  tps_cpu=t_final-t_initial
+      !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
+
+   END DO
+
+   IF(isoCheck) WRITE(*,*) 'advtrac 402'
+   CALL check_isotopes_seq(q,ip1jmp1,'advtrac 397')
+
+   !-------------------------------------------------------------------------
+   !   on reinitialise a zero les flux de masse cumules
+   !-------------------------------------------------------------------------
+   iadvtr=0
 
 END SUBROUTINE advtrac
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/check_isotopes.F
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/check_isotopes.F	(revision 4013)
+++ 	(revision )
@@ -1,201 +1,0 @@
-        subroutine check_isotopes_seq(q,ip1jmp1,err_msg)
-        USE infotrac
-        implicit none
-
-#include "dimensions.h"
-
-        ! inputs
-        integer ip1jmp1
-        real q(ip1jmp1,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)
-            do k=1,llm
-              DO i = 1,ip1jmp1
-                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 = 1,ip1jmp1
-            enddo !do k=1,llm
-          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)
-            do k=1,llm
-            DO i = 1,ip1jmp1  
-              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=',ixt,phase
-                  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 = 1,ip1jmp1
-            enddo !do k=1,llm
-          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)
-            do k=1,llm
-            DO i = 1,ip1jmp1
-            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 = 1,ip1jmp1
-            enddo !do k=1,llm
-          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)
-            do k=1,llm
-            DO i = 1,ip1jmp1
-            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 = 1,ip1jmp1
-            enddo !do k=1,llm
-          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)
-               do k=1,llm
-                DO i = 1,ip1jmp1
-                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 = 1,ip1jmp1
-                enddo  ! do k=1,llm
-              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)
-              do k=1,llm
-                DO i = 1,ip1jmp1
-                   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 = 1,ip1jmp1
-              enddo !do k=1,llm
-           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: LMDZ6/branches/Ocean_skin/libf/dyn3d/check_isotopes.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/check_isotopes.F90	(revision 4368)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/check_isotopes.F90	(revision 4368)
@@ -0,0 +1,174 @@
+SUBROUTINE check_isotopes_seq(q, ip1jmp1, err_msg)
+   USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
+   USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
+                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, getKey
+   IMPLICIT NONE
+   include "dimensions.h"
+   REAL,             INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
+   INTEGER,          INTENT(IN)    :: ip1jmp1
+   CHARACTER(LEN=*), INTENT(IN)    :: err_msg    !--- Error message to display
+   CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
+   INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
+   INTEGER, ALLOCATABLE ::   ix(:)
+   REAL,    ALLOCATABLE :: tnat(:)
+   REAL    :: xtractot, xiiso, deltaD, q1, q2
+   REAL, PARAMETER :: borne     = 1e19,  &
+                      errmax    = 1e-8,  &       !--- Max. absolute error
+                      errmaxrel = 1e-3,  &       !--- Max. relative error
+                      qmin      = 1e-11, &
+                      deltaDmax =1000.0, &
+                      deltaDmin =-999.0, &
+                      ridicule  = 1e-12
+   INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, &
+                             iso_O17, iso_HTO
+   LOGICAL, SAVE :: first=.TRUE.
+
+   modname='check_isotopes'
+   IF(.NOT.isoCheck)    RETURN                   !--- No need to check => finished
+   IF(isoSelect('H2O')) RETURN                   !--- No H2O isotopes group found
+   IF(niso == 0)        RETURN                   !--- No isotopes => finished
+   IF(first) THEN
+      iso_eau = strIdx(isoName,'H2[16]O')
+      iso_HDO = strIdx(isoName,'H[2]HO')
+      iso_O18 = strIdx(isoName,'H2[18]O')
+      iso_O17 = strIdx(isoName,'H2[17]O')
+      iso_HTO = strIdx(isoName,'H[3]HO')
+      IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1)
+      first = .FALSE.
+   END IF
+   CALL msg('31: err_msg='//TRIM(err_msg), modname)
+
+   !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
+   modname = 'check_isotopes:iso_verif_noNaN'
+   DO ixt = 1, ntiso
+      DO ipha = 1, nphas
+         iq = iqIsoPha(ixt,ipha)
+         DO k = 1, llm
+            DO i = 1, ip1jmp1
+               IF(ABS(q(i,k,iq)) < borne) CYCLE
+               WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
+               CALL msg(msg1, modname)
+               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+            END DO
+         END DO
+      END DO
+   END DO
+
+   !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
+   modname = 'check_isotopes:iso_verif_egalite'
+   ixt = iso_eau
+   IF(ixt /= 0) THEN
+      DO ipha = 1, nphas
+         iq = iqIsoPha(ixt,ipha)
+         iqpar = tracers(iq)%iqParent
+         DO k = 1, llm
+            DO i = 1, ip1jmp1
+               q1 = q(i,k,iqpar)
+               q2 = q(i,k,iq)
+!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
+!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
+!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
+!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
+               IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) <= errmaxrel) THEN
+                  q(i,k,iq) = q1                 !--- Bidouille pour convergence
+!                 q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2
+                  CYCLE
+               END IF
+               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
+               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
+               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
+               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
+               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+            END DO
+         END DO
+      END DO
+   END IF
+
+   !--- CHECK DELTA ANOMALIES
+   modname = 'check_isotopes:iso_verif_aberrant'
+   ix = [ iso_HDO  ,   iso_O18 ]
+   nm = ['deltaD  ', 'deltaO18']
+   DO iiso = 1, SIZE(ix)
+      ixt = ix(iiso)
+      IF(ixt  == 0) CYCLE
+      DO ipha = 1, nphas
+         iq = iqIsoPha(ixt,ipha)
+         iqpar = tracers(iq)%iqParent
+         DO k = 1, llm
+            DO i = 1, ip1jmp1
+               q1 = q(i,k,iqpar)
+               q2 = q(i,k,iq)
+!--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
+!    This would be probably required to sum from smallest to highest concentrations ; the corresponding
+!    indices vector can be computed once only (in the initializations stage), using mean concentrations.
+!              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
+               IF(q2 <= qmin) CYCLE
+               deltaD = (q2/q1/tnat(ixt)-1.)*1000.
+               IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
+               CALL msg('ixt, iq = '//TRIM(strStack(int2str([ixt,iq]))), modname)
+               msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
+               CALL msg(TRIM(tracers(iqpar)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
+               CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
+               CALL msg(TRIM(nm(iiso))//TRIM(real2str(deltaD)), modname)
+               CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+            END DO
+         END DO
+      END DO
+   END DO
+
+   IF(nzone == 0) RETURN
+
+   !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
+   modname = 'check_isotopes:iso_verif_aberrant'
+   IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN
+      DO izon = 1, nzone
+         ixt  = itZonIso(izon, iso_HDO)
+         ieau = itZonIso(izon, iso_eau)
+         DO ipha = 1, nphas
+            iq    = iqIsoPha(ixt,  ipha)
+            iqeau = iqIsoPha(ieau, ipha)
+            DO k = 1, llm
+               DO i = 1, ip1jmp1
+                  q1 = q(i,k,iqeau)
+                  q2 = q(i,k,iq)
+                  IF(q2<=qmin) CYCLE
+                  deltaD = (q2/q1/tnat(iso_HDO)-1.)*1000.
+                  IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
+                  CALL msg('izon, ipha = '//TRIM(strStack(int2str([izon, ipha]))), modname)
+                  CALL msg( 'ixt, ieau = '//TRIM(strStack(int2str([ ixt, ieau]))), modname)
+                  msg1 = '('//TRIM(strStack(int2str([i,k])))//')'
+                  CALL msg(TRIM(tracers(iqeau)%name)//TRIM(msg1)//' = '//TRIM(real2str(q1)), modname)
+                  CALL msg(TRIM(tracers(iq   )%name)//TRIM(msg1)//' = '//TRIM(real2str(q2)), modname)
+                  CALL msg('deltaD = '//TRIM(real2str(deltaD)), modname)
+                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+               END DO
+            END DO
+         END DO
+      END DO
+   END IF
+
+   !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
+   DO iiso = 1, niso
+      DO ipha = 1, nphas
+         iq = iqIsoPha(iiso, ipha)
+         DO k = 1, llm
+            DO i = 1, ip1jmp1
+               xiiso = q(i,k,iq)
+               xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone,iiso), ipha)))
+               IF(ABS(xtractot-xiiso) > errmax .AND. ABS(xtractot-xiiso)/MAX(MAX(ABS(xtractot),ABS(xiiso)),1e-18) > errmaxrel) THEN
+                  CALL msg('Error in iso_verif_aberrant trac: '//TRIM(err_msg))
+                  CALL msg('iiso, ipha = '//TRIM(strStack(int2str([iiso, ipha]))), modname)
+                  CALL msg('q('//TRIM(strStack(int2str([i,k])))//',:) = '//TRIM(strStack(real2str(q(i,k,:)))), modname)
+                  CALL abort_gcm(modname, 'Error with isotopes: '//TRIM(err_msg), 1)
+               END IF
+               IF(ABS(xtractot) <= ridicule) CYCLE
+               DO izon = 1, nzone
+                  q(i,k,iq) = q(i,k,iq) / xtractot * xiiso !--- Bidouille pour convergence
+               END DO
+            END DO
+         END DO
+      END DO
+   END DO
+
+END SUBROUTINE check_isotopes_seq
+
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/conf_gcm.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/conf_gcm.F90	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/conf_gcm.F90	(revision 4368)
@@ -116,13 +116,13 @@
   CALL getin('calend', calend)
 ! initialize year_len for aquaplanets and 1D
-      if (calend == 'earth_360d') then
-        year_len=360
-      else if (calend == 'earth_365d') then
-        year_len=365
-      else if (calend == 'earth_366d') then
-        year_len=366
-      else 
-        year_len=1
-      endif
+  IF (calend == 'earth_360d') THEN
+     year_len=360
+  ELSE IF (calend == 'earth_365d') THEN
+     year_len=365
+  ELSE IF (calend == 'earth_366d') THEN
+     year_len=366
+  ELSE 
+     year_len=1
+  ENDIF
 
   !Config  Key  = dayref
@@ -315,5 +315,4 @@
   CALL getin('ngroup',ngroup)
 
-
   ! top_bound sponge: only active if ok_strato=.true. and iflag_top_bound!=0
   !                   iflag_top_bound=0 for no sponge
@@ -322,4 +321,6 @@
   iflag_top_bound=1
   CALL getin('iflag_top_bound',iflag_top_bound)
+  IF (iflag_top_bound < 0 .or. iflag_top_bound > 2) &
+       call abort_gcm("conf_gcm", 'iflag_top_bound must be 0, 1 or 2', 1)
 
   ! mode_top_bound : fields towards which sponge relaxation will be done:
@@ -394,5 +395,5 @@
   !     .........   (  modif  le 17/04/96 )   .........
 
-  test_etatinit: IF (.not. etatinit) then
+  test_etatinit: IF (.not. etatinit) THEN
      !Config  Key  = clon
      !Config  Desc = centre du zoom, longitude
@@ -598,14 +599,4 @@
      type_trac = 'lmdz'
      CALL getin('type_trac',type_trac)
-
-     !Config  Key  = config_inca
-     !Config  Desc = Choix de configuration de INCA
-     !Config  Def  = none
-     !Config  Help = Choix de configuration de INCA :
-     !Config         'none' = sans INCA
-     !Config         'chem' = INCA avec calcul de chemie
-     !Config         'aero' = INCA avec calcul des aerosols 
-     config_inca = 'none'
-     CALL getin('config_inca',config_inca)
 
      !Config  Key  = ok_dynzon 
@@ -671,9 +662,8 @@
      write(lunout,*)' offline = ', offline
      write(lunout,*)' type_trac = ', type_trac
-     write(lunout,*)' config_inca = ', config_inca
      write(lunout,*)' ok_dynzon = ', ok_dynzon 
      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
      write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave 
-  else
+  ELSE
      !Config  Key  = clon
      !Config  Desc = centre du zoom, longitude
@@ -795,14 +785,4 @@
      CALL getin('type_trac',type_trac)
 
-     !Config  Key  = config_inca
-     !Config  Desc = Choix de configuration de INCA
-     !Config  Def  = none
-     !Config  Help = Choix de configuration de INCA :
-     !Config         'none' = sans INCA
-     !Config         'chem' = INCA avec calcul de chemie
-     !Config         'aero' = INCA avec calcul des aerosols 
-     config_inca = 'none'
-     CALL getin('config_inca',config_inca)
-
      !Config  Key  = ok_dynzon 
      !Config  Desc = sortie des transports zonaux dans la dynamique
@@ -875,6 +855,5 @@
 
      write(lunout,*)' #########################################'
-     write(lunout,*)' Configuration des parametres de cel0' &
-          //'_limit: '
+     write(lunout,*)' Configuration des parametres de cel0_limit: '
      write(lunout,*)' planet_type = ', planet_type
      write(lunout,*)' calend = ', calend
@@ -912,5 +891,4 @@
      write(lunout,*)' offline = ', offline
      write(lunout,*)' type_trac = ', type_trac
-     write(lunout,*)' config_inca = ', config_inca
      write(lunout,*)' ok_dynzon = ', ok_dynzon
      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
@@ -920,6 +898,7 @@
      write(lunout,*)' ok_limit = ', ok_limit
      write(lunout,*)' ok_etat0 = ', ok_etat0
+     write(lunout,*)' ok_guide = ', ok_guide
      write(lunout,*)' read_orop = ', read_orop
-  end IF test_etatinit
+  ENDIF test_etatinit
 
 END SUBROUTINE conf_gcm
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/dynetat0.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/dynetat0.F90	(revision 4368)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/dynetat0.F90	(revision 4368)
@@ -0,0 +1,228 @@
+SUBROUTINE dynetat0(fichnom,vcov,ucov,teta,q,masse,ps,phis,time)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , L.Fairhead
+!-------------------------------------------------------------------------------
+! Purpose: Initial state reading.
+!-------------------------------------------------------------------------------
+  USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName
+  USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str
+  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, &
+                         NF90_CLOSE, NF90_GET_VAR, NF90_NoErr
+  USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3, getKey
+  USE control_mod, ONLY: planet_type
+  USE assert_eq_m, ONLY: assert_eq
+  USE comvert_mod, ONLY: pa,preff
+  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, omeg, rad
+  USE logic_mod, ONLY: fxyhypb, ysinus
+  USE serre_mod, ONLY: clon, clat, grossismx, grossismy
+  USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn, start_time
+  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "description.h"
+  include "iniprint.h"
+!===============================================================================
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
+  REAL, INTENT(OUT) ::  vcov(iip1,jjm, llm)        !--- V COVARIANT WIND
+  REAL, INTENT(OUT) ::  ucov(iip1,jjp1,llm)        !--- U COVARIANT WIND
+  REAL, INTENT(OUT) ::  teta(iip1,jjp1,llm)        !--- POTENTIAL TEMP.
+  REAL, INTENT(OUT) ::     q(iip1,jjp1,llm,nqtot)  !--- TRACERS
+  REAL, INTENT(OUT) :: masse(iip1,jjp1,llm)        !--- MASS PER CELL
+  REAL, INTENT(OUT) ::    ps(iip1,jjp1)            !--- GROUND PRESSURE
+  REAL, INTENT(OUT) ::  phis(iip1,jjp1)            !--- GEOPOTENTIAL
+!===============================================================================
+! Local variables:
+  CHARACTER(LEN=maxlen) :: mesg, var, modname, oldVar
+  INTEGER, PARAMETER :: length=100
+  INTEGER :: iq, fID, vID, idecal, iqParent, iName, iZone, iPhase
+  REAL    :: time, tnat, alpha_ideal, tab_cntrl(length)    !--- RUN PARAMS TABLE
+  LOGICAL :: lSkip, ll
+!-------------------------------------------------------------------------------
+  modname="dynetat0"
+
+!--- Initial state file opening
+  var=fichnom
+  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+  CALL get_var1("controle",tab_cntrl)
+
+!!! AS: idecal is a hack to be able to read planeto starts...
+!!!     .... while keeping everything OK for LMDZ EARTH
+  IF(planet_type=="generic") THEN
+    CALL msg('NOTE NOTE NOTE : Planeto-like start files', modname)
+    idecal = 4
+    annee_ref  = 2000
+  ELSE
+    CALL msg('NOTE NOTE NOTE : Earth-like start files', modname)
+    idecal = 5
+    annee_ref  = tab_cntrl(5)
+  END IF
+  im         = tab_cntrl(1)
+  jm         = tab_cntrl(2)
+  lllm       = tab_cntrl(3)
+  day_ref    = tab_cntrl(4)
+  rad        = tab_cntrl(idecal+1)
+  omeg       = tab_cntrl(idecal+2)
+  g          = tab_cntrl(idecal+3)
+  cpp        = tab_cntrl(idecal+4)
+  kappa      = tab_cntrl(idecal+5)
+  daysec     = tab_cntrl(idecal+6)
+  dtvr       = tab_cntrl(idecal+7)
+  etot0      = tab_cntrl(idecal+8)
+  ptot0      = tab_cntrl(idecal+9)
+  ztot0      = tab_cntrl(idecal+10)
+  stot0      = tab_cntrl(idecal+11)
+  ang0       = tab_cntrl(idecal+12)
+  pa         = tab_cntrl(idecal+13)
+  preff      = tab_cntrl(idecal+14)
+!
+  clon       = tab_cntrl(idecal+15)
+  clat       = tab_cntrl(idecal+16)
+  grossismx  = tab_cntrl(idecal+17)
+  grossismy  = tab_cntrl(idecal+18)
+!
+  IF ( tab_cntrl(idecal+19)==1. )  THEN
+    fxyhypb  = .TRUE.
+!   dzoomx   = tab_cntrl(25)
+!   dzoomy   = tab_cntrl(26)
+!   taux     = tab_cntrl(28)
+!   tauy     = tab_cntrl(29)
+  ELSE
+    fxyhypb = .FALSE.
+    ysinus  = tab_cntrl(idecal+22)==1.
+  END IF
+
+  day_ini    = tab_cntrl(30)
+  itau_dyn   = tab_cntrl(31)
+  start_time = tab_cntrl(32)
+
+!-------------------------------------------------------------------------------
+  CALL msg('rad, omeg, g, cpp, kappa = '//TRIM(strStack(real2str([rad,omeg,g,cpp,kappa]))), modname)
+  CALL check_dim(im,iim,'im','im')
+  CALL check_dim(jm,jjm,'jm','jm')
+  CALL check_dim(lllm,llm,'lm','lllm')
+  CALL get_var1("rlonu",rlonu)
+  CALL get_var1("rlatu",rlatu)
+  CALL get_var1("rlonv",rlonv)
+  CALL get_var1("rlatv",rlatv)
+  CALL get_var2("cu"   ,cu)
+  CALL get_var2("cv"   ,cv)
+  CALL get_var2("aire" ,aire)
+  var="temps"
+  IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
+    CALL msg('missing field <temps> ; trying with <Time>', modname)
+    var="Time"
+    CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  END IF
+  CALL err(NF90_GET_VAR(fID,vID,time),"get",var)
+  CALL get_var2("phisinit",phis)
+  CALL get_var3("ucov",ucov)
+  CALL get_var3("vcov",vcov)
+  CALL get_var3("teta",teta)
+  CALL get_var3("masse",masse)
+  CALL get_var2("ps",ps)
+
+!--- Tracers
+  ll=.FALSE.
+#ifdef REPROBUS
+  ll = NF90_INQ_VARID(fID, 'HNO3tot', vID) /= NF90_NoErr                                 !--- DETECT OLD REPRO start.nc FILE
+#endif
+  DO iq=1,nqtot
+    var = tracers(iq)%name
+    oldVar = new2oldH2O(var)
+    lSkip = ll .AND. var == 'HNO3'                                                       !--- FORCE "HNO3_g" READING FOR "HNO3"
+#ifdef REPROBUS
+    ix = strIdx(newHNO3, var); IF(ix /= 0) oldVar = oldHNO3(ix)                          !--- REPROBUS HNO3 exceptions
+#endif
+#ifdef INCA
+    IF(var == 'O3') oldVar = 'OX'                                                        !--- DEAL WITH INCA OZONE EXCEPTION
+#endif
+    !--------------------------------------------------------------------------------------------------------------------------
+    IF(NF90_INQ_VARID(fID, var, vID) == NF90_NoErr .AND. .NOT.lSkip) THEN                !=== REGULAR CASE: AVAILABLE VARIABLE
+      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",var)
+    !--------------------------------------------------------------------------------------------------------------------------
+    ELSE IF(NF90_INQ_VARID(fID, oldVar, vID) == NF90_NoErr) THEN                         !=== TRY WITH ALTERNATE NAME
+      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to <'//TRIM(oldVar)//'>', modname)
+      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",oldVar)
+    !--------------------------------------------------------------------------------------------------------------------------
+    ELSE IF(tracers(iq)%iso_iGroup == iH2O .AND. niso > 0) THEN                          !=== WATER ISOTOPES
+      iName    = tracers(iq)%iso_iName
+      iPhase   = tracers(iq)%iso_iPhase
+      iqParent = tracers(iq)%iqParent
+      IF(tracers(iq)%iso_iZone == 0) THEN
+         IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
+            CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized with a simplified Rayleigh distillation law.', modname)
+         q(:,:,:,iq) = q(:,:,:,iqParent)*tnat*(q(:,:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
+      ELSE
+         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to its parent isotope concentration.', modname)
+         q(:,:,:,iq) = q(:,:,:,iqIsoPha(iName,iPhase))
+      END IF
+    !--------------------------------------------------------------------------------------------------------------------------
+    ELSE                                                                                 !=== MISSING: SET TO 0
+      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to zero', modname)
+      q(:,:,:,iq)=0.
+    !--------------------------------------------------------------------------------------------------------------------------
+    END IF
+  END DO
+
+  CALL err(NF90_CLOSE(fID),"close",fichnom)
+  day_ini=day_ini+INT(time)
+  time=time-INT(time)
+
+
+  CONTAINS
+
+
+SUBROUTINE check_dim(n1,n2,str1,str2)
+  INTEGER,          INTENT(IN) :: n1, n2
+  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
+  CHARACTER(LEN=maxlen) :: s1, s2
+  IF(n1/=n2) CALL abort_gcm(TRIM(modname), 'value of "'//TRIM(str1)//'" = '//TRIM(int2str(n1))// &
+   ' read in starting file differs from gcm value of "'//TRIM(str2)//'" = '//TRIM(int2str(n2)), 1)
+END SUBROUTINE check_dim
+
+
+SUBROUTINE get_var1(var,v)
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var1
+
+
+SUBROUTINE get_var2(var,v)
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:,:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var2
+
+
+SUBROUTINE get_var3(var,v)
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:,:,:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var3
+
+
+SUBROUTINE err(ierr,typ,nam)
+  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
+  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
+  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
+  IF(ierr==NF90_NoERR) RETURN
+  SELECT CASE(typ)
+    CASE('inq');   mesg="Field <"//TRIM(nam)//"> is missing"
+    CASE('get');   mesg="Reading failed for <"//TRIM(nam)//">"
+    CASE('open');  mesg="File opening failed for <"//TRIM(nam)//">"
+    CASE('close'); mesg="File closing failed for <"//TRIM(nam)//">"
+  END SELECT
+  CALL ABORT_gcm(TRIM(modname),TRIM(mesg),ierr)
+END SUBROUTINE err
+
+END SUBROUTINE dynetat0
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/dynetat0.f90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/dynetat0.f90	(revision 4013)
+++ 	(revision )
@@ -1,203 +1,0 @@
-SUBROUTINE dynetat0(fichnom,vcov,ucov,teta,q,masse,ps,phis,time)
-!
-!-------------------------------------------------------------------------------
-! Authors: P. Le Van , L.Fairhead
-!-------------------------------------------------------------------------------
-! Purpose: Initial state reading.
-!-------------------------------------------------------------------------------
-  USE infotrac
-  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, NF90_NoErr, &
-                         NF90_CLOSE, NF90_GET_VAR
-  USE control_mod, ONLY: planet_type
-  USE assert_eq_m, ONLY: assert_eq
-  USE comvert_mod, ONLY: pa,preff
-  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, omeg, rad
-  USE logic_mod, ONLY: fxyhypb, ysinus
-  USE serre_mod, ONLY: clon, clat, grossismx, grossismy
-  USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn, start_time
-  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
-
-  IMPLICIT NONE
-  include "dimensions.h"
-  include "paramet.h"
-  include "comgeom2.h"
-  include "description.h"
-  include "iniprint.h"
-!===============================================================================
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
-  REAL, INTENT(OUT) ::  vcov(iip1,jjm, llm)        !--- V COVARIANT WIND
-  REAL, INTENT(OUT) ::  ucov(iip1,jjp1,llm)        !--- U COVARIANT WIND
-  REAL, INTENT(OUT) ::  teta(iip1,jjp1,llm)        !--- POTENTIAL TEMP.
-  REAL, INTENT(OUT) ::     q(iip1,jjp1,llm,nqtot)  !--- TRACERS
-  REAL, INTENT(OUT) :: masse(iip1,jjp1,llm)        !--- MASS PER CELL
-  REAL, INTENT(OUT) ::    ps(iip1,jjp1)            !--- GROUND PRESSURE
-  REAL, INTENT(OUT) ::  phis(iip1,jjp1)            !--- GEOPOTENTIAL
-!===============================================================================
-! Local variables:
-  CHARACTER(LEN=256) :: msg, var, modname
-  INTEGER, PARAMETER :: length=100
-  INTEGER :: iq, fID, vID, idecal!, iml, jml, lml, nqt
-  REAL    :: time, tab_cntrl(length)               !--- RUN PARAMS TABLE
-!-------------------------------------------------------------------------------
-  modname="dynetat0"
-
-!--- Initial state file opening
-  var=fichnom
-  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
-  CALL get_var1("controle",tab_cntrl)
-
-!!! AS: idecal is a hack to be able to read planeto starts...
-!!!     .... while keeping everything OK for LMDZ EARTH
-  IF(planet_type=="generic") THEN
-    WRITE(lunout,*)'NOTE NOTE NOTE : Planeto-like start files'
-    idecal = 4
-    annee_ref  = 2000
-  ELSE
-    WRITE(lunout,*)'NOTE NOTE NOTE : Earth-like start files'
-    idecal = 5
-    annee_ref  = tab_cntrl(5)
-  END IF
-  im         = tab_cntrl(1)
-  jm         = tab_cntrl(2)
-  lllm       = tab_cntrl(3)
-  day_ref    = tab_cntrl(4)
-  rad        = tab_cntrl(idecal+1)
-  omeg       = tab_cntrl(idecal+2)
-  g          = tab_cntrl(idecal+3)
-  cpp        = tab_cntrl(idecal+4)
-  kappa      = tab_cntrl(idecal+5)
-  daysec     = tab_cntrl(idecal+6)
-  dtvr       = tab_cntrl(idecal+7)
-  etot0      = tab_cntrl(idecal+8)
-  ptot0      = tab_cntrl(idecal+9)
-  ztot0      = tab_cntrl(idecal+10)
-  stot0      = tab_cntrl(idecal+11)
-  ang0       = tab_cntrl(idecal+12)
-  pa         = tab_cntrl(idecal+13)
-  preff      = tab_cntrl(idecal+14)
-!
-  clon       = tab_cntrl(idecal+15)
-  clat       = tab_cntrl(idecal+16)
-  grossismx  = tab_cntrl(idecal+17)
-  grossismy  = tab_cntrl(idecal+18)
-!
-  IF ( tab_cntrl(idecal+19)==1. )  THEN
-    fxyhypb  = .TRUE.
-!   dzoomx   = tab_cntrl(25)
-!   dzoomy   = tab_cntrl(26)
-!   taux     = tab_cntrl(28)
-!   tauy     = tab_cntrl(29)
-  ELSE
-    fxyhypb = .FALSE.
-    ysinus  = tab_cntrl(idecal+22)==1.
-  END IF
-
-  day_ini    = tab_cntrl(30)
-  itau_dyn   = tab_cntrl(31)
-  start_time = tab_cntrl(32)
-
-!-------------------------------------------------------------------------------
-  WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
-  CALL check_dim(im,iim,'im','im')
-  CALL check_dim(jm,jjm,'jm','jm')
-  CALL check_dim(lllm,llm,'lm','lllm')
-  CALL get_var1("rlonu",rlonu)
-  CALL get_var1("rlatu",rlatu)
-  CALL get_var1("rlonv",rlonv)
-  CALL get_var1("rlatv",rlatv)
-  CALL get_var2("cu"   ,cu)
-  CALL get_var2("cv"   ,cv)
-  CALL get_var2("aire" ,aire)
-  var="temps"
-  IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
-    WRITE(lunout,*)TRIM(modname)//": missing field <temps>"
-    WRITE(lunout,*)TRIM(modname)//": trying with <Time>"; var="Time"
-    CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  END IF
-  CALL err(NF90_GET_VAR(fID,vID,time),"get",var)
-  CALL get_var2("phisinit",phis)
-  CALL get_var3("ucov",ucov)
-  CALL get_var3("vcov",vcov)
-  CALL get_var3("teta",teta)
-  CALL get_var3("masse",masse)
-  CALL get_var2("ps",ps)
-
-!--- Tracers
-  DO iq=1,nqtot
-    var=tname(iq)
-    IF(NF90_INQ_VARID(fID,var,vID)==NF90_NoErr) THEN
-      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",var); CYCLE
-    END IF
-    WRITE(lunout,*)TRIM(modname)//": Tracer <"//TRIM(var)//"> is missing"
-    WRITE(lunout,*)"         It is hence initialized to zero"
-    q(:,:,:,iq)=0.
-   !--- CRisi: for isotops, theoretical initialization using very simplified
-   !           Rayleigh distillation las.
-    IF(ok_isotopes.AND.iso_num(iq)>0) THEN
-      IF(zone_num(iq)==0) q(:,:,:,iq)=q(:,:,:,iqpere(iq))*tnat(iso_num(iq))    &
-     &             *(q(:,:,:,iqpere(iq))/30.e-3)**(alpha_ideal(iso_num(iq))-1)
-      IF(zone_num(iq)==1) q(:,:,:,iq)=q(:,:,:,iqiso(iso_indnum(iq),phase_num(iq)))
-    END IF
-  END DO
-
-  CALL err(NF90_CLOSE(fID),"close",fichnom)
-  day_ini=day_ini+INT(time)
-  time=time-INT(time)
-
-
-  CONTAINS
-
-
-SUBROUTINE check_dim(n1,n2,str1,str2)
-  INTEGER,          INTENT(IN) :: n1, n2
-  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
-  CHARACTER(LEN=100) :: s1, s2
-  IF(n1/=n2) THEN
-    s1='value of '//TRIM(str1)//' ='
-    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
-    WRITE(msg,'(10x,a,i4,2x,a,i4)'),TRIM(ADJUSTL(s1)),n1,TRIM(ADJUSTL(s2)),n2
-    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
-  END IF
-END SUBROUTINE check_dim
-
-
-SUBROUTINE get_var1(var,v)
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var1
-
-
-SUBROUTINE get_var2(var,v)
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:,:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var2
-
-
-SUBROUTINE get_var3(var,v)
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:,:,:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var3
-
-
-SUBROUTINE err(ierr,typ,nam)
-  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
-  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
-  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
-  IF(ierr==NF90_NoERR) RETURN
-  SELECT CASE(typ)
-    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
-    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
-    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
-    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
-  END SELECT
-  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
-END SUBROUTINE err
-
-END SUBROUTINE dynetat0
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/dynredem.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/dynredem.F90	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/dynredem.F90	(revision 4368)
@@ -7,11 +7,11 @@
   USE IOIPSL
 #endif
-  USE infotrac
+  USE strings_mod, ONLY: maxlen
+  USE infotrac, ONLY: nqtot, tracers
   USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL,    &
                     NF90_CLOSE,  NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER,   &
                     NF90_64BIT_OFFSET
   USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, err, modname, fil
-  USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt,pa,preff, &
-                              nivsig,nivsigs
+  USE comvert_mod,  ONLY: ap, bp, presnivs, pa, preff, nivsig, nivsigs
   USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad
   USE logic_mod, ONLY: fxyhypb, ysinus
@@ -34,14 +34,14 @@
 !===============================================================================
 ! Local variables:
-  INTEGER :: iq, l
+  INTEGER :: iq
   INTEGER, PARAMETER :: length=100
   REAL    :: tab_cntrl(length)                     !--- RUN PARAMETERS TABLE
 !   For NetCDF:
-  CHARACTER(LEN=30) :: unites
+  CHARACTER(LEN=maxlen) :: unites
   INTEGER :: indexID
   INTEGER :: rlonuID, rlonvID, rlatuID, rlatvID
-  INTEGER :: sID, sigID, nID, vID, timID
+  INTEGER :: sID, sigID, nID, timID
   INTEGER :: yyears0, jjour0, mmois0
-  REAL    :: zan0, zjulian, hours
+  REAL    :: zjulian, hours
 !===============================================================================
   modname='dynredem0'; fil=fichnom
@@ -138,5 +138,5 @@
 
 !--- Define fields saved later
-  WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')"),&
+  WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')") &
                yyears0,mmois0,jjour0
   CALL cre_var(nid,"temps","Temps de simulation",[timID],unites)
@@ -145,5 +145,5 @@
   CALL cre_var(nid,"teta" ,"Temperature",[rlonvID,rlatuID,sID,timID])
   DO iq=1,nqtot
-    CALL cre_var(nid,tname(iq),ttext(iq),[rlonvID,rlatuID,sID,timID])
+    CALL cre_var(nid,tracers(iq)%name,tracers(iq)%longName,[rlonvID,rlatuID,sID,timID])
   END DO
   CALL cre_var(nid,"masse","Masse d air"    ,[rlonvID,rlatuID,sID,timID])
@@ -166,5 +166,6 @@
 ! Purpose: Write the NetCDF restart file (append).
 !-------------------------------------------------------------------------------
-  USE infotrac
+  USE strings_mod, ONLY: maxlen
+  USE infotrac, ONLY: nqtot, tracers, types_trac
   USE control_mod
   USE netcdf,   ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_GET_VAR, NF90_INQ_VARID,  &
@@ -192,9 +193,9 @@
 !===============================================================================
 ! Local variables:
-  INTEGER :: l, iq, nid, vID, ierr, nid_trac, vID_trac
+  INTEGER :: iq, nid, vID, ierr, nid_trac, vID_trac
   INTEGER, SAVE :: nb=0
   INTEGER, PARAMETER :: length=100
   REAL               :: tab_cntrl(length) ! tableau des parametres du run
-  CHARACTER(LEN=256) :: var, dum
+  CHARACTER(LEN=maxlen) :: var, dum
   LOGICAL            :: lread_inca
 !===============================================================================
@@ -227,9 +228,9 @@
 !--- Tracers in file "start_trac.nc" (added by Anne)
   lread_inca=.FALSE.; fil="start_trac.nc"
-  IF(type_trac=='inca' .OR. type_trac=='inco') INQUIRE(FILE=fil,EXIST=lread_inca)
+  IF(ANY(types_trac=='inca') .OR. ANY(types_trac=='inco')) INQUIRE(FILE=fil,EXIST=lread_inca)
   IF(lread_inca) CALL err(NF90_OPEN(fil,NF90_NOWRITE,nid_trac),"open")
 
 !--- Save tracers
-  DO iq=1,nqtot; var=tname(iq); ierr=-1
+  DO iq=1,nqtot; var=TRIM(tracers(iq)%name); ierr=-1
     IF(lread_inca) THEN                  !--- Possibly read from "start_trac.nc"
       fil="start_trac.nc"
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/gcm.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/gcm.F90	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/gcm.F90	(revision 4368)
@@ -20,5 +20,5 @@
 
   USE filtreg_mod
-  USE infotrac
+  USE infotrac, ONLY: nqtot, init_infotrac
   USE control_mod
   USE mod_const_mpi, ONLY: COMM_LMDZ
@@ -178,5 +178,5 @@
 #ifdef CPP_IOIPSL
   if (calend == 'earth_360d') then
-     call ioconf_calendar('360d')
+     call ioconf_calendar('360_day')
      write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
   else if (calend == 'earth_365d') then
@@ -205,5 +205,5 @@
   !  Choix du nombre de traceurs et du schema pour l'advection
   !  dans fichier traceur.def, par default ou via INCA
-  call infotrac_init
+  call init_infotrac
 
   ! Allocation de la tableau q : champs advectes   
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/guide_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/guide_mod.F90	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/guide_mod.F90	(revision 4368)
@@ -39,4 +39,6 @@
   REAL, PRIVATE, SAVE     :: lon_min_g,lon_max_g
   REAL, PRIVATE, SAVE     :: tau_lon,tau_lat
+
+  REAL, PRIVATE, SAVE     :: plim_guide_BL
 
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_u,alpha_v 
@@ -67,4 +69,5 @@
   SUBROUTINE guide_init
 
+    use netcdf, only: nf90_noerr
     USE control_mod, ONLY: day_step
     USE serre_mod, ONLY: grossismx
@@ -113,4 +116,6 @@
     CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
     CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
+    CALL getpar('plim_guide_BL',85000.,plim_guide_BL,'BL top presnivs value')
+
 
 ! Sauvegarde du for�age
@@ -169,5 +174,5 @@
        if (ncidpl.eq.-99) then
           rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
-          if (rcod.NE.NF_NOERR) THEN
+          if (rcod.NE.NF90_NOERR) THEN
              abort_message=' Nudging error -> no file apbp.nc'
              CALL abort_gcm(modname,abort_message,1)
@@ -177,5 +182,5 @@
        if (ncidpl.EQ.-99) then
           rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
-          if (rcod.NE.NF_NOERR) THEN
+          if (rcod.NE.NF90_NOERR) THEN
              abort_message=' Nudging error -> no file P.nc'
              CALL abort_gcm(modname,abort_message,1)
@@ -186,5 +191,5 @@
            if (ncidpl.eq.-99) then
                rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
-               if (rcod.NE.NF_NOERR) THEN
+               if (rcod.NE.NF90_NOERR) THEN
                   CALL abort_gcm(modname, &
                        ' Nudging error -> no file u.nc',1)
@@ -195,5 +200,5 @@
            if (ncidpl.eq.-99) then
                rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
-               if (rcod.NE.NF_NOERR) THEN
+               if (rcod.NE.NF90_NOERR) THEN
                   CALL abort_gcm(modname, &
                        ' Nudging error -> no file v.nc',1)
@@ -203,5 +208,5 @@
            if (ncidpl.eq.-99) then
                rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
-               if (rcod.NE.NF_NOERR) THEN
+               if (rcod.NE.NF90_NOERR) THEN
                   CALL abort_gcm(modname, &
                        ' Nudging error -> no file T.nc',1)
@@ -211,5 +216,5 @@
            if (ncidpl.eq.-99) then
                rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
-               if (rcod.NE.NF_NOERR) THEN
+               if (rcod.NE.NF90_NOERR) THEN
                   CALL abort_gcm(modname, &
                        ' Nudging error -> no file hur.nc',1)
@@ -220,6 +225,6 @@
     endif 
     error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
-    IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
-    IF (error.NE.NF_NOERR) THEN
+    IF (error.NE.NF90_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
+    IF (error.NE.NF90_NOERR) THEN
         CALL abort_gcm(modname,'Nudging: error reading pressure levels',1)
     ENDIF
@@ -399,5 +404,5 @@
         else
             do l=1,llm
-                alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
+               alpha_pcor(l)=(1.+tanh(((plim_guide_BL-presnivs(l))/preff)/0.05))/2.
             enddo
         endif
@@ -1078,7 +1083,8 @@
   SUBROUTINE guide_read(timestep)
 
+    use netcdf, only: NF90_GET_VAR, nf90_noerr
+
     IMPLICIT NONE
 
-    include "netcdf.inc"
     include "dimensions.h"
     include "paramet.h"
@@ -1103,24 +1109,24 @@
     if (first) then
          ncidpl=-99
-         write(*,*),trim(modname)//': opening nudging files '
+         write(*,*) trim(modname)//': opening nudging files '
 ! Niveaux de pression si non constants
          if (guide_plevs.EQ.1) then
-             write(*,*),trim(modname)//' Reading nudging on model levels'
+             write(*,*) trim(modname)//' Reading nudging on model levels'
              rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no file apbp.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
              rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no AP variable in file apbp.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
              rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no BP variable in file apbp.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidpl,varidap',ncidpl,varidap
+             write(*,*) trim(modname)//' ncidpl,varidap',ncidpl,varidap
          endif
 
@@ -1128,14 +1134,14 @@
          if (guide_plevs.EQ.2) then
              rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no file P.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
              rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
-             IF (rcode.NE.NF_NOERR) THEN 
+             IF (rcode.NE.NF90_NOERR) THEN 
               abort_message='Nudging: error -> no PRES variable in file P.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidp,varidp',ncidp,varidp
+             write(*,*) trim(modname)//' ncidp,varidp',ncidp,varidp
              if (ncidpl.eq.-99) ncidpl=ncidp
          endif
@@ -1144,14 +1150,14 @@
          if (guide_u) then
              rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no file u.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
              rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no UWND variable in file u.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidu,varidu',ncidu,varidu
+             write(*,*) trim(modname)//' ncidu,varidu',ncidu,varidu
              if (ncidpl.eq.-99) ncidpl=ncidu
 
@@ -1175,14 +1181,14 @@
          if (guide_v) then
              rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no file v.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
              rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no VWND variable in file v.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidv,varidv',ncidv,varidv
+             write(*,*) trim(modname)//' ncidv,varidv',ncidv,varidv
              if (ncidpl.eq.-99) ncidpl=ncidv
              
@@ -1208,14 +1214,14 @@
          if (guide_T) then
              rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no file T.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
              rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no AIR variable in file T.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidT,varidT',ncidt,varidt
+             write(*,*) trim(modname)//' ncidT,varidT',ncidt,varidt
              if (ncidpl.eq.-99) ncidpl=ncidt
 
@@ -1239,14 +1245,14 @@
          if (guide_Q) then
              rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no file hur.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
              rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no RH variable in file hur.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
+             write(*,*) trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
              if (ncidpl.eq.-99) ncidpl=ncidQ
 
@@ -1270,14 +1276,14 @@
          if ((guide_P).OR.(guide_modele)) then
              rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no file ps.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
              rcode = nf90_inq_varid(ncidps, 'SP', varidps)
-             IF (rcode.NE.NF_NOERR) THEN
+             IF (rcode.NE.NF90_NOERR) THEN
               abort_message='Nudging: error -> no SP variable in file ps.nc'
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidps,varidps',ncidps,varidps
+             write(*,*) trim(modname)//' ncidps,varidps',ncidps,varidps
          endif
 ! Coordonnee verticale
@@ -1285,21 +1291,12 @@
               rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
               IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
-              write(*,*),trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
+              write(*,*) trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
          endif
 ! Coefs ap, bp pour calcul de la pression aux differents niveaux
          if (guide_plevs.EQ.1) then
-#ifdef NC_DOUBLE
-             status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
-             status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
-#else
-             status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
-             status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
-#endif
+             status=NF90_GET_VAR(ncidpl,varidap,apnc,[1],[nlevnc])
+             status=NF90_GET_VAR(ncidpl,varidbp,bpnc,[1],[nlevnc])
          ELSEIF (guide_plevs.EQ.0) THEN
-#ifdef NC_DOUBLE
-             status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
-#else
-             status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
-#endif
+             status=NF90_GET_VAR(ncidpl,varidpl,apnc,[1],[nlevnc])
 !FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous
              IF(convert_Pa) apnc=apnc*100.! conversion en Pascals
@@ -1326,9 +1323,5 @@
 ! Pression 
      if (guide_plevs.EQ.2) then
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,pnat2)
-#else
-         status=NF_GET_VARA_REAL(ncidp,varidp,start,count,pnat2)
-#endif
+         status=NF90_GET_VAR(ncidp,varidp,pnat2,start,count)
          IF (invert_y) THEN
 !           PRINT*,"Invertion impossible actuellement"
@@ -1340,9 +1333,5 @@
 !  Vent zonal
      if (guide_u) then
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2)
-#else
-         status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unat2)
-#endif
+         status=NF90_GET_VAR(ncidu,varidu,unat2,start,count)
          IF (invert_y) THEN
            CALL invert_lat(iip1,jjp1,nlevnc,unat2)
@@ -1352,9 +1341,5 @@
 !  Temperature
      if (guide_T) then
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2)
-#else
-         status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnat2)
-#endif
+         status=NF90_GET_VAR(ncidt,varidt,tnat2,start,count)
          IF (invert_y) THEN
            CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
@@ -1364,9 +1349,5 @@
 !  Humidite
      if (guide_Q) then
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2)
-#else
-         status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,qnat2)
-#endif
+         status=NF90_GET_VAR(ncidQ,varidQ,qnat2,start,count)
          IF (invert_y) THEN
            CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
@@ -1378,9 +1359,5 @@
      if (guide_v) then
          count(2)=jjm
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2)
-#else
-         status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnat2)
-#endif
+         status=NF90_GET_VAR(ncidv,varidv,vnat2,start,count)
          IF (invert_y) THEN
            CALL invert_lat(iip1,jjm,nlevnc,vnat2)
@@ -1395,9 +1372,5 @@
          count(3)=1
          count(4)=0
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2)
-#else
-         status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnat2)
-#endif
+         status=NF90_GET_VAR(ncidps,varidps,psnat2,start,count)
          IF (invert_y) THEN
            CALL invert_lat(iip1,jjp1,1,psnat2)
@@ -1410,7 +1383,8 @@
   SUBROUTINE guide_read2D(timestep)
 
+    use netcdf, only: nf90_get_var, nf90_noerr
+
     IMPLICIT NONE
 
-    include "netcdf.inc"
     include "dimensions.h"
     include "paramet.h"
@@ -1443,15 +1417,15 @@
            write(*,*)trim(modname)//' Reading nudging on model levels'
            rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no file apbp.nc'
            CALL abort_gcm(modname,abort_message,1)
            ENDIF
            rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no AP variable in file apbp.nc'
            CALL abort_gcm(modname,abort_message,1)
            ENDIF
            rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no BP variable in file apbp.nc'
              CALL abort_gcm(modname,abort_message,1)
@@ -1462,10 +1436,10 @@
          if (guide_plevs.EQ.2) then
            rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no file P.nc'
              CALL abort_gcm(modname,abort_message,1)
            ENDIF
            rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no PRES variable in file P.nc'
              CALL abort_gcm(modname,abort_message,1)
@@ -1477,10 +1451,10 @@
          if (guide_u) then
            rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no file u.nc'
              CALL abort_gcm(modname,abort_message,1)
            ENDIF
            rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no UWND variable in file u.nc'
              CALL abort_gcm(modname,abort_message,1)
@@ -1492,10 +1466,10 @@
          if (guide_v) then
            rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no file v.nc'
              CALL abort_gcm(modname,abort_message,1)
            ENDIF
            rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no VWND variable in file v.nc'
              CALL abort_gcm(modname,abort_message,1)
@@ -1507,10 +1481,10 @@
          if (guide_T) then
            rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no file T.nc'
              CALL abort_gcm(modname,abort_message,1)
            ENDIF
            rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no AIR variable in file T.nc'
              CALL abort_gcm(modname,abort_message,1)
@@ -1522,10 +1496,10 @@
          if (guide_Q) then
            rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no file hur.nc'
              CALL abort_gcm(modname,abort_message,1)
            ENDIF
            rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no RH,variable in file hur.nc'
              CALL abort_gcm(modname,abort_message,1)
@@ -1537,10 +1511,10 @@
          if ((guide_P).OR.(guide_modele)) then
            rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no file ps.nc'
              CALL abort_gcm(modname,abort_message,1)
            ENDIF
            rcode = nf90_inq_varid(ncidps, 'SP', varidps)
-           IF (rcode.NE.NF_NOERR) THEN
+           IF (rcode.NE.NF90_NOERR) THEN
              abort_message='Nudging: error -> no SP variable in file ps.nc'
              CALL abort_gcm(modname,abort_message,1)
@@ -1556,17 +1530,8 @@
 ! Coefs ap, bp pour calcul de la pression aux differents niveaux
          if (guide_plevs.EQ.1) then
-#ifdef NC_DOUBLE
-             status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
-             status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
-#else
-             status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
-             status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
-#endif
+             status=NF90_GET_VAR(ncidpl,varidap,apnc,[1],[nlevnc])
+             status=NF90_GET_VAR(ncidpl,varidbp,bpnc,[1],[nlevnc])
          elseif (guide_plevs.EQ.0) THEN
-#ifdef NC_DOUBLE
-             status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
-#else
-             status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
-#endif
+             status=NF90_GET_VAR(ncidpl,varidpl,apnc,[1],[nlevnc])
              apnc=apnc*100.! conversion en Pascals
              bpnc(:)=0.
@@ -1592,9 +1557,5 @@
 !  Pression
      if (guide_plevs.EQ.2) then
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,zu)
-#else
-         status=NF_GET_VARA_REAL(ncidp,varidp,start,count,zu)
-#endif
+         status=NF90_GET_VAR(ncidp,varidp,zu,start,count)
          DO i=1,iip1
              pnat2(i,:,:)=zu(:,:)
@@ -1609,9 +1570,5 @@
 !  Vent zonal
      if (guide_u) then
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu)
-#else
-         status=NF_GET_VARA_REAL(ncidu,varidu,start,count,zu)
-#endif
+         status=NF90_GET_VAR(ncidu,varidu,zu,start,count)
          DO i=1,iip1
              unat2(i,:,:)=zu(:,:)
@@ -1626,9 +1583,5 @@
 !  Temperature
      if (guide_T) then
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu)
-#else
-         status=NF_GET_VARA_REAL(ncidt,varidt,start,count,zu)
-#endif
+         status=NF90_GET_VAR(ncidt,varidt,zu,start,count)
          DO i=1,iip1
              tnat2(i,:,:)=zu(:,:)
@@ -1643,9 +1596,5 @@
 !  Humidite
      if (guide_Q) then
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu)
-#else
-         status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,zu)
-#endif
+         status=NF90_GET_VAR(ncidQ,varidQ,zu,start,count)
          DO i=1,iip1
              qnat2(i,:,:)=zu(:,:)
@@ -1661,9 +1610,5 @@
      if (guide_v) then
          count(2)=jjm
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv)
-#else
-         status=NF_GET_VARA_REAL(ncidv,varidv,start,count,zv)
-#endif
+         status=NF90_GET_VAR(ncidv,varidv,zv,start,count)
          DO i=1,iip1
              vnat2(i,:,:)=zv(:,:)
@@ -1683,9 +1628,5 @@
          count(3)=1
          count(4)=0
-#ifdef NC_DOUBLE
-         status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1))
-#else
-         status=NF_GET_VARA_REAL(ncidps,varidps,start,count,zu(:,1))
-#endif
+         status=NF90_GET_VAR(ncidps,varidps,zu(:,1),start,count)
          DO i=1,iip1
              psnat2(i,:)=zu(:,1)
@@ -1706,5 +1647,5 @@
     USE comvert_mod, ONLY: presnivs
     use netcdf95, only: nf95_def_var, nf95_put_var
-    use netcdf, only: nf90_float
+    use netcdf, only: nf90_float, nf90_def_var
     
     IMPLICIT NONE
@@ -1748,13 +1689,13 @@
 
 ! Creation des variables dimensions
-        ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
-        ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
-        ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
-        ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
-        ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
-        ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
-        ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
-        ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
-        ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
+        ierr=NF90_DEF_VAR(nid,"LONU",NF90_FLOAT,id_lonu,vid_lonu)
+        ierr=NF90_DEF_VAR(nid,"LONV",NF90_FLOAT,id_lonv,vid_lonv)
+        ierr=NF90_DEF_VAR(nid,"LATU",NF90_FLOAT,id_latu,vid_latu)
+        ierr=NF90_DEF_VAR(nid,"LATV",NF90_FLOAT,id_latv,vid_latv)
+        ierr=NF90_DEF_VAR(nid,"LEVEL",NF90_FLOAT,id_lev,vid_lev)
+        ierr=NF90_DEF_VAR(nid,"cu",NF90_FLOAT,(/id_lonu,id_latu/),vid_cu)
+        ierr=NF90_DEF_VAR(nid,"cv",NF90_FLOAT,(/id_lonv,id_latv/),vid_cv)
+        ierr=NF90_DEF_VAR(nid,"au",NF90_FLOAT,(/id_lonu,id_latu/),vid_au)
+        ierr=NF90_DEF_VAR(nid,"av",NF90_FLOAT,(/id_lonv,id_latv/),vid_av)
         call nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
              varid_alpha_t)
@@ -1794,33 +1735,33 @@
 ! Pressure (GCM)
         dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-        ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,4,dim4,varid)
+        ierr = NF90_DEF_VAR(nid,"SP",NF90_FLOAT,dim4,varid)
 ! Surface pressure (guidage)
         IF (guide_P) THEN
             dim3=(/id_lonv,id_latu,id_tim/)
-            ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
+            ierr = NF90_DEF_VAR(nid,"ps",NF90_FLOAT,dim3,varid)
         ENDIF
 ! Zonal wind
         IF (guide_u) THEN
             dim4=(/id_lonu,id_latu,id_lev,id_tim/)
-            ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
-            ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
-            ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
+            ierr = NF90_DEF_VAR(nid,"u",NF90_FLOAT,dim4,varid)
+            ierr = NF90_DEF_VAR(nid,"ua",NF90_FLOAT,dim4,varid)
+            ierr = NF90_DEF_VAR(nid,"ucov",NF90_FLOAT,dim4,varid)
         ENDIF
 ! Merid. wind
         IF (guide_v) THEN
             dim4=(/id_lonv,id_latv,id_lev,id_tim/)
-            ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
-            ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
-            ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
+            ierr = NF90_DEF_VAR(nid,"v",NF90_FLOAT,dim4,varid)
+            ierr = NF90_DEF_VAR(nid,"va",NF90_FLOAT,dim4,varid)
+            ierr = NF90_DEF_VAR(nid,"vcov",NF90_FLOAT,dim4,varid)
         ENDIF
 ! Pot. Temperature
         IF (guide_T) THEN
             dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-            ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
+            ierr = NF90_DEF_VAR(nid,"teta",NF90_FLOAT,dim4,varid)
         ENDIF
 ! Specific Humidity
         IF (guide_Q) THEN
             dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-            ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
+            ierr = NF90_DEF_VAR(nid,"q",NF90_FLOAT,dim4,varid)
         ENDIF
         
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/iniacademic.F90	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/iniacademic.F90	(revision 4368)
@@ -5,6 +5,8 @@
 
   USE filtreg_mod, ONLY: inifilr
-  USE infotrac
+  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
   USE control_mod, ONLY: day_step,planet_type
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
 #ifdef CPP_IOIPSL
   USE IOIPSL, ONLY: getin
@@ -14,6 +16,4 @@
 #endif
   USE Write_Field
-  use exner_hyb_m, only: exner_hyb
-  use exner_milieu_m, only: exner_milieu
   USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm
   USE logic_mod, ONLY: iflag_phys, read_start
@@ -21,4 +21,5 @@
   USE temps_mod, ONLY: annee_ref, day_ini, day_ref
   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+  USE readTracFiles_mod, ONLY: addPhase
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -61,5 +62,5 @@
   real tetastrat ! potential temperature in the stratosphere, in K
   real tetajl(jjp1,llm)
-  INTEGER i,j,l,lsup,ij
+  INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
 
   REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
@@ -72,9 +73,8 @@
   integer idum
 
-  REAL zdtvr
+  REAL zdtvr, tnat, alpha_ideal
   
   character(len=*),parameter :: modname="iniacademic"
   character(len=80) :: abort_message
-
 
   ! Sanity check: verify that options selected by user are not incompatible
@@ -130,7 +130,9 @@
         ps(:)=preff
      endif
+
      ! ground geopotential
      phis(:)=0.
      CALL pression ( ip1jmp1, ap, bp, ps, p       )
+
      if (pressure_exner) then
        CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
@@ -193,5 +195,4 @@
      CALL getin('tetanoise',tetanoise)
 
-
      ! 2. Initialize fields towards which to relax
      ! Friction
@@ -246,5 +247,4 @@
      IF (.NOT. read_start) THEN
         ! bulk initialization of temperature
-
         IF (iflag_phys>10000) THEN
         ! Particular case to impose a constant temperature T0=0.01*iflag_physx
@@ -253,5 +253,4 @@
            teta(:,:)=tetarappel(:,:)
         ENDIF
-
         ! geopotential
         CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
@@ -275,23 +274,22 @@
         if (planet_type=="earth") then
            ! Earth: first two tracers will be water
-           do i=1,nqtot
-              if (i == 1) q(:,:,i)=1.e-10
-              if (i == 2) q(:,:,i)=1.e-15
-              if (i.gt.2) q(:,:,i)=0.
+           do iq=1,nqtot
+              q(:,:,iq)=0.
+              IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10
+              IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15
 
               ! 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(:,:,i)=q(:,:,iqpere(i))             &
-      &                  *tnat(iso_num(i))               &
-      &                  *(q(:,:,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(:,:,i)=q(:,:,iqiso(iso_indnum(i),phase_num(i)))
-                endif
-              endif !if (ok_isotopes) then
-
+              iName    = tracers(iq)%iso_iName
+              if (niso <= 0 .OR. iName <= 0) CYCLE
+              iPhase   = tracers(iq)%iso_iPhase
+              iqParent = tracers(iq)%iqParent
+              IF(tracers(iq)%iso_iZone == 0) THEN
+                 IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
+                    CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+                 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
+              ELSE
+                 q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
+              END IF
            enddo
         else
@@ -299,7 +297,5 @@
         endif ! of if (planet_type=="earth")
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
-        endif !if (ok_iso_verif) then
+        call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
 
         ! add random perturbation to temperature
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/leapfrog.F
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/leapfrog.F	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/leapfrog.F	(revision 4368)
@@ -11,5 +11,5 @@
       use IOIPSL
 #endif
-      USE infotrac, ONLY: nqtot,ok_iso_verif
+      USE infotrac, ONLY: nqtot, isoCheck
       USE guide_mod, ONLY : guide_main
       USE write_field, ONLY: writefield
@@ -26,4 +26,5 @@
       USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
      &                        start_time,dt
+      USE strings_mod, ONLY: msg
 
       IMPLICIT NONE
@@ -237,7 +238,5 @@
       jH_cur = jH_cur - int(jH_cur)
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
-        endif !if (ok_iso_verif) then
+      call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
 
 #ifdef CPP_IOIPSL
@@ -271,7 +270,5 @@
 !      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
-        endif !if (ok_iso_verif) then
+      call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
 
    2  CONTINUE ! Matsuno backward or leapfrog step begins here
@@ -324,7 +321,5 @@
 
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
-        endif !if (ok_iso_verif) then
+      call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
 
 c-----------------------------------------------------------------------
@@ -345,8 +340,6 @@
 c   -------------------------------------------------------------
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,
+      call check_isotopes_seq(q,ip1jmp1,
      &           'leapfrog 686: avant caladvtrac')
-        endif !if (ok_iso_verif) then
 
       IF( forward. OR . leapf )  THEN
@@ -376,8 +369,6 @@
 c   ----------------------------------
 
-        if (ok_iso_verif) then
-           write(*,*) 'leapfrog 720' 
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
-        endif !if (ok_iso_verif) then
+       CALL msg('720', modname, isoCheck)
+       call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
         
        CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
@@ -385,8 +376,6 @@
 !     $              finvmaold                                    )
 
-       if (ok_iso_verif) then
-          write(*,*) 'leapfrog 724'
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
-        endif !if (ok_iso_verif) then
+       CALL msg('724', modname, isoCheck)
+       call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
 
 c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
@@ -451,5 +440,5 @@
 c+jld
 
-c  Diagnostique de conservation de l'énergie : initialisation
+c  Diagnostique de conservation de l'energie : initialisation
          IF (ip_ebil_dyn.ge.1 ) THEN 
           ztit='bil dyn'
@@ -498,5 +487,5 @@
        
 c
-c  Diagnostique de conservation de l'énergie : difference
+c  Diagnostique de conservation de l'energie : difference
          IF (ip_ebil_dyn.ge.1 ) THEN 
           ztit='bil phys'
@@ -552,7 +541,5 @@
         CALL massdair(p,masse)
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
-        endif !if (ok_iso_verif) then
+        call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
 
 c-----------------------------------------------------------------------
@@ -639,7 +626,5 @@
 c   preparation du pas d'integration suivant  ......
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
-        endif !if (ok_iso_verif) then
+      call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
 
       IF ( .NOT.purmats ) THEN
@@ -703,7 +688,5 @@
             ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
-        endif !if (ok_iso_verif) then
+            call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
 
 c-----------------------------------------------------------------------
@@ -790,7 +773,5 @@
       ELSE ! of IF (.not.purmats)
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
-        endif !if (ok_iso_verif) then
+            call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
 
 c       ........................................................
@@ -817,7 +798,5 @@
             ELSE ! of IF(forward) i.e. backward step
  
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
-        endif !if (ok_iso_verif) then  
+              call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
 
               IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/qminimum.F
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/qminimum.F	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/qminimum.F	(revision 4368)
@@ -4,5 +4,7 @@
       SUBROUTINE qminimum( q,nqtot,deltap )
 
-      USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif
+      USE infotrac, ONLY: niso, ntiso,iqIsoPha, tracers
+      USE strings_mod, ONLY: strIdx
+      USE readTracFiles_mod, ONLY: addPhase
       IMPLICIT none
 c
@@ -16,10 +18,8 @@
       REAL q(ip1jmp1,llm,nqtot), deltap(ip1jmp1,llm)
 c
-      INTEGER iq_vap, iq_liq
-      PARAMETER ( iq_vap = 1 ) ! indice pour l'eau vapeur
-      PARAMETER ( iq_liq = 2 ) ! indice pour l'eau liquide
-      REAL seuil_vap, seuil_liq
-      PARAMETER ( seuil_vap = 1.0e-10 ) ! seuil pour l'eau vapeur
-      PARAMETER ( seuil_liq = 1.0e-11 ) ! seuil pour l'eau liquide
+      LOGICAL, SAVE :: first=.TRUE.
+      INTEGER, SAVE :: iq_vap, iq_liq        ! indices pour l'eau vapeur/liquide
+      REAL, PARAMETER :: seuil_vap = 1.0e-10 ! seuil pour l'eau vapeur
+      REAL, PARAMETER :: seuil_liq = 1.0e-11 ! seuil pour l'eau liquide
 c
 c  NB. ....( Il est souhaitable mais non obligatoire que les valeurs des
@@ -43,4 +43,10 @@
       !INTEGER nb_pump
       INTEGER ixt
+
+      IF(first) THEN
+         iq_vap = strIdx(tracers(:)%name, addPhase('H2O', 'g'))
+         iq_liq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
+         first = .FALSE.
+      END IF
 c
 c Quand l'eau liquide est trop petite (ou negative), on prend
@@ -49,7 +55,5 @@
 c
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'qminimum 52')   
-        endif !if (ok_iso_verif) then     
+      call check_isotopes_seq(q,ip1jmp1,'qminimum 52')   
 
       zx_defau_diag(:,:,:)=0.0
@@ -59,8 +63,6 @@
           if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
 
-              if (ok_isotopes) then
-                 zx_defau_diag(i,k,iq_liq)=AMAX1
+              if (niso > 0) 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
@@ -80,7 +82,6 @@
           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
+            if (niso > 0)
+     &        zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 )
 
             q(i,k-1,iq) =  q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) *
@@ -110,5 +111,5 @@
 
       !write(*,*) 'qminimum 128'
-      if (ok_isotopes) then
+      if (niso > 0) then
       ! CRisi: traiter de même les traceurs d'eau
       ! Mais il faut les prendre à l'envers pour essayer de conserver la
@@ -130,14 +131,16 @@
           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)
+              do ixt=1,ntiso
+               q(i,k,iqIsoPha(ixt,iq_vap))=q(i,k,iqIsoPha(ixt,iq_vap))
+     :           +zx_defau_diag(i,k,iq_vap)
+     :           *q(i,k-1,iqIsoPha(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
                 
               ! et on la retranche en k-1
-               q(i,k-1,iqiso(ixt,iq_vap))=q(i,k-1,iqiso(ixt,iq_vap))
+               q(i,k-1,iqIsoPha(ixt,iq_vap))=
+     :            q(i,k-1,iqIsoPha(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)
+     :              *q(i,k-1,iqIsoPha(ixt,iq_vap))
+     :              /q_follow(i,k-1,iq_vap)
 
               enddo !do ixt=1,niso
@@ -151,7 +154,5 @@
        enddo !do k=2,llm
 
-        if (ok_iso_verif) then     
-           call check_isotopes_seq(q,ip1jmp1,'qminimum 168')
-        endif !if (ok_iso_verif) then
+       call check_isotopes_seq(q,ip1jmp1,'qminimum 168')
         
       
@@ -163,12 +164,12 @@
 
               ! 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))
+              do ixt=1,ntiso
+               q(i,k,iqIsoPha(ixt,iq_liq))=q(i,k,iqIsoPha(ixt,iq_liq))
      :              +zx_defau_diag(i,k,iq_liq)
-     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)
+     :              *q(i,k,iqIsoPha(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))
+               q(i,k,iqIsoPha(ixt,iq_vap))=q(i,k,iqIsoPha(ixt,iq_vap))
      :              -zx_defau_diag(i,k,iq_liq)
-     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)   
+     :              *q(i,k,iqIsoPha(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)
@@ -180,9 +181,7 @@
        enddo !do k=2,llm  
 
-        if (ok_iso_verif) then
-           call check_isotopes_seq(q,ip1jmp1,'qminimum 197')
-        endif !if (ok_iso_verif) then
+       call check_isotopes_seq(q,ip1jmp1,'qminimum 197')
 
-      endif !if (ok_isotopes) then
+      endif !if (niso > 0) then
       !write(*,*) 'qminimum 188'
       
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/replay3d.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/replay3d.F90	(revision 4368)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/replay3d.F90	(revision 4368)
@@ -0,0 +1,195 @@
+!
+! $Id: gcm.F90 3579 2019-10-09 13:11:07Z fairhead $
+!
+!
+!
+PROGRAM replay3d
+
+
+
+USE comvert_mod, only :  preff, pa
+USE inigeomphy_mod, ONLY: inigeomphy
+
+
+  USE control_mod
+  USE temps_mod, ONLY: calend,start_time,annee_ref,day_ref, &
+                     itau_dyn,itau_phy,day_ini,jD_ref,jH_ref,day_end
+  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, g, r, rad
+  USE logic_mod, ONLY: ecripar, iflag_phys, read_start
+
+  USE serre_mod, ONLY : clon,clat,transx,transy,alphax,alphay,pxo,pyo,              &
+        grossismx, grossismy, dzoomx, dzoomy,taux,tauy
+  USE mod_const_mpi, ONLY: comm_lmdz
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
+  ! A nettoyer. On ne veut qu'une ou deux routines d'interface 
+  ! dynamique -> physique pour l'initialisation
+#ifdef CPP_PHYS
+  USE iniphysiq_mod, ONLY: iniphysiq
+#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 "comdissnew.h"
+  include "comgeom2.h"
+
+  REAL zdtvr
+
+
+  ! variables dynamiques intermediaire pour le transport
+
+  !   variables pour le fichier histoire
+  REAL dtav      ! intervalle de temps elementaire
+
+  REAL time_0
+
+  LOGICAL lafin
+
+  integer :: ntime=10000,it,klon,klev
+
+
+  ! locales pour gestion du temps
+  INTEGER :: an, mois, jour
+  REAL :: heure
+
+  !   Initialisations:
+  !   ----------------
+
+
+  !----------------------------------------------------------------------
+  !  lecture des fichiers gcm.def ou run.def
+  !  ---------------------------------------
+  !
+preff=101325.
+ pa=50000.
+ clon=0.
+ clat=0.
+ taux=3.
+ tauy=3.
+ dzoomx=0.1
+ dzoomy=0.1
+ grossismx=1.
+ grossismx=1.
+ transx=0.
+ transy=0.
+
+  CALL conf_gcm( 99, .TRUE.)
+
+  if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm", &
+       "iphysiq must be a multiple of iperiod", 1)
+
+ rad=6400000
+g=9.81
+
+
+
+  !-----------------------------------------------------------------------
+  !   Choix du calendrier
+  !   -------------------
+
+  !      calend = 'earth_365d'
+
+
+  !-----------------------------------------------------------------------
+  !   Lecture de l'etat initial :
+  !   ---------------------------
+
+  !-----------------------------------------------------------------------
+  !   Lecture des parametres de controle pour la simulation :
+  !   -------------------------------------------------------
+  !  on recalcule eventuellement le pas de temps
+
+
+  zdtvr    = daysec/REAL(day_step)
+
+  !
+  ! on remet le calendrier \`a zero si demande
+  !
+     annee_ref = anneeref
+     day_ref = dayref
+     day_ini = dayref
+     itau_dyn = 0
+     itau_phy = 0
+     time_0 = 0.
+
+  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)
+
+
+
+     dtvr = zdtvr
+     CALL iniconst
+     print*,'APRES inisconst'
+     CALL inigeom
+
+
+  CALL inigeomphy(iim,jjm,llm, &
+               1, comm_lmdz, &
+               rlatu,rlatv, &
+               rlonu,rlonv, &
+               aire,cu,cv)
+
+CALL suphel
+!open(82,file='dump_param.bin',form='unformatted',status='old')
+
+
+
+     CALL iophys_ini(900.)
+print*,'Rlatu=',rlatu
+klon=2+iim*(jjm-1)
+klev=llm
+
+!---------------------------------------------------------------------
+! Initialisation de la parametrisation
+!---------------------------------------------------------------------
+      call get_ini_module
+
+!---------------------------------------------------------------------
+! Boucle en temps sur l'appel à la parametrisation
+!---------------------------------------------------------------------
+      DO it=1,ntime
+         print*,'Pas de temps ',it,klon,klev
+         call call_param_replay(klon,klev)
+      ENDDO
+
+
+END PROGRAM replay3d
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/vlsplt.F
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/vlsplt.F	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/vlsplt.F	(revision 4368)
@@ -4,5 +4,5 @@
 
       SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
-      USE infotrac, ONLY: nqtot,nqdesc,iqfils
+      USE infotrac, ONLY: nqtot,tracers
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -38,27 +38,16 @@
 c   ---------
 c
-      INTEGER i,ij,l,j,ii
-      INTEGER ijlqmin,iqmin,jqmin,lqmin
-c
-      REAL zm(ip1jmp1,llm,nqtot),newmasse
+      INTEGER ij,l
+c
+      REAL zm(ip1jmp1,llm,nqtot)
       REAL mu(ip1jmp1,llm)
       REAL mv(ip1jm,llm)
       REAL mw(ip1jmp1,llm+1)
-      REAL zq(ip1jmp1,llm,nqtot),zz
-      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
-      REAL second,temps0,temps1,temps2,temps3
-      REAL ztemps1,ztemps2,ztemps3
+      REAL zq(ip1jmp1,llm,nqtot)
       REAL zzpbar, zzw
-      LOGICAL testcpu
-      SAVE testcpu
-      SAVE temps1,temps2,temps3
-      INTEGER iminn,imaxx
       INTEGER ifils,iq2 ! CRisi
 
       REAL qmin,qmax
       DATA qmin,qmax/0.,1.e33/
-      DATA testcpu/.false./
-      DATA temps1,temps2,temps3/0.,0.,0./
-
 
         zzpbar = 0.5 * pdt
@@ -83,10 +72,8 @@
       CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
         
-      if (nqdesc(iq).gt.0) then  
-        do ifils=1,nqdesc(iq)
-          iq2=iqfils(ifils,iq)
-          CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
-        enddo  
-      endif !if (nqfils(iq).gt.0) then
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
+      enddo  
 
 cprint*,'Entree vlx1'
@@ -122,23 +109,21 @@
       ENDDO
       ! CRisi: aussi pour les fils
-      if (nqdesc(iq).gt.0) then
-      do ifils=1,nqdesc(iq)
-        iq2=iqfils(ifils,iq)
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
         DO l=1,llm
-         DO ij=1,ip1jmp1
-           q(ij,l,iq2)=zq(ij,l,iq2)
-         ENDDO
-         DO ij=1,ip1jm+1,iip1
+          DO ij=1,ip1jmp1
+            q(ij,l,iq2)=zq(ij,l,iq2)
+          ENDDO
+          DO ij=1,ip1jm+1,iip1
             q(ij+iim,l,iq2)=q(ij,l,iq2)
-         ENDDO
+          ENDDO
         ENDDO
-      enddo !do ifils=1,nqdesc(iq)    
-      endif ! if (nqdesc(iq).gt.0) then    
+      enddo
 
       RETURN
       END
       RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
-      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils, ! CRisi
-     &                     qperemin,masseqmin,ratiomin ! MVals et CRisi
+      USE infotrac, ONLY : nqtot,tracers, ! CRisi
+     &                     min_qParent,min_qMass,min_ratio ! MVals et CRisi
 
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -161,7 +146,6 @@
 c   ----------
       REAL masse(ip1jmp1,llm,nqtot),pente_max
-      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
+      REAL u_m( ip1jmp1,llm )
       REAL q(ip1jmp1,llm,nqtot)
-      REAL w(ip1jmp1,llm)
       INTEGER iq ! CRisi
 c
@@ -173,5 +157,6 @@
 c
       REAL new_m,zu_m,zdum(ip1jmp1,llm)
-      REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
+c      REAL sigu(ip1jmp1)
+      REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
       REAL zz(ip1jmp1)
       REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
@@ -182,23 +167,7 @@
       INTEGER ifils,iq2 ! CRisi
 
-      Logical extremum,first,testcpu
-      SAVE first,testcpu
-
-      REAL      SSUM
-      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
-      SAVE temps0,temps1,temps2,temps3,temps4,temps5
-
-      REAL z1,z2,z3
-
-      DATA first,testcpu/.true.,.false./
-
-      IF(first) THEN
-         temps1=0.
-         temps2=0.
-         temps3=0.
-         temps4=0.
-         temps5=0.
-         first=.false.
-      ENDIF
+      Logical first
+      SAVE first
+      DATA first/.true./
 
 c   calcul de la pente a droite et a gauche de la maille
@@ -436,5 +405,5 @@
          ENDDO
       ENDIF  ! n0.gt.0 
-9999    continue
+c9999    continue
 
 
@@ -449,29 +418,27 @@
 ! 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,nqfils(iq)=',iq,nqfils(iq) 
+      !write(*,*) 'vlsplt 326: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen
       
-      if (nqdesc(iq).gt.0) then  
-       do ifils=1,nqdesc(iq)
-         iq2=iqfils(ifils,iq)
-         DO l=1,llm
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
           DO ij=iip2,ip1jm
-           ! On a besoin de q et masse seulement entre iip2 et ip1jm 
-           !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
-           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
-           !Mvals: veiller a ce qu'on n'ait pas de denominateur nul
-           masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)
-           if (q(ij,l,iq).gt.qperemin) then
-             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
-           else
-             Ratio(ij,l,iq2)=ratiomin
-           endif
+            ! On a besoin de q et masse seulement entre iip2 et ip1jm 
+            !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+            !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+            !Mvals: veiller a ce qu'on n'ait pas de denominateur nul
+            masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
+            if (q(ij,l,iq).gt.min_qParent) then
+              Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+            else
+              Ratio(ij,l,iq2)=min_ratio
+            endif
           enddo   
-         enddo
-        enddo !do ifils=1,nqdesc(iq)
-        do ifils=1,nqfils(iq)
-         iq2=iqfils(ifils,iq)
-         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
-        enddo !do ifils=1,nqfils(iq)
-      endif !if (nqfils(iq).gt.0) then
+        enddo
+      enddo
+      do ifils=1,tracers(iq)%nqChildren
+        iq2=tracers(iq)%iqDescen(ifils)
+        call vlx(Ratio,pente_max,masseq,u_mq,iq2)
+      enddo
 ! end CRisi
 
@@ -482,5 +449,5 @@
          DO ij=iip2+1,ip1jm
             !MVals: veiller a ce qu'on ait pas de denominateur nul
-            new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),masseqmin)
+            new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),min_qMass)
             q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
      &      u_mq(ij-1,l)-u_mq(ij,l))
@@ -498,8 +465,7 @@
       ! On calcule q entre iip2+1,ip1jm -> 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)  
-         DO l=1,llm
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
           DO ij=iip2+1,ip1jm
             q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
@@ -508,7 +474,6 @@
              q(ij-iim,l,iq2)=q(ij,l,iq2)
           enddo ! DO ij=ijb+iip1-1,ije,iip1
-         enddo !DO l=1,llm
-        enddo !do ifils=1,nqdesc(iq)
-      endif !if (nqfils(iq).gt.0) then
+        enddo !DO l=1,llm
+      enddo
 
 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
@@ -519,6 +484,6 @@
       END
       RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
-      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils, ! CRisi
-     &                     qperemin,masseqmin,ratiomin ! MVals et CRisi
+      USE infotrac, ONLY : nqtot,tracers, ! CRisi
+     &                     min_qParent,min_qMass,min_ratio ! MVals et CRisi
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -544,5 +509,5 @@
       REAL masse(ip1jmp1,llm,nqtot),pente_max
       REAL masse_adv_v( ip1jm,llm)
-      REAL q(ip1jmp1,llm,nqtot), dq( ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot)
       INTEGER iq ! CRisi
 c
@@ -553,14 +518,13 @@
 c
       REAL airej2,airejjm,airescb(iim),airesch(iim)
-      REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
+      REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
       REAL adyqv(ip1jm),dyqmax(ip1jmp1)
       REAL qbyv(ip1jm,llm)
 
-      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
+      REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
+c     REAL appn apps
 c     REAL newq,oldmasse
-      Logical extremum,first,testcpu
-      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
-      SAVE temps0,temps1,temps2,temps3,temps4,temps5
-      SAVE first,testcpu
+      LOGICAL first
+      SAVE first
 
       REAL convpn,convps,convmpn,convmps
@@ -578,6 +542,5 @@
       REAL      SSUM
 
-      DATA first,testcpu/.true.,.false./
-      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
+      DATA first/.true./
 
       !write(*,*) 'vly 578: entree, iq=',iq
@@ -778,31 +741,29 @@
 ! 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) 
+      !write(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen
    
-      if (nqfils(iq).gt.0) then  
-       do ifils=1,nqdesc(iq)
-         iq2=iqfils(ifils,iq)
-         DO l=1,llm
-         DO ij=1,ip1jmp1
-           ! attention, chaque fils doit avoir son masseq, sinon, le 1er
-           ! fils ecrase le masseq de ses freres.
-           !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
-           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
-           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
-           masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)
-           if (q(ij,l,iq).gt.qperemin) then
-             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
-           else
-             Ratio(ij,l,iq2)=ratiomin
-           endif 
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
+          DO ij=1,ip1jmp1
+            ! attention, chaque fils doit avoir son masseq, sinon, le 1er
+            ! fils ecrase le masseq de ses freres.
+            !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+            !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
+            !MVals: veiller a ce qu'on n'ait pas de denominateur nul
+            masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
+            if (q(ij,l,iq).gt.min_qParent) then
+              Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+            else
+              Ratio(ij,l,iq2)=min_ratio
+            endif 
           enddo   
-         enddo
-        enddo !do ifils=1,nqdesc(iq)
-
-        do ifils=1,nqfils(iq)
-         iq2=iqfils(ifils,iq)
-         call vly(Ratio,pente_max,masseq,qbyv,iq2)
-        enddo !do ifils=1,nqfils(iq)
-      endif !if (nqfils(iq).gt.0) then
+        enddo
+      enddo
+
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        call vly(Ratio,pente_max,masseq,qbyv,iq2)
+      enddo
 
       DO l=1,llm
@@ -872,14 +833,12 @@
   
 ! 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)  
-         DO l=1,llm
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
           DO ij=1,ip1jmp1
             q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
           enddo
-         enddo
-        enddo !do ifils=1,nqdesc(iq)
-      endif !if (nqfils(iq).gt.0) then
+        enddo
+      enddo
 
       !write(*,*) 'vly 853: sortie'
@@ -888,6 +847,6 @@
       END
       RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq)
-      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils, ! CRisi 
-     &                     qperemin,masseqmin,ratiomin ! MVals et CRisi
+      USE infotrac, ONLY : nqtot,tracers, ! CRisi 
+     &                     min_qParent,min_qMass,min_ratio ! MVals et CRisi
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -917,5 +876,5 @@
 c   ---------
 c
-      INTEGER i,ij,l,j,ii
+      INTEGER ij,l
 c
       REAL wq(ip1jmp1,llm+1),newmasse
@@ -930,10 +889,11 @@
       SAVE testcpu
 
-      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
-      SAVE temps0,temps1,temps2,temps3,temps4,temps5
-      REAL      SSUM
+#ifdef BIDON
+      REAL temps0,temps1,second
+      SAVE temps0,temps1
 
       DATA testcpu/.false./
-      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
+      DATA temps0,temps1/0.,0./
+#endif
 
 c    On oriente tout dans le sens de la pression c'est a dire dans le
@@ -1009,28 +969,26 @@
 ! 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)
-         DO l=1,llm
+      !write(*,*) 'vlsplt 942: iq,nqChildren(iq)=',iq,nqChildren(iq)
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
           DO ij=1,ip1jmp1
-           !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
-           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)        
-           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
-           masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)
-           if (q(ij,l,iq).gt.qperemin) then
-             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
-           else
-             Ratio(ij,l,iq2)=ratiomin
-           endif      
+            !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+            !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)        
+            !MVals: veiller a ce qu'on n'ait pas de denominateur nul
+            masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
+            if (q(ij,l,iq).gt.min_qParent) then
+              Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+            else
+              Ratio(ij,l,iq2)=min_ratio
+            endif      
           enddo   
-         enddo
-        enddo !do ifils=1,nqdesc(iq)
+        enddo
+      enddo
         
-        do ifils=1,nqfils(iq)
-         iq2=iqfils(ifils,iq)         
-         call vlz(Ratio,pente_max,masseq,wq,iq2)
-        enddo !do ifils=1,nqfils(iq)
-      endif !if (nqfils(iq).gt.0) then
+      do ifils=1,tracers(iq)%nqChildren
+        iq2=tracers(iq)%iqDescen(ifils)
+        call vlz(Ratio,pente_max,masseq,wq,iq2)
+      enddo
 ! end CRisi  
 
@@ -1045,14 +1003,12 @@
 
 ! 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)  
-         DO l=1,llm
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
           DO ij=1,ip1jmp1
             q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
           enddo
-         enddo
-        enddo !do ifils=1,nqdesc(iq)
-      endif !if (nqfils(iq).gt.0) then
+        enddo
+      enddo
       !write(*,*) 'vlsplt 1032'
 
@@ -1099,4 +1055,5 @@
       real zzq(iip1,jjp1,llm)
 
+#ifdef isminmax
       integer imin,jmin,lmin,ijlmin
       integer imax,jmax,lmax,ijlmax
@@ -1104,5 +1061,4 @@
       integer ismin,ismax
 
-#ifdef isminismax
       call scopy (ip1jmp1*llm,zq,1,zzq,1)
 
@@ -1132,5 +1088,5 @@
 #endif
       return
-9999  format(a20,'  q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
+c9999  format(a20,'  q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
       end
 
Index: LMDZ6/branches/Ocean_skin/libf/dyn3d/vlspltqs.F
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/dyn3d/vlspltqs.F	(revision 4013)
+++ LMDZ6/branches/Ocean_skin/libf/dyn3d/vlspltqs.F	(revision 4368)
@@ -4,5 +4,5 @@
        SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
      ,                                  p,pk,teta,iq             )
-       USE infotrac, ONLY: nqtot,nqdesc,iqfils
+       USE infotrac, ONLY: nqtot,tracers
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron 
@@ -121,10 +121,8 @@
       CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
       CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
-      if (nqdesc(iq).gt.0) then  
-       do ifils=1,nqdesc(iq)
-        iq2=iqfils(ifils,iq)
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
         CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
-       enddo  
-      endif !if (nqfils(iq).gt.0) then
+      enddo  
 
 c      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
@@ -162,17 +160,15 @@
       ENDDO
       ! CRisi: aussi pour les fils
-      if (nqdesc(iq).gt.0) then
-      do ifils=1,nqdesc(iq)
-        iq2=iqfils(ifils,iq)
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
         DO l=1,llm
-         DO ij=1,ip1jmp1
-           q(ij,l,iq2)=zq(ij,l,iq2)
-         ENDDO
-         DO ij=1,ip1jm+1,iip1
+          DO ij=1,ip1jmp1
+            q(ij,l,iq2)=zq(ij,l,iq2)
+          ENDDO
+          DO ij=1,ip1jm+1,iip1
             q(ij+iim,l,iq2)=q(ij,l,iq2)
-         ENDDO
+          ENDDO
         ENDDO
-      enddo !do ifils=1,nqdesc(iq)  
-      endif ! if (nqfils(iq).gt.0) then
+      enddo 
       !write(*,*) 'vlspltqs 183: fin de la routine'
 
@@ -180,5 +176,5 @@
       END
       SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq)
-      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
+      USE infotrac, ONLY : nqtot,tracers ! CRisi
 
 c
@@ -483,22 +479,21 @@
 ! CRisi: appel récursif de l'advection sur les fils.
 ! Il faut faire ça avant d'avoir mis à jour q et masse
-      !write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq) 
+      !write(*,*) 'vlspltqs 326: iq,nqChildren(iq)=',iq,
+!     &                 tracers(iq)%nqChildren
       
-      if (nqfils(iq).gt.0) then  
-       do ifils=1,nqdesc(iq)
-         iq2=iqfils(ifils,iq)
-         DO l=1,llm
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
           DO ij=iip2,ip1jm
-           ! On a besoin de q et masse seulement entre iip2 et ip1jm 
-           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
-           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+          ! On a besoin de q et masse seulement entre iip2 et ip1jm 
+            masseq(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
-        enddo !do ifils=1,nqdesc(iq)
-        do ifils=1,nqfils(iq)
-         iq2=iqfils(ifils,iq)
-         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
-        enddo !do ifils=1,nqfils(iq)
-      endif !if (nqfils(iq).gt.0) then
+        enddo
+      enddo
+      do ifils=1,tracers(iq)%nqChildren
+        iq2=tracers(iq)%iqDescen(ifils)
+        call vlx(Ratio,pente_max,masseq,u_mq,iq2)
+      enddo
 ! end CRisi
 
@@ -523,17 +518,15 @@
       ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 
       ! puis on boucle en longitude
-      if (nqdesc(iq).gt.0) then  
-       do ifils=1,nqdesc(iq)
-         iq2=iqfils(ifils,iq)  
-         DO l=1,llm
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
           DO ij=iip2+1,ip1jm
             q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
           enddo
           DO ij=iip1+iip1,ip1jm,iip1
-             q(ij-iim,l,iq2)=q(ij,l,iq2)
-          enddo ! DO ij=ijb+iip1-1,ije,iip1
-         enddo !DO l=1,llm
-        enddo !do ifils=1,nqdesc(iq)
-      endif !if (nqfils(iq).gt.0) then
+            q(ij-iim,l,iq2)=q(ij,l,iq2)
+          enddo
+        enddo
+      enddo
 
 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
@@ -544,5 +537,5 @@
       END
       SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
-      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
+      USE infotrac, ONLY : nqtot,tracers ! CRisi
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -794,23 +787,21 @@
 ! 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) 
+      !write(*,*) 'vlyqs 689: iq,nqChildren(iq)=',iq,
+!     &              tracers(iq)%nqChildren
    
-      if (nqfils(iq).gt.0) then  
-       do ifils=1,nqdesc(iq)
-         iq2=iqfils(ifils,iq)
-         DO l=1,llm
-         DO ij=1,ip1jmp1
-           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
-           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
+          DO ij=1,ip1jmp1
+            masseq(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
-        enddo !do ifils=1,nqdesc(iq)
-
-        do ifils=1,nqfils(iq)
-         iq2=iqfils(ifils,iq)
-         !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
-         call vly(Ratio,pente_max,masseq,qbyv,iq2)
-        enddo !do ifils=1,nqfils(iq)
-      endif !if (nqfils(iq).gt.0) then
+        enddo
+      enddo
+      do ifils=1,tracers(iq)%nqChildren
+        iq2=tracers(iq)%iqDescen(ifils)
+        !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
+        call vly(Ratio,pente_max,masseq,qbyv,iq2)
+      enddo
 
       DO l=1,llm
@@ -868,14 +859,12 @@
 
 ! retablir les fils en rapport de melange par rapport a l'air:
-      if (nqdesc(iq).gt.0) then  
-       do ifils=1,nqdesc(iq)
-         iq2=iqfils(ifils,iq)  
-         DO l=1,llm
+      do ifils=1,tracers(iq)%nqDescen
+        iq2=tracers(iq)%iqDescen(ifils)
+        DO l=1,llm
           DO ij=1,ip1jmp1
             q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
           enddo
-         enddo
-        enddo !do ifils=1,nqdesc(iq)
-      endif !if (nqfils(iq).gt.0) then
+        enddo
+      enddo
       !write(*,*) 'vly 879'
 
