source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advz.f90 @ 5186

! $Header$

SUBROUTINE advz(limit,dtz,w,sm,s0,sx,sy,sz)
  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
  USE lmdz_paramet
  IMPLICIT NONE

  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !                                                                C
  !  first-order moments (FOM) advection of tracer in Z direction  C
  !                                                                C
  !  Source : Pascal Simon (Meteo,CNRM)                            C
  !  Adaptation : A.Armengaud (LGGE) juin 94                       C
  !                                                                C
  !                                                                C
  !  sont des arguments d'entree pour le s-pg...                   C
  !                                                                C
  !  dq est l'argument de sortie pour le s-pg                      C
  !                                                                C
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

  !  parametres principaux du modele
  !



  ! INCLUDE "traceur.h"

  !  Arguments :
  !  -----------
  !  dtz : frequence fictive d'appel du transport
  !  w : flux de masse en z en Pa.m2.s-1

  INTEGER :: ntra
  PARAMETER (ntra = 1)

  REAL :: dtz
  REAL :: w ( 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) &
        ,sz(iip1,jjp1,llm,ntra)


  !  Local :
  !  -------

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

  REAL :: WGRI(iip1,jjp1,0:llm)


  !  the moments F are used as temporary  storage for
  !  portions of grid boxes in transit at the current latitude

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

  !  work arrays

  REAL :: ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
  REAL :: TEMPTM            ! Just temporal variable
  REAL :: sqi,sqf

  LOGICAL :: LIMIT
  INTEGER :: lon,lat,niv
  INTEGER :: i,j,jv,k,l,lp

  lon = iim
  lat = jjp1
  niv = llm

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

  ! DO 399 l = 1, llm
  ! DO 399 j = 1, jjp1
  ! DO 399 i = 1, iip1
  !    IF (S0(i,j,l,ntra) .lt. 0. ) THEN
  !       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)
  !       PRINT*, 'AIE !! debut ADVZ - pbl arg. passage dans ADVZ'
  !        STOP
  !    ENDIF
  399   CONTINUE

  !-----------------------------------------------------------------
  ! *** Test : diag de la qqtite totale de traceur
         ! dans l'atmosphere avant l'advection en z
  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 ADVZ - ENTREE ---------'
  PRINT*,'sqi=',sqi

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

  !  Conversion du flux de masse en kg.s-1

  DO l = 1,llm
     DO j = 1,jjp1
        DO i = 1,iip1
         ! wgri (i,j,llm+1-l) =  w (i,j,l) / g
           wgri (i,j,llm+1-l) =  w (i,j,l)
          ! wgri (i,j,0) = 0.                ! a detruire ult.
          ! wgri (i,j,l) = 0.1               !    w (i,j,l)
          ! wgri (i,j,llm) = 0.              ! a detruire ult.
  END DO
  END DO
  END DO
     DO  j = 1,jjp1
        DO i = 1,iip1
           wgri(i,j,0)=0.
        enddo
     enddo

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

  !  start here
  !  boucle sur les latitudes

  DO K=1,LAT

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

  IF(.NOT.LIMIT) GO TO 101

  DO JV=1,NTRA
  DO L=1,NIV
     DO I=1,LON
        sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.), &
              ABS(sz(I,K,L,JV))),sz(I,K,L,JV))
  END DO
  END DO
  END DO

 101   CONTINUE

  !  boucle sur les niveaux intercouches de 1 a NIV-1
  !   (flux nul au sommet L=0 et a la base L=NIV)

  !  calculate flux and moments between adjacent boxes
  ! (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
  !  1- create temporary moments/masses for partial boxes in transit
  !  2- reajusts moments remaining in the box

  DO L=1,NIV-1
  LP=L+1

  DO I=1,LON

     IF(WGRI(I,K,L)<0.) THEN
       FM(I,L)=-WGRI(I,K,L)*DTZ
       ALF(I)=FM(I,L)/SM(I,K,LP)
       SM(I,K,LP)=SM(I,K,LP)-FM(I,L)
     ELSE
       FM(I,L)=WGRI(I,K,L)*DTZ
       ALF(I)=FM(I,L)/SM(I,K,L)
       SM(I,K,L)=SM(I,K,L)-FM(I,L)
     ENDIF

     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,LON

     IF(WGRI(I,K,L)<0.) THEN

       F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) )
       FZ(I,L,JV)=ALFQ(I)*sz(I,K,LP,JV)
       FX(I,L,JV)=ALF (I)*sx(I,K,LP,JV)
       FY(I,L,JV)=ALF (I)*sy(I,K,LP,JV)

       S0(I,K,LP,JV)=S0(I,K,LP,JV)-F0(I,L,JV)
       sz(I,K,LP,JV)=ALF1Q(I)*sz(I,K,LP,JV)
       sx(I,K,LP,JV)=sx(I,K,LP,JV)-FX(I,L,JV)
       sy(I,K,LP,JV)=sy(I,K,LP,JV)-FY(I,L,JV)

     ELSE

       F0(I,L,JV)=ALF (I)*(S0(I,K,L,JV)+ALF1(I)*sz(I,K,L,JV) )
       FZ(I,L,JV)=ALFQ(I)*sz(I,K,L,JV)
       FX(I,L,JV)=ALF (I)*sx(I,K,L,JV)
       FY(I,L,JV)=ALF (I)*sy(I,K,L,JV)

       S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,L,JV)
       sz(I,K,L,JV)=ALF1Q(I)*sz(I,K,L,JV)
       sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,L,JV)
       sy(I,K,L,JV)=sy(I,K,L,JV)-FY(I,L,JV)

     ENDIF

  END DO
  END DO

  END DO

  !  puts the temporary moments Fi into appropriate neighboring boxes

  DO L=1,NIV-1
  LP=L+1

  DO I=1,LON

     IF(WGRI(I,K,L)<0.) THEN
       SM(I,K,L)=SM(I,K,L)+FM(I,L)
       ALF(I)=FM(I,L)/SM(I,K,L)
     ELSE
       SM(I,K,LP)=SM(I,K,LP)+FM(I,L)
       ALF(I)=FM(I,L)/SM(I,K,LP)
     ENDIF

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

  END DO

  DO JV=1,NTRA
  DO I=1,LON

     IF(WGRI(I,K,L)<0.) THEN

       TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
       S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)
       sz(I,K,L,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,L,JV)+3.*TEMPTM
       sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,L,JV)
       sy(I,K,L,JV)=sy(I,K,L,JV)+FY(I,L,JV)

     ELSE

       TEMPTM=ALF(I)*S0(I,K,LP,JV)-ALF1(I)*F0(I,L,JV)
       S0(I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)
       sz(I,K,LP,JV)=ALF(I)*FZ(I,L,JV)+ALF1(I)*sz(I,K,LP,JV) &
             +3.*TEMPTM
       sx(I,K,LP,JV)=sx(I,K,LP,JV)+FX(I,L,JV)
       sy(I,K,LP,JV)=sy(I,K,LP,JV)+FY(I,L,JV)

     ENDIF

  END DO
  END DO

  END DO

  !  fin de la boucle principale sur les latitudes

  END DO

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

  ! ----------- AA Test en fin de ADVX ------ Controle des S*

  ! DO 9999 l = 1, llm
  ! DO 9999 j = 1, jjp1
  ! DO 9999 i = 1, iip1
  !    IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
  !       PRINT*, '-------------------'
  !       PRINT*, 'En fin de ADVZ'
  !       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 ADVZ1'
  !        STOP
  !    ENDIF
 9999   CONTINUE

  ! *** ------------------- bouclage cyclique  en X ------------

   ! DO l = 1,llm
   !    DO j = 1,jjp1
   !       SM(iip1,j,l) = SM(1,j,l)
   !       S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
   !       sx(iip1,j,l,ntra) = sx(1,j,l,ntra)
   !       sy(iip1,j,l,ntra) = sy(1,j,l,ntra)
   !       sz(iip1,j,l,ntra) = sz(1,j,l,ntra)
   !    ENDDO
   ! ENDDO

  !-------------------------------------------------------------
  ! *** Test : diag de la qqtite totale de traceur
   !       dans l'atmosphere avant l'advection en z
  DO l = 1,llm
     DO j = 1,jjp1
        DO i = 1,iim
  !IM 240305            sqf = sqf + S0(i,j,l,9)
           sqf = sqf + S0(i,j,l,ntra)
        ENDDO
     ENDDO
  ENDDO
  PRINT*,'-------- DIAG DANS ADVZ - SORTIE ---------'
  PRINT*,'sqf=', sqf

  !-------------------------------------------------------------
  RETURN
END SUBROUTINE advz
!_______________________________________________________________
!_______________________________________________________________