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

! $Header$

SUBROUTINE advy(limit, dty, pbarv, sm, s0, sx, sy, sz)
  USE lmdz_comgeom2

  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
  USE lmdz_paramet
  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

  !



  !  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 :: 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