!
! $Id: advtrac_p.F 1403 2010-07-01 09:02:53Z asima $
!
c
c
      SUBROUTINE advtrac_p(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 parallel
      USE Write_Field_p
      USE Bands
      USE mod_hallo
      USE Vampir
      USE times
      USE infotrac
      USE control_mod
      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 "ener.h"
#include "description.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,SAVE::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
      SAVE iadvtr, massem, pbaruc, pbarvc
      DATA iadvtr/0/
c$OMP THREADPRIVATE(iadvtr)
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./
      REAL,SAVE :: finmasse(ip1jmp1,llm)
      integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
      type(Request) :: Request_vanleer
      REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
      REAL,SAVE :: teta_tmp(ip1jmp1,llm)
      REAL,SAVE :: pk_tmp(ip1jmp1,llm)

      ijb_u=ij_begin
      ije_u=ij_end
      
      ijb_v=ij_begin-iip1
      ije_v=ij_end
      if (pole_nord) ijb_v=ij_begin
      if (pole_sud)  ije_v=ij_end-iip1

      IF(iadvtr.EQ.0) THEN
c         CALL initial0(ijp1llm,pbaruc)
c         CALL initial0(ijmllm,pbarvc)
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)        
	DO l=1,llm   
          pbaruc(ijb_u:ije_u,l)=0.
          pbarvc(ijb_v:ije_v,l)=0.
        ENDDO
c$OMP END DO NOWAIT  
      ENDIF

c   accumulation des flux de masse horizontaux
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
         DO ij = ijb_u,ije_u
            pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
         ENDDO
         DO ij = ijb_v,ije_v
            pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
         ENDDO
      ENDDO
c$OMP END DO NOWAIT

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

c         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
          ijb=ij_begin
          ije=ij_end

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
          massem(ijb:ije,l)=masse(ijb:ije,l)
       ENDDO
c$OMP END DO NOWAIT

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

      iadvtr   = iadvtr+1

c$OMP MASTER
      iapptrac = iadvtr
c$OMP END MASTER

c   Test pour savoir si on advecte a ce pas de temps

      IF ( iadvtr.EQ.iapp_tracvl ) THEN
c$OMP MASTER
        call suspend_timer(timer_caldyn)
c$OMP END MASTER
      
      ijb=ij_begin
      ije=ij_end
      

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_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )

c$OMP BARRIER

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llmp1
        p_tmp(ijb:ije,l)=p(ijb:ije,l)
      ENDDO
c$OMP END DO NOWAIT
      
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l=1,llm
        pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
        teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
      ENDDO
c$OMP END DO NOWAIT

c$OMP MASTER
      call VTb(VTHallo)
c$OMP END MASTER

      call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm,
     *                             jj_Nb_vanleer,0,0,Request_vanleer)
      call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm,
     *                             jj_Nb_vanleer,1,0,Request_vanleer)
      call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm,
     *                             jj_Nb_vanleer,0,0,Request_vanleer)
      call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm,
     *                             jj_Nb_vanleer,0,0,Request_vanleer)
      call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm,
     *                             jj_Nb_vanleer,1,1,Request_vanleer)
      call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1,
     *                             jj_Nb_vanleer,1,1,Request_vanleer)
      call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,
     *                             jj_Nb_vanleer,1,1,Request_vanleer)
      do j=1,nqtot
        call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
     *                             jj_nb_vanleer,0,0,Request_vanleer)
      enddo

      call SendRequest(Request_vanleer)
c$OMP BARRIER
      call WaitRequest(Request_vanleer)


c$OMP BARRIER
c$OMP MASTER      
      call SetDistrib(jj_nb_vanleer)
      call VTe(VTHallo)
      call VTb(VTadvection)
      call start_timer(timer_vanleer)
c$OMP END MASTER
c$OMP BARRIER
      
      ! ... Flux de masse diaganostiques traceurs
         ijb=ij_begin
         ije=ij_end
         flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)

c  test sur l'eventuelle creation de valeurs negatives de la masse
         ijb=ij_begin
         ije=ij_end
         if (pole_nord) ijb=ij_begin+iip1
         if (pole_sud) ije=ij_end-iip1
         
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)          
         DO l=1,llm-1
            DO ij = ijb+1,ije
              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
            
c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
            
            do ij=ijb,ije-iip1+1,iip1
              zdp(ij)=zdp(ij+iip1-1)
            enddo
            
            DO ij = ijb,ije
               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 
            ENDDO 


c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
c  ym ---> eventuellement a revoir
            CALL minmax ( ije-ijb+1, zdp(ijb), 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$OMP END DO NOWAIT

c-------------------------------------------------------------------
c   Advection proprement dite (Modification Le Croller (07/2001)
c-------------------------------------------------------------------

c----------------------------------------------------
c        Calcul des moyennes basées sur la masse
c----------------------------------------------------

cym      ----> Normalement, inutile pour les schémas classiques
cym      ----> Revérifier lors de la parallélisation des autres schemas
   
cym          call massbar_p(massem,massebx,masseby)          

          call vlspltgen_p( q,iadv, 2., massem, wg ,
     *                    pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )

         
	  GOTO 1234     
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_p(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
cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
           CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
     *                 pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
c   ----------------------------------------------------------------
c   Schema de Frederic Hourdin
c   ----------------------------------------------------------------
        else if(iadv(iq).eq.12) then
          stop 'advtrac : schema non parallelise'
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
          stop 'advtrac : schema non parallelise'
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
          stop 'advtrac : schema non parallelise'

            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
          stop 'advtrac : schema non parallelise'
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

           stop 'advtrac : schema non parallelise'

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

1234  CONTINUE
c$OMP BARRIER

      ijb=ij_begin
      ije=ij_end

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
       DO l = 1, llm
         DO ij = ijb, ije
           finmasse(ij,l) =  p(ij,l) - p(ij,l+1) 
         ENDDO
       ENDDO
c$OMP END DO

       CALL qminimum_p( q, 2, finmasse )

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

c$OMP MASTER
	call VTe(VTadvection)
        call stop_timer(timer_vanleer)
        call VTb(VThallo)
c$OMP END MASTER

	do j=1,nqtot
          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
     *                             jj_nb_caldyn,0,0,Request_vanleer)
        enddo

        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
     *       jj_nb_caldyn,0,0,Request_vanleer)

        call SendRequest(Request_vanleer)
c$OMP BARRIER
        call WaitRequest(Request_vanleer)      

c$OMP BARRIER
c$OMP MASTER
        call SetDistrib(jj_nb_caldyn)
	call VTe(VThallo)
	call resume_timer(timer_caldyn)
c$OMP END MASTER
c$OMP BARRIER	
          iadvtr=0
       ENDIF ! if iadvtr.EQ.iapp_tracvl

       RETURN
       END