source: LMDZ6/branches/Amaury_dev/libf/dyn3d/tetaleveli1j.F90 @ 5116

!================================================================
!================================================================
SUBROUTINE tetaleveli1j(ilon, ilev, lnew, pgcm, pres, Qgcm, Qpres)
  !================================================================
  !================================================================

  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
  IMPLICIT none

  include "dimensions.h"

  !================================================================
  !
  ! 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
  logical :: lnew

  REAL :: pgcm(ilon, ilev)
  REAL :: Qgcm(ilon, ilev)
  REAL :: pres
  REAL :: Qpres(ilon)

  !   local :
  !   -------

  !IM 211004
  ! INTEGER lt(klon), lb(klon)
  ! REAL ptop, pbot, aist(klon), aisb(klon)
  !
  include "paramet.h"
  !
  INTEGER :: lt(ip1jm), lb(ip1jm)
  REAL :: ptop, pbot, aist(ip1jm), aisb(ip1jm)
  !MI 211004
  save lt, lb, ptop, pbot, aist, aisb

  INTEGER :: i, k
  !
  ! PRINT*,'tetalevel pres=',pres
  !=====================================================================
  if (lnew) THEN
    !   on réinitialise les réindicages 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 i = 1, ilon
      !IM      IF ( ABS(pres-pgcm(i,ilev) ) .LT.
      IF (ABS(pres - pgcm(i, ilev)) > &
              ABS(pres - pgcm(i, 1))) THEN
        lt(i) = ilev     ! 2
        lb(i) = ilev - 1   ! 1
      ELSE
        lt(i) = 2
        lb(i) = 1
      ENDIF
      !IM   PRINT*,'i, ABS(pres-pgcm),ABS(pres-pgcm)',
      !IM  .i, ABS(pres-pgcm(i,ilev)),ABS(pres-pgcm(i,1))
    END DO
    DO k = 1, ilev - 1
      DO i = 1, ilon
        pbot = pgcm(i, k)
        ptop = pgcm(i, k + 1)
        !IM         IF (ptop.LE.pres .AND. pbot.GE.pres) THEN
        IF (ptop>=pres .AND. pbot<=pres) THEN
          lt(i) = k + 1
          lb(i) = k
        ENDIF
      END DO
    END DO
    !
    ! Interpolation lineaire:
    !
    DO i = 1, ilon
      ! interpolation en logarithme de pression:
      !
      ! ...   Modif . P. Le Van    ( 20/01/98) ....
      !   Modif Frédéric Hourdin (3/01/02)

      IF(pgcm(i, lb(i))==0.OR. &
              pgcm(i, lt(i))==0.) THEN
        !
        PRINT*, 'i,lb,lt,2pgcm,pres', i, lb(i), &
                lt(i), pgcm(i, lb(i)), pgcm(i, lt(i)), pres
        !
      ENDIF
      !
      aist(i) = LOG(pgcm(i, lb(i)) / pres) &
              / LOG(pgcm(i, lb(i)) / pgcm(i, lt(i)))
      aisb(i) = LOG(pres / pgcm(i, lt(i))) &
              / LOG(pgcm(i, lb(i)) / pgcm(i, lt(i)))
    enddo

  endif ! lnew

  !======================================================================
  !    inteprollation
  !======================================================================

  do i = 1, ilon
    Qpres(i) = Qgcm(i, lb(i)) * aisb(i) + Qgcm(i, lt(i)) * aist(i)
    !IM      PRINT*,'i,Qgcm,Qpres',i,Qgcm(i,lb(i)),aisb(i),
    !IM  $   Qgcm(i,lt(i)),aist(i),Qpres(i)
  enddo
  !
  ! Je mets les vents a zero quand je rencontre une montagne
  do i = 1, ilon
    !IM      if (pgcm(i,1).LT.pres) THEN
    if (pgcm(i, 1)>pres) THEN
      ! Qpres(i)=1e33
      Qpres(i) = 1e+20
      !IM         PRINT*,'i,pgcm(i,1),pres =',i,pgcm(i,1),pres
    endif
  enddo

  !

END SUBROUTINE tetaleveli1j