source: LMDZ6/trunk/libf/dyn3d/advtrac.F90 @ 5161

! $Id: advtrac.F90 4143 2022-05-09 10:35:40Z aborella $

#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)
   ! 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