Index: /LMDZ5/trunk/libf/phylmd/fisrtilp.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 2806)
+++ /LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 2807)
@@ -17,4 +17,7 @@
   USE cloudth_mod 
   USE ioipsl_getin_p_mod, ONLY : getin_p
+  USE phys_local_var_mod, ONLY: ql_seri,qs_seri
+  ! flag to include modifications to ensure energy conservation (if flag >0)
+  USE add_phys_tend_mod, only : fl_cor_ebil 
   IMPLICIT none
   !======================================================================
@@ -26,8 +29,20 @@
   !             et ilp (il pleut, L. Li)
   ! Principales parties:
+  ! P0> Thermalisation des precipitations venant de la couche du dessus
   ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
   ! P2> Formation du nuage (en k) 
+  ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
+  ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour 
+  ! les valeurs de T et Q initiales
+  ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
+  ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
+  ! P2.B> Nuage "tout ou rien"
+  ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
   ! P3> Formation de la precipitation (en k)
   !======================================================================
+  ! JLD: 
+  ! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
+  ! * fl_cor_ebil doit etre > 0 ; 
+  !   fl_cor_ebil= 0 pour reproduire anciens bugs
   !======================================================================
   include "YOMCST.h"
@@ -115,4 +130,5 @@
   INTEGER i, k, n, kk
   REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5  
+  REAL zdqsdT_raw(klon)
   REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
   LOGICAL convergence(klon)
@@ -140,6 +156,8 @@
 !!!!
 !  Variables pour Bergeron
-  REAL zcp, coef1, DeltaT
+  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
   REAL zqpreci(klon), zqprecl(klon)
+! Variable pour conservation enegie des precipitations
+  REAL zmqc(klon)
   !
   LOGICAL appel1er
@@ -170,5 +188,5 @@
   REAL   beta(klon,klev) ! taux de conversion de l'eau cond
 ! RomP <<<
-  REAL zmair, zcpair, zcpeau
+  REAL zmair(klon), zcpair, zcpeau
   !     Pour la conversion eau-neige
   REAL zlh_solid(klon), zm_solid
@@ -180,4 +198,5 @@
                      ! (Heymsfield & Donner, 1990)
   REAL zzz
+
   include "YOETHF.h"
   include "FCTTRE.h"
@@ -312,31 +331,55 @@
      ENDDO
      !
+     ! ----------------------------------------------------------------
+     ! P0> Thermalisation des precipitations venant de la couche du dessus
+     ! ----------------------------------------------------------------
      ! Calculer la varition de temp. de l'air du a la chaleur sensible
-     ! transporter par la pluie.
-     ! Il resterait a rajouter cet effet de la chaleur sensible sur les
-     ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
-     ! surface.
+     ! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
+     ! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
+     ! des differentes transformations thermodynamiques. Cette masse d'eau doit
+     ! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
+     ! de l'enthalpie  de la couche avec la temperature
+     ! Variables calculees ou modifiees:
+     !   -  zt: temperature de la cocuhe
+     !   - zmqc: masse de precip qui doit etre thermalisee
      !
      IF(k.LE.klevm1) THEN         
         DO i = 1, klon
            !IM
-           zmair=(paprs(i,k)-paprs(i,k+1))/RG
+           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+           ! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
            zcpair=RCPD*(1.0+RVTMP2*zq(i))
            zcpeau=RCPD*RVTMP2
+         if (fl_cor_ebil .GT. 0) then
+           ! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm 
+           ! pour s'assurer que la precip arrivant au sol aura bien la temperature de la 
+           ! derniere couche
+           zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
+           ! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
+           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
+                 / (zcpair + zmqc(i)*zcpeau)
+         else ! si on maintient les anciennes erreurs
            zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
-                + zmair*zcpair*zt(i) ) &
-                / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
-           !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
+                + zmair(i)*zcpair*zt(i) ) &
+                / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
+         end if
         ENDDO
-     ENDIF
-     ! ----------------------------------------------------------------
-     ! P1> Debut evaporation de la precipitation
-     ! ----------------------------------------------------------------
+     ENDIF ! end IF(k.LE.klevm1)
+     !
+     ! ----------------------------------------------------------------
+     ! P1> Calcul de l'evaporation de la precipitation
+     ! ----------------------------------------------------------------
+     ! On evapore une partie des precipitations venant de la maille du dessus.
+     ! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
+     ! ce que la fraction de cette couche qui est sous le nuage soit saturee.
+     ! Variables calculees ou modifiees:
+     !   - zrfl et zifl: flux de precip liquide et glace
+     !   - zq, zt: humidite et temperature de la cocuhe
+     !   - zmqc: masse de precip qui doit etre thermalisee
+     !
      IF (evap_prec) THEN
         DO i = 1, klon
-!AJ<
-!!           IF (zrfl(i) .GT.0.) THEN
+!          S'il y a des precipitations
            IF (zrfl(i)+zifl(i).GT.0.) THEN
-!>AJ
               ! Calcul du qsat
               IF (thermcep) THEN
@@ -356,12 +399,41 @@
         ENDDO
 !AJ<
+
        IF (.NOT. ice_thermo) THEN
         DO i = 1, klon
-!AJ<
-!!           IF (zrfl(i) .GT.0.) THEN
+!          S'il y a des precipitations
            IF (zrfl(i)+zifl(i).GT.0.) THEN
-!>AJ
                 ! Evap max pour ne pas saturer la fraction sous le nuage
+                ! Evap max jusqu'à atteindre la saturation dans la partie
+                ! de la maille qui est sous le nuage de la couche du dessus
+                !!! On ne tient compte de cette fraction que sous une seule
+                !!! couche sous le nuage
                 zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
+             ! Ajout de la prise en compte des precip a thermiser
+             ! avec petite reecriture 
+             if  (fl_cor_ebil .GT. 0) then ! nouveau
+                ! Calcul de l'evaporation du flux de precip herite
+                !   d'au-dessus
+                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
+                     * zmair(i)/pplay(i,k)*zt(i)*RD
+                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
+
+                ! Seuil pour ne pas saturer la fraction sous le nuage
+                zqev = MIN (zqev, zqevt)
+                ! Nouveau flux de precip
+                zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
+                ! Aucun flux liquide pour T < t_coup, on reevapore tout.
+                IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
+                  zrfln(i)=0.
+                  zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
+                END IF
+                ! Nouvelle vapeur 
+                zq(i) = zq(i) + zqev
+                zmqc(i) = zmqc(i)-zqev
+                ! Nouvelle temperature (chaleur latente)
+                zt(i) = zt(i) - zqev &
+                     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+!!JLD debut de partie a supprimer a therme
+            else ! if  (fl_cor_ebil .GT. 0)
                 ! Calcul de l'evaporation du flux de precip herite
                 !   d'au-dessus
@@ -384,4 +456,6 @@
                      * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
                      * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
+              end if ! if  (fl_cor_ebil .GT. 0)
+!!JLD fin de partie a supprimer a therme
                 zrfl(i) = zrfln(i)
                 zifl(i) = 0.
@@ -390,24 +464,19 @@
 !
        ELSE ! (.NOT. ice_thermo) 
-!
+!      ================================
+!      Avec thermodynamique de la glace
+!      ================================
         DO i = 1, klon
 !AJ<
-!!           IF (zrfl(i) .GT.0.) THEN
-           IF (zrfl(i)+zifl(i).GT.0.) THEN
-!>AJ
-     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     ! Modification de l'evaporation avec la glace
-     ! Differentiation entre precipitation liquide et solide 
-     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
+!        S'il y a des precipitations
+         IF (zrfl(i)+zifl(i).GT.0.) THEN
      
-     ! Evap max pour ne pas saturer la fraction sous le nuage
+        ! Evap max pour ne pas saturer la fraction sous le nuage
          zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
-     !    zqev0 = MAX (0.0, zqs(i)-zq(i) ) 
-
-     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     ! On differencie qsat pour l'eau et la glace
-     ! Si zdelta=1. --> glace
-     ! Si zdelta=0. --> eau liquide
-     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+         !JAM
+         ! On differencie qsat pour l'eau et la glace
+         ! Si zdelta=1. --> glace
+         ! Si zdelta=0. --> eau liquide
         
          ! Calcul du qsat par rapport a l'eau liquide
@@ -417,6 +486,5 @@
          qsl= qsl*zcor
          
-         ! Calcul de l'evaporation du flux de precip herite
-         !   d'au-dessus
+         ! Calcul de l'evaporation du flux de precip venant du dessus
          ! Formulation en racine du flux de precip
          ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
@@ -425,5 +493,4 @@
          zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
               *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 
-
          
          ! Calcul du qsat par rapport a la glace
@@ -440,21 +507,21 @@
               *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
 
-              
-     !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-     ! Verification sur l'evaporation
-     ! On s'assure qu'on ne sature pas 
-     ! la fraction sous le nuage sinon on
-     ! repartit zqev0 en gardant la proportion
-     ! liquide / glace
-     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        !JAM
+        ! Limitation de l'evaporation. On s'assure qu'on ne sature pas 
+        ! la fraction de la couche sous le nuage sinon on repartit zqev0 
+        ! en conservant la proportion liquide / glace
      
          IF (zqevt+zqevti.GT.zqev0) THEN
-	          zqev=zqev0*zqevt/(zqevt+zqevti)
-	          zqevi=zqev0*zqevti/(zqevt+zqevti)
-             
+            zqev=zqev0*zqevt/(zqevt+zqevti)
+            zqevi=zqev0*zqevti/(zqevt+zqevti)
          ELSE
+!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
+!       liquides et solides meme si on ne sature pas la couche. 
+!       A mon avis, le test est inutile, et il faudrait tout remplacer par:
+!            zqev=zqevt
+!            zqevi=zqevti
              IF (zqevt+zqevti.GT.0.) THEN
-	          zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
-	          zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
+                zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
+                zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
              ELSE
              zqev=0.
@@ -471,4 +538,14 @@
          zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+       if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
+         zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                  * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
+                  + (zifln(i)-zifl(i)) &
+                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                  * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+       else ! sans correction thermalisation des precips
          zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
@@ -477,5 +554,6 @@
                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
                   * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
-     
+       end if
+        ! Nouvelles vaeleurs des precips liquides et solides 
          zrfl(i) = zrfln(i)
          zifl(i) = zifln(i)
@@ -486,16 +564,25 @@
 !jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
            zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
+           ! fraction de la precip solide qui est fondue
            zmelt = MIN(MAX(zmelt,0.),1.)
            ! Fusion de la glace
            zrfl(i)=zrfl(i)+zmelt*zifl(i)
-           zifl(i)=zifl(i)*(1.-zmelt)
-!           print*,zt(i),'octavio1'
+           if (fl_cor_ebil .LE. 0) then
+             ! the following line should not be here. Indeed, if zifl is modified
+             ! now, zifl(i)*zmelt is no more the amount of ice that has melt
+             ! and therefore the change in temperature computed below is wrong
+             zifl(i)=zifl(i)*(1.-zmelt)
+           end if
            ! Chaleur latente de fusion
+        if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
+           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
+                      *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+        else ! sans correction thermalisation des precips
            zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
                       *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
-!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
-!fin CR
-
-
+        end if
+           if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
+             zifl(i)=zifl(i)*(1.-zmelt)
+           end if
 
            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
@@ -515,5 +602,9 @@
            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
            zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+       if  (fl_cor_ebil .GT. 0) then ! nouveau
+           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+       else    
            zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
+       endif
            zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
            zqs(i) = MIN(0.5,zqs(i))
@@ -521,4 +612,6 @@
            zqs(i) = zqs(i)*zcor
            zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
+           zdqsdT_raw(i) = zdqs(i)*  &
+         &         RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
         ENDDO
      ELSE
@@ -547,7 +640,15 @@
      ! P2> Formation du nuage
      ! ----------------------------------------------------------------
+     ! Variables calculees:
+     !   rneb  : fraction nuageuse
+     !   zcond : eau condensee moyenne dans la maille.
+     !   rhcl: humidite relative ciel-clair
+     !   zt : temperature de la maille
+     ! ----------------------------------------------------------------
+     !
      IF (cpartiel) THEN
-
-        !        print*,'Dans partiel k=',k
+        ! -------------------------
+        ! P2.A> Nuage fractionnaire
+        ! -------------------------
         !
         !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
@@ -561,4 +662,7 @@
         !   Version avec les raqts
 
+        ! ----------------------------------------------------------------
+        ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
+        ! ----------------------------------------------------------------
         if (iflag_pdf.eq.0) then
 
@@ -570,7 +674,9 @@
            enddo
 
-        else
-           !
-           !   Version avec les nouvelles PDFs.
+        else !  if (iflag_pdf.eq.0)
+           ! ----------------------------------------------------------------
+           ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour 
+           ! les valeurs de T et Q initiales
+           ! ----------------------------------------------------------------
            do i=1,klon
               if(zq(i).lt.1.e-15) then
@@ -620,4 +726,17 @@
            enddo
 
+        ! ----------------------------------------------------------------
+        ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
+        ! Calcul des grandeurs nuageuses en tenant compte de l'effet de
+        ! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
+        ! qui en dependent. Ce changement de temperature est provisoire, et
+        ! la valeur definitive sera calcule plus tard.
+        ! Variables calculees: 
+        !   rneb : nebulosite
+        !   zcond: eau condensee en moyenne dans la maille
+        ! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
+        ! pas clair, il n'y a probablement pas de prise en compte de l'effet de
+        ! T sur qsat
+        ! ----------------------------------------------------------------
 
 !Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
@@ -629,7 +748,19 @@
                do i=1,klon
 !                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
+                 ! !! convergence = .true. tant que l'on n'a pas converge !!
+                 ! ------------------------------
                  convergence(i)=abs(DT(i)).gt.DDT0
                  if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
-                 Tbef(i)=Tbef(i)+DT(i) 
+                 ! si on n'a pas converge
+                 !
+                 ! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
+                 ! ---------------------------------------------------------------
+                 ! Variables calculees: 
+                 ! rneb : nebulosite
+                 ! zqn : eau condensee, dans le nuage (in cloud water content)
+                 ! zcond: eau condensee en moyenne dans la maille
+                 ! rhcl: humidite relative ciel-clair
+                 !
+                 Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
                  if (.not.ice_thermo) then   
                     zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
@@ -641,12 +772,16 @@
                           zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
                        else
-!avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
+                          !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
                           zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
                        endif
                     endif
                  endif
-                 ! Calcul des PDF lognormales
+                 ! Calcul de rneb, qzn et zcond pour les PDF lognormales
                  zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+               if (fl_cor_ebil .GT. 0) then
+                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+               else
                  zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
+               end if
                  zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
                  zqs(i) = MIN(0.5,zqs(i))
@@ -677,11 +812,21 @@
                enddo ! boucle en i
 
+                 ! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT 
+                 !         due a la condensation.
+                 ! ---------------------------------------------------------------
+                 ! Variables calculees: 
+                 ! DT : variation de temperature due a la condensation 
+
                  if (.not. ice_thermo) then
-
+                 ! --------------------------
                  do i=1,klon
                  if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
 
                  qlbef(i)=max(0.,zqn(i)-zqs(i))
+               if (fl_cor_ebil .GT. 0) then
+                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
+               else
                  num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
+               end if
                  denom=1.+rneb(i,k)*zdqs(i)
                  DT(i)=num/denom
@@ -690,6 +835,6 @@
                  enddo
 
-                 else
-                 ! Iteration pour convergence avec qsat(T)
+                 else ! if (.not. ice_thermo)
+                 ! --------------------------
                  if (iflag_t_glace.ge.1) then
                  CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
@@ -703,9 +848,11 @@
                     zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
                     zfice(i) = zfice(i)**exposant_glace_old
-                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - RTT)
+                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
+          &                     / (t_glace_min_old - RTT)
                  endif
                  
                  if (iflag_t_glace.ge.1) then
-                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) / (t_glace_min - t_glace_max)
+                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
+          &                    / (t_glace_min - t_glace_max)
                  endif
                
@@ -721,7 +868,16 @@
 
                  qlbef(i)=max(0.,zqn(i)-zqs(i))
-                 num=-Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
-                 denom=1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
+               if (fl_cor_ebil .GT. 0) then
+                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
+           &          +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
+                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
+                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)    &
+           &               *qlbef(i)*dzfice(i)
+               else
+                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
+           &         +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
+                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
                          -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
+               end if
                  DT(i)=num/denom
                  n_i(i)=n_i(i)+1
@@ -732,16 +888,13 @@
                  endif !ice_thermo
 
-!                 endif
-!               enddo
-  
-          
              enddo ! iter=1,iflag_fisrtilp_qsat+1
              ! Fin d'iteration pour condensation avec variation de qsat(T) 
              ! -----------------------------------------------------------
-           endif
-
+           endif !  if (iflag_fisrtilp_qsat.ge.0)
+     ! ----------------------------------------------------------------
+     ! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
+     ! ----------------------------------------------------------------
 
         endif ! iflag_pdf
-
 
 !        if (iflag_fisrtilp_qsat.eq.-1) then 
@@ -771,9 +924,17 @@
 
 !        ELSE
-
-        ! Calcul de l'eau in-cloud (zqn),
-        ! moyenne dans la maille (zcond),
-        ! fraction nuageuse (rneb) et
-        ! humidite relative ciel-clair (rhcl)
+        ! ----------------------------------------------------------------
+        ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
+        ! Variables calculees: 
+        !   rneb : nebulosite
+        !   zcond: eau condensee en moyenne dans la maille
+        !   zq : eau vapeur dans la maille
+        !   zt : temperature de la maille
+        !   rhcl: humidite relative ciel-clair
+        ! ----------------------------------------------------------------
+        !
+        ! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb) 
+        ! Calcule de l'eau condensee moyenne dans la maille (zcond),
+        ! et de l'humidite relative ciel-clair (rhcl)
         DO i=1,klon
            IF (rneb(i,k) .LE. 0.0) THEN
@@ -793,9 +954,11 @@
         ENDDO
 
-
 !        ENDIF
 
      ELSE ! de IF (cpartiel) 
-        ! Cas "tout ou rien"
+        ! -------------------------
+        ! P2.B> Nuage "tout ou rien"
+        ! -------------------------
+        ! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
         DO i = 1, klon
            IF (zq(i).GT.zqs(i)) THEN
@@ -806,10 +969,8 @@
            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
         ENDDO
-     ENDIF
-     ! ----------------------------------------------------------------
-     ! Fin de formation du nuage
-     ! ----------------------------------------------------------------
+     ENDIF ! de IF (cpartiel) 
      !
      ! Mise a jour vapeur d'eau
+     ! -------------------------
      DO i = 1, klon
         zq(i) = zq(i) - zcond(i)
@@ -817,6 +978,10 @@
      ENDDO
 !AJ<
-     ! Chaleur latente apres formation nuage
+     ! ------------------------------------
+     ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
      ! -------------------------------------
+     ! Variable calcule:
+     !   zt : temperature de la maille
+     !
      IF (.NOT. ice_thermo) THEN
         if (iflag_fisrtilp_qsat.lt.1) then
@@ -826,5 +991,9 @@
         else if (iflag_fisrtilp_qsat.gt.0) then
            DO i= 1, klon
+    if (fl_cor_ebil .GT. 0) then
+             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
+    else
              zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
+    end if
            ENDDO
         endif
@@ -856,6 +1025,12 @@
                     zfice(i) = zfice(i)**exposant_glace_old
               endif
+        if (fl_cor_ebil .GT. 0) then
+              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & 
+           &             * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
+                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
+        else 
               zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
-                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
+                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
+        end if
            ENDDO
          endif
@@ -863,4 +1038,5 @@
        ENDIF
 !>AJ
+
      ! ----------------------------------------------------------------
      ! P3> Formation des precipitations
@@ -958,4 +1134,8 @@
      !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
      !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
+!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
+!      temporaires et ne doivent pas etre calcule (alors qu'elles le sont
+!      si iflag_bergeron <> 2
+!      A SUPPRIMER A TERME
               zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
               zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
@@ -966,4 +1146,5 @@
         ENDDO  ! i = 1,klon
      ENDDO     ! n = 1,ninter
+
      ! ----------------------------------------------------------------
      !
@@ -994,4 +1175,26 @@
              zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
              zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
+           if (fl_cor_ebil .GT. 0) then
+             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
+             coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
+!            Calcul de DT si toute les precips liquides congelent
+             DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
+!            On ne veut pas que T devienne superieur a la temp. de congelation.
+!            donc que Delta > RTT-zt(i
+             DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
+             zt(i)      = zt(i)      + DeltaT
+!            Eau vaporisee du fait de l'augmentation de T
+             Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
+!            on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
+             zq(i) = zq(i) +  Deltaq
+!            Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
+             zcond(i)   = max( zcond(i)- Deltaq, 0. )
+!            precip liquide qui congele ou qui s'evapore
+             Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
+             zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
+!            bilan eau glacee
+             zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
+           else ! if (fl_cor_ebil .GT. 0)
+!            ancien calcul
              zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
              coef1 = RLMLT*zdqs(i)/RLVTT
@@ -1002,4 +1205,5 @@
              zq(i)      = zq(i)      + zcp/RLVTT*zdqs(i)*DeltaT
              zt(i)      = zt(i)      + DeltaT
+           end if ! if (fl_cor_ebil .GT. 0)
            ENDIF  ! rneb(i,k) .GT. 0.0
          ENDDO
@@ -1021,11 +1225,12 @@
          IF (rneb(i,k).GT.0.0) THEN
 !CR on prend en compte la phase glace
-           if (.not.ice_thermo) then
-           d_ql(i,k) = zoliq(i)
-           d_qi(i,k) = 0.
-           else
+!JLD inutile car on ne passe jamais ici si .not.ice_thermo
+!           if (.not.ice_thermo) then
+!           d_ql(i,k) = zoliq(i)
+!           d_qi(i,k) = 0.
+!           else
            d_ql(i,k) = (1-zfice(i))*zoliq(i)
            d_qi(i,k) = zfice(i)*zoliq(i)
-           endif
+!           endif
 !AJ<
            zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
@@ -1043,6 +1248,11 @@
               zifl(i) = zifl(i)+zrfl(i)
               zrfl(i) = 0.
+           if (fl_cor_ebil .GT. 0) then
+              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
+                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+           else
               zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
                       *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
+           end if
            ENDIF  ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
 !RC   
@@ -1094,5 +1304,4 @@
      ! Fin de formation des precipitations
      ! ----------------------------------------------------------------
-     !
      !
      ! Calculer les tendances de q et de t:
@@ -1202,7 +1411,7 @@
      DO i = 1, klon
         zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
-        zmair=(paprs(i,k)-paprs(i,k+1))/RG
+        zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
         zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
-        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
+        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i))
      END DO
   END DO
Index: /LMDZ5/trunk/libf/phylmd/reevap.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/reevap.F90	(revision 2806)
+++ /LMDZ5/trunk/libf/phylmd/reevap.F90	(revision 2807)
@@ -24,5 +24,5 @@
     ! Re-evaporer l'eau liquide nuageuse
     !
-
+print *,'rrevap ; fl_cor_ebil:',fl_cor_ebil,' iflag_ice_thermo:',iflag_ice_thermo,' RVTMP2',RVTMP2
     DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
        DO i = 1, klon
@@ -46,4 +46,5 @@
 
              zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
+  zdelta = 0.
              zb = MAX(0.0,ql_seri(i,k))
              za = - MAX(0.0,ql_seri(i,k)) &
