MODULE lmdz_plevel_new
  IMPLICIT NONE; PRIVATE
  PUBLIC plevel_new
CONTAINS

  SUBROUTINE plevel_new(ilon, ilev, klevstd, lnew, pgcm, pres, qgcm, qpres)
    USE dimphy
    USE phys_state_var_mod, ONLY: missing_val_nf90
    USE lmdz_wxios, ONLY: missing_val_xios => missing_val, using_xios

    IMPLICIT NONE

    ! ================================================================

    ! Interpoler des champs 3-D u, v et g du modele a un niveau de
    ! pression donnee (pres)

    ! INPUT:  ilon ----- nombre de points
    ! ilev ----- nombre de couches
    ! lnew ----- true si on doit reinitialiser les poids
    ! pgcm ----- pressions modeles
    ! pres ----- pression vers laquelle on interpolle
    ! Qgcm ----- champ GCM
    ! Qpres ---- champ interpolle au niveau pres

    ! ================================================================

    ! arguments :
    ! -----------

    INTEGER ilon, ilev, klevstd
    LOGICAL lnew

    REAL pgcm(ilon, ilev)
    REAL qgcm(ilon, ilev)
    REAL pres(klevstd)
    REAL qpres(ilon, klevstd)

    ! local :
    ! -------

    ! ym      INTEGER lt(klon), lb(klon)
    ! ym      REAL ptop, pbot, aist(klon), aisb(klon)

    ! ym      save lt,lb,ptop,pbot,aist,aisb
    INTEGER, ALLOCATABLE, SAVE, DIMENSION (:, :) :: lt, lb
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :) :: aist, aisb
    !$OMP THREADPRIVATE(lt,lb,aist,aisb)
    REAL, SAVE :: ptop, pbot
    !$OMP THREADPRIVATE(ptop, pbot)
    LOGICAL, SAVE :: first = .TRUE.
    INTEGER :: nlev
    !$OMP THREADPRIVATE(first)
    INTEGER i, k

    REAL :: missing_val

    IF (using_xios) THEN
      missing_val = missing_val_xios
    ELSE
      missing_val = missing_val_nf90
    ENDIF

    IF (first) THEN
      ALLOCATE (lt(klon, klevstd), lb(klon, klevstd))
      ALLOCATE (aist(klon, klevstd), aisb(klon, klevstd))
      first = .FALSE.
    END IF

    ! =====================================================================
    IF (lnew) THEN
      ! on reinitialise les reindicages et les poids
      ! =====================================================================


      ! Chercher les 2 couches les plus proches du niveau a obtenir

      ! Eventuellement, faire l'extrapolation a partir des deux couches
      ! les plus basses ou les deux couches les plus hautes:

      DO nlev = 1, klevstd
        DO i = 1, klon
          IF (abs(pres(nlev) - pgcm(i, ilev))<abs(pres(nlev) - pgcm(i, 1))) THEN
            lt(i, nlev) = ilev ! 2
            lb(i, nlev) = ilev - 1 ! 1
          ELSE
            lt(i, nlev) = 2
            lb(i, nlev) = 1
          END IF
        END DO
        DO k = 1, ilev - 1
          DO i = 1, klon
            pbot = pgcm(i, k)
            ptop = pgcm(i, k + 1)
            IF (ptop<=pres(nlev) .AND. pbot>=pres(nlev)) THEN
              lt(i, nlev) = k + 1
              lb(i, nlev) = k
            END IF
          END DO
        END DO

        ! Interpolation lineaire:
        DO i = 1, klon
          ! interpolation en logarithme de pression:

          ! ...   Modif . P. Le Van    ( 20/01/98) ....
          ! Modif Frederic Hourdin (3/01/02)

          aist(i, nlev) = log(pgcm(i, lb(i, nlev)) / pres(nlev)) / log(pgcm(i, lb(i, &
                  nlev)) / pgcm(i, lt(i, nlev)))
          aisb(i, nlev) = log(pres(nlev) / pgcm(i, lt(i, nlev))) / log(pgcm(i, lb(i, &
                  nlev)) / pgcm(i, lt(i, nlev)))
        END DO
      END DO

    END IF ! lnew

    ! ======================================================================
    ! inteprollation
    ! ET je mets les vents a zero quand je rencontre une montagne
    ! ======================================================================

    DO nlev = 1, klevstd
      DO i = 1, klon
        IF (pgcm(i, 1)<pres(nlev)) THEN
          qpres(i, nlev) = missing_val
        ELSE
          qpres(i, nlev) = qgcm(i, lb(i, nlev)) * aisb(i, nlev) + &
                  qgcm(i, lt(i, nlev)) * aist(i, nlev)
        END IF
      END DO
    END DO

  END SUBROUTINE plevel_new
END MODULE lmdz_plevel_new