source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advx.f90 @ 5209

! $Header$

SUBROUTINE  advx(limit,dtx,pbaru,sm,s0, &
        sx,sy,sz,lati,latf)
  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
  USE lmdz_paramet
  IMPLICIT NONE

  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !                                                                C
  !  first-order moments (FOM) advection of tracer in X direction  C
  !                                                                C
  !  Source : Pascal Simon (Meteo,CNRM)                            C
  !  Adaptation : A.Armengaud (LGGE) juin 94                       C
  !                                                                C
  !  limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz                       C
  !  sont des arguments d'entree pour le s-pg...                   C
  !                                                                C
  !  sm,s0,sx,sy,sz                                                C
  !  sont les arguments de sortie pour le s-pg                     C
  !                                                                C
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

  !  parametres principaux du modele
  !



  !  Arguments :
  !  -----------
  !  dtx : frequence fictive d'appel du transport
  !  pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1

   INTEGER :: ntra
   PARAMETER (ntra = 1)

  ! ATTENTION partout ou on trouve ntra, insertion de boucle
        ! possible dans l'avenir.

  REAL :: dtx
  REAL :: pbaru ( iip1,jjp1,llm )

  !  moments: SM  total mass in each grid box
        ! S0  mass of tracer in each grid box
        ! Si  1rst order moment in i direction

  REAL :: SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
  REAL :: sx(iip1,jjp1,llm,ntra) &
        ,sy(iip1,jjp1,llm,ntra)
  REAL :: sz(iip1,jjp1,llm,ntra)

  !  Local :
  !  -------

  !  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
  !  mass fluxes in kg
  !  declaration :

  REAL :: UGRI(iip1,jjp1,llm)

  !  Rem : VGRI et WGRI ne sont pas utilises dans
  !  cette SUBROUTINE ( advection en x uniquement )

  !  Ti are the moments for the current latitude and level

  REAL :: TM(iim)
  REAL :: T0(iim,ntra),TX(iim,ntra)
  REAL :: TY(iim,ntra),TZ(iim,ntra)
  REAL :: TEMPTM                ! just a temporary variable

  !  the moments F are similarly defined and used as temporary
  !  storage for portions of the grid boxes in transit

  REAL :: FM(iim)
  REAL :: F0(iim,ntra),FX(iim,ntra)
  REAL :: FY(iim,ntra),FZ(iim,ntra)

  !  work arrays

  REAL :: ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)

  REAL :: SMNEW(iim),UEXT(iim)

  REAL :: sqi,sqf

  LOGICAL :: LIMIT
  INTEGER :: NUM(jjp1),LONK,NUMK
  INTEGER :: lon,lati,latf,niv
  INTEGER :: i,i2,i3,j,jv,l,k,itrac

  lon = iim
  niv = llm

  ! *** Test de passage d'arguments ******


  !  -------------------------------------
  DO j = 1,jjp1
     NUM(j) = 1
  END DO
  sqi = 0.
  sqf = 0.

  DO l = 1,llm
     DO j = 1,jjp1
        DO i = 1,iim
  !IM 240305            sqi = sqi + S0(i,j,l,9)
           sqi = sqi + S0(i,j,l,ntra)
        ENDDO
     ENDDO
  ENDDO
  PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
  PRINT*,'sqi=',sqi


  !  Interface : adaptation nouveau modele
  !  -------------------------------------

  !  ---------------------------------------------------------
  !  Conversion des flux de masses en kg/s
  !  pbaru est en N/s d'ou :
  !  ugri est en kg/s

  DO l = 1,llm
     DO j = 1,jjm+1
        DO i = 1,iip1
         ! ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
         ugri (i,j,llm+1-l) = pbaru (i,j,l)
  END DO
  END DO
  END DO


  !  ---------------------------------------------------------
  !  ---------------------------------------------------------
  !  ---------------------------------------------------------

  !  start here

  !  boucle principale sur les niveaux et les latitudes

  DO L=1,NIV
  DO K=lati,latf

  !  initialisation

  !  program assumes periodic boundaries in X

  DO I=2,LON
     SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
  END DO
  SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX

  !  modifications for extended polar zones

  NUMK=NUM(K)
  LONK=LON/NUMK

  IF(NUMK>1) THEN

  DO I=1,LON
     TM(I)=0.
  END DO
  DO JV=1,NTRA
  DO I=1,LON
     T0(I,JV)=0.
     TX(I,JV)=0.
     TY(I,JV)=0.
     TZ(I,JV)=0.
  END DO
  END DO

  DO I2=1,NUMK

     DO I=1,LONK
        I3=(I-1)*NUMK+I2
        TM(I)=TM(I)+SM(I3,K,L)
        ALF(I)=SM(I3,K,L)/TM(I)
        ALF1(I)=1.-ALF(I)
  END DO

     DO  JV=1,NTRA
     DO  I=1,LONK
        I3=(I-1)*NUMK+I2
        TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I) &
              *S0(I3,K,L,JV)
        T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
        TX(I,JV)=ALF(I)  *sx(I3,K,L,JV)+ &
              ALF1(I)*TX(I,JV) +3.*TEMPTM
        TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
        TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
     ENDDO
     ENDDO

  END DO

  ELSE

  DO I=1,LON
     TM(I)=SM(I,K,L)
  END DO
  DO JV=1,NTRA
  DO I=1,LON
     T0(I,JV)=S0(I,K,L,JV)
     TX(I,JV)=sx(I,K,L,JV)
     TY(I,JV)=sy(I,K,L,JV)
     TZ(I,JV)=sz(I,K,L,JV)
  END DO
  END DO

  ENDIF

  DO I=1,LONK
     UEXT(I)=UGRI(I*NUMK,K,L)
  END DO

  !  place limits on appropriate moments before transport
  !  (if flux-limiting is to be applied)

  IF(.NOT.LIMIT) GO TO 13

  DO JV=1,NTRA
  DO I=1,LONK
    TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
  END DO
  END DO

 13   CONTINUE

  !  calculate flux and moments between adjacent boxes
  !  1- create temporary moments/masses for partial boxes in transit
  !  2- reajusts moments remaining in the box

  !  flux from IP to I if U(I).lt.0

  DO I=1,LONK-1
     IF(UEXT(I)<0.) THEN
       FM(I)=-UEXT(I)*DTX
       ALF(I)=FM(I)/TM(I+1)
       TM(I+1)=TM(I+1)-FM(I)
     ENDIF
  END DO

  I=LONK
  IF(UEXT(I)<0.) THEN
    FM(I)=-UEXT(I)*DTX
    ALF(I)=FM(I)/TM(1)
    TM(1)=TM(1)-FM(I)
  ENDIF

  !  flux from I to IP if U(I).gt.0

  DO I=1,LONK
     IF(UEXT(I)>=0.) THEN
       FM(I)=UEXT(I)*DTX
       ALF(I)=FM(I)/TM(I)
       TM(I)=TM(I)-FM(I)
     ENDIF
  END DO

  DO I=1,LONK
     ALFQ(I)=ALF(I)*ALF(I)
     ALF1(I)=1.-ALF(I)
     ALF1Q(I)=ALF1(I)*ALF1(I)
  END DO

  DO JV=1,NTRA
  DO I=1,LONK-1

     IF(UEXT(I)<0.) THEN

       F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
       FX(I,JV)=ALFQ(I)*TX(I+1,JV)
       FY(I,JV)=ALF (I)*TY(I+1,JV)
       FZ(I,JV)=ALF (I)*TZ(I+1,JV)

       T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
       TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
       TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
       TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)

     ENDIF

  END DO
  END DO

  I=LONK
  IF(UEXT(I)<0.) THEN

    DO JV=1,NTRA

       F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
       FX (I,JV)=ALFQ(I)*TX(1,JV)
       FY (I,JV)=ALF (I)*TY(1,JV)
       FZ (I,JV)=ALF (I)*TZ(1,JV)

       T0(1,JV)=T0(1,JV)-F0(I,JV)
       TX(1,JV)=ALF1Q(I)*TX(1,JV)
       TY(1,JV)=TY(1,JV)-FY(I,JV)
       TZ(1,JV)=TZ(1,JV)-FZ(I,JV)

  END DO

  ENDIF

  DO JV=1,NTRA
  DO I=1,LONK

     IF(UEXT(I)>=0.) THEN

       F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
       FX(I,JV)=ALFQ(I)*TX(I,JV)
       FY(I,JV)=ALF (I)*TY(I,JV)
       FZ(I,JV)=ALF (I)*TZ(I,JV)

       T0(I,JV)=T0(I,JV)-F0(I,JV)
       TX(I,JV)=ALF1Q(I)*TX(I,JV)
       TY(I,JV)=TY(I,JV)-FY(I,JV)
       TZ(I,JV)=TZ(I,JV)-FZ(I,JV)

     ENDIF

  END DO
  END DO

  !  puts the temporary moments Fi into appropriate neighboring boxes

  DO I=1,LONK
     IF(UEXT(I)<0.) THEN
       TM(I)=TM(I)+FM(I)
       ALF(I)=FM(I)/TM(I)
     ENDIF
  END DO

  DO I=1,LONK-1
     IF(UEXT(I)>=0.) THEN
       TM(I+1)=TM(I+1)+FM(I)
       ALF(I)=FM(I)/TM(I+1)
     ENDIF
  END DO

  I=LONK
  IF(UEXT(I)>=0.) THEN
    TM(1)=TM(1)+FM(I)
    ALF(I)=FM(I)/TM(1)
  ENDIF

  DO I=1,LONK
     ALF1(I)=1.-ALF(I)
  END DO

  DO JV=1,NTRA
  DO I=1,LONK

     IF(UEXT(I)<0.) THEN

       TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
       T0(I,JV)=T0(I,JV)+F0(I,JV)
       TX(I,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM
       TY(I,JV)=TY(I,JV)+FY(I,JV)
       TZ(I,JV)=TZ(I,JV)+FZ(I,JV)

     ENDIF

  END DO
  END DO

  DO JV=1,NTRA
  DO I=1,LONK-1

     IF(UEXT(I)>=0.) THEN

       TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
       T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
       TX(I+1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I+1,JV)+3.*TEMPTM
       TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
       TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)

     ENDIF

  END DO
  END DO

  I=LONK
  IF(UEXT(I)>=0.) THEN
    DO JV=1,NTRA
       TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
       T0(1,JV)=T0(1,JV)+F0(I,JV)
       TX(1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM
       TY(1,JV)=TY(1,JV)+FY(I,JV)
       TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
  END DO
  ENDIF

  !  retour aux mailles d'origine (passage des Tij aux Sij)

  IF(NUMK>1) THEN

  DO I2=1,NUMK

     DO I=1,LONK

        I3=I2+(I-1)*NUMK
        SM(I3,K,L)=SMNEW(I3)
        ALF(I)=SMNEW(I3)/TM(I)
        TM(I)=TM(I)-SMNEW(I3)

        ALFQ(I)=ALF(I)*ALF(I)
        ALF1(I)=1.-ALF(I)
        ALF1Q(I)=ALF1(I)*ALF1(I)

  END DO
  END DO

     DO  JV=1,NTRA
     DO  I=1,LONK

        I3=I2+(I-1)*NUMK
        S0(I3,K,L,JV)=ALF (I) &
              * (T0(I,JV)-ALF1(I)*TX(I,JV))
        sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
        sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
        sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)

  !   reajusts moments remaining in the box

        T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
        TX(I,JV)=ALF1Q(I)*TX(I,JV)
        TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
        TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
      ENDDO
      ENDDO


  ELSE

  DO I=1,LON
     SM(I,K,L)=TM(I)
  END DO
  DO JV=1,NTRA
  DO I=1,LON
     S0(I,K,L,JV)=T0(I,JV)
     sx(I,K,L,JV)=TX(I,JV)
     sy(I,K,L,JV)=TY(I,JV)
     sz(I,K,L,JV)=TZ(I,JV)
  END DO
  END DO

  ENDIF

  END DO
  END DO

  ! ----------- AA Test en fin de ADVX ------ Controle des S*
  ! OK
  !  DO 9998 l = 1, llm
  !  DO 9998 j = 1, jjp1
  !  DO 9998 i = 1, iip1
  !     IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
  !        PRINT*, '-------------------'
  !        PRINT*, 'En fin de ADVX'
  !        PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
  !        PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
  !        PRINT*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
  !        PRINT*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
  !        PRINT*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
  !        WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
  !c            STOP
  !     ENDIF
  ! 9998 CONTINUE

  ! ---------- bouclage cyclique
  DO itrac=1,ntra
  DO l = 1,llm
    DO j = lati,latf
       SM(iip1,j,l) = SM(1,j,l)
       S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
       sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
       sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
       sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
    END DO
  END DO
  ENDDO

  ! ----------- qqtite totale de traceur dans tte l'atmosphere
  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iim
  !IM 240405          sqf = sqf + S0(i,j,l,9)
         sqf = sqf + S0(i,j,l,ntra)
      END DO
    END DO
  END DO

  PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
  PRINT*,'sqf=',sqf
  !-------------

  RETURN
END SUBROUTINE advx
!_________________________________________________________________
!_________________________________________________________________