Index: /LMDZ6/trunk/libf/phylmd/fisrtilp.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/fisrtilp.F90	(revision 3874)
+++ /LMDZ6/trunk/libf/phylmd/fisrtilp.F90	(revision 3875)
@@ -1,3 +1,2 @@
-!
 ! $Id$
 !
@@ -105,4 +104,10 @@
   !
   REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela
+  !<LTP
+  REAL smallestreal
+  REAL, SAVE :: rain_int_min=0.001 !intensité locale minimum pour la pluie avant diminution de la fraction précipitante associée = 0.001 mm/s
+  !>LTP  
+
+
   !$OMP THREADPRIVATE(seuil_neb)
 
@@ -149,6 +154,15 @@
   REAL qcloud(klon)
  
-  REAL zrfl(klon), zrfln(klon), zqev, zqevt 
+  REAL zrfl(klon), zrfln(klon), zqev, zqevt
+!<LTP
+  REAL zrflclr(klon), zrflcld(klon)
+  REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)
+  REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon) 
+!>LTP
+
   REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
+!<LTP
+  REAL ziflclr(klon), ziflcld(klon) 
+!>LTP
   REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq 
   REAL zoliqp(klon), zoliqi(klon)
@@ -161,4 +175,10 @@
   REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
   REAL zchau      ,zfroi      ,zfice(klon),zneb(klon),znebprecip(klon)
+!<LTP
+  REAL znebprecipclr(klon), znebprecipcld(klon)
+  REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)
+  REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
+!>LTP
+
   REAL zmelt, zpluie, zice
   REAL dzfice(klon)
@@ -219,6 +239,13 @@
 !  ice_thermo = iflag_ice_thermo .GE. 1
 
+  
   itap=itap+1
   znebprecip(:)=0.
+
+!<LTP
+  smallestreal=1.e-9
+  znebprecipclr(:)=0.
+  znebprecipcld(:)=0.
+!>LTP
 
   ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
@@ -232,9 +259,16 @@
      CALL getin_p('iflag_evap_prec',iflag_evap_prec)
      CALL getin_p('seuil_neb',seuil_neb)
+!<LTP   
+     CALL getin_p('rain_int_min',rain_int_min)
+!>LTP
      write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp 
      !
      WRITE(lunout,*) 'fisrtilp, ninter:', ninter
      WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec
+!<LTP   
+     WRITE(lunout,*) 'fisrtilp, rain_int_min:', rain_int_min
+!>LTP   
      WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
+     WRITE(lunout,*) 'FISRTILP VERSION LUDO'
      
      IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
@@ -303,4 +337,5 @@
 
   !cdir collapse
+
   DO k = 1, klev
      DO i = 1, klon
@@ -326,4 +361,14 @@
      zrfl(i) = 0.0
      zifl(i) = 0.0
+!<LTP
+     zrflclr(i) = 0.0
+     ziflclr(i) = 0.0
+     zrflcld(i) = 0.0
+     ziflcld(i) = 0.0
+     tot_zneb(i) = 0.0
+     tot_znebn(i) = 0.0
+     d_tot_zneb(i) = 0.0
+!>LTP
+
      zneb(i) = seuil_neb
   ENDDO
@@ -492,7 +537,18 @@
 !      ================================
         DO i = 1, klon
+
+
 !AJ<
 !        S'il y a des precipitations
          IF (zrfl(i)+zifl(i).GT.0.) THEN
+
+        !LTP<
+        !On ne tient compte que du flux de précipitation en ciel clair dans le calcul de l'évaporation.
+                IF (iflag_evap_prec==4) THEN
+                        zrfl(i) = zrflclr(i)
+                        zifl(i) = ziflclr(i)
+                ENDIF
+        
+        !>LTP
 
          IF (iflag_evap_prec==1) THEN
@@ -501,7 +557,12 @@
             znebprecip(i)=MAX(zneb(i),znebprecip(i))
          ENDIF
-     
+         
+         IF (iflag_evap_prec==4) THEN
+        ! Evap max pour ne pas saturer toute la maille
+         zqev0 = MAX (0.0, zqs(i)-zq(i))
+         ELSE 
         ! Evap max pour ne pas saturer la fraction sous le nuage
          zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
+         ENDIF
 
          !JAM
@@ -523,5 +584,11 @@
               *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
               *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-         ELSE
+!<LTP
+         ELSE IF (iflag_evap_prec==4) THEN
+         zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl) &
+              *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
+              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+!>LTP
+        ELSE 
          zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
               *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
@@ -544,4 +611,10 @@
               *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
               *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+!<LTP
+         ELSE IF (iflag_evap_prec==4) THEN
+         zqevti = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsi) &
+              *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
+              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+!>LTP
          ELSE
          zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
@@ -551,4 +624,5 @@
               *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
 
+       
         !JAM
         ! Limitation de l'evaporation. On s'assure qu'on ne sature pas 
@@ -573,4 +647,5 @@
              ENDIF
          ENDIF
+
          ! Nouveaux flux de precip liquide et solide 
          zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
@@ -602,4 +677,18 @@
          zrfl(i) = zrfln(i)
          zifl(i) = zifln(i)
+
+!<LTP
+        IF (iflag_evap_prec==4) THEN
+                zrflclr(i) = zrfl(i)
+                ziflclr(i) = zifl(i)    
+                IF(zrflclr(i) + ziflclr(i) .LE. 0) THEN
+                        znebprecipclr(i) = 0.
+                ENDIF   
+                zrfl(i) = zrflclr(i) + zrflcld(i)
+                zifl(i) = ziflclr(i) + ziflcld(i) 
+        ENDIF
+!>LTP        
+
+
 !        print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1)
 
@@ -612,5 +701,13 @@
            zmelt = MIN(MAX(zmelt,0.),1.)
            ! Fusion de la glace
-           zrfl(i)=zrfl(i)+zmelt*zifl(i)
+!<LTP
+           IF (iflag_evap_prec==4) THEN
+                   zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
+                   zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
+                   zrfl(i)=zrflclr(i)+zrflcld(i)
+!>LTP        
+           ELSE
+                   zrfl(i)=zrfl(i)+zmelt*zifl(i)
+           ENDIF
            if (fl_cor_ebil .LE. 0) then
              ! the following line should not be here. Indeed, if zifl is modified
@@ -628,5 +725,13 @@
         end if
            if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
-             zifl(i)=zifl(i)*(1.-zmelt)
+!<LTP
+             IF (iflag_evap_prec==4) THEN
+                   ziflclr(i)=ziflclr(i)*(1.-zmelt)
+                   ziflcld(i)=ziflcld(i)*(1.-zmelt)
+                   zifl(i)=ziflclr(i)+ziflcld(i)
+!>LTP        
+             ELSE
+                   zifl(i)=zifl(i)*(1.-zmelt)
+             ENDIF
            end if
 
@@ -817,5 +922,5 @@
                  ! Variables calculees: 
                  ! rneb : nebulosite
-                 ! zqn : eau totale dans le nuage (in cloud total water content)
+                 ! zqn : eau condensee, dans le nuage (in cloud water content)
                  ! zcond: eau condensee en moyenne dans la maille
                  ! rhcl: humidite relative ciel-clair
@@ -1019,4 +1124,6 @@
            ENDIF
         ENDDO
+
+       
         ! If vertical heterogeneity, change fraction by volume as well
         if (iflag_cloudth_vert>=3) then
@@ -1116,4 +1223,56 @@
      ! Partager l'eau condensee en precipitation et eau liquide nuageuse
      !
+
+!<LTP
+
+IF (iflag_evap_prec==4) THEN
+        !Partitionnement des precipitations venant du dessus en précipitations nuageuses
+        !et précipitations ciel clair
+
+        !0) Calculate tot_zneb, la fraction nuageuse totale au-dessus du nuage
+        !en supposant un recouvrement maximum aléatoire (voir Jakob and Klein, 2000)
+        
+        DO i=1, klon
+                tot_znebn(i) = 1 - (1-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
+                        /(1-min(zneb(i),1-smallestreal))
+                d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
+                tot_zneb(i) = tot_znebn(i)
+
+
+                !1) Cloudy to clear air
+                d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i)) 
+                IF (znebprecipcld(i) .GT. 0) THEN
+                        d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i) 
+                        d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i) 
+                ELSE 
+                        d_zrfl_cld_clr(i) = 0. 
+                        d_zifl_cld_clr(i) = 0. 
+                ENDIF
+
+                !2) Clear to cloudy air
+                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) &
+                        - d_tot_zneb(i) - zneb(i)))
+                IF (znebprecipclr(i) .GT. 0) THEN
+                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i) 
+                        d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
+                ELSE
+                        d_zrfl_clr_cld(i) = 0. 
+                        d_zifl_clr_cld(i) = 0. 
+                ENDIF 
+
+                !Update variables
+                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)  
+                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i) 
+                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i) 
+                ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
+                zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
+                ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
+
+        ENDDO
+ENDIF
+
+!>LTP
+
+
 
      ! Initialisation de zoliq (eau condensee moyenne dans la maille)
@@ -1293,8 +1452,22 @@
              d_ql(i,k) = (1-zfice(i))*zoliq(i)
              d_qi(i,k) = zfice(i)*zoliq(i)
-             zrfl(i) = zrfl(i)+ zqprecl(i) &
+!<LTP
+             IF (iflag_evap_prec == 4) THEN
+                zrflcld(i) = zrflcld(i)+zqprecl(i) &
+                 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+                ziflcld(i) = ziflcld(i)+ zqpreci(i) &
+                      *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+                znebprecipcld(i) = rneb(i,k)
+                zrfl(i) = zrflcld(i) + zrflclr(i)
+                zifl(i) = ziflcld(i) + ziflclr(i)        
+!>LTP
+             ELSE
+                zrfl(i) = zrfl(i)+ zqprecl(i) &
                  *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
-             zifl(i) = zifl(i)+ zqpreci(i) &
+                zifl(i) = zifl(i)+ zqpreci(i) &
                       *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)  
+             
+             ENDIF !iflag_evap_prec==4
+
            ENDIF                      
          ENDDO
@@ -1314,19 +1487,44 @@
            d_qi(i,k) = zfice(i)*zoliq(i)
 !           endif
+!<LTP
+             IF (iflag_evap_prec == 4) THEN
+                zrflcld(i) = zrflcld(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
+                       *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
+                ziflcld(i) = ziflcld(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
+                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)  
+                znebprecipcld(i) = rneb(i,k)
+                zrfl(i) = zrflcld(i) + zrflclr(i)
+                zifl(i) = ziflcld(i) + ziflclr(i)        
+!>LTP
+             ELSE
 !AJ<
-           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
-               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
-           zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
-                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)  
+                   zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
+                       *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
+                        zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
+                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)  
      !      zrfl(i) = zrfl(i)+  zpluie                         &
      !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
      !      zifl(i) = zifl(i)+  zice                    &
      !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                    
+             ENDIF !iflag_evap_prec == 4             
 
 !CR : on prend en compte l'effet Bergeron dans les flux de precipitation
            IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
-              zsolid = zrfl(i)
-              zifl(i) = zifl(i)+zrfl(i)
-              zrfl(i) = 0.
+!<LTP
+                IF (iflag_evap_prec == 4) THEN
+                     zsolid = zrfl(i)
+                     ziflclr(i) = ziflclr(i) +zrflclr(i)
+                     ziflcld(i) = ziflcld(i) +zrflcld(i)
+                     zifl(i) = ziflclr(i)+ziflcld(i)
+                     zrflcld(i)=0.
+                     zrflclr(i)=0.   
+                     zrfl(i) = zrflclr(i)+zrflcld(i)
+!>LTP
+                ELSE
+                     zsolid = zrfl(i)
+                     zifl(i) = zifl(i)+zrfl(i)
+                     zrfl(i) = 0.
+                 ENDIF!iflag_evap_prec==4
+
            if (fl_cor_ebil .GT. 0) then
               zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
@@ -1358,4 +1556,30 @@
 !       ENDDO
 !     ENDIF
+
+
+!<LTP
+
+!Limitation de la fraction surfacique couverte par les précipitations lorsque l'intensité locale du flux de précipitation descend en
+!dessous de rain_int_min
+    IF (iflag_evap_prec==4) THEN
+        DO i=1, klon
+            IF (zrflclr(i) + ziflclr(i) .GT. 0 ) THEN
+                znebprecipclr(i) = min(znebprecipclr(i), max(zrflclr(i)/(znebprecipclr(i)*rain_int_min), ziflclr(i)/(znebprecipclr(i)*rain_int_min)))
+            ELSE
+                znebprecipclr(i)=0.
+            ENDIF
+
+            IF (zrflcld(i) + ziflcld(i) .GT. 0 ) THEN
+                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/(znebprecipcld(i)*rain_int_min), ziflcld(i)/(znebprecipcld(i)*rain_int_min)))
+            ELSE
+                znebprecipcld(i)=0.
+            ENDIF
+       ENDDO
+    ENDIf
+
+!>LTP
+
+
+
 
        
