! $Header$

SUBROUTINE conlmd(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, rain, snow, &
    ibas, itop)
  USE dimphy
  USE yomcst_mod_h
  USE yoethf_mod_h
IMPLICIT NONE
  ! ======================================================================
  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
  ! Objet: Schema de convection utilis'e dans le modele du LMD
  ! Ajustement humide (Manabe) + Ajustement convectif (Kuo)
  ! ======================================================================

  ! Arguments:

  REAL dtime ! pas d'integration (s)
  REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
  REAL t(klon, klev) ! temperature (K)
  REAL q(klon, klev) ! humidite specifique (kg/kg)
  REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)

  REAL d_t(klon, klev) ! incrementation temperature
  REAL d_q(klon, klev) ! incrementation humidite
  REAL rain(klon) ! pluies (mm/s)
  REAL snow(klon) ! neige (mm/s)
  INTEGER ibas(klon) ! niveau du bas
  INTEGER itop(klon) ! niveau du haut

  LOGICAL usekuo ! utiliser convection profonde (schema Kuo)
  PARAMETER (usekuo=.TRUE.)

  REAL d_t_bis(klon, klev)
  REAL d_q_bis(klon, klev)
  REAL rain_bis(klon)
  REAL snow_bis(klon)
  INTEGER ibas_bis(klon)
  INTEGER itop_bis(klon)
  REAL d_ql(klon, klev), d_ql_bis(klon, klev)
  REAL rneb(klon, klev), rneb_bis(klon, klev)

  INTEGER i, k
  REAL zlvdcp, zlsdcp, zdelta, zz, za, zb

  ! cc      CALL fiajh ! ancienne version de Convection Manabe
  CALL conman &                    ! nouvelle version de Convection
                                   ! Manabe
    (dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, ibas, itop)

  IF (usekuo) THEN
    ! cc      CALL fiajc ! ancienne version de Convection Kuo
    CALL conkuo &                  ! nouvelle version de Convection
                                   ! Kuo
      (dtime, paprs, pplay, t, q, conv_q, d_t_bis, d_q_bis, d_ql_bis, &
      rneb_bis, rain_bis, snow_bis, ibas_bis, itop_bis)
    DO k = 1, klev
      DO i = 1, klon
        d_t(i, k) = d_t(i, k) + d_t_bis(i, k)
        d_q(i, k) = d_q(i, k) + d_q_bis(i, k)
        d_ql(i, k) = d_ql(i, k) + d_ql_bis(i, k)
      END DO
    END DO
    DO i = 1, klon
      rain(i) = rain(i) + rain_bis(i)
      snow(i) = snow(i) + snow_bis(i)
      ibas(i) = min(ibas(i), ibas_bis(i))
      itop(i) = max(itop(i), itop_bis(i))
    END DO
  END IF

  ! L'eau liquide convective est dispersee dans l'air:

  DO k = 1, klev
    DO i = 1, klon
      zlvdcp = rlvtt/rcpd/(1.0+rvtmp2*q(i,k))
      zlsdcp = rlstt/rcpd/(1.0+rvtmp2*q(i,k))
      zdelta = max(0., sign(1.,rtt-t(i,k)))
      zz = d_ql(i, k) ! re-evap. de l'eau liquide
      zb = max(0.0, zz)
      za = -max(0.0, zz)*(zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
      d_t(i, k) = d_t(i, k) + za
      d_q(i, k) = d_q(i, k) + zb
    END DO
  END DO

  RETURN
END SUBROUTINE conlmd
SUBROUTINE conman(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, &
    snow, ibas, itop)
  USE dimphy
  USE yomcst_mod_h
  USE yoethf_mod_h
  IMPLICIT NONE
  ! ======================================================================
  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19970324
  ! Objet: ajustement humide convectif avec la possibilite de faire
  ! l'ajustement sur une fraction de la maille.
  ! Methode: On impose une distribution uniforme pour la vapeur d'eau
  ! au sein d'une maille. On applique la procedure d'ajustement
  ! successivement a la totalite, 75%, 50%, 25% et 5% de la maille
  ! jusqu'a ce que l'ajustement a lieu. J'espere que ceci augmente
  ! les activites convectives et corrige le biais "trop froid et sec"
  ! du modele.
  ! ======================================================================
  REAL dtime ! pas d'integration (s)
  REAL t(klon, klev) ! temperature (K)
  REAL q(klon, klev) ! humidite specifique (kg/kg)
  REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)

  REAL d_t(klon, klev) ! incrementation temperature
  REAL d_q(klon, klev) ! incrementation humidite
  REAL d_ql(klon, klev) ! incrementation eau liquide
  REAL rneb(klon, klev) ! nebulosite
  REAL rain(klon) ! pluies (mm/s)
  REAL snow(klon) ! neige (mm/s)
  INTEGER ibas(klon) ! niveau du bas
  INTEGER itop(klon) ! niveau du haut

  LOGICAL afaire(klon) ! .TRUE. implique l'ajustement
  LOGICAL accompli(klon) ! .TRUE. si l'ajustement est effectif

  INTEGER nb ! nombre de sous-fractions a considere
  PARAMETER (nb=1)
  ! cc      PARAMETER (nb=3)

  REAL ratqs ! largeur de la distribution pour vapeur d'eau
  PARAMETER (ratqs=0.05)

  REAL w_q(klon, klev)
  REAL w_d_t(klon, klev), w_d_q(klon, klev), w_d_ql(klon, klev)
  REAL w_rneb(klon, klev)
  REAL w_rain(klon), w_snow(klon)
  INTEGER w_ibas(klon), w_itop(klon)
  REAL zq1, zq2
  INTEGER i, k, n

  REAL t_coup
  PARAMETER (t_coup=234.0)
  REAL zdp1, zdp2
  REAL zqs1, zqs2, zdqs1, zdqs2
  REAL zgamdz
  REAL zflo ! flotabilite
  REAL zsat ! sur-saturation
  REAL zdelta, zcor, zcvm5
  LOGICAL imprim

  INTEGER ncpt
  SAVE ncpt
  !$OMP THREADPRIVATE(ncpt)
  REAL frac(nb) ! valeur de la maille fractionnelle
  SAVE frac
  !$OMP THREADPRIVATE(frac)
  INTEGER opt_cld(nb) ! option pour le modele nuageux
  SAVE opt_cld
  !$OMP THREADPRIVATE(opt_cld)
  LOGICAL appel1er
  SAVE appel1er
  !$OMP THREADPRIVATE(appel1er)

  ! Fonctions thermodynamiques:

  include "FCTTRE.h"

  DATA frac/1.0/
  DATA opt_cld/4/
  ! cc      DATA frac    / 1.0, 0.50, 0.25/
  ! cc      DATA opt_cld / 4,   4,    4/

  DATA appel1er/.TRUE./
  DATA ncpt/0/

  IF (appel1er) THEN
    PRINT *, 'conman, nb:', nb
    PRINT *, 'conman, frac:', frac
    PRINT *, 'conman, opt_cld:', opt_cld
    appel1er = .FALSE.
  END IF

  ! Initialiser les sorties a zero:

  DO k = 1, klev
    DO i = 1, klon
      d_t(i, k) = 0.0
      d_q(i, k) = 0.0
      d_ql(i, k) = 0.0
      rneb(i, k) = 0.0
    END DO
  END DO
  DO i = 1, klon
    ibas(i) = klev
    itop(i) = 1
    rain(i) = 0.0
    snow(i) = 0.0
  END DO

  ! S'il n'y a pas d'instabilite conditionnelle,
  ! pas la penne de se fatiguer:

  DO i = 1, klon
    afaire(i) = .FALSE.
  END DO
  DO k = 1, klev - 1
    DO i = 1, klon
      IF (thermcep) THEN
        zdelta = max(0., sign(1.,rtt-t(i,k)))
        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
        zqs1 = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
        zqs1 = min(0.5, zqs1)
        zcor = 1./(1.-retv*zqs1)
        zqs1 = zqs1*zcor
        zdqs1 = foede(t(i,k), zdelta, zcvm5, zqs1, zcor)

        zdelta = max(0., sign(1.,rtt-t(i,k+1)))
        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k+1))
        zqs2 = r2es*foeew(t(i,k+1), zdelta)/pplay(i, k+1)
        zqs2 = min(0.5, zqs2)
        zcor = 1./(1.-retv*zqs2)
        zqs2 = zqs2*zcor
        zdqs2 = foede(t(i,k+1), zdelta, zcvm5, zqs2, zcor)
      ELSE
        IF (t(i,k)<t_coup) THEN
          zqs1 = qsats(t(i,k))/pplay(i, k)
          zdqs1 = dqsats(t(i,k), zqs1)

          zqs2 = qsats(t(i,k+1))/pplay(i, k+1)
          zdqs2 = dqsats(t(i,k+1), zqs2)
        ELSE
          zqs1 = qsatl(t(i,k))/pplay(i, k)
          zdqs1 = dqsatl(t(i,k), zqs1)

          zqs2 = qsatl(t(i,k+1))/pplay(i, k+1)
          zdqs2 = dqsatl(t(i,k+1), zqs2)
        END IF
      END IF
      zdp1 = paprs(i, k) - paprs(i, k+1)
      zdp2 = paprs(i, k+1) - paprs(i, k+2)
      zgamdz = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd*(rd*(t(i, &
        k)*zdp1+t(i,k+1)*zdp2)/(zdp1+zdp2)+rlvtt*(zqs1*zdp1+zqs2*zdp2)/(zdp1+ &
        zdp2))/(1.0+(zdqs1*zdp1+zdqs2*zdp2)/(zdp1+zdp2))
      zflo = t(i, k) + zgamdz - t(i, k+1)
      zsat = (q(i,k)-zqs1)*zdp1 + (q(i,k+1)-zqs2)*zdp2
      IF (zflo>0.0) afaire(i) = .TRUE.
      ! erreur         IF (zflo.GT.0.0 .AND. zsat.GT.0.0) afaire(i) = .TRUE.
    END DO
  END DO

  imprim = mod(ncpt, 48) == 0
  DO n = 1, nb

    DO k = 1, klev
      DO i = 1, klon
        IF (afaire(i)) THEN
          zq1 = q(i, k)*(1.0-ratqs)
          zq2 = q(i, k)*(1.0+ratqs)
          w_q(i, k) = zq2 - frac(n)/2.0*(zq2-zq1)
        END IF
      END DO
    END DO

    CALL conmanv(dtime, paprs, pplay, t, w_q, afaire, opt_cld(n), w_d_t, &
      w_d_q, w_d_ql, w_rneb, w_rain, w_snow, w_ibas, w_itop, accompli, &
      imprim)
    DO k = 1, klev
      DO i = 1, klon
        IF (afaire(i) .AND. accompli(i)) THEN
          d_t(i, k) = w_d_t(i, k)*frac(n)
          d_q(i, k) = w_d_q(i, k)*frac(n)
          d_ql(i, k) = w_d_ql(i, k)*frac(n)
          IF (nint(w_rneb(i,k))==1) rneb(i, k) = frac(n)
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (afaire(i) .AND. accompli(i)) THEN
        rain(i) = w_rain(i)*frac(n)
        snow(i) = w_snow(i)*frac(n)
        ibas(i) = min(ibas(i), w_ibas(i))
        itop(i) = max(itop(i), w_itop(i))
      END IF
    END DO
    DO i = 1, klon
      IF (afaire(i) .AND. accompli(i)) afaire(i) = .FALSE.
    END DO

  END DO

  ncpt = ncpt + 1

  RETURN
END SUBROUTINE conman
SUBROUTINE conmanv(dtime, paprs, pplay, t, q, afaire, opt_cld, d_t, d_q, &
    d_ql, rneb, rain, snow, ibas, itop, accompli, imprim)
  USE dimphy
  USE yomcst_mod_h
  USE yoethf_mod_h
IMPLICIT NONE
  ! ======================================================================
  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
  ! Objet: ajustement humide (convection proposee par Manabe).
  ! Pour une colonne verticale, il peut avoir plusieurs blocs
  ! necessitant l'ajustement. ibas est le bas du plus bas bloc
  ! et itop est le haut du plus haut bloc
  ! ======================================================================


  ! Arguments:

  REAL dtime ! pas d'integration (s)
  REAL t(klon, klev) ! temperature (K)
  REAL q(klon, klev) ! humidite specifique (kg/kg)
  REAL paprs(klon, klev+1) ! pression inter-couche (Pa)
  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
  INTEGER opt_cld ! comment traiter l'eau liquide
  LOGICAL afaire(klon) ! .TRUE. si le point est a faire (Input)
  LOGICAL imprim ! .T. pour imprimer quelques diagnostiques

  REAL d_t(klon, klev) ! incrementation temperature
  REAL d_q(klon, klev) ! incrementation humidite
  REAL d_ql(klon, klev) ! incrementation eau liquide
  REAL rneb(klon, klev) ! nebulosite
  REAL rain(klon) ! pluies (mm/s)
  REAL snow(klon) ! neige (mm/s)
  INTEGER ibas(klon) ! niveau du bas
  INTEGER itop(klon) ! niveau du haut
  LOGICAL accompli(klon) ! .TRUE. si l'ajustement a eu lieu (Output)

  ! Quelques options:

  LOGICAL new_top ! re-calculer sommet quand re-ajustement est fait
  PARAMETER (new_top=.FALSE.)
  LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
  PARAMETER (evap_prec=.TRUE.)
  REAL coef_eva
  PARAMETER (coef_eva=1.0E-05)
  REAL t_coup
  PARAMETER (t_coup=234.0)
  REAL seuil_vap
  PARAMETER (seuil_vap=1.0E-10)
  LOGICAL old_tau ! implique precip nulle, si vrai.
  PARAMETER (old_tau=.FALSE.)
  REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
  REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
  PARAMETER (dpmin=0.15, tomax=0.97)
  REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
  PARAMETER (dpmax=0.30, tomin=0.05)
  REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
  PARAMETER (deep_sig=0.50, deep_to=0.05)
  LOGICAL exigent ! implique un calcul supplementaire pour Qs
  PARAMETER (exigent=.FALSE.)

  INTEGER kbase
  PARAMETER (kbase=0)

  ! Variables locales:

  INTEGER nexpo
  INTEGER i, k, k1min, k1max, k2min, k2max, is
  REAL zgamdz(klon, klev-1)
  REAL zt(klon, klev), zq(klon, klev)
  REAL zqs(klon, klev), zdqs(klon, klev)
  REAL zqmqsdp(klon, klev)
  REAL ztnew(klon, klev), zqnew(klon, klev)
  REAL zcond(klon), zvapo(klon), zrapp(klon)
  REAL zrfl(klon), zrfln, zqev, zqevt
  REAL zsat(klon) ! sur-saturation
  REAL zflo(klon) ! flotabilite
  REAL za(klon), zb(klon), zc(klon)
  INTEGER k1(klon), k2(klon)
  REAL zdelta, zcor, zcvm5
  REAL delp(klon, klev)
  LOGICAL possible(klon), todo(klon), etendre(klon)
  LOGICAL aller(klon), todobis(klon)
  REAL zalfa
  INTEGER nbtodo, nbdone

  ! Fonctions thermodynamiques:

  include "FCTTRE.h"

  DO k = 1, klev
    DO i = 1, klon
      delp(i, k) = paprs(i, k) - paprs(i, k+1)
    END DO
  END DO

  ! Initialiser les sorties a zero

  DO k = 1, klev
    DO i = 1, klon
      d_t(i, k) = 0.0
      d_q(i, k) = 0.0
      d_ql(i, k) = 0.0
      rneb(i, k) = 0.0
    END DO
  END DO
  DO i = 1, klon
    ibas(i) = klev
    itop(i) = 1
    rain(i) = 0.0
    snow(i) = 0.0
    accompli(i) = .FALSE.
  END DO

  ! Preparations

  DO k = 1, klev
    DO i = 1, klon
      IF (afaire(i)) THEN
        zt(i, k) = t(i, k)
        zq(i, k) = q(i, k)

        ! Calculer Qs et L/Cp*dQs/dT

        IF (thermcep) THEN
          zdelta = max(0., sign(1.,rtt-zt(i,k)))
          zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
          zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i,k))
          zqs(i, k) = r2es*foeew(zt(i,k), zdelta)/pplay(i, k)
          zqs(i, k) = min(0.5, zqs(i,k))
          zcor = 1./(1.-retv*zqs(i,k))
          zqs(i, k) = zqs(i, k)*zcor
          zdqs(i, k) = foede(zt(i,k), zdelta, zcvm5, zqs(i,k), zcor)
        ELSE
          IF (zt(i,k)<t_coup) THEN
            zqs(i, k) = qsats(zt(i,k))/pplay(i, k)
            zdqs(i, k) = dqsats(zt(i,k), zqs(i,k))
          ELSE
            zqs(i, k) = qsatl(zt(i,k))/pplay(i, k)
            zdqs(i, k) = dqsatl(zt(i,k), zqs(i,k))
          END IF
        END IF

        ! Calculer (q-qs)*dp
        zqmqsdp(i, k) = (zq(i,k)-zqs(i,k))*delp(i, k)
      END IF
    END DO
  END DO

  ! -----zgama is the moist convective lapse rate (-dT/dz).
  ! -----zgamdz(*,k) est la difference minimale autorisee des temperatures
  ! -----entre deux couches (k et k+1), c.a.d. si T(k+1)-T(k) est inferieur
  ! -----a zgamdz(*,k), alors ces 2 couches sont instables conditionnellement

  DO k = 1, klev - 1
    DO i = 1, klon
      IF (afaire(i)) THEN
        zgamdz(i, k) = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd*(rd*(zt( &
          i,k)*delp(i,k)+zt(i,k+1)*delp(i,k+1))/(delp(i,k)+delp(i, &
          k+1))+rlvtt*(zqs(i,k)*delp(i,k)+zqs(i,k+1)*delp(i,k+1))/(delp(i, &
          k)+delp(i,k+1)))/(1.0+(zdqs(i,k)*delp(i,k)+zdqs(i,k+1)*delp(i, &
          k+1))/(delp(i,k)+delp(i,k+1)))
      END IF
    END DO
  END DO

  ! On cherche la presence simultanee d'instabilite conditionnelle
  ! et de sur-saturation. Sinon, pas la penne de se fatiguer:

  DO i = 1, klon
    possible(i) = .FALSE.
  END DO
  DO k = 2, klev
    DO i = 1, klon
      IF (afaire(i)) THEN
        zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
        zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k-1)
        IF (zflo(i)>0.0 .AND. zsat(i)>0.0) possible(i) = .TRUE.
      END IF
    END DO
  END DO

  DO i = 1, klon
    IF (possible(i)) THEN
      k1(i) = kbase
      k2(i) = k1(i) + 1
    END IF
  END DO

810 CONTINUE ! chercher le bas de la colonne a ajuster

  k2min = klev
  DO i = 1, klon
    todo(i) = .FALSE.
    aller(i) = .TRUE.
    IF (possible(i)) k2min = min(k2min, k2(i))
  END DO
  IF (k2min==klev) GO TO 860
  DO k = k2min, klev - 1
    DO i = 1, klon
      IF (possible(i) .AND. k>=k2(i) .AND. aller(i)) THEN
        zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k+1)
        zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k+1)
        IF (zflo(i)>0.0 .AND. zsat(i)>0.0) THEN
          k1(i) = k
          k2(i) = k + 1
          todo(i) = .TRUE.
          aller(i) = .FALSE.
        END IF
      END IF
    END DO
  END DO
  DO i = 1, klon
    IF (possible(i) .AND. aller(i)) THEN
      todo(i) = .FALSE.
      k1(i) = klev
      k2(i) = klev
    END IF
  END DO

  ! CC      DO i = 1, klon
  ! CC      IF (possible(i)) THEN
  ! CC  811    k2(i) = k2(i) + 1
  ! CC         IF (k2(i) .GT. klev) THEN
  ! CC            todo(i) = .FALSE.
  ! CC            GOTO 812
  ! CC         ENDIF
  ! CC         k = k2(i)
  ! CC         zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
  ! CC         zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1)
  ! CC         IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 811
  ! CC         k1(i) = k2(i) - 1
  ! CC         todo(i) = .TRUE.
  ! CC      ENDIF
  ! CC  812 CONTINUE
  ! CC      ENDDO

820 CONTINUE ! chercher le haut de la colonne

  k2min = klev
  DO i = 1, klon
    aller(i) = .TRUE.
    IF (todo(i)) k2min = min(k2min, k2(i))
  END DO
  IF (k2min<klev) THEN
    DO k = k2min, klev
      DO i = 1, klon
        IF (todo(i) .AND. k>k2(i) .AND. aller(i)) THEN
          zsat(i) = zsat(i) + zqmqsdp(i, k)
          zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
          IF (zflo(i)<=0.0 .OR. zsat(i)<=0.0) THEN
            aller(i) = .FALSE.
          ELSE
            k2(i) = k
          END IF
        END IF
      END DO
    END DO
    ! error      is = 0
    ! error      DO i = 1, klon
    ! error      IF(todo(i).AND.aller(i)) THEN
    ! error         is = is + 1
    ! error         todo(i) = .FALSE.
    ! error         k2(i) = klev
    ! error      ENDIF
    ! error      ENDDO
    ! error      IF (is.GT.0) THEN
    ! error         PRINT*, "Bizard. je pourrais continuer mais j arrete"
    ! error         CALL abort
    ! error      ENDIF
  END IF

  ! CC      DO i = 1, klon
  ! CC      IF (todo(i)) THEN
  ! CC  821    CONTINUE
  ! CC         IF (k2(i) .EQ. klev) GOTO 822
  ! CC         k = k2(i) + 1
  ! CC         zsat(i) = zsat(i) + zqmqsdp(i,k)
  ! CC         zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
  ! CC         IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 822
  ! CC         k2(i) = k
  ! CC         GOTO 821
  ! CC      ENDIF
  ! CC  822 CONTINUE
  ! CC      ENDDO

830 CONTINUE ! faire l'ajustement en sachant k1 et k2

  is = 0
  DO i = 1, klon
    IF (todo(i)) THEN
      IF (k2(i)<=k1(i)) is = is + 1
    END IF
  END DO
  IF (is>0) THEN
    PRINT *, 'Impossible: k1 trop grand ou k2 trop petit'
    PRINT *, 'is=', is
    CALL abort
  END IF

  k1min = klev
  k1max = 1
  k2max = 1
  DO i = 1, klon
    IF (todo(i)) THEN
      k1min = min(k1min, k1(i))
      k1max = max(k1max, k1(i))
      k2max = max(k2max, k2(i))
    END IF
  END DO

  DO i = 1, klon
    IF (todo(i)) THEN
      k = k1(i)
      za(i) = 0.
      zb(i) = (rcpd*(1.+zdqs(i,k))*(zt(i,k)-za(i))-rlvtt*(zqs(i,k)-zq(i, &
        k)))*delp(i, k)
      zc(i) = delp(i, k)*rcpd*(1.+zdqs(i,k))
    END IF
  END DO

  DO k = k1min, k2max
    DO i = 1, klon
      IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) THEN
        za(i) = za(i) + zgamdz(i, k-1)
        zb(i) = zb(i) + (rcpd*(1.+zdqs(i,k))*(zt(i,k)-za(i))-rlvtt*(zqs(i, &
          k)-zq(i,k)))*delp(i, k)
        zc(i) = zc(i) + delp(i, k)*rcpd*(1.+zdqs(i,k))
      END IF
    END DO
  END DO

  DO i = 1, klon
    IF (todo(i)) THEN
      k = k1(i)
      ztnew(i, k) = zb(i)/zc(i)
      zqnew(i, k) = zqs(i, k) + (ztnew(i,k)-zt(i,k))*rcpd/rlvtt*zdqs(i, k)
    END IF
  END DO

  DO k = k1min, k2max
    DO i = 1, klon
      IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) THEN
        ztnew(i, k) = ztnew(i, k-1) + zgamdz(i, k-1)
        zqnew(i, k) = zqs(i, k) + (ztnew(i,k)-zt(i,k))*rcpd/rlvtt*zdqs(i, k)
      END IF
    END DO
  END DO

  ! Quantite de condensation produite pendant l'ajustement:

  DO i = 1, klon
    zcond(i) = 0.0
  END DO
  DO k = k1min, k2max
    DO i = 1, klon
      IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
        rneb(i, k) = 1.0
        zcond(i) = zcond(i) + (zq(i,k)-zqnew(i,k))*delp(i, k)/rg
      END IF
    END DO
  END DO

  ! Si condensation negative, effort completement perdu:

  DO i = 1, klon
    IF (todo(i) .AND. zcond(i)<=0.) todo(i) = .FALSE.
  END DO

  ! L'ajustement a ete accompli, meme les calculs accessoires
  ! ne sont pas encore faits:

  DO i = 1, klon
    IF (todo(i)) accompli(i) = .TRUE.
  END DO

  ! =====
  ! Une fois que la condensation a lieu, on doit construire un
  ! "modele nuageux" pour partager la condensation entre l'eau
  ! liquide nuageuse et la precipitation (leur rapport toliq
  ! est calcule selon l'epaisseur nuageuse). Je suppose que
  ! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
  ! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
  ! lineaire entre dpmin et dpmax).
  ! =====
  DO i = 1, klon
    IF (todo(i)) THEN
      toliq(i) = tomax - ((paprs(i,k1(i))-paprs(i,k2(i)+1))/paprs(i,1)-dpmin) &
        *(tomax-tomin)/(dpmax-dpmin)
      toliq(i) = max(tomin, min(tomax,toliq(i)))
      IF (pplay(i,k2(i))/paprs(i,1)<=deep_sig) toliq(i) = deep_to
      IF (old_tau) toliq(i) = 1.0
    END IF
  END DO
  ! =====
  ! On doit aussi determiner la distribution verticale de
  ! l'eau nuageuse. Plusieurs options sont proposees:

  ! (0) La condensation precipite integralement (toliq ne sera
  ! pas utilise).
  ! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
  ! a la vapeur d'eau locale.
  ! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
  ! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
  ! est effectivement diminuee pendant le processus d'ajustement.
  ! (4) Elle est en fonction (lineaire ou exponentielle) de la
  ! distance (epaisseur en pression) avec le niveau k1 (la couche
  ! k1 n'aura donc pas d'eau liquide).
  ! =====

  IF (opt_cld==0) THEN

    DO i = 1, klon
      IF (todo(i)) zrfl(i) = zcond(i)/dtime
    END DO

  ELSE IF (opt_cld==1) THEN

    DO i = 1, klon
      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
    END DO
    DO k = k1min, k2max
      DO i = 1, klon
        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
          zqnew(i, k)*delp(i, k)/rg
      END DO
    END DO
    DO i = 1, klon
      IF (todo(i)) THEN
        zrapp(i) = toliq(i)*zcond(i)/zvapo(i)
        zrapp(i) = max(0., min(1.,zrapp(i)))
        zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
      END IF
    END DO
    DO k = k1min, k2max
      DO i = 1, klon
        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
          d_ql(i, k) = d_ql(i, k) + zrapp(i)*zqnew(i, k)
        END IF
      END DO
    END DO

  ELSE IF (opt_cld==2) THEN

    DO i = 1, klon
      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
    END DO
    DO k = k1min, k2max
      DO i = 1, klon
        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
          delp(i, k)/rg
      END DO
    END DO
    DO k = k1min, k2max
      DO i = 1, klon
        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
          d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
    END DO

  ELSE IF (opt_cld==3) THEN

    DO i = 1, klon
      IF (todo(i)) zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
    END DO
    DO k = k1min, k2max
      DO i = 1, klon
        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
          max(0.0, zq(i,k)-zqnew(i,k))*delp(i, k)/rg
      END DO
    END DO
    DO k = k1min, k2max
      DO i = 1, klon
        IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i) .AND. zvapo(i)>0.0) d_ql(i, &
          k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)*max(0.0, zq(i,k)-zqnew &
          (i,k))
      END DO
    END DO
    DO i = 1, klon
      IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
    END DO

  ELSE IF (opt_cld==4) THEN

    nexpo = 3
    ! cc         nexpo = 1 ! distribution lineaire

    DO i = 1, klon
      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
    END DO ! (avec ponderation)
    DO k = k1min, k2max
      DO i = 1, klon
        IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) zvapo(i) = zvapo(i) + &
          delp(i, k)/rg*(pplay(i,k1(i))-pplay(i,k))**nexpo
      END DO
    END DO
    DO k = k1min, k2max
      DO i = 1, klon
        IF (todo(i) .AND. k>=(k1(i)+1) .AND. k<=k2(i)) d_ql(i, k) = d_ql(i, &
          k) + toliq(i)*zcond(i)/zvapo(i)*(pplay(i,k1(i))-pplay(i,k))**nexpo
      END DO
    END DO
    DO i = 1, klon
      IF (todo(i)) zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
    END DO

  ELSE ! valeur non-prevue pour opt_cld

    PRINT *, 'opt_cld est faux:', opt_cld
    CALL abort

  END IF ! fin de opt_cld

  ! L'eau precipitante peut etre evaporee:

  zalfa = 0.05
  IF (evap_prec .AND. (k1max>=2)) THEN
    DO k = k1max - 1, 1, -1
      DO i = 1, klon
        IF (todo(i) .AND. k<k1(i) .AND. zrfl(i)>0.0) THEN
          zqev = max(0.0, (zqs(i,k)-zq(i,k))*zalfa)
          zqevt = coef_eva*(1.0-zq(i,k)/zqs(i,k))*sqrt(zrfl(i))*delp(i, k)/ &
            pplay(i, k)*zt(i, k)*rd/rg
          zqevt = max(0.0, min(zqevt,zrfl(i)))*rg*dtime/delp(i, k)
          zqev = min(zqev, zqevt)
          zrfln = zrfl(i) - zqev*(delp(i,k))/rg/dtime
          zq(i, k) = zq(i, k) - (zrfln-zrfl(i))*(rg/(delp(i,k)))*dtime
          zt(i, k) = zt(i, k) + (zrfln-zrfl(i))*(rg/(delp(i, &
            k)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i,k))
          zrfl(i) = zrfln
        END IF
      END DO
    END DO
  END IF

  ! La temperature de la premiere couche determine la pluie ou la neige:

  DO i = 1, klon
    IF (todo(i)) THEN
      IF (zt(i,1)>rtt) THEN
        rain(i) = rain(i) + zrfl(i)
      ELSE
        snow(i) = snow(i) + zrfl(i)
      END IF
    END IF
  END DO

  ! Mise a jour de la temperature et de l'humidite

  DO k = k1min, k2max
    DO i = 1, klon
      IF (todo(i) .AND. k>=k1(i) .AND. k<=k2(i)) THEN
        zt(i, k) = ztnew(i, k)
        zq(i, k) = zqnew(i, k)
      END IF
    END DO
  END DO

  ! Re-calculer certaines variables pour etendre et re-ajuster la colonne

  IF (exigent) THEN
    DO k = 1, klev
      DO i = 1, klon
        IF (todo(i)) THEN
          IF (thermcep) THEN
            zdelta = max(0., sign(1.,rtt-zt(i,k)))
            zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
            zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i,k))
            zqs(i, k) = r2es*foeew(zt(i,k), zdelta)/pplay(i, k)
            zqs(i, k) = min(0.5, zqs(i,k))
            zcor = 1./(1.-retv*zqs(i,k))
            zqs(i, k) = zqs(i, k)*zcor
            zdqs(i, k) = foede(zt(i,k), zdelta, zcvm5, zqs(i,k), zcor)
          ELSE
            IF (zt(i,k)<t_coup) THEN
              zqs(i, k) = qsats(zt(i,k))/pplay(i, k)
              zdqs(i, k) = dqsats(zt(i,k), zqs(i,k))
            ELSE
              zqs(i, k) = qsatl(zt(i,k))/pplay(i, k)
              zdqs(i, k) = dqsatl(zt(i,k), zqs(i,k))
            END IF
          END IF
        END IF
      END DO
    END DO
  END IF

  IF (exigent) THEN
    DO k = 1, klev - 1
      DO i = 1, klon
        IF (todo(i)) THEN
          zgamdz(i, k) = -(pplay(i,k)-pplay(i,k+1))/paprs(i, k+1)/rcpd*(rd*( &
            zt(i,k)*delp(i,k)+zt(i,k+1)*delp(i,k+1))/(delp(i,k)+delp(i, &
            k+1))+rlvtt*(zqs(i,k)*delp(i,k)+zqs(i,k+1)*delp(i,k+1))/(delp(i, &
            k)+delp(i,k+1)))/(1.0+(zdqs(i,k)*delp(i,k)+zdqs(i,k+1)*delp(i, &
            k+1))/(delp(i,k)+delp(i,k+1)))
        END IF
      END DO
    END DO
  END IF

  ! Puisque l'humidite a ete modifiee, on re-fait (q-qs)*dp

  DO k = 1, klev
    DO i = 1, klon
      IF (todo(i)) THEN
        zqmqsdp(i, k) = (zq(i,k)-zqs(i,k))*delp(i, k)
      END IF
    END DO
  END DO

  ! Verifier si l'on peut etendre le bas de la colonne

  DO i = 1, klon
    etendre(i) = .FALSE.
  END DO

  k1max = 1
  DO i = 1, klon
    IF (todo(i) .AND. k1(i)>(kbase+1)) THEN
      k = k1(i)
      zflo(i) = zt(i, k-1) + zgamdz(i, k-1) - zt(i, k)
      zsat(i) = zqmqsdp(i, k) + zqmqsdp(i, k-1)
      ! sc voici l'ancienne ligne:
      ! sc         IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) THEN
      ! sc sylvain: il faut RESPECTER les 2 criteres:
      IF (zflo(i)>0.0 .AND. zsat(i)>0.0) THEN
        etendre(i) = .TRUE.
        k1(i) = k1(i) - 1
        k1max = max(k1max, k1(i))
        aller(i) = .TRUE.
      END IF
    END IF
  END DO

  IF (k1max>(kbase+1)) THEN
    DO k = k1max, kbase + 1, -1
      DO i = 1, klon
        IF (etendre(i) .AND. k<k1(i) .AND. aller(i)) THEN
          zsat(i) = zsat(i) + zqmqsdp(i, k)
          zflo(i) = zt(i, k) + zgamdz(i, k) - zt(i, k+1)
          IF (zsat(i)<=0.0 .OR. zflo(i)<=0.0) THEN
            aller(i) = .FALSE.
          ELSE
            k1(i) = k
          END IF
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (etendre(i) .AND. aller(i)) THEN
        k1(i) = 1
      END IF
    END DO
  END IF

  ! CC      DO i = 1, klon
  ! CC      IF (etendre(i)) THEN
  ! CC  840    k = k1(i)
  ! CC         IF (k.GT.1) THEN
  ! CC            zsat(i) = zsat(i) + zqmqsdp(i,k-1)
  ! CC            zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k)
  ! CC            IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN
  ! CC               k1(i) = k - 1
  ! CC               GOTO 840
  ! CC            ENDIF
  ! CC         ENDIF
  ! CC      ENDIF
  ! CC      ENDDO

  DO i = 1, klon
    todobis(i) = todo(i)
    todo(i) = .FALSE.
  END DO
  is = 0
  DO i = 1, klon
    IF (etendre(i)) THEN
      todo(i) = .TRUE.
      is = is + 1
    END IF
  END DO
  IF (is>0) THEN
    IF (new_top) THEN
      GO TO 820 ! chercher de nouveau le sommet k2
    ELSE
      GO TO 830 ! supposer que le sommet est celui deja trouve
    END IF
  END IF

  DO i = 1, klon
    possible(i) = .FALSE.
  END DO
  is = 0
  DO i = 1, klon
    IF (todobis(i) .AND. k2(i)<klev) THEN
      is = is + 1
      possible(i) = .TRUE.
    END IF
  END DO
  IF (is>0) GO TO 810 !on cherche en haut d'autres blocks
  ! a ajuster a partir du sommet de la colonne precedente

860 CONTINUE ! Calculer les tendances et diagnostiques
  ! cc      print*, "Apres 860"

  DO k = 1, klev
    DO i = 1, klon
      IF (accompli(i)) THEN
        d_t(i, k) = zt(i, k) - t(i, k)
        zq(i, k) = max(zq(i,k), seuil_vap)
        d_q(i, k) = zq(i, k) - q(i, k)
      END IF
    END DO
  END DO

  DO i = 1, klon
    IF (accompli(i)) THEN
      DO k = 1, klev
        IF (rneb(i,k)>0.0) THEN
          ibas(i) = k
          GO TO 807
        END IF
      END DO
807   CONTINUE
      DO k = klev, 1, -1
        IF (rneb(i,k)>0.0) THEN
          itop(i) = k
          GO TO 808
        END IF
      END DO
808   CONTINUE
    END IF
  END DO

  IF (imprim) THEN
    nbtodo = 0
    nbdone = 0
    DO i = 1, klon
      IF (afaire(i)) nbtodo = nbtodo + 1
      IF (accompli(i)) nbdone = nbdone + 1
    END DO
    PRINT *, 'nbTodo, nbDone=', nbtodo, nbdone
  END IF

  RETURN
END SUBROUTINE conmanv
SUBROUTINE conkuo(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
    rain, snow, ibas, itop)
  USE dimphy
  USE yomcst_mod_h
  USE yoethf_mod_h
IMPLICIT NONE
  ! ======================================================================
  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
  ! Objet: Schema de convection de type Kuo (1965).
  ! Cette version du code peut calculer le niveau de depart
  ! N.B. version vectorielle (le 6 oct. 1997)
  ! ======================================================================


  ! Arguments:

  REAL dtime ! intervalle du temps (s)
  REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
  REAL t(klon, klev) ! temperature (K)
  REAL q(klon, klev) ! humidite specifique
  REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s)

  REAL d_t(klon, klev) ! incrementation temperature
  REAL d_q(klon, klev) ! incrementation humidite
  REAL d_ql(klon, klev) ! incrementation eau liquide
  REAL rneb(klon, klev) ! nebulosite
  REAL rain(klon) ! pluies (mm/s)
  REAL snow(klon) ! neige (mm/s)
  INTEGER itop(klon) ! niveau du sommet
  INTEGER ibas(klon) ! niveau du bas

  LOGICAL ldcum(klon) ! convection existe
  LOGICAL todo(klon)

  ! Quelsques options:

  LOGICAL calcfcl ! calculer le niveau de convection libre
  PARAMETER (calcfcl=.TRUE.)
  INTEGER ldepar ! niveau fixe de convection libre
  PARAMETER (ldepar=4)
  INTEGER opt_cld ! comment traiter l'eau liquide
  PARAMETER (opt_cld=4) ! valeur possible: 0, 1, 2, 3 ou 4
  LOGICAL evap_prec ! evaporation de pluie au-dessous de convection
  PARAMETER (evap_prec=.TRUE.)
  REAL coef_eva
  PARAMETER (coef_eva=1.0E-05)
  LOGICAL new_deh ! nouvelle facon de calculer dH
  PARAMETER (new_deh=.FALSE.)
  REAL t_coup
  PARAMETER (t_coup=234.0)
  LOGICAL old_tau ! implique precipitation nulle
  PARAMETER (old_tau=.FALSE.)
  REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante
  REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande
  PARAMETER (dpmin=0.15, tomax=0.97)
  REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible
  PARAMETER (dpmax=0.30, tomin=0.05)
  REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to
  PARAMETER (deep_sig=0.50, deep_to=0.05)

  ! Variables locales:

  INTEGER nexpo
  LOGICAL nuage(klon)
  INTEGER i, k, kbmin, kbmax, khmax
  REAL ztotal(klon, klev), zdeh(klon, klev)
  REAL zgz(klon, klev)
  REAL zqs(klon, klev)
  REAL zdqs(klon, klev)
  REAL ztemp(klon, klev)
  REAL zpres(klon, klev)
  REAL zconv(klon) ! convergence d'humidite
  REAL zvirt(klon) ! convergence virtuelle d'humidite
  REAL zfrac(klon) ! fraction convective
  INTEGER kb(klon), kh(klon)

  REAL zcond(klon), zvapo(klon), zrapp(klon)
  REAL zrfl(klon), zrfln, zqev, zqevt
  REAL zdelta, zcvm5, zcor
  REAL zvar

  LOGICAL appel1er
  SAVE appel1er
  !$OMP THREADPRIVATE(appel1er)

  ! Fonctions thermodynamiques

  include "FCTTRE.h"

  DATA appel1er/.TRUE./

  IF (appel1er) THEN
    PRINT *, 'conkuo, calcfcl:', calcfcl
    IF (.NOT. calcfcl) PRINT *, 'conkuo, ldepar:', ldepar
    PRINT *, 'conkuo, opt_cld:', opt_cld
    PRINT *, 'conkuo, evap_prec:', evap_prec
    PRINT *, 'conkuo, new_deh:', new_deh
    appel1er = .FALSE.
  END IF

  ! Initialiser les sorties a zero

  DO k = 1, klev
    DO i = 1, klon
      d_q(i, k) = 0.0
      d_t(i, k) = 0.0
      d_ql(i, k) = 0.0
      rneb(i, k) = 0.0
    END DO
  END DO
  DO i = 1, klon
    rain(i) = 0.0
    snow(i) = 0.0
    ibas(i) = 0
    itop(i) = 0
  END DO

  ! Calculer la vapeur d'eau saturante Qs et sa derive L/Cp * dQs/dT

  DO k = 1, klev
    DO i = 1, klon
      IF (thermcep) THEN
        zdelta = max(0., sign(1.,rtt-t(i,k)))
        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
        zqs(i, k) = r2es*foeew(t(i,k), zdelta)/pplay(i, k)
        zqs(i, k) = min(0.5, zqs(i,k))
        zcor = 1./(1.-retv*zqs(i,k))
        zqs(i, k) = zqs(i, k)*zcor
        zdqs(i, k) = foede(t(i,k), zdelta, zcvm5, zqs(i,k), zcor)
      ELSE
        IF (t(i,k)<t_coup) THEN
          zqs(i, k) = qsats(t(i,k))/pplay(i, k)
          zdqs(i, k) = dqsats(t(i,k), zqs(i,k))
        ELSE
          zqs(i, k) = qsatl(t(i,k))/pplay(i, k)
          zdqs(i, k) = dqsatl(t(i,k), zqs(i,k))
        END IF
      END IF
    END DO
  END DO

  ! Calculer gz (energie potentielle)

  DO i = 1, klon
    zgz(i, 1) = rd*t(i, 1)/(0.5*(paprs(i,1)+pplay(i, &
      1)))*(paprs(i,1)-pplay(i,1))
  END DO
  DO k = 2, klev
    DO i = 1, klon
      zgz(i, k) = zgz(i, k-1) + rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i &
        ,k-1)-pplay(i,k))
    END DO
  END DO

  ! Calculer l'energie statique humide saturee (Cp*T + gz + L*Qs)

  DO k = 1, klev
    DO i = 1, klon
      ztotal(i, k) = rcpd*t(i, k) + rlvtt*zqs(i, k) + zgz(i, k)
    END DO
  END DO

  ! Determiner le niveau de depart et calculer la difference de
  ! l'energie statique humide saturee (ztotal) entre la couche
  ! de depart et chaque couche au-dessus.

  IF (calcfcl) THEN
    DO k = 1, klev
      DO i = 1, klon
        zpres(i, k) = pplay(i, k)
        ztemp(i, k) = t(i, k)
      END DO
    END DO
    CALL kuofcl(ztemp, q, zgz, zpres, ldcum, kb)
    DO i = 1, klon
      IF (ldcum(i)) THEN
        k = kb(i)
        IF (new_deh) THEN
          zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
        ELSE
          zdeh(i, k) = rcpd*(t(i,k-1)-t(i,k)) - rd*0.5*(t(i,k-1)+t(i,k))/ &
            paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
            rlvtt*(zqs(i,k-1)-zqs(i,k))
        END IF
        zdeh(i, k) = zdeh(i, k)*0.5
      END IF
    END DO
    DO k = 1, klev
      DO i = 1, klon
        IF (ldcum(i) .AND. k>=(kb(i)+1)) THEN
          IF (new_deh) THEN
            zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
          ELSE
            zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
              rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)* &
              (pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
          END IF
        END IF
      END DO
    END DO
  ELSE
    DO i = 1, klon
      k = ldepar
      kb(i) = ldepar
      ldcum(i) = .TRUE.
      IF (new_deh) THEN
        zdeh(i, k) = ztotal(i, k-1) - ztotal(i, k)
      ELSE
        zdeh(i, k) = rcpd*(t(i,k-1)-t(i,k)) - rd*0.5*(t(i,k-1)+t(i,k))/paprs( &
          i, k)*(pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
      END IF
      zdeh(i, k) = zdeh(i, k)*0.5
    END DO
    DO k = ldepar + 1, klev
      DO i = 1, klon
        IF (new_deh) THEN
          zdeh(i, k) = zdeh(i, k-1) + (ztotal(i,k-1)-ztotal(i,k))
        ELSE
          zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
            rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
            rlvtt*(zqs(i,k-1)-zqs(i,k))
        END IF
      END DO
    END DO
  END IF

  ! -----Chercher le sommet du nuage
  ! -----Calculer la convergence de l'humidite (en kg/m**2 a un facteur
  ! -----psolpa/RG pres) du bas jusqu'au sommet du nuage.
  ! -----Calculer la convergence virtuelle pour que toute la maille
  ! -----deviennt nuageuse (du bas jusqu'au sommet du nuage)

  DO i = 1, klon
    nuage(i) = .TRUE.
    zconv(i) = 0.0
    zvirt(i) = 0.0
    kh(i) = -999
  END DO
  DO k = 1, klev
    DO i = 1, klon
      IF (k>=kb(i) .AND. ldcum(i)) THEN
        nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
        IF (nuage(i)) THEN
          kh(i) = k
          zconv(i) = zconv(i) + conv_q(i, k)*dtime*(paprs(i,k)-paprs(i,k+1))
          zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k) &
            -paprs(i,k+1))
        END IF
      END IF
    END DO
  END DO

  DO i = 1, klon
    todo(i) = ldcum(i) .AND. kh(i) > kb(i) .AND. zconv(i) > 0.0
  END DO

  kbmin = klev
  kbmax = 0
  khmax = 0
  DO i = 1, klon
    IF (todo(i)) THEN
      kbmin = min(kbmin, kb(i))
      kbmax = max(kbmax, kb(i))
      khmax = max(khmax, kh(i))
    END IF
  END DO

  ! -----Calculer la surface couverte par le nuage

  DO i = 1, klon
    IF (todo(i)) THEN
      zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
    END IF
  END DO

  ! -----Calculs essentiels:

  DO i = 1, klon
    IF (todo(i)) THEN
      zcond(i) = 0.0
    END IF
  END DO
  DO k = kbmin, khmax
    DO i = 1, klon
      IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
        zvar = zdeh(i, k)/(1.+zdqs(i,k))
        d_t(i, k) = zvar*zfrac(i)/rcpd
        d_q(i, k) = (zvar*zdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))*zfrac(i) - &
          conv_q(i, k)*dtime
        zcond(i) = zcond(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
        rneb(i, k) = zfrac(i)
      END IF
    END DO
  END DO

  DO i = 1, klon
    IF (todo(i) .AND. zcond(i)<0.0) THEN
      PRINT *, 'WARNING: cond. negative (Kuo) ', i, kb(i), kh(i), zcond(i)
      zcond(i) = 0.0
      DO k = kb(i), kh(i)
        d_t(i, k) = 0.0
        d_q(i, k) = 0.0
      END DO
      todo(i) = .FALSE. ! effort totalement perdu
    END IF
  END DO

  ! =====
  ! Une fois que la condensation a lieu, on doit construire un
  ! "modele nuageux" pour partager la condensation entre l'eau
  ! liquide nuageuse et la precipitation (leur rapport toliq
  ! est calcule selon l'epaisseur nuageuse). Je suppose que
  ! toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin,
  ! et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation
  ! lineaire entre dpmin et dpmax).
  ! =====
  DO i = 1, klon
    IF (todo(i)) THEN
      toliq(i) = tomax - ((paprs(i,kb(i))-paprs(i,kh(i)+1))/paprs(i,1)-dpmin) &
        *(tomax-tomin)/(dpmax-dpmin)
      toliq(i) = max(tomin, min(tomax,toliq(i)))
      IF (pplay(i,kh(i))/paprs(i,1)<=deep_sig) toliq(i) = deep_to
      IF (old_tau) toliq(i) = 1.0
    END IF
  END DO
  ! =====
  ! On doit aussi determiner la distribution verticale de
  ! l'eau nuageuse. Plusieurs options sont proposees:

  ! (0) La condensation precipite integralement (toliq ne sera
  ! pas utilise).
  ! (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle
  ! a la vapeur d'eau locale.
  ! (2) Elle est distribuee entre k1 et k2 avec une valeur constante.
  ! (3) Elle est seulement distribuee aux couches ou la vapeur d'eau
  ! est effectivement diminuee pendant le processus d'ajustement.
  ! (4) Elle est en fonction (lineaire ou exponentielle) de la
  ! distance (epaisseur en pression) avec le niveau k1 (la couche
  ! k1 n'aura donc pas d'eau liquide).
  ! =====

  IF (opt_cld==0) THEN

    DO i = 1, klon
      IF (todo(i)) zrfl(i) = zcond(i)/dtime
    END DO

  ELSE IF (opt_cld==1) THEN

    DO i = 1, klon
      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau
    END DO
    DO k = kbmin, khmax
      DO i = 1, klon
        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
          zvapo(i) = zvapo(i) + (q(i,k)+d_q(i,k))*(paprs(i,k)-paprs(i,k+1))/ &
            rg
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (todo(i)) THEN
        zrapp(i) = toliq(i)*zcond(i)/zvapo(i)
        zrapp(i) = max(0., min(1.,zrapp(i)))
      END IF
    END DO
    DO k = kbmin, khmax
      DO i = 1, klon
        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
          d_ql(i, k) = zrapp(i)*(q(i,k)+d_q(i,k))
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (todo(i)) THEN
        zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
      END IF
    END DO

  ELSE IF (opt_cld==2) THEN

    DO i = 1, klon
      IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse
    END DO
    DO k = kbmin, khmax
      DO i = 1, klon
        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
          zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg
        END IF
      END DO
    END DO
    DO k = kbmin, khmax
      DO i = 1, klon
        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
          d_ql(i, k) = toliq(i)*zcond(i)/zvapo(i)
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (todo(i)) THEN
        zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
      END IF
    END DO

  ELSE IF (opt_cld==3) THEN

    DO i = 1, klon
      IF (todo(i)) THEN
        zvapo(i) = 0.0 ! quantite de l'eau strictement condensee
      END IF
    END DO
    DO k = kbmin, khmax
      DO i = 1, klon
        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i)) THEN
          zvapo(i) = zvapo(i) + max(0.0, -d_q(i,k))*(paprs(i,k)-paprs(i,k+1)) &
            /rg
        END IF
      END DO
    END DO
    DO k = kbmin, khmax
      DO i = 1, klon
        IF (todo(i) .AND. k>=kb(i) .AND. k<=kh(i) .AND. zvapo(i)>0.0) THEN
          d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)*max(0.0, -d_q( &
            i,k))
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (todo(i)) THEN
        zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
      END IF
    END DO

  ELSE IF (opt_cld==4) THEN

    nexpo = 3
    ! cc         nexpo = 1 ! distribution lineaire

    DO i = 1, klon
      IF (todo(i)) THEN
        zvapo(i) = 0.0 ! quantite integrale de masse (avec ponderation)
      END IF
    END DO
    DO k = kbmin, khmax
      DO i = 1, klon
        IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
          zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/rg*(pplay(i,kb(i))- &
            pplay(i,k))**nexpo
        END IF
      END DO
    END DO
    DO k = kbmin, khmax
      DO i = 1, klon
        IF (todo(i) .AND. k>=(kb(i)+1) .AND. k<=kh(i)) THEN
          d_ql(i, k) = d_ql(i, k) + toliq(i)*zcond(i)/zvapo(i)*(pplay(i,kb(i) &
            )-pplay(i,k))**nexpo
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (todo(i)) THEN
        zrfl(i) = (1.0-toliq(i))*zcond(i)/dtime
      END IF
    END DO

  ELSE ! valeur non-prevue pour opt_cld

    PRINT *, 'opt_cld est faux:', opt_cld
    CALL abort

  END IF ! fin de opt_cld

  ! L'eau precipitante peut etre re-evaporee:

  IF (evap_prec .AND. kbmax>=2) THEN
    DO k = kbmax, 1, -1
      DO i = 1, klon
        IF (todo(i) .AND. k<=(kb(i)-1) .AND. zrfl(i)>0.0) THEN
          zqev = max(0.0, (zqs(i,k)-q(i,k))*zfrac(i))
          zqevt = coef_eva*(1.0-q(i,k)/zqs(i,k))*sqrt(zrfl(i))* &
            (paprs(i,k)-paprs(i,k+1))/pplay(i, k)*t(i, k)*rd/rg
          zqevt = max(0.0, min(zqevt,zrfl(i)))*rg*dtime/ &
            (paprs(i,k)-paprs(i,k+1))
          zqev = min(zqev, zqevt)
          zrfln = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime
          d_q(i, k) = -(zrfln-zrfl(i))*(rg/(paprs(i,k)-paprs(i,k+1)))*dtime
          d_t(i, k) = (zrfln-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
            k+1)))*dtime*rlvtt/rcpd
          zrfl(i) = zrfln
        END IF
      END DO
    END DO
  END IF

  ! La temperature de la premiere couche determine la pluie ou la neige:

  DO i = 1, klon
    IF (todo(i)) THEN
      IF (t(i,1)>rtt) THEN
        rain(i) = rain(i) + zrfl(i)
      ELSE
        snow(i) = snow(i) + zrfl(i)
      END IF
    END IF
  END DO

  RETURN
END SUBROUTINE conkuo
SUBROUTINE kuofcl(pt, pq, pg, pp, ldcum, kcbot)
  USE yomcst_mod_h
  USE dimphy
  USE yoethf_mod_h
  IMPLICIT NONE
  ! ======================================================================
  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
  ! adaptation du code de Tiedtke du ECMWF
  ! Objet: calculer le niveau de convection libre
  ! (FCL: Free Convection Level)
  ! ======================================================================
  ! Arguments:
  ! pt---input-R- temperature (K)
  ! pq---input-R- vapeur d'eau (kg/kg)
  ! pg---input-R- geopotentiel (g*z ou z est en metre)
  ! pp---input-R- pression (Pa)

  ! LDCUM---output-L- Y-t-il la convection
  ! kcbot---output-I- Niveau du bas de la convection
  ! ======================================================================

  REAL pt(klon, klev), pq(klon, klev), pg(klon, klev), pp(klon, klev)
  INTEGER kcbot(klon)
  LOGICAL ldcum(klon)

  REAL ztu(klon, klev), zqu(klon, klev), zlu(klon, klev)
  REAL zqold(klon), zbuo
  INTEGER is, i, k

  ! klab=1: on est sous le nuage convectif
  ! klab=2: le bas du nuage convectif
  ! klab=0: autres couches
  INTEGER klab(klon, klev)

  ! quand lflag=.true., on est sous le nuage, il faut donc appliquer
  ! le processus d'elevation.
  LOGICAL lflag(klon)

  DO k = 1, klev
    DO i = 1, klon
      ztu(i, k) = pt(i, k)
      zqu(i, k) = pq(i, k)
      zlu(i, k) = 0.0
      klab(i, k) = 0
    END DO
  END DO
  ! ----------------------------------------------------------------------
  DO i = 1, klon
    klab(i, 1) = 1
    kcbot(i) = 2
    ldcum(i) = .FALSE.
  END DO

  DO k = 2, klev - 1

    is = 0
    DO i = 1, klon
      IF (klab(i,k-1)==1) is = is + 1
      lflag(i) = .FALSE.
      IF (klab(i,k-1)==1) lflag(i) = .TRUE.
    END DO
    IF (is==0) GO TO 290

    ! on eleve le parcel d'air selon l'adiabatique sec

    DO i = 1, klon
      IF (lflag(i)) THEN
        zqu(i, k) = zqu(i, k-1)
        ztu(i, k) = ztu(i, k-1) + (pg(i,k-1)-pg(i,k))/rcpd
        zbuo = ztu(i, k)*(1.+retv*zqu(i,k)) - pt(i, k)*(1.+retv*pq(i,k)) + &
          0.5
        IF (zbuo>0.) klab(i, k) = 1
        zqold(i) = zqu(i, k)
      END IF
    END DO

    ! on calcule la condensation eventuelle

    CALL adjtq(pp(1,k), ztu(1,k), zqu(1,k), lflag, 1)

    ! s'il y a la condensation et la "buoyancy" force est positive
    ! c'est bien le bas de la tour de convection

    DO i = 1, klon
      IF (lflag(i) .AND. zqu(i,k)/=zqold(i)) THEN
        klab(i, k) = 2
        zlu(i, k) = zlu(i, k) + zqold(i) - zqu(i, k)
        zbuo = ztu(i, k)*(1.+retv*zqu(i,k)) - pt(i, k)*(1.+retv*pq(i,k)) + &
          0.5
        IF (zbuo>0.) THEN
          kcbot(i) = k
          ldcum(i) = .TRUE.
        END IF
      END IF
    END DO

290 END DO

  RETURN
END SUBROUTINE kuofcl
SUBROUTINE adjtq(pp, pt, pq, ldflag, kcall)
  USE yomcst_mod_h
  USE dimphy
  USE yoethf_mod_h
  IMPLICIT NONE
  ! ======================================================================
  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927
  ! adaptation du code de Tiedtke du ECMWF
  ! Objet: ajustement entre T et Q
  ! ======================================================================
  ! Arguments:
  ! pp---input-R- pression (Pa)
  ! pt---input/output-R- temperature (K)
  ! pq---input/output-R- vapeur d'eau (kg/kg)
  ! ======================================================================
  ! TO PRODUCE T,Q AND L VALUES FOR CLOUD ASCENT

  ! NOTE: INPUT PARAMETER KCALL DEFINES CALCULATION AS
  ! KCALL=0    ENV. T AND QS IN*CUINI*
  ! KCALL=1  CONDENSATION IN UPDRAFTS  (E.G. CUBASE, CUASC)
  ! KCALL=2  EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF)


  REAL pt(klon), pq(klon), pp(klon)
  LOGICAL ldflag(klon)
  INTEGER kcall

  REAL t_coup
  PARAMETER (t_coup=234.0)

  REAL zcond(klon), zcond1
  REAL zdelta, zcvm5, zldcp, zqsat, zcor, zdqsat
  INTEGER is, i
  include "FCTTRE.h"

  DO i = 1, klon
    zcond(i) = 0.0
  END DO

  DO i = 1, klon
    IF (ldflag(i)) THEN
      zdelta = max(0., sign(1.,rtt-pt(i)))
      zldcp = rlvtt*(1.-zdelta) + zdelta*rlstt
      zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
      IF (thermcep) THEN
        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
        zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
        zqsat = min(0.5, zqsat)
        zcor = 1./(1.-retv*zqsat)
        zqsat = zqsat*zcor
        zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
      ELSE
        IF (pt(i)<t_coup) THEN
          zqsat = qsats(pt(i))/pp(i)
          zdqsat = dqsats(pt(i), zqsat)
        ELSE
          zqsat = qsatl(pt(i))/pp(i)
          zdqsat = dqsatl(pt(i), zqsat)
        END IF
      END IF
      zcond(i) = (pq(i)-zqsat)/(1.+zdqsat)
      IF (kcall==1) zcond(i) = max(zcond(i), 0.)
      IF (kcall==2) zcond(i) = min(zcond(i), 0.)
      pt(i) = pt(i) + zldcp*zcond(i)
      pq(i) = pq(i) - zcond(i)
    END IF
  END DO

  is = 0
  DO i = 1, klon
    IF (zcond(i)/=0.) is = is + 1
  END DO
  IF (is==0) GO TO 230

  DO i = 1, klon
    IF (ldflag(i) .AND. zcond(i)/=0.) THEN
      zdelta = max(0., sign(1.,rtt-pt(i)))
      zldcp = rlvtt*(1.-zdelta) + zdelta*rlstt
      zldcp = zldcp/rcpd/(1.0+rvtmp2*pq(i))
      IF (thermcep) THEN
        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*pq(i))
        zqsat = r2es*foeew(pt(i), zdelta)/pp(i)
        zqsat = min(0.5, zqsat)
        zcor = 1./(1.-retv*zqsat)
        zqsat = zqsat*zcor
        zdqsat = foede(pt(i), zdelta, zcvm5, zqsat, zcor)
      ELSE
        IF (pt(i)<t_coup) THEN
          zqsat = qsats(pt(i))/pp(i)
          zdqsat = dqsats(pt(i), zqsat)
        ELSE
          zqsat = qsatl(pt(i))/pp(i)
          zdqsat = dqsatl(pt(i), zqsat)
        END IF
      END IF
      zcond1 = (pq(i)-zqsat)/(1.+zdqsat)
      pt(i) = pt(i) + zldcp*zcond1
      pq(i) = pq(i) - zcond1
    END IF
  END DO

230 CONTINUE
  RETURN
END SUBROUTINE adjtq
SUBROUTINE fiajh(dtime, paprs, pplay, t, q, d_t, d_q, d_ql, rneb, rain, snow, &
    ibas, itop)
  USE dimphy
  USE yomcst_mod_h
  USE yoethf_mod_h
IMPLICIT NONE

  ! Ajustement humide (Schema de convection de Manabe)
  ! .


  ! Arguments:

  REAL dtime ! intervalle du temps (s)
  REAL t(klon, klev) ! temperature (K)
  REAL q(klon, klev) ! humidite specifique (kg/kg)
  REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)

  REAL d_t(klon, klev) ! incrementation pour la temperature
  REAL d_q(klon, klev) ! incrementation pour vapeur d'eau
  REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
  REAL rneb(klon, klev) ! fraction nuageuse

  REAL rain(klon) ! variable non utilisee
  REAL snow(klon) ! variable non utilisee
  INTEGER ibas(klon) ! variable non utilisee
  INTEGER itop(klon) ! variable non utilisee

  REAL t_coup
  PARAMETER (t_coup=234.0)
  REAL seuil_vap
  PARAMETER (seuil_vap=1.0E-10)

  ! Variables locales:

  INTEGER i, k
  INTEGER k1, k1p, k2, k2p
  LOGICAL itest(klon)
  REAL delta_q(klon, klev)
  REAL cp_new_t(klev)
  REAL cp_delta_t(klev)
  REAL new_qb(klev)
  REAL v_cptj(klev), v_cptjk1, v_ssig
  REAL v_cptt(klon, klev), v_p, v_t
  REAL v_qs(klon, klev), v_qsd(klon, klev)
  REAL zq1(klon), zq2(klon)
  REAL gamcpdz(klon, 2:klev)
  REAL zdp, zdpm

  REAL zsat ! sur-saturation
  REAL zflo ! flotabilite

  REAL local_q(klon, klev), local_t(klon, klev)

  REAL zdelta, zcor, zcvm5

  include "FCTTRE.h"

  DO k = 1, klev
    DO i = 1, klon
      local_q(i, k) = q(i, k)
      local_t(i, k) = t(i, k)
      rneb(i, k) = 0.0
      d_ql(i, k) = 0.0
      d_t(i, k) = 0.0
      d_q(i, k) = 0.0
    END DO
  END DO
  DO i = 1, klon
    rain(i) = 0.0
    snow(i) = 0.0
    ibas(i) = 0
    itop(i) = 0
  END DO

  ! Calculer v_qs et v_qsd:

  DO k = 1, klev
    DO i = 1, klon
      v_cptt(i, k) = rcpd*local_t(i, k)
      v_t = local_t(i, k)
      v_p = pplay(i, k)

      IF (thermcep) THEN
        zdelta = max(0., sign(1.,rtt-v_t))
        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
        v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
        v_qs(i, k) = min(0.5, v_qs(i,k))
        zcor = 1./(1.-retv*v_qs(i,k))
        v_qs(i, k) = v_qs(i, k)*zcor
        v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
      ELSE
        IF (v_t<t_coup) THEN
          v_qs(i, k) = qsats(v_t)/v_p
          v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
        ELSE
          v_qs(i, k) = qsatl(v_t)/v_p
          v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
        END IF
      END IF
    END DO
  END DO

  ! Calculer Gamma * Cp * dz: (gamm est le gradient critique)

  DO k = 2, klev
    DO i = 1, klon
      zdp = paprs(i, k) - paprs(i, k+1)
      zdpm = paprs(i, k-1) - paprs(i, k)
      gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)*(v_cptt(i,k-1)*zdpm+ &
        v_cptt(i,k)*zdp)+rlvtt/(zdpm+zdp)*(v_qs(i,k-1)*zdpm+ &
        v_qs(i,k)*zdp))*(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
        k-1)*zdpm+v_qsd(i,k)*zdp)/(zdpm+zdp))
    END DO
  END DO

  ! ------------------------------------ modification des profils instables
  DO i = 1, klon
    itest(i) = .FALSE.

    k1 = 0
    k2 = 1

810 CONTINUE ! chercher k1, le bas de la colonne
    k2 = k2 + 1
    IF (k2>klev) GO TO 9999
    zflo = v_cptt(i, k2-1) - v_cptt(i, k2) - gamcpdz(i, k2)
    zsat = (local_q(i,k2-1)-v_qs(i,k2-1))*(paprs(i,k2-1)-paprs(i,k2)) + &
      (local_q(i,k2)-v_qs(i,k2))*(paprs(i,k2)-paprs(i,k2+1))
    IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 810
    k1 = k2 - 1
    itest(i) = .TRUE.

820 CONTINUE ! chercher k2, le haut de la colonne
    IF (k2==klev) GO TO 821
    k2p = k2 + 1
    zsat = zsat + (paprs(i,k2p)-paprs(i,k2p+1))*(local_q(i,k2p)-v_qs(i,k2p))
    zflo = v_cptt(i, k2p-1) - v_cptt(i, k2p) - gamcpdz(i, k2p)
    IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 821
    k2 = k2p
    GO TO 820
821 CONTINUE

    ! ------------------------------------------------------ ajustement local
830 CONTINUE ! ajustement proprement dit
    v_cptj(k1) = 0.0
    zdp = paprs(i, k1) - paprs(i, k1+1)
    v_cptjk1 = ((1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1))+rlvtt*(local_q(i, &
      k1)-v_qs(i,k1)))*zdp
    v_ssig = zdp*(1.0+v_qsd(i,k1))

    k1p = k1 + 1
    DO k = k1p, k2
      zdp = paprs(i, k) - paprs(i, k+1)
      v_cptj(k) = v_cptj(k-1) + gamcpdz(i, k)
      v_cptjk1 = v_cptjk1 + zdp*((1.0+v_qsd(i,k))*(v_cptt(i, &
        k)+v_cptj(k))+rlvtt*(local_q(i,k)-v_qs(i,k)))
      v_ssig = v_ssig + zdp*(1.0+v_qsd(i,k))
    END DO

    DO k = k1, k2
      cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
      cp_delta_t(k) = cp_new_t(k) - v_cptt(i, k)
      new_qb(k) = v_qs(i, k) + v_qsd(i, k)*cp_delta_t(k)/rlvtt
      local_q(i, k) = new_qb(k)
      local_t(i, k) = cp_new_t(k)/rcpd
    END DO

    ! --------------------------------------------------- sondage vers le bas
    ! -- on redefinit les variables prognostiques dans
    ! -- la colonne qui vient d'etre ajustee

    DO k = k1, k2
      v_cptt(i, k) = rcpd*local_t(i, k)
      v_t = local_t(i, k)
      v_p = pplay(i, k)

      IF (thermcep) THEN
        zdelta = max(0., sign(1.,rtt-v_t))
        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*local_q(i,k))
        v_qs(i, k) = r2es*foeew(v_t, zdelta)/v_p
        v_qs(i, k) = min(0.5, v_qs(i,k))
        zcor = 1./(1.-retv*v_qs(i,k))
        v_qs(i, k) = v_qs(i, k)*zcor
        v_qsd(i, k) = foede(v_t, zdelta, zcvm5, v_qs(i,k), zcor)
      ELSE
        IF (v_t<t_coup) THEN
          v_qs(i, k) = qsats(v_t)/v_p
          v_qsd(i, k) = dqsats(v_t, v_qs(i,k))
        ELSE
          v_qs(i, k) = qsatl(v_t)/v_p
          v_qsd(i, k) = dqsatl(v_t, v_qs(i,k))
        END IF
      END IF
    END DO
    DO k = 2, klev
      zdpm = paprs(i, k-1) - paprs(i, k)
      zdp = paprs(i, k) - paprs(i, k+1)
      gamcpdz(i, k) = ((rd/rcpd/(zdpm+zdp)*(v_cptt(i,k-1)*zdpm+ &
        v_cptt(i,k)*zdp)+rlvtt/(zdpm+zdp)*(v_qs(i,k-1)*zdpm+ &
        v_qs(i,k)*zdp))*(pplay(i,k-1)-pplay(i,k))/paprs(i,k))/(1.0+(v_qsd(i, &
        k-1)*zdpm+v_qsd(i,k)*zdp)/(zdpm+zdp))
    END DO

    ! Verifier si l'on peut etendre la colonne vers le bas

    IF (k1==1) GO TO 841 ! extension echouee
    zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
    zsat = (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1)) + &
      (local_q(i,k1)-v_qs(i,k1))*(paprs(i,k1)-paprs(i,k1+1))
    IF (zflo<=0.0 .OR. zsat<=0.0) GO TO 841 ! extension echouee

840 CONTINUE
    k1 = k1 - 1
    IF (k1==1) GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
    zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1))
    zflo = v_cptt(i, k1-1) - v_cptt(i, k1) - gamcpdz(i, k1)
    IF (zflo>0.0 .AND. zsat>0.0) THEN
      GO TO 840
    ELSE
      GO TO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
    END IF
841 CONTINUE

    GO TO 810 ! chercher d'autres blocks en haut

9999 END DO ! boucle sur tous les points
  ! -----------------------------------------------------------------------

  ! Determiner la fraction nuageuse (hypothese: la nebulosite a lieu
  ! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):

  DO k = 1, klev
    DO i = 1, klon
      IF (itest(i)) THEN
        delta_q(i, k) = local_q(i, k) - q(i, k)
        IF (delta_q(i,k)<0.) rneb(i, k) = 1.0
      END IF
    END DO
  END DO

  ! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
  ! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
  ! diminue et d'une maniere proportionnelle a cet diminution):

  DO i = 1, klon
    IF (itest(i)) THEN
      zq1(i) = 0.0
      zq2(i) = 0.0
    END IF
  END DO
  DO k = 1, klev
    DO i = 1, klon
      IF (itest(i)) THEN
        zdp = paprs(i, k) - paprs(i, k+1)
        zq1(i) = zq1(i) - delta_q(i, k)*zdp
        zq2(i) = zq2(i) - min(0.0, delta_q(i,k))*zdp
      END IF
    END DO
  END DO
  DO k = 1, klev
    DO i = 1, klon
      IF (itest(i)) THEN
        IF (zq2(i)/=0.0) d_ql(i, k) = -min(0.0, delta_q(i,k))*zq1(i)/zq2(i)
      END IF
    END DO
  END DO

  DO k = 1, klev
    DO i = 1, klon
      local_q(i, k) = max(local_q(i,k), seuil_vap)
    END DO
  END DO

  DO k = 1, klev
    DO i = 1, klon
      d_t(i, k) = local_t(i, k) - t(i, k)
      d_q(i, k) = local_q(i, k) - q(i, k)
    END DO
  END DO

  RETURN
END SUBROUTINE fiajh
SUBROUTINE fiajc(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, d_ql, rneb, &
    rain, snow, ibas, itop)
  USE dimphy
  USE yomcst_mod_h
  USE yoethf_mod_h
IMPLICIT NONE



  ! Options:

  INTEGER plb ! niveau de depart pour la convection
  PARAMETER (plb=4)

  ! Mystere: cette option n'est pas innocente pour les resultats !
  ! Qui peut resoudre ce mystere ? (Z.X.Li mars 1995)
  LOGICAL vector ! calcul vectorise
  PARAMETER (vector=.FALSE.)

  REAL t_coup
  PARAMETER (t_coup=234.0)

  ! Arguments:

  REAL q(klon, klev) ! humidite specifique (kg/kg)
  REAL t(klon, klev) ! temperature (K)
  REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
  REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
  REAL dtime ! intervalle du temps (s)
  REAL conv_q(klon, klev) ! taux de convergence de l'humidite
  REAL rneb(klon, klev) ! fraction nuageuse
  REAL d_q(klon, klev) ! incrementaion pour la vapeur d'eau
  REAL d_ql(klon, klev) ! incrementation pour l'eau liquide
  REAL d_t(klon, klev) ! incrementation pour la temperature
  REAL rain(klon) ! variable non-utilisee
  REAL snow(klon) ! variable non-utilisee
  INTEGER itop(klon) ! variable non-utilisee
  INTEGER ibas(klon) ! variable non-utilisee

  INTEGER kh(klon), i, k
  LOGICAL nuage(klon), test(klon, klev)
  REAL zconv(klon), zdeh(klon, klev), zvirt(klon)
  REAL zdqs(klon, klev), zqs(klon, klev)
  REAL ztt, zvar, zfrac(klon)
  REAL zq1(klon), zq2(klon)
  REAL zdelta, zcor, zcvm5

  include "FCTTRE.h"

  ! Initialiser les sorties:

  DO k = 1, klev
    DO i = 1, klon
      rneb(i, k) = 0.0
      d_ql(i, k) = 0.0
      d_t(i, k) = 0.0
      d_q(i, k) = 0.0
    END DO
  END DO
  DO i = 1, klon
    itop(i) = 0
    ibas(i) = 0
    rain(i) = 0.0
    snow(i) = 0.0
  END DO

  ! Calculer Qs et L/Cp * dQs/dT:

  DO k = 1, klev
    DO i = 1, klon
      ztt = t(i, k)
      IF (thermcep) THEN
        zdelta = max(0., sign(1.,rtt-ztt))
        zcvm5 = r5les*rlvtt*(1.-zdelta) + zdelta*r5ies*rlstt
        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*q(i,k))
        zqs(i, k) = r2es*foeew(ztt, zdelta)/pplay(i, k)
        zqs(i, k) = min(0.5, zqs(i,k))
        zcor = 1./(1.-retv*zqs(i,k))
        zqs(i, k) = zqs(i, k)*zcor
        zdqs(i, k) = foede(ztt, zdelta, zcvm5, zqs(i,k), zcor)
      ELSE
        IF (ztt<t_coup) THEN
          zqs(i, k) = qsats(ztt)/pplay(i, k)
          zdqs(i, k) = dqsats(ztt, zqs(i,k))
        ELSE
          zqs(i, k) = qsatl(ztt)/pplay(i, k)
          zdqs(i, k) = dqsatl(ztt, zqs(i,k))
        END IF
      END IF
    END DO
  END DO

  ! Determiner la difference de l'energie totale saturee:

  DO i = 1, klon
    k = plb
    zdeh(i, k) = rcpd*(t(i,k-1)-t(i,k)) - rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k &
      )*(pplay(i,k-1)-pplay(i,k)) + rlvtt*(zqs(i,k-1)-zqs(i,k))
    zdeh(i, k) = zdeh(i, k)*0.5 ! on prend la moitie
  END DO
  DO k = plb + 1, klev
    DO i = 1, klon
      zdeh(i, k) = zdeh(i, k-1) + rcpd*(t(i,k-1)-t(i,k)) - &
        rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) + &
        rlvtt*(zqs(i,k-1)-zqs(i,k))
    END DO
  END DO

  ! Determiner le sommet du nuage selon l'instabilite
  ! Calculer les convergences d'humidite (reelle et virtuelle)

  DO i = 1, klon
    nuage(i) = .TRUE.
    zconv(i) = 0.0
    zvirt(i) = 0.0
    kh(i) = -999
  END DO
  DO k = plb, klev
    DO i = 1, klon
      nuage(i) = nuage(i) .AND. zdeh(i, k) > 0.0
      IF (nuage(i)) THEN
        kh(i) = k
        zconv(i) = zconv(i) + conv_q(i, k)*dtime*(paprs(i,k)-paprs(i,k+1))
        zvirt(i) = zvirt(i) + (zdeh(i,k)/rlvtt+zqs(i,k)-q(i,k))*(paprs(i,k)- &
          paprs(i,k+1))
      END IF
    END DO
  END DO

  IF (vector) THEN


    DO k = plb, klev
      DO i = 1, klon
        IF (k<=kh(i) .AND. kh(i)>plb .AND. zconv(i)>0.0) THEN
          test(i, k) = .TRUE.
          zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
        ELSE
          test(i, k) = .FALSE.
        END IF
      END DO
    END DO

    DO k = plb, klev
      DO i = 1, klon
        IF (test(i,k)) THEN
          zvar = zdeh(i, k)/(1.0+zdqs(i,k))
          d_q(i, k) = (zvar*zdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))*zfrac(i) - &
            conv_q(i, k)*dtime
          d_t(i, k) = zvar*zfrac(i)/rcpd
        END IF
      END DO
    END DO

    DO i = 1, klon
      zq1(i) = 0.0
      zq2(i) = 0.0
    END DO
    DO k = plb, klev
      DO i = 1, klon
        IF (test(i,k)) THEN
          IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
          zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
          zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
        END IF
      END DO
    END DO

    DO k = plb, klev
      DO i = 1, klon
        IF (test(i,k)) THEN
          IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
        END IF
      END DO
    END DO

  ELSE ! (.NOT. vector)

    DO i = 1, klon
      IF (kh(i)>plb .AND. zconv(i)>0.0) THEN
        ! cc         IF (kh(i).LE.plb) GOTO 999 ! il n'y a pas d'instabilite
        ! cc         IF (zconv(i).LE.0.0) GOTO 999 ! convergence insuffisante
        zfrac(i) = max(0.0, min(zconv(i)/zvirt(i),1.0))
        DO k = plb, kh(i)
          zvar = zdeh(i, k)/(1.0+zdqs(i,k))
          d_q(i, k) = (zvar*zdqs(i,k)/rlvtt+zqs(i,k)-q(i,k))*zfrac(i) - &
            conv_q(i, k)*dtime
          d_t(i, k) = zvar*zfrac(i)/rcpd
        END DO

        zq1(i) = 0.0
        zq2(i) = 0.0
        DO k = plb, kh(i)
          IF (d_q(i,k)<0.0) rneb(i, k) = zfrac(i)
          zq1(i) = zq1(i) - d_q(i, k)*(paprs(i,k)-paprs(i,k+1))
          zq2(i) = zq2(i) - min(0.0, d_q(i,k))*(paprs(i,k)-paprs(i,k+1))
        END DO
        DO k = plb, kh(i)
          IF (zq2(i)/=0.) d_ql(i, k) = -min(0.0, d_q(i,k))*zq1(i)/zq2(i)
        END DO
      END IF
    END DO

  END IF ! fin de teste sur vector

  RETURN
END SUBROUTINE fiajc