module exner_hyb_loc_m

  IMPLICIT NONE

contains

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

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

    !                                 -------- z
    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
    !    ( voir note de Fr.Hourdin )  ,

    !    on determine successivement , du haut vers le bas des couches, les
    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,
    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .

    USE parallel_lmdz
    USE lmdz_filtreg_p
    USE write_field_loc
    USE comconst_mod, ONLY: cpp, kappa, r, jmp1
    USE comvert_mod, ONLY: preff
    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO

    IMPLICIT NONE

    include "dimensions.h"
    include "paramet.h"
    include "comgeom.h"

    INTEGER  ngrid
    REAL p(ijb_u:ije_u, llmp1), pk(ijb_u:ije_u, llm)
    REAL, optional :: pkf(ijb_u:ije_u, llm)
    REAL ps(ijb_u:ije_u), pks(ijb_u:ije_u)
    REAL alpha(ijb_u:ije_u, llm), beta(ijb_u:ije_u, llm)

    !    .... variables locales   ...

    INTEGER l, ij
    REAL unpl2k, dellta

    INTEGER ije, ijb, jje, jjb
    logical, save :: firstcall = .TRUE.
    !$OMP THREADPRIVATE(firstcall)
    CHARACTER(LEN = *), parameter :: modname = "exner_hyb_loc"

    !$OMP BARRIER

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

    !$OMP BARRIER

    ! Specific behaviour for Shallow Water (1 vertical layer) case:
    if (llm==1) THEN
      ! Compute pks(:),pk(:),pkf(:)
      ijb = ij_begin
      ije = ij_end
      !$OMP DO SCHEDULE(STATIC)
      DO ij = ijb, ije
        pks(ij) = (cpp / preff) * ps(ij)
        pk(ij, 1) = .5 * pks(ij)
        if (present(pkf)) pkf(ij, 1) = pk(ij, 1)
      ENDDO
      !$OMP ENDDO

      !$OMP BARRIER
      if (present(pkf)) THEN
        jjb = jj_begin
        jje = jj_end
        CALL filtreg_p (pkf, jjb_u, jje_u, jjb, jje, jmp1, llm, &
                2, 1, .TRUE., 1)
      end if

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

    ! General case:

    unpl2k = 1. + 2. * kappa

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

    ijb = ij_begin
    ije = ij_end

    !$OMP DO SCHEDULE(STATIC)
    DO   ij = ijb, ije
      pks(ij) = cpp * (ps(ij) / preff) ** kappa
    ENDDO
    !$OMP ENDDO
    ! Synchro OPENMP ici

    !$OMP BARRIER


    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..

    !$OMP DO SCHEDULE(STATIC)
    DO     ij = ijb, ije
      alpha(ij, llm) = 0.
      beta (ij, llm) = 1. / unpl2k
    ENDDO
    !$OMP ENDDO NOWAIT

    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...

    DO l = llm - 1, 2, -1

      !$OMP DO SCHEDULE(STATIC)
      DO ij = ijb, ije
        dellta = p(ij, l) * unpl2k + p(ij, l + 1) * (beta(ij, l + 1) - unpl2k)
        alpha(ij, l) = - p(ij, l + 1) / dellta * alpha(ij, l + 1)
        beta (ij, l) = p(ij, l) / dellta
      ENDDO
      !$OMP ENDDO NOWAIT
    ENDDO

    !  ***********************************************************************
    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....

    !$OMP DO SCHEDULE(STATIC)
    DO   ij = ijb, ije
      pk(ij, 1) = (p(ij, 1) * pks(ij) - 0.5 * alpha(ij, 2) * p(ij, 2)) / &
              (p(ij, 1) * (1. + kappa) + 0.5 * (beta(ij, 2) - unpl2k) * p(ij, 2))
    ENDDO
    !$OMP ENDDO NOWAIT

    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........

    DO l = 2, llm
      !$OMP DO SCHEDULE(STATIC)
      DO   ij = ijb, ije
        pk(ij, l) = alpha(ij, l) + beta(ij, l) * pk(ij, l - 1)
      ENDDO
      !$OMP ENDDO NOWAIT
    ENDDO

    if (present(pkf)) THEN
      !    calcul de pkf

      DO l = 1, llm
        !$OMP DO SCHEDULE(STATIC)
        DO   ij = ijb, ije
          pkf(ij, l) = pk(ij, l)
        ENDDO
        !$OMP ENDDO NOWAIT
      ENDDO

      !$OMP BARRIER

      jjb = jj_begin
      jje = jj_end
      IF (CPPKEY_DEBUGIO) THEN
        CALL WriteField_u('pkf', pkf)
      END IF
      CALL filtreg_p (pkf, jjb_u, jje_u, jjb, jje, jmp1, llm, &
              2, 1, .TRUE., 1)
      IF (CPPKEY_DEBUGIO) THEN
        CALL WriteField_u('pkf', pkf)
      END IF
    end if

  END SUBROUTINE exner_hyb_loc

end module exner_hyb_loc_m