source: LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_advect.f90 @ 5227

MODULE lmdz_advect
  IMPLICIT NONE; PRIVATE
  PUBLIC advect

CONTAINS

  SUBROUTINE advect(ucov, vcov, teta, w, massebx, masseby, du, dv, dteta)

    USE comconst_mod, ONLY: daysec
    USE logic_mod, ONLY: conser
    USE ener_mod, ONLY: gtot
    USE lmdz_ssum_scopy, ONLY: ssum
    USE lmdz_comgeom

    USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
    USE lmdz_paramet
    IMPLICIT NONE
    !=======================================================================

    !   Auteurs:  P. Le Van , Fr. Hourdin  .
    !   -------

    !   Objet:
    !   ------

    !   *************************************************************
    !   .... calcul des termes d'advection vertic.pour u,v,teta,q ...
    !   *************************************************************
    !    ces termes sont ajoutes a du,dv,dteta et dq .
    !  Modif F.Forget 03/94 : on retire q de advect

    !=======================================================================
    !-----------------------------------------------------------------------
    !   Declarations:
    !   -------------




    !   Arguments:
    !   ----------

    REAL :: vcov(ip1jm, llm), ucov(ip1jmp1, llm), teta(ip1jmp1, llm)
    REAL :: massebx(ip1jmp1, llm), masseby(ip1jm, llm), w(ip1jmp1, llm)
    REAL :: dv(ip1jm, llm), du(ip1jmp1, llm), dteta(ip1jmp1, llm)

    !   Local:
    !   ------

    REAL :: uav(ip1jmp1, llm), vav(ip1jm, llm), wsur2(ip1jmp1)
    REAL :: unsaire2(ip1jmp1), ge(ip1jmp1)
    REAL :: deuxjour, ww, gt, uu, vv

    INTEGER :: ij, l

    !-----------------------------------------------------------------------
    !   2. Calculs preliminaires:
    !   -------------------------

    IF (conser)  THEN
      deuxjour = 2. * daysec

      DO   ij = 1, ip1jmp1
        unsaire2(ij) = unsaire(ij) * unsaire(ij)
      END DO
    END IF


    !------------------  -yy ----------------------------------------------
    !   .  Calcul de     u

    DO  l = 1, llm
      DO    ij = iip2, ip1jmp1
        uav(ij, l) = 0.25 * (ucov(ij, l) + ucov(ij - iip1, l))
      ENDDO
      DO    ij = iip2, ip1jm
        uav(ij, l) = uav(ij, l) + uav(ij + iip1, l)
      ENDDO
      DO      ij = 1, iip1
        uav(ij, l) = 0.
        uav(ip1jm + ij, l) = 0.
      ENDDO
    ENDDO

    !------------------  -xx ----------------------------------------------
    !   .  Calcul de     v

    DO  l = 1, llm
      DO    ij = 2, ip1jm
        vav(ij, l) = 0.25 * (vcov(ij, l) + vcov(ij - 1, l))
      ENDDO
      DO    ij = 1, ip1jm, iip1
        vav(ij, l) = vav(ij + iim, l)
      ENDDO
      DO    ij = 1, ip1jm - 1
        vav(ij, l) = vav(ij, l) + vav(ij + 1, l)
      ENDDO
      DO    ij = 1, ip1jm, iip1
        vav(ij + iim, l) = vav(ij, l)
      ENDDO
    ENDDO

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

    DO l = 1, llmm1


      ! ......   calcul de  - w/2.    au niveau  l+1   .......

      DO ij = 1, ip1jmp1
        wsur2(ij) = - 0.5 * w(ij, l + 1)
      END DO


      ! .....................     calcul pour  du     ..................

      DO ij = iip2, ip1jm - 1
        ww = wsur2 (ij) + wsur2(ij + 1)
        uu = 0.5 * (ucov(ij, l) + ucov(ij, l + 1))
        du(ij, l) = du(ij, l) - ww * (uu - uav(ij, l)) / massebx(ij, l)
        du(ij, l + 1) = du(ij, l + 1) + ww * (uu - uav(ij, l + 1)) / massebx(ij, l + 1)
      END DO

      ! .....  correction pour  du(iip1,j,l)  ........
      ! .....     du(iip1,j,l)= du(1,j,l)   .....

      !DIR$ IVDEP
      DO   ij = iip1 + iip1, ip1jm, iip1
        du(ij, l) = du(ij - iim, l)
        du(ij, l + 1) = du(ij - iim, l + 1)
      END DO

      ! .................    calcul pour   dv      .....................

      DO ij = 1, ip1jm
        ww = wsur2(ij + iip1) + wsur2(ij)
        vv = 0.5 * (vcov(ij, l) + vcov(ij, l + 1))
        dv(ij, l) = dv(ij, l) - ww * (vv - vav(ij, l)) / masseby(ij, l)
        dv(ij, l + 1) = dv(ij, l + 1) + ww * (vv - vav(ij, l + 1)) / masseby(ij, l + 1)
      END DO

      !

      ! ............................................................
      ! ...............    calcul pour   dh      ...................
      ! ............................................................

      !                   ---z
      !   calcul de  - d( teta  * w )      qu'on ajoute a   dh
      !               ...............

      DO ij = 1, ip1jmp1
        ww = wsur2(ij) * (teta(ij, l) + teta(ij, l + 1))
        dteta(ij, l) = dteta(ij, l) - ww
        dteta(ij, l + 1) = dteta(ij, l + 1) + ww
      END DO

      IF(conser)  THEN
        DO ij = 1, ip1jmp1
          ge(ij) = wsur2(ij) * wsur2(ij) * unsaire2(ij)
        END DO
        gt = SSUM(ip1jmp1, ge, 1)
        gtot(l) = deuxjour * SQRT(gt / ip1jmp1)
      END IF

    END DO

  END SUBROUTINE advect

END MODULE lmdz_advect