source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/exner_milieu_m.F90 @ 5136

module exner_milieu_m

  IMPLICIT NONE

CONTAINS

  SUBROUTINE exner_milieu(ngrid, ps, p, pks, pk, pkf)

    !     Auteurs :  F. Forget , Y. Wanherdrick
    ! P.Le Van  , Fr. Hourdin  .
    !    ..........

    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...

    !   ************************************************************************
    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    !   ************************************************************************
    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    !    la pression et la fonction d'Exner  au  sol  .

    !     WARNING : CECI est une version speciale de exner_hyb originale
    !               Utilise dans la version martienne pour pouvoir 
    !               tourner avec des coordonnees verticales complexe
    !              => Il ne verifie PAS la condition la proportionalite en 
    !              energie totale/ interne / potentielle (F.Forget 2001)
    !    ( voir note de Fr.Hourdin )  ,

    USE comconst_mod, ONLY: jmp1, cpp, kappa, r
    USE comvert_mod, ONLY: preff
    USE lmdz_filtreg, ONLY: filtreg
    USE lmdz_comgeom

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "paramet.h"

    INTEGER  ngrid
    REAL p(ngrid, llmp1), pk(ngrid, llm)
    REAL, optional :: pkf(ngrid, llm)
    REAL ps(ngrid), pks(ngrid)

    !    .... variables locales   ...

    INTEGER l, ij
    REAL dum1

    logical, save :: firstcall = .TRUE.
    CHARACTER(LEN = *), parameter :: modname = "exner_milieu"

    ! Sanity check
    IF (firstcall) THEN
      ! sanity checks for Shallow Water case (1 vertical layer)
      IF (llm==1) THEN
        IF (kappa/=1) THEN
          CALL abort_gcm(modname, &
                  "kappa!=1 , but running in Shallow Water mode!!", 42)
        endif
        IF (cpp/=r) THEN
          CALL abort_gcm(modname, &
                  "cpp!=r , but running in Shallow Water mode!!", 42)
        endif
      endif ! of if (llm.EQ.1)

      firstcall = .FALSE.
    endif ! of if (firstcall)

    ! Specific behaviour for Shallow Water (1 vertical layer) case:
    IF (llm==1) THEN
      ! Compute pks(:),pk(:),pkf(:)

      DO   ij = 1, ngrid
        pks(ij) = (cpp / preff) * ps(ij)
        pk(ij, 1) = .5 * pks(ij)
      ENDDO

      IF (present(pkf)) THEN
        pkf = pk
        CALL filtreg (pkf, jmp1, llm, 2, 1, .TRUE., 1)
      end if

      ! our work is done, exit routine
      RETURN
    endif ! of if (llm.EQ.1)

    ! General case:

    !     -------------
    !     Calcul de pks
    !     -------------

    DO   ij = 1, ngrid
      pks(ij) = cpp * (ps(ij) / preff) ** kappa
    ENDDO

    !    .... Calcul de pk  pour la couche l 
    !    --------------------------------------------

    dum1 = cpp * (2 * preff)**(-kappa)
    DO l = 1, llm - 1
      DO   ij = 1, ngrid
        pk(ij, l) = dum1 * (p(ij, l) + p(ij, l + 1))**kappa
      ENDDO
    ENDDO

    !    .... Calcul de pk  pour la couche l = llm ..
    !    (on met la meme distance (en log pression)  entre Pk(llm)
    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)

    DO   ij = 1, ngrid
      pk(ij, llm) = pk(ij, llm - 1)**2 / pk(ij, llm - 2)
    ENDDO

    IF (present(pkf)) THEN
      !    calcul de pkf
      pkf = pk
      CALL filtreg (pkf, jmp1, llm, 2, 1, .TRUE., 1)
    end if

  END SUBROUTINE exner_milieu

END MODULE exner_milieu_m