source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advy.f90 @ 5119

! $Header$

SUBROUTINE advy(limit,dty,pbarv,sm,s0,sx,sy,sz)
  IMPLICIT NONE

  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !                                                                C
  !  first-order moments (SOM) advection of tracer in Y direction  C
  !                                                                C
  !  Source : Pascal Simon ( Meteo, CNRM )                         C
  !  Adaptation : A.A. (LGGE)                                      C
  !  Derniere Modif : 15/12/94 LAST
                                                             ! C
  !  sont les arguments d'entree pour le s-pg                      C
  !                                                                C
  !  argument de sortie du s-pg                                    C
  !                                                                C
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !
  !  Rem : Probleme aux poles il faut reecrire ce cas specifique
  !    Attention au sens de l'indexation
  !
  !  parametres principaux du modele
  !
  !
  include "dimensions.h"
  include "paramet.h"
  include "comgeom2.h"

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

  INTEGER :: lon,lat,niv
  INTEGER :: i,j,jv,k,kp,l
  INTEGER :: ntra
  PARAMETER (ntra = 1)

  REAL :: dty
  REAL :: pbarv ( iip1,jjm, 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 :: VGRI(iip1,0:jjp1,llm)

  !  Rem : UGRI et WGRI ne sont pas utilises dans
  !  cette SUBROUTINE ( advection en y uniquement )
  !  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
  !
  !  the moments F are similarly defined and used as temporary
  !  storage for portions of the grid boxes in transit
  !
  REAL :: F0(iim,0:jjp1,ntra),FM(iim,0:jjp1)
  REAL :: FX(iim,jjm,ntra),FY(iim,jjm,ntra)
  REAL :: FZ(iim,jjm,ntra)
  REAL :: S00(ntra)
  REAL :: SM0             ! Just temporal variable
  !
  !  work arrays
  !
  REAL :: ALF(iim,0:jjp1),ALF1(iim,0:jjp1)
  REAL :: ALFQ(iim,0:jjp1),ALF1Q(iim,0:jjp1)
  REAL :: TEMPTM          ! Just temporal variable
  !
  !  Special pour poles
  !
  REAL :: sbms,sfms,sfzs,sbmn,sfmn,sfzn
  REAL :: sns0(ntra),snsz(ntra),snsm
  REAL :: s1v(llm),slatv(llm)
  REAL :: qy1(iim,llm,ntra),qylat(iim,llm,ntra)
  REAL :: cx1(llm,ntra), cxLAT(llm,ntra)
  REAL :: cy1(llm,ntra), cyLAT(llm,ntra)
  REAL :: z1(iim), zcos(iim), zsin(iim)
  REAL :: smpn,smps,s0pn,s0ps
  REAL :: SSUM
  EXTERNAL SSUM
  !
  REAL :: sqi,sqf
  LOGICAL :: LIMIT

  lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
  lat = jjp1        ! a cause des dim. differentes entre les
  niv=llm

  !
  !  the moments Fi are used as temporary storage for
  !  portions of the grid boxes in transit at the current level
  !
  !  work arrays
  !

  DO l = 1,llm
     DO j = 1,jjm
        DO i = 1,iip1
        vgri (i,j,llm+1-l)=-1.*pbarv(i,j,l)
        enddo
     enddo
     do i=1,iip1
         vgri(i,0,l) = 0.
         vgri(i,jjp1,l) = 0.
     enddo
  enddo

  DO L=1,NIV
  !
  !  place limits on appropriate moments before transport
  !  (if flux-limiting is to be applied)
  !
  IF(.NOT.LIMIT) GO TO 11
  !
  DO JV=1,NTRA
  DO K=1,LAT
  DO I=1,LON
     sy(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.), &
           ABS(sy(I,K,L,JV))),sy(I,K,L,JV))
  END DO
  END DO
  END DO
  !
 11   CONTINUE
  !
  !  le flux a travers le pole Nord est traite separement
  !
  SM0=0.
  DO JV=1,NTRA
     S00(JV)=0.
  END DO
  !
  DO I=1,LON
  !
     IF(VGRI(I,0,L)<=0.) THEN
       FM(I,0)=-VGRI(I,0,L)*DTY
       ALF(I,0)=FM(I,0)/SM(I,1,L)
       SM(I,1,L)=SM(I,1,L)-FM(I,0)
       SM0=SM0+FM(I,0)
     ENDIF
  !
     ALFQ(I,0)=ALF(I,0)*ALF(I,0)
     ALF1(I,0)=1.-ALF(I,0)
     ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
  !
  END DO
  !
  DO JV=1,NTRA
  DO I=1,LON
  !
     IF(VGRI(I,0,L)<=0.) THEN
  !
       F0(I,0,JV)=ALF(I,0)* &
             ( S0(I,1,L,JV)-ALF1(I,0)*sy(I,1,L,JV) )
  !
       S00(JV)=S00(JV)+F0(I,0,JV)
       S0(I,1,L,JV)=S0(I,1,L,JV)-F0(I,0,JV)
       sy(I,1,L,JV)=ALF1Q(I,0)*sy(I,1,L,JV)
       sx(I,1,L,JV)=ALF1 (I,0)*sx(I,1,L,JV)
       sz(I,1,L,JV)=ALF1 (I,0)*sz(I,1,L,JV)
  !
     ENDIF
  !
  END DO
  END DO
  !
  DO I=1,LON
     IF(VGRI(I,0,L)>0.) THEN
       FM(I,0)=VGRI(I,0,L)*DTY
       ALF(I,0)=FM(I,0)/SM0
     ENDIF
  END DO
  !
  DO JV=1,NTRA
  DO I=1,LON
     IF(VGRI(I,0,L)>0.) THEN
       F0(I,0,JV)=ALF(I,0)*S00(JV)
     ENDIF
  END DO
  END DO
  !
  !  puts the temporary moments Fi into appropriate neighboring boxes
  !
  DO I=1,LON
  !
     IF(VGRI(I,0,L)>0.) THEN
       SM(I,1,L)=SM(I,1,L)+FM(I,0)
       ALF(I,0)=FM(I,0)/SM(I,1,L)
     ENDIF
  !
     ALF1(I,0)=1.-ALF(I,0)
  !
  END DO
  !
  DO JV=1,NTRA
  DO I=1,LON
  !
     IF(VGRI(I,0,L)>0.) THEN
  !
     TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
     S0(I,1,L,JV)=S0(I,1,L,JV)+F0(I,0,JV)
     sy(I,1,L,JV)=ALF1(I,0)*sy(I,1,L,JV)+3.*TEMPTM
  !
     ENDIF
  !
  END DO
  END DO
  !
  !  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 KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
  !
  DO K=1,LAT-1
  KP=K+1
  DO I=1,LON
  !
     IF(VGRI(I,K,L)<0.) THEN
       FM(I,K)=-VGRI(I,K,L)*DTY
       ALF(I,K)=FM(I,K)/SM(I,KP,L)
       SM(I,KP,L)=SM(I,KP,L)-FM(I,K)
     ELSE
       FM(I,K)=VGRI(I,K,L)*DTY
       ALF(I,K)=FM(I,K)/SM(I,K,L)
       SM(I,K,L)=SM(I,K,L)-FM(I,K)
     ENDIF
  !
     ALFQ(I,K)=ALF(I,K)*ALF(I,K)
     ALF1(I,K)=1.-ALF(I,K)
     ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
  !
  END DO
  END DO
  !
  DO JV=1,NTRA
  DO K=1,LAT-1
  KP=K+1
  DO I=1,LON
  !
     IF(VGRI(I,K,L)<0.) THEN
  !
       F0(I,K,JV)=ALF (I,K)* &
             ( S0(I,KP,L,JV)-ALF1(I,K)*sy(I,KP,L,JV) )
       FY(I,K,JV)=ALFQ(I,K)*sy(I,KP,L,JV)
       FX(I,K,JV)=ALF (I,K)*sx(I,KP,L,JV)
       FZ(I,K,JV)=ALF (I,K)*sz(I,KP,L,JV)
  !
       S0(I,KP,L,JV)=S0(I,KP,L,JV)-F0(I,K,JV)
       sy(I,KP,L,JV)=ALF1Q(I,K)*sy(I,KP,L,JV)
       sx(I,KP,L,JV)=sx(I,KP,L,JV)-FX(I,K,JV)
       sz(I,KP,L,JV)=sz(I,KP,L,JV)-FZ(I,K,JV)
  !
     ELSE
  !
       F0(I,K,JV)=ALF (I,K)* &
             ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
       FY(I,K,JV)=ALFQ(I,K)*sy(I,K,L,JV)
       FX(I,K,JV)=ALF(I,K)*sx(I,K,L,JV)
       FZ(I,K,JV)=ALF(I,K)*sz(I,K,L,JV)
  !
       S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,K,JV)
       sy(I,K,L,JV)=ALF1Q(I,K)*sy(I,K,L,JV)
       sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,K,JV)
       sz(I,K,L,JV)=sz(I,K,L,JV)-FZ(I,K,JV)
  !
     ENDIF
  !
  END DO
  END DO
  END DO
  !
  !  puts the temporary moments Fi into appropriate neighboring boxes
  !
  DO K=1,LAT-1
  KP=K+1
  DO I=1,LON
  !
     IF(VGRI(I,K,L)<0.) THEN
       SM(I,K,L)=SM(I,K,L)+FM(I,K)
       ALF(I,K)=FM(I,K)/SM(I,K,L)
     ELSE
       SM(I,KP,L)=SM(I,KP,L)+FM(I,K)
       ALF(I,K)=FM(I,K)/SM(I,KP,L)
     ENDIF
  !
     ALF1(I,K)=1.-ALF(I,K)
  !
  END DO
  END DO
  !
  DO JV=1,NTRA
  DO K=1,LAT-1
  KP=K+1
  DO I=1,LON
  !
     IF(VGRI(I,K,L)<0.) THEN
  !
     TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
     S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
     sy(I,K,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*sy(I,K,L,JV) &
           +3.*TEMPTM
     sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,K,JV)
     sz(I,K,L,JV)=sz(I,K,L,JV)+FZ(I,K,JV)
  !
     ELSE
  !
     TEMPTM=ALF(I,K)*S0(I,KP,L,JV)-ALF1(I,K)*F0(I,K,JV)
     S0(I,KP,L,JV)=S0(I,KP,L,JV)+F0(I,K,JV)
     sy(I,KP,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*sy(I,KP,L,JV) &
           +3.*TEMPTM
     sx(I,KP,L,JV)=sx(I,KP,L,JV)+FX(I,K,JV)
     sz(I,KP,L,JV)=sz(I,KP,L,JV)+FZ(I,K,JV)
  !
     ENDIF
  !
  END DO
  END DO
  END DO
  !
  !  traitement special pour le pole Sud (idem pole Nord)
  !
  K=LAT
  !
  SM0=0.
  DO JV=1,NTRA
     S00(JV)=0.
  END DO
  !
  DO I=1,LON
  !
     IF(VGRI(I,K,L)>=0.) THEN
       FM(I,K)=VGRI(I,K,L)*DTY
       ALF(I,K)=FM(I,K)/SM(I,K,L)
       SM(I,K,L)=SM(I,K,L)-FM(I,K)
       SM0=SM0+FM(I,K)
     ENDIF
  !
     ALFQ(I,K)=ALF(I,K)*ALF(I,K)
     ALF1(I,K)=1.-ALF(I,K)
     ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
  !
  END DO
  !
  DO JV=1,NTRA
  DO I=1,LON
  !
     IF(VGRI(I,K,L)>=0.) THEN
       F0 (I,K,JV)=ALF(I,K)* &
             ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
       S00(JV)=S00(JV)+F0(I,K,JV)
  !
       S0(I,K,L,JV)=S0 (I,K,L,JV)-F0 (I,K,JV)
       sy(I,K,L,JV)=ALF1Q(I,K)*sy(I,K,L,JV)
       sx(I,K,L,JV)=ALF1(I,K)*sx(I,K,L,JV)
       sz(I,K,L,JV)=ALF1(I,K)*sz(I,K,L,JV)
     ENDIF
  !
  END DO
  END DO
  !
  DO I=1,LON
     IF(VGRI(I,K,L)<0.) THEN
       FM(I,K)=-VGRI(I,K,L)*DTY
       ALF(I,K)=FM(I,K)/SM0
     ENDIF
  END DO
  !
  DO JV=1,NTRA
  DO I=1,LON
     IF(VGRI(I,K,L)<0.) THEN
       F0(I,K,JV)=ALF(I,K)*S00(JV)
     ENDIF
  END DO
  END DO
  !
  !  puts the temporary moments Fi into appropriate neighboring boxes
  !
  DO I=1,LON
  !
     IF(VGRI(I,K,L)<0.) THEN
       SM(I,K,L)=SM(I,K,L)+FM(I,K)
       ALF(I,K)=FM(I,K)/SM(I,K,L)
     ENDIF
  !
     ALF1(I,K)=1.-ALF(I,K)
  !
  END DO
  !
  DO JV=1,NTRA
  DO I=1,LON
  !
     IF(VGRI(I,K,L)<0.) THEN
  !
     TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
     S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
     sy(I,K,L,JV)=ALF1(I,K)*sy(I,K,L,JV)+3.*TEMPTM
  !
     ENDIF
  !
  END DO
  END DO
  !
  END DO
  !
  RETURN
END SUBROUTINE advy