!
! $Id: advtrac.F 1176 2009-06-11 08:54:10Z oboucher $
!
c
c
      SUBROUTINE advtrac(pbaru,pbarv ,
     *                   p,  masse,q,iapptrac,teta,
     *                  flxw,
     *                  pk)
c     Auteur :  F. Hourdin
c
c     Modif. P. Le Van     (20/12/97)
c            F. Codron     (10/99)
c            D. Le Croller (07/2001)
c            M.A Filiberti (04/2002)
c
      USE infotrac

      IMPLICIT NONE
c
#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
#include "comvert.h"
#include "comdissip.h"
#include "comgeom2.h"
#include "logic.h"
#include "temps.h"
#include "control.h"
#include "ener.h"
#include "description.h"
#include "iniprint.h"

c-------------------------------------------------------------------
c     Arguments
c-------------------------------------------------------------------
c     Ajout PPM
c--------------------------------------------------------
      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
c--------------------------------------------------------
      INTEGER iapptrac
      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
      REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
      REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
      REAL pk(ip1jmp1,llm)
      REAL flxw(ip1jmp1,llm)

c-------------------------------------------------------------
c     Variables locales
c-------------------------------------------------------------

      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/
c----------------------------------------------------------
c     Rajouts pour PPM
c----------------------------------------------------------
      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
         CALL initial0(ijp1llm,pbaruc)
         CALL initial0(ijmllm,pbarvc)
      ENDIF

c   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

c   selection de la masse instantannee des mailles avant le transport.
      IF(iadvtr.EQ.0) THEN

         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
c
      ENDIF

      iadvtr   = iadvtr+1
      iapptrac = iadvtr


c   Test pour savoir si on advecte a ce pas de temps
      IF ( iadvtr.EQ.iapp_tracvl ) THEN

cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
cc

c   traitement des flux de masse avant advection.
c     1. calcul de w
c     2. groupement des mailles pres du pole.

        CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )

      ! ... Flux de masse diaganostiques traceurs
      flxw = wg / FLOAT(iapp_tracvl)

c  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)
     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
     s                  +       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,
     s        '   MAX:', zdpmax
            ENDIF

         ENDDO


c-------------------------------------------------------------------
! Calcul des criteres CFL en X, Y et Z
c-------------------------------------------------------------------

      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, CFLmax '
     s   ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
         enddo
         countcfl=0
      endif
    
c-------------------------------------------------------------------
c   Advection proprement dite (Modification Le Croller (07/2001)
c-------------------------------------------------------------------

c----------------------------------------------------
c        Calcul des moyennes basées sur la masse
c----------------------------------------------------
          call massbar(massem,massebx,masseby)          

c-----------------------------------------------------------
c     Appel des sous programmes d'advection
c-----------------------------------------------------------
      do iq=1,nqtot
c        call clock(t_initial)
        if(iadv(iq) == 0) cycle 
c   ----------------------------------------------------------------
c   Schema de Van Leer I MUSCL
c   ----------------------------------------------------------------
        if(iadv(iq).eq.10) THEN
            call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)


c   ----------------------------------------------------------------
c   Schema "pseudo amont" + test sur humidite specifique
C    pour la vapeur d'eau. F. Codron
c   ----------------------------------------------------------------
        else if(iadv(iq).eq.14) then
c
           CALL vlspltqs( q(1,1,1), 2., massem, wg ,
     *                 pbarug,pbarvg,dtvr,p,pk,teta )
c   ----------------------------------------------------------------
c   Schema de Frederic Hourdin
c   ----------------------------------------------------------------
        else if(iadv(iq).eq.12) then
c            Pas de temps adaptatif
           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
           if (n.GT.1) then
           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
     s             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
c            Pas de temps adaptatif
           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
           if (n.GT.1) then
           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
     s             dtvr,'n=',n
           endif
          do indice=1,n
            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
          end do
c   ----------------------------------------------------------------
c   Schema de pente SLOPES
c   ----------------------------------------------------------------
        else if (iadv(iq).eq.20) then
            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)

c   ----------------------------------------------------------------
c   Schema de Prather
c   ----------------------------------------------------------------
        else if (iadv(iq).eq.30) then
c            Pas de temps adaptatif
           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
           if (n.GT.1) then
           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
     s             dtvr,'n=',n
           endif
           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
     s                     n,dtbon)

c   ----------------------------------------------------------------
c   Schemas PPM Lin et Rood
c   ----------------------------------------------------------------
         else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
     s                     iadv(iq).LE.18)) then

c        Test sur le flux horizontal
c        Pas de temps adaptatif
         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
         if (n.GT.1) then
         write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
     s             dtvr,'n=',n
         endif
c        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

c-----------------------------------------------------------
c        Ss-prg interface LMDZ.4->PPM3d
c-----------------------------------------------------------

          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
     s                 unatppm,vnatppm,psppm)

          do indice=1,n
c---------------------------------------------------------------------
c                         VL (version PPM) horiz. et PPM vert.
c---------------------------------------------------------------------
                if (iadv(iq).eq.11) then
c                  Ss-prg PPM3d de Lin
                  call ppm3d(1,qppm(1,1,iq),
     s                       psppm,psppm,
     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
     s                       fill,dum,220.)

c----------------------------------------------------------------------
c                           Monotonic PPM
c----------------------------------------------------------------------
               else if (iadv(iq).eq.16) then
c                  Ss-prg PPM3d de Lin
                  call ppm3d(1,qppm(1,1,iq),
     s                       psppm,psppm,
     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
     s                       fill,dum,220.)
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c                           Semi Monotonic PPM
c---------------------------------------------------------------------
               else if (iadv(iq).eq.17) then
c                  Ss-prg PPM3d de Lin
                  call ppm3d(1,qppm(1,1,iq),
     s                       psppm,psppm,
     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
     s                       fill,dum,220.)
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c                         Positive Definite PPM
c---------------------------------------------------------------------
                else if (iadv(iq).eq.18) then
c                  Ss-prg PPM3d de Lin
                  call ppm3d(1,qppm(1,1,iq),
     s                       psppm,psppm,
     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
     s                       fill,dum,220.)
c---------------------------------------------------------------------
                endif
            enddo
c-----------------------------------------------------------------
c               Ss-prg interface PPM3d-LMDZ.4
c-----------------------------------------------------------------
                  call interpost(q(1,1,iq),qppm(1,1,iq))
            endif
c----------------------------------------------------------------------

c-----------------------------------------------------------------
c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
c et Nord j=1
c-----------------------------------------------------------------

c                  call traceurpole(q(1,1,iq),massem)

c calcul du temps cpu pour un schema donne

c                  call clock(t_final)
cym                  tps_cpu=t_final-t_initial
cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu

       end DO


c------------------------------------------------------------------
c   on reinitialise a zero les flux de masse cumules
c---------------------------------------------------
          iadvtr=0

       ENDIF ! if iadvtr.EQ.iapp_tracvl

       RETURN
       END