Index: LMDZ6/trunk/libf/dyn3d/advtrac.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d/advtrac.F90	(revision 4050)
+++ LMDZ6/trunk/libf/dyn3d/advtrac.F90	(revision 4056)
@@ -1,407 +1,358 @@
 ! $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, tracers, 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)
-  INTEGER :: iadv
-
-  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,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
+   !---------------------------------------------------------------------------
+   INTEGER :: ij, l, iq, iiq, 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)
-       iadv = tracers(iq)%iadv
-       SELECT CASE(iadv)
+   ! 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(ok_iso_verif) THEN
+      WRITE(*,*) 'advtrac 227'
+      CALL check_isotopes_seq(q,ip1jmp1,'advtrac 162')
+   END IF
+
+   !-------------------------------------------------------------------------
+   !       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
-        !   ----------------------------------------------------------------
-           ! 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)
-           
-         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
-           !   ----------------------------------------------------------------
-           !            Pas de temps adaptatif
-           call adaptdt(iadv,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
-         CASE(13)
-           !            Pas de temps adaptatif
-           call adaptdt(iadv,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
-         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.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)
-
-         CASE(11,16,17,18)
-           !   ----------------------------------------------------------------
-           !   Schemas PPM Lin et Rood
-           !   ----------------------------------------------------------------
-
-           !        Test sur le flux horizontal
-           !        Pas de temps adaptatif
-           call adaptdt(iadv,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 (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
-           enddo
-           !-----------------------------------------------------------------
-           !               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 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
+         !--------------------------------------------------------------------
+         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(ok_iso_verif) then
+      WRITE(*,*) 'advtrac 402'
+      CALL check_isotopes_seq(q,ip1jmp1,'advtrac 397')
+   END IF
+
+   !-------------------------------------------------------------------------
+   !   on reinitialise a zero les flux de masse cumules
+   !-------------------------------------------------------------------------
+   iadvtr=0
 
 END SUBROUTINE advtrac
