source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_lscp_old.F90 @ 5087

Last change on this file since 5087 was 5087, checked in by abarral, 4 months ago

remove fixed-form \s+& remaining in .f90,.F90

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 64.2 KB
Line 
1! $Id: lmdz_lscp_old.F90 5087 2024-07-20 10:00:23Z abarral $
2!
3!
4MODULE lmdz_lscp_old
5CONTAINS
6SUBROUTINE fisrtilp(klon,klev,dtime,paprs,pplay,t,q,ptconv,ratqs, &
7     d_t, d_q, d_ql, d_qi, rneb,rneblsvol,radliq, rain, snow,          &
8     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
9     frac_impa, frac_nucl, beta,                        &
10     prfl, psfl, rhcl, zqta, fraca,                     &
11     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
12     iflag_ice_thermo,                                  &
13     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
14
15
16  !
17  USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
18  USE lmdz_cloudth, only : cloudth, cloudth_v3, cloudth_v6
19
20  USE lmdz_lscp_ini, ONLY: prt_level, lunout
21  USE lmdz_lscp_ini, ONLY : fl_cor_ebil
22  USE lmdz_lscp_ini, ONLY: iflag_t_glace,t_glace_min, t_glace_max, exposant_glace
23  USE lmdz_lscp_ini, ONLY : seuil_neb, rain_int_min, iflag_evap_prec, iflag_oldbug_fisrtilp,a_tr_sca
24  USE lmdz_lscp_ini, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol
25  USE lmdz_lscp_ini, ONLY: coef_eva, ffallv_lsc, ffallv_con
26  USE lmdz_lscp_ini, ONLY: cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
27  USE lmdz_lscp_ini, ONLY: reevap_ice, iflag_bergeron, iflag_fisrtilp_qsat, iflag_pdf
28
29
30
31  IMPLICIT none
32  !======================================================================
33  ! Auteur(s): Z.X. Li (LMD/CNRS)
34  ! Date: le 20 mars 1995
35  ! Objet: condensation et precipitation stratiforme.
36  !        schema de nuage
37  ! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
38  !             et ilp (il pleut, L. Li)
39  ! Principales parties:
40  ! P0> Thermalisation des precipitations venant de la couche du dessus
41  ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
42  ! P2> Formation du nuage (en k)
43  ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
44  ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
45  ! les valeurs de T et Q initiales
46  ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
47  ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
48  ! P2.B> Nuage "tout ou rien"
49  ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
50  ! P3> Formation de la precipitation (en k)
51  !======================================================================
52  ! JLD:
53  ! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
54  ! * fl_cor_ebil doit etre > 0 ;
55  !   fl_cor_ebil= 0 pour reproduire anciens bugs
56  !======================================================================
57  include "YOMCST.h"
58  include "YOETHF.h"
59  include "FCTTRE.h"
60  !
61  ! Principaux inputs:
62  !
63  REAL, INTENT(IN)                              :: dtime  ! intervalle du temps (s)
64  INTEGER, INTENT(IN)                           :: klon, klev
65  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs  ! pression a inter-couche
66  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay  ! pression au milieu de couche
67  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t      ! temperature (K)
68  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q      ! humidite specifique (kg/kg)
69  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv ! points ou le schema de conv. prof. est actif
70  INTEGER,                         INTENT(IN)   :: iflag_cld_th
71  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo
72  !
73  ! Inputs lies aux thermiques
74  !
75  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv
76  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta, fraca
77  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk, ztla
78  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: zthl
79  !
80  !  Input/output
81  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs  ! determine la largeur de distribution de vapeur
82  !
83  ! Principaux outputs:
84  !
85  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t  ! incrementation de la temperature (K)
86  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q  ! incrementation de la vapeur d'eau
87  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql ! incrementation de l'eau liquide
88  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi ! incrementation de l'eau glace
89  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb, rneblsvol ! fraction nuageuse
90  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta ! taux de conversion de l'eau cond
91  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radliq ! eau liquide utilisee dans rayonnements
92  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl ! humidite relative en ciel clair
93  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain
94  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow
95  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl
96  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl
97
98  !AA
99  ! Coeffients de fraction lessivee : pour OFF-LINE
100  !
101  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_nucl
102  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_1nucl
103  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_impa
104  !
105  ! Fraction d'aerosols lessivee par impaction et par nucleation
106  ! POur ON-LINE
107  !
108  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa
109  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl
110  REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth,cloudth_senv
111  REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath,cloudth_sigmaenv
112  !AA
113  ! --------------------------------------------------------------------------------
114  !
115  ! Options du programme:
116  !
117
118  REAL :: smallestreal
119
120  INTEGER, PARAMETER :: ninter=5 ! sous-intervals pour la precipitation
121  LOGICAL, PARAMETER :: cpartiel=.TRUE. ! condensation partielle
122  REAL, PARAMETER :: t_coup=234.0
123  REAL, PARAMETER :: DDT0=.01
124  REAL, PARAMETER :: ztfondue=278.15
125  ! --------------------------------------------------------------------------------
126  !
127  ! Variables locales:
128  !
129  INTEGER :: i, k, n, kk
130  REAL :: qsl, qsi
131  REAL :: zct      ,zcl
132  INTEGER :: ncoreczq 
133  REAL, DIMENSION(klon,klev) :: ctot,ctot_vol
134  REAL, DIMENSION(klon) :: zqs, zdqs, zdqsdT_raw, Tbef,qlbef,DT
135  REAL :: zdelta, zcor, zcvm5 
136  REAL ::num,denom
137
138  LOGICAL, DIMENSION(klon) :: lognormale,convergence
139  LOGICAL :: ice_thermo
140  INTEGER, DIMENSION(klon) :: n_i
141  INTEGER :: iter
142  REAL :: cste
143
144  REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta, Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2, qcloud
145  REAL :: erf   
146 
147  REAL :: zqev, zqevt, zqev0,zqevi, zqevti, zdelq
148  REAL, DIMENSION(klon) :: zrfl(klon), zrfln(klon), zrflclr(klon), zrflcld(klon), d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon), d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon)
149
150  REAL, DIMENSION(klon) :: zifl, zifln, ziflclr, ziflcld, zoliq, zcond, zq, zqn, zoliqp, zoliqi, zt
151! JBM (3/14) nexpo is replaced by exposant_glace
152! REAL nexpo ! exponentiel pour glace/eau
153! INTEGER, PARAMETER :: nexpo=6
154  INTEGER :: exposant_glace_old
155  REAL :: t_glace_min_old, ztot
156  REAL, DIMENSION(klon) ::  zdz,zrho , zrhol, zfice,zneb,znebprecip
157  REAL :: zchau      ,zfroi     
158  REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld, tot_zneb, tot_znebn, d_tot_zneb, d_znebprecip_clr_cld, d_znebprecip_cld_clr, dzfice
159  REAL :: zmelt, zpluie, zice
160  REAL :: zsolid
161!!!!
162!  Variables pour Bergeron
163  REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
164  REAL, DIMENSION(klon) :: zqpreci, zqprecl
165! Variable pour conservation enegie des precipitations
166  REAL, DIMENSION(klon) :: zmqc
167  !
168  LOGICAL, SAVE :: appel1er=.TRUE.
169  !$OMP THREADPRIVATE(appel1er)
170  !
171! iflag_oldbug_fisrtilp=0 enleve le BUG par JYG : tglace_min -> tglace_max
172! iflag_oldbug_fisrtilp=1 ajoute le BUG
173  !---------------------------------------------------------------
174  !
175  ! Fonctions en ligne:
176  !
177  REAL ::  fallvs,fallvc, zzz ! Vitesse de chute pour cristaux de glace
178                     ! (Heymsfield & Donner, 1990)
179  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
180  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
181  !
182  !---------------------------------------------------------------
183  !AA Variables traceurs:
184  !AA  Provisoire !!! Parametres alpha du lessivage
185  !AA  A priori on a 4 scavenging # possibles
186  !
187  ! Variables intermediaires
188  !
189  REAL :: zalpha_tr, zfrac_lessi
190  REAL, DIMENSION(klon) :: zprec_cond
191  !AA
192! RomP >>> 15 nov 2012
193! RomP <<<
194  REAL, DIMENSION(klon) :: zmair
195  REAL :: zcpair, zcpeau
196  !     Pour la conversion eau-neige
197  REAL, DIMENSION(klon) :: zlh_solid
198  REAL :: zm_solid
199
200
201  !
202  !ym
203!CR: pour iflag_ice_thermo=2, on active que la convection
204!  ice_thermo = iflag_ice_thermo .GE. 1
205
206 
207  znebprecip(:)=0.
208
209!<LTP
210  smallestreal=1.e-9
211  znebprecipclr(:)=0.
212  znebprecipcld(:)=0.
213!>LTP
214
215  ice_thermo = (iflag_ice_thermo == 1).OR.(iflag_ice_thermo >= 3)
216  zdelq=0.0
217  ctot_vol(1:klon,1:klev)=0.0
218  rneblsvol(1:klon,1:klev)=0.0
219
220  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
221  IF (appel1er) THEN
222     WRITE(lunout,*) 'fisrtilp, ninter:', ninter
223     WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
224     WRITE(lunout,*) 'FISRTILP VERSION LUDO'
225     
226     IF (ABS(dtime/REAL(ninter)-360.0)>0.001) THEN
227        WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
228        WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
229        !         CALL abort
230     ENDIF
231     appel1er = .FALSE.
232     !
233     !cdir collapse
234     DO k = 1, klev
235        DO i = 1, klon
236           pfrac_nucl(i,k)=1.
237           pfrac_1nucl(i,k)=1.
238           pfrac_impa(i,k)=1.
239           beta(i,k)=0.  !RomP initialisation
240        ENDDO
241     ENDDO
242
243  ENDIF          !  test sur appel1er
244  !
245  !MAf Initialisation a 0 de zoliq
246  !      DO i = 1, klon
247  !         zoliq(i)=0.
248  !      ENDDO
249  ! Determiner les nuages froids par leur temperature
250  !  nexpo regle la raideur de la transition eau liquide / eau glace.
251  !
252!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
253  IF (iflag_t_glace==0) THEN
254!   ztglace = RTT - 15.0
255    t_glace_min_old = RTT - 15.0
256    !AJ<
257    IF (ice_thermo) THEN
258!     nexpo = 2
259      exposant_glace_old = 2
260    ELSE
261!     nexpo = 6
262      exposant_glace_old = 6
263    ENDIF
264   
265  ENDIF
266 
267!!  RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
268!!  RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
269!>AJ
270  !cc      nexpo = 1
271  !
272  ! Initialiser les sorties:
273  !
274  !cdir collapse
275  DO k = 1, klev+1
276     DO i = 1, klon
277        prfl(i,k) = 0.0
278        psfl(i,k) = 0.0
279     ENDDO
280  ENDDO
281
282  !cdir collapse
283
284  DO k = 1, klev
285     DO i = 1, klon
286        d_t(i,k) = 0.0
287        d_q(i,k) = 0.0
288        d_ql(i,k) = 0.0
289        d_qi(i,k) = 0.0
290        rneb(i,k) = 0.0
291        radliq(i,k) = 0.0
292        frac_nucl(i,k) = 1.
293        frac_impa(i,k) = 1.
294     ENDDO
295  ENDDO
296  DO i = 1, klon
297     rain(i) = 0.0
298     snow(i) = 0.0
299     zoliq(i)=0.
300     !     ENDDO
301     !
302     ! Initialiser le flux de precipitation a zero
303     !
304     !     DO i = 1, klon
305     zrfl(i) = 0.0
306     zifl(i) = 0.0
307!<LTP
308     zrflclr(i) = 0.0
309     ziflclr(i) = 0.0
310     zrflcld(i) = 0.0
311     ziflcld(i) = 0.0
312     tot_zneb(i) = 0.0
313     tot_znebn(i) = 0.0
314     d_tot_zneb(i) = 0.0
315!>LTP
316
317     zneb(i) = seuil_neb
318  ENDDO
319  !
320  !
321  !AA Pour plus de securite
322
323  zalpha_tr   = 0.
324  zfrac_lessi = 0.
325
326  !AA==================================================================
327  !
328  ncoreczq=0
329  ! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
330  !
331  DO k = klev, 1, -1
332     !
333     !AA===============================================================
334     !
335     ! Initialisation temperature et vapeur
336     DO i = 1, klon
337        zt(i)=t(i,k)
338        zq(i)=q(i,k)
339     ENDDO
340     !
341     ! ----------------------------------------------------------------
342     ! P0> Thermalisation des precipitations venant de la couche du dessus
343     ! ----------------------------------------------------------------
344     ! Calculer la varition de temp. de l'air du a la chaleur sensible
345     ! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
346     ! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
347     ! des differentes transformations thermodynamiques. Cette masse d'eau doit
348     ! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
349     ! de l'enthalpie  de la couche avec la temperature
350     ! Variables calculees ou modifiees:
351     !   -  zt: temperature de la cocuhe
352     !   - zmqc: masse de precip qui doit etre thermalisee
353     !
354     IF(k<=klev-1) THEN
355        DO i = 1, klon
356           !IM
357           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
358           ! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
359           zcpair=RCPD*(1.0+RVTMP2*zq(i))
360           zcpeau=RCPD*RVTMP2
361         if (fl_cor_ebil > 0) then
362           ! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm
363           ! pour s'assurer que la precip arrivant au sol aura bien la temperature de la
364           ! derniere couche
365           zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
366           ! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
367           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
368                 / (zcpair + zmqc(i)*zcpeau)
369         else ! si on maintient les anciennes erreurs
370           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
371                + zmair(i)*zcpair*zt(i) ) &
372                / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
373         end if
374        ENDDO
375     ELSE  ! IF(k.LE.klev-1)
376        DO i = 1, klon
377           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
378           zmqc(i) = 0.
379        ENDDO
380     ENDIF ! end IF(k.LE.klev-1)
381     !
382     ! ----------------------------------------------------------------
383     ! P1> Calcul de l'evaporation de la precipitation
384     ! ----------------------------------------------------------------
385     ! On evapore une partie des precipitations venant de la maille du dessus.
386     ! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
387     ! ce que la fraction de cette couche qui est sous le nuage soit saturee.
388     ! Variables calculees ou modifiees:
389     !   - zrfl et zifl: flux de precip liquide et glace
390     !   - zq, zt: humidite et temperature de la cocuhe
391     !   - zmqc: masse de precip qui doit etre thermalisee
392     !
393     IF (iflag_evap_prec>=1) THEN
394        DO i = 1, klon
395!          S'il y a des precipitations
396           IF (zrfl(i)+zifl(i)>0.) THEN
397              ! Calcul du qsat
398              IF (thermcep) THEN
399                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
400                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
401                 zqs(i)=MIN(0.5,zqs(i))
402                 zcor=1./(1.-RETV*zqs(i))
403                 zqs(i)=zqs(i)*zcor
404              ELSE
405                 IF (zt(i) < t_coup) THEN
406                    zqs(i) = qsats(zt(i)) / pplay(i,k)
407                 ELSE
408                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
409                 ENDIF
410              ENDIF
411           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
412        ENDDO
413!AJ<
414
415       IF (.NOT. ice_thermo) THEN
416        DO i = 1, klon
417!          S'il y a des precipitations
418           IF (zrfl(i)+zifl(i)>0.) THEN
419                ! Evap max pour ne pas saturer la fraction sous le nuage
420                ! Evap max jusqu'à atteindre la saturation dans la partie
421                ! de la maille qui est sous le nuage de la couche du dessus
422                !!! On ne tient compte de cette fraction que sous une seule
423                !!! couche sous le nuage
424                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
425             ! Ajout de la prise en compte des precip a thermiser
426             ! avec petite reecriture
427             if  (fl_cor_ebil > 0) then ! nouveau
428                ! Calcul de l'evaporation du flux de precip herite
429                !   d'au-dessus
430                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
431                     * zmair(i)/pplay(i,k)*zt(i)*RD
432                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
433
434                ! Seuil pour ne pas saturer la fraction sous le nuage
435                zqev = MIN (zqev, zqevt)
436                ! Nouveau flux de precip
437                zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
438                ! Aucun flux liquide pour T < t_coup, on reevapore tout.
439                IF (zt(i) < t_coup.and.reevap_ice) THEN
440                  zrfln(i)=0.
441                  zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
442                END IF
443                ! Nouvelle vapeur
444                zq(i) = zq(i) + zqev
445                zmqc(i) = zmqc(i)-zqev
446                ! Nouvelle temperature (chaleur latente)
447                zt(i) = zt(i) - zqev &
448                     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
449!!JLD debut de partie a supprimer a terme
450            else ! if  (fl_cor_ebil .GT. 0)
451                ! Calcul de l'evaporation du flux de precip herite
452                !   d'au-dessus
453                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
454                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
455                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
456                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
457                ! Seuil pour ne pas saturer la fraction sous le nuage
458                zqev = MIN (zqev, zqevt)
459                ! Nouveau flux de precip
460                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
461                     /RG/dtime
462                ! Aucun flux liquide pour T < t_coup
463                IF (zt(i) < t_coup.and.reevap_ice) zrfln(i)=0.
464                ! Nouvelle vapeur
465                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
466                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
467                ! Nouvelle temperature (chaleur latente)
468                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
469                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
470                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
471              end if ! if  (fl_cor_ebil .GT. 0)
472!!JLD fin de partie a supprimer a terme
473                zrfl(i) = zrfln(i)
474                zifl(i) = 0.
475           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
476        ENDDO
477!
478       ELSE ! (.NOT. ice_thermo)
479!      ================================
480!      Avec thermodynamique de la glace
481!      ================================
482        DO i = 1, klon
483
484
485!AJ<
486!        S'il y a des precipitations
487         IF (zrfl(i)+zifl(i)>0.) THEN
488
489        !LTP<
490        !On ne tient compte que du flux de précipitation en ciel clair dans le calcul de l'évaporation.
491                IF (iflag_evap_prec==4) THEN
492                        zrfl(i) = zrflclr(i)
493                        zifl(i) = ziflclr(i)
494                ENDIF
495       
496        !>LTP
497
498         IF (iflag_evap_prec==1) THEN
499            znebprecip(i)=zneb(i)
500         ELSE
501            znebprecip(i)=MAX(zneb(i),znebprecip(i))
502         ENDIF
503         
504         IF (iflag_evap_prec==4) THEN
505        ! Evap max pour ne pas saturer toute la maille
506         zqev0 = MAX (0.0, zqs(i)-zq(i))
507         ELSE
508        ! Evap max pour ne pas saturer la fraction sous le nuage
509         zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) )
510         ENDIF
511
512         !JAM
513         ! On differencie qsat pour l'eau et la glace
514         ! Si zdelta=1. --> glace
515         ! Si zdelta=0. --> eau liquide
516       
517         ! Calcul du qsat par rapport a l'eau liquide
518         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
519         qsl= MIN(0.5,qsl)
520         zcor= 1./(1.-RETV*qsl)
521         qsl= qsl*zcor
522         
523         ! Calcul de l'evaporation du flux de precip venant du dessus
524         ! Formulation en racine du flux de precip
525         ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
526         IF (iflag_evap_prec==3) THEN
527         zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) &
528              *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
529              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
530!<LTP
531         ELSE IF (iflag_evap_prec==4) THEN
532         zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl) &
533              *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
534              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
535!>LTP
536        ELSE
537         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
538              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
539         ENDIF
540
541
542         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
543              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
544         
545         ! Calcul du qsat par rapport a la glace
546         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
547         qsi= MIN(0.5,qsi)
548         zcor= 1./(1.-RETV*qsi)
549         qsi= qsi*zcor
550
551         ! Calcul de la sublimation du flux de precip solide herite
552         !   d'au-dessus
553         IF (iflag_evap_prec==3) THEN
554         zqevti = znebprecip(i)*coef_eva*(1.0-zq(i)/qsi) &
555              *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
556              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
557!<LTP
558         ELSE IF (iflag_evap_prec==4) THEN
559         zqevti = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsi) &
560              *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
561              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
562!>LTP
563         ELSE
564         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
565              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
566         ENDIF
567         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
568              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
569
570       
571        !JAM
572        ! Limitation de l'evaporation. On s'assure qu'on ne sature pas
573        ! la fraction de la couche sous le nuage sinon on repartit zqev0
574        ! en conservant la proportion liquide / glace
575     
576         IF (zqevt+zqevti>zqev0) THEN
577            zqev=zqev0*zqevt/(zqevt+zqevti)
578            zqevi=zqev0*zqevti/(zqevt+zqevti)
579         ELSE
580!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
581!       liquides et solides meme si on ne sature pas la couche.
582!       A mon avis, le test est inutile, et il faudrait tout remplacer par:
583!            zqev=zqevt
584!            zqevi=zqevti
585             IF (zqevt+zqevti>0.) THEN
586                zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
587                zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
588             ELSE
589             zqev=0.
590             zqevi=0.
591             ENDIF
592         ENDIF
593
594         ! Nouveaux flux de precip liquide et solide
595         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
596                                 /RG/dtime)
597         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
598                                 /RG/dtime)
599
600         ! Mise a jour de la vapeur, temperature et flux de precip
601         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
602                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
603       if (fl_cor_ebil > 0) then ! avec correction thermalisation des precips
604         zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
605                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
606         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
607                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
608                  * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
609                  + (zifln(i)-zifl(i)) &
610                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
611                  * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
612       else ! sans correction thermalisation des precips
613         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
614                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
615                  * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
616                  + (zifln(i)-zifl(i)) &
617                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
618                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
619       end if
620        ! Nouvelles vaeleurs des precips liquides et solides
621         zrfl(i) = zrfln(i)
622         zifl(i) = zifln(i)
623
624!<LTP
625        IF (iflag_evap_prec==4) THEN
626                zrflclr(i) = zrfl(i)
627                ziflclr(i) = zifl(i)   
628                IF(zrflclr(i) + ziflclr(i) <= 0) THEN
629                        znebprecipclr(i) = 0.
630                ENDIF   
631                zrfl(i) = zrflclr(i) + zrflcld(i)
632                zifl(i) = ziflclr(i) + ziflcld(i)
633        ENDIF
634!>LTP       
635
636
637!CR ATTENTION: deplacement de la fonte de la glace
638!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
639!!!        zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2  !!!!!!!!! jyg
640!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
641           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
642           ! fraction de la precip solide qui est fondue
643           zmelt = MIN(MAX(zmelt,0.),1.)
644           ! Fusion de la glace
645!<LTP
646           IF (iflag_evap_prec==4) THEN
647                   zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
648                   zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
649                   zrfl(i)=zrflclr(i)+zrflcld(i)
650!>LTP       
651           ELSE
652                   zrfl(i)=zrfl(i)+zmelt*zifl(i)
653           ENDIF
654           if (fl_cor_ebil <= 0) then
655             ! the following line should not be here. Indeed, if zifl is modified
656             ! now, zifl(i)*zmelt is no more the amount of ice that has melt
657             ! and therefore the change in temperature computed below is wrong
658             zifl(i)=zifl(i)*(1.-zmelt)
659           end if
660           ! Chaleur latente de fusion
661        if (fl_cor_ebil > 0) then ! avec correction thermalisation des precips
662           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
663                      *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
664        else ! sans correction thermalisation des precips
665           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
666                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
667        end if
668           if (fl_cor_ebil > 0) then ! correction bug, deplacement ligne precedente
669!<LTP
670             IF (iflag_evap_prec==4) THEN
671                   ziflclr(i)=ziflclr(i)*(1.-zmelt)
672                   ziflcld(i)=ziflcld(i)*(1.-zmelt)
673                   zifl(i)=ziflclr(i)+ziflcld(i)
674!>LTP       
675             ELSE
676                   zifl(i)=zifl(i)*(1.-zmelt)
677             ENDIF
678           end if
679
680           ELSE
681              ! Si on n'a plus de pluies, on reinitialise a 0 la farcion
682              ! sous nuageuse utilisee pour la pluie.
683              znebprecip(i)=0.
684           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
685        ENDDO
686
687      ENDIF ! (.NOT. ice_thermo)
688     
689     ! ----------------------------------------------------------------
690     ! Fin evaporation de la precipitation
691     ! ----------------------------------------------------------------
692     ENDIF ! (iflag_evap_prec>=1)
693     !
694     ! Calculer Qs et L/Cp*dQs/dT:
695     !
696     IF (thermcep) THEN
697        DO i = 1, klon
698           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
699           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
700       if  (fl_cor_ebil > 0) then ! nouveau
701           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
702       else   
703           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
704       endif
705           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
706           zqs(i) = MIN(0.5,zqs(i))
707           zcor = 1./(1.-RETV*zqs(i))
708           zqs(i) = zqs(i)*zcor
709           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
710           zdqsdT_raw(i) = zdqs(i)*  &
711           RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
712        ENDDO
713     ELSE
714        DO i = 1, klon
715           IF (zt(i)<t_coup) THEN
716              zqs(i) = qsats(zt(i))/pplay(i,k)
717              zdqs(i) = dqsats(zt(i),zqs(i))
718           ELSE
719              zqs(i) = qsatl(zt(i))/pplay(i,k)
720              zdqs(i) = dqsatl(zt(i),zqs(i))
721           ENDIF
722        ENDDO
723     ENDIF
724     !
725     ! Determiner la condensation partielle et calculer la quantite
726     ! de l'eau condensee:
727     !
728!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
729!       if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
730!         write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
731!        " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
732!         stop
733!       endif
734
735     ! ----------------------------------------------------------------
736     ! P2> Formation du nuage
737     ! ----------------------------------------------------------------
738     ! Variables calculees:
739     !   rneb  : fraction nuageuse
740     !   zcond : eau condensee moyenne dans la maille.
741     !   rhcl: humidite relative ciel-clair
742     !   zt : temperature de la maille
743     ! ----------------------------------------------------------------
744     !
745     IF (cpartiel) THEN
746        ! -------------------------
747        ! P2.A> Nuage fractionnaire
748        ! -------------------------
749        !
750        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
751        !   nuageuse a partir des PDF de Sandrine Bony.
752        !   rneb  : fraction nuageuse
753        !   zqn   : eau totale dans le nuage
754        !   zcond : eau condensee moyenne dans la maille.
755        !  on prend en compte le réchauffement qui diminue la partie
756        ! condensee
757        !
758        !   Version avec les raqts
759
760        ! ----------------------------------------------------------------
761        ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
762        ! ----------------------------------------------------------------
763        if (iflag_pdf==0) then
764
765           ! version creneau de (Li, 1998)
766           do i=1,klon
767              zdelq = min(ratqs(i,k),0.99) * zq(i)
768              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
769              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
770           enddo
771
772        else !  if (iflag_pdf.eq.0)
773           ! ----------------------------------------------------------------
774           ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
775           ! les valeurs de T et Q initiales
776           ! ----------------------------------------------------------------
777           do i=1,klon
778              if(zq(i)<1.e-15) then
779                 ncoreczq=ncoreczq+1
780                 zq(i)=1.e-15
781              endif
782           enddo
783
784           if (iflag_cld_th>=5) then
785
786              if (iflag_cloudth_vert<=2) then
787               call cloudth(klon,klev,k,ztv, &
788                   zq,zqta,fraca, &
789                   qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
790                   ratqs,zqs,t, &
791                   cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
792
793              elseif (iflag_cloudth_vert>=3 .and. iflag_cloudth_vert<=5) then
794               call cloudth_v3(klon,klev,k,ztv, &
795                   zq,zqta,fraca, &
796                   qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
797                   ratqs,zqs,t, &
798                   cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
799              !----------------------------------
800              !Version these Jean Jouhaud, Decembre 2018
801              !----------------------------------             
802              elseif (iflag_cloudth_vert==6) then
803               call cloudth_v6(klon,klev,k,ztv, &
804                   zq,zqta,fraca, &
805                   qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
806                   ratqs,zqs,t, &
807                   cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
808
809              endif
810              do i=1,klon
811                 rneb(i,k)=ctot(i,k)
812                 rneblsvol(i,k)=ctot_vol(i,k)
813                 zqn(i)=qcloud(i)
814              enddo
815
816           endif
817
818           if (iflag_cld_th <= 4) then
819              lognormale = .true.
820           elseif (iflag_cld_th >= 6) then
821              ! lognormale en l'absence des thermiques
822              lognormale = fraca(:,k) < 1e-10
823           else
824              ! Dans le cas iflag_cld_th=5, on prend systématiquement la
825              ! bi-gaussienne
826              lognormale = .false.
827           end if
828
829!CR: variation de qsat avec T en presence de glace ou non
830!initialisations
831           do i=1,klon
832              DT(i) = 0.
833              n_i(i)=0
834              Tbef(i)=zt(i)
835              qlbef(i)=0.
836           enddo
837
838        ! ----------------------------------------------------------------
839        ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
840        ! Calcul des grandeurs nuageuses en tenant compte de l'effet de
841        ! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
842        ! qui en dependent. Ce changement de temperature est provisoire, et
843        ! la valeur definitive sera calcule plus tard.
844        ! Variables calculees:
845        !   rneb : nebulosite
846        !   zcond: eau condensee en moyenne dans la maille
847        ! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
848        ! pas clair, il n'y a probablement pas de prise en compte de l'effet de
849        ! T sur qsat
850        ! ----------------------------------------------------------------
851
852!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
853           if (iflag_fisrtilp_qsat>=0) then
854             ! Iteration pour condensation avec variation de qsat(T)
855             ! -----------------------------------------------------
856             do iter=1,iflag_fisrtilp_qsat+1
857               
858               do i=1,klon
859!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
860                 ! !! convergence = .true. tant que l'on n'a pas converge !!
861                 ! ------------------------------
862                 convergence(i)=abs(DT(i))>DDT0
863                 if ((convergence(i).or.(n_i(i)==0)).and.lognormale(i)) then
864                 ! si on n'a pas converge
865                 !
866                 ! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
867                 ! ---------------------------------------------------------------
868                 ! Variables calculees:
869                 ! rneb : nebulosite
870                 ! zqn : eau condensee, dans le nuage (in cloud water content)
871                 ! zcond: eau condensee en moyenne dans la maille
872                 ! rhcl: humidite relative ciel-clair
873                 !
874                 Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
875                 if (.not.ice_thermo) then   
876                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
877                 else
878                    if (iflag_t_glace==0) then
879                    zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
880                    else if (iflag_t_glace>=1) then
881                       if (iflag_oldbug_fisrtilp==0) then
882                          zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
883                       else
884                          !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
885                          zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
886                       endif
887                    endif
888                 endif
889                 ! Calcul de rneb, qzn et zcond pour les PDF lognormales
890                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
891               if (fl_cor_ebil > 0) then
892                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
893               else
894                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
895               end if
896                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
897                 zqs(i) = MIN(0.5,zqs(i))
898                 zcor = 1./(1.-RETV*zqs(i))
899                 zqs(i) = zqs(i)*zcor
900                 zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
901                 zpdf_sig(i)=ratqs(i,k)*zq(i)
902                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
903                 zpdf_delta(i)=log(zq(i)/zqs(i))
904                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
905                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
906                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
907                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
908                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
909                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
910                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
911                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
912             
913                 if (zpdf_e1(i)<1.e-10) then
914                    rneb(i,k)=0.
915                    zqn(i)=zqs(i)
916                 else
917                    rneb(i,k)=0.5*zpdf_e1(i)
918                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
919                 endif
920
921                 ! If vertical heterogeneity, change fraction by volume as well
922                 if (iflag_cloudth_vert>=3) then
923                   ctot_vol(i,k)=rneb(i,k)
924                   rneblsvol(i,k)=ctot_vol(i,k)
925                 endif
926
927                 endif !convergence
928
929               enddo ! boucle en i
930
931                 ! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT
932                 !         due a la condensation.
933                 ! ---------------------------------------------------------------
934                 ! Variables calculees:
935                 ! DT : variation de temperature due a la condensation
936
937                 if (.not. ice_thermo) then
938                 ! --------------------------
939                 do i=1,klon
940                 if ((convergence(i).or.(n_i(i)==0)).and.lognormale(i)) then
941
942                 qlbef(i)=max(0.,zqn(i)-zqs(i))
943               if (fl_cor_ebil > 0) then
944                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
945               else
946                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
947               end if
948                 denom=1.+rneb(i,k)*zdqs(i)
949                 DT(i)=num/denom
950                 n_i(i)=n_i(i)+1
951                 endif
952                 enddo
953
954                 else ! if (.not. ice_thermo)
955                 ! --------------------------
956                 if (iflag_t_glace>=1) then
957                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
958                 endif
959
960                 do i=1,klon
961                 if ((convergence(i).or.(n_i(i)==0)).and.lognormale(i)) then
962                 
963                 if (iflag_t_glace==0) then
964                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
965                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
966                    zfice(i) = zfice(i)**exposant_glace_old
967                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
968                       / (t_glace_min_old - RTT)
969                 endif
970                 
971                 if (iflag_t_glace>=1.and.zfice(i)>0.) then
972                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
973                      / (t_glace_min - t_glace_max)
974                 endif
975               
976                 if ((zfice(i)==0).or.(zfice(i)==1)) then
977                    dzfice(i)=0.
978                 endif
979
980                 if (zfice(i)<1) then
981                    cste=RLVTT
982                 else
983                    cste=RLSTT
984                 endif
985
986                 qlbef(i)=max(0.,zqn(i)-zqs(i))
987               if (fl_cor_ebil > 0) then
988                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
989            +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
990                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
991                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)    &
992                 *qlbef(i)*dzfice(i)
993               else
994                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
995           +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
996                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
997                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
998               end if
999                 DT(i)=num/denom
1000                 n_i(i)=n_i(i)+1
1001
1002                 endif ! fin convergence
1003                 enddo ! fin boucle i
1004
1005                 endif !ice_thermo
1006
1007             enddo ! iter=1,iflag_fisrtilp_qsat+1
1008             ! Fin d'iteration pour condensation avec variation de qsat(T)
1009             ! -----------------------------------------------------------
1010           endif !  if (iflag_fisrtilp_qsat.ge.0)
1011     ! ----------------------------------------------------------------
1012     ! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
1013     ! ----------------------------------------------------------------
1014
1015        endif ! iflag_pdf
1016
1017!        if (iflag_fisrtilp_qsat.eq.-1) then
1018!------------------------------------------
1019!CR: ATTENTION option fausse mais a existe:
1020! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
1021! activer les lignes suivantes:
1022       IF (1==0) THEN
1023       DO i=1,klon
1024           IF (rneb(i,k) <= 0.0) THEN
1025              zqn(i) = 0.0
1026              rneb(i,k) = 0.0
1027              zcond(i) = 0.0
1028              rhcl(i,k)=zq(i)/zqs(i)
1029           ELSE IF (rneb(i,k) >= 1.0) THEN
1030              zqn(i) = zq(i)
1031              rneb(i,k) = 1.0                 
1032              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
1033              rhcl(i,k)=1.0
1034           ELSE
1035              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
1036              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1037           ENDIF
1038       ENDDO
1039       ENDIF
1040!------------------------------------------
1041
1042!        ELSE
1043        ! ----------------------------------------------------------------
1044        ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
1045        ! Variables calculees:
1046        !   rneb : nebulosite
1047        !   zcond: eau condensee en moyenne dans la maille
1048        !   zq : eau vapeur dans la maille
1049        !   zt : temperature de la maille
1050        !   rhcl: humidite relative ciel-clair
1051        ! ----------------------------------------------------------------
1052        !
1053        ! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb)
1054        ! Calcule de l'eau condensee moyenne dans la maille (zcond),
1055        ! et de l'humidite relative ciel-clair (rhcl)
1056        DO i=1,klon
1057           IF (rneb(i,k) <= 0.0) THEN
1058              zqn(i) = 0.0
1059              rneb(i,k) = 0.0
1060              zcond(i) = 0.0
1061              rhcl(i,k)=zq(i)/zqs(i)
1062           ELSE IF (rneb(i,k) >= 1.0) THEN
1063              zqn(i) = zq(i)
1064              rneb(i,k) = 1.0
1065              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1066              rhcl(i,k)=1.0
1067           ELSE
1068              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1069              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
1070           ENDIF
1071        ENDDO
1072
1073       
1074        ! If vertical heterogeneity, change fraction by volume as well
1075        if (iflag_cloudth_vert>=3) then
1076          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1077          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1078        endif
1079
1080!        ENDIF
1081
1082     ELSE ! de IF (cpartiel)
1083        ! -------------------------
1084        ! P2.B> Nuage "tout ou rien"
1085        ! -------------------------
1086        ! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
1087        DO i = 1, klon
1088           IF (zq(i)>zqs(i)) THEN
1089              rneb(i,k) = 1.0
1090           ELSE
1091              rneb(i,k) = 0.0
1092           ENDIF
1093           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
1094        ENDDO
1095     ENDIF ! de IF (cpartiel)
1096     !
1097     ! Mise a jour vapeur d'eau
1098     ! -------------------------
1099     DO i = 1, klon
1100        zq(i) = zq(i) - zcond(i)
1101        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
1102     ENDDO
1103!AJ<
1104     ! ------------------------------------
1105     ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
1106     ! -------------------------------------
1107     ! Variable calcule:
1108     !   zt : temperature de la maille
1109     !
1110     IF (.NOT. ice_thermo) THEN
1111        if (iflag_fisrtilp_qsat<1) then
1112           DO i = 1, klon
1113             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
1114           ENDDO
1115        else if (iflag_fisrtilp_qsat>0) then
1116           DO i= 1, klon
1117    if (fl_cor_ebil > 0) then
1118             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1119    else
1120             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1121    end if
1122           ENDDO
1123        endif
1124     ELSE
1125         if (iflag_t_glace>=1) then
1126            CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1127         endif
1128         if (iflag_fisrtilp_qsat<1) then
1129           DO i = 1, klon
1130! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1131!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1132!                                     t_glace_max, exposant_glace)
1133              if (iflag_t_glace==0) then
1134                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1135                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1136                    zfice(i) = zfice(i)**exposant_glace_old
1137              endif
1138              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
1139                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
1140           ENDDO
1141         else
1142           DO i=1, klon
1143! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1144!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1145!                                     t_glace_max, exposant_glace)
1146              if (iflag_t_glace==0) then
1147                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
1148                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1149                    zfice(i) = zfice(i)**exposant_glace_old
1150              endif
1151        if (fl_cor_ebil > 0) then
1152              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1153               * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1154                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1155        else
1156              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
1157                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
1158        end if
1159           ENDDO
1160         endif
1161!         print*,zt(i),zrfl(i),zifl(i),'temp1'
1162       ENDIF
1163!>AJ
1164
1165     ! ----------------------------------------------------------------
1166     ! P3> Formation des precipitations
1167     ! ----------------------------------------------------------------
1168     !
1169     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
1170     !
1171
1172!<LTP
1173
1174IF (iflag_evap_prec==4) THEN
1175        !Partitionnement des precipitations venant du dessus en précipitations nuageuses
1176        !et précipitations ciel clair
1177
1178        !0) Calculate tot_zneb, la fraction nuageuse totale au-dessus du nuage
1179        !en supposant un recouvrement maximum aléatoire (voir Jakob and Klein, 2000)
1180       
1181        DO i=1, klon
1182                tot_znebn(i) = 1 - (1-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
1183                        /(1-min(zneb(i),1-smallestreal))
1184                d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
1185                tot_zneb(i) = tot_znebn(i)
1186
1187
1188                !1) Cloudy to clear air
1189                d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
1190                IF (znebprecipcld(i) > 0) THEN
1191                        d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
1192                        d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
1193                ELSE
1194                        d_zrfl_cld_clr(i) = 0.
1195                        d_zifl_cld_clr(i) = 0.
1196                ENDIF
1197
1198                !2) Clear to cloudy air
1199                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) &
1200                        - d_tot_zneb(i) - zneb(i)))
1201                IF (znebprecipclr(i) > 0) THEN
1202                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
1203                        d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
1204                ELSE
1205                        d_zrfl_clr_cld(i) = 0.
1206                        d_zifl_clr_cld(i) = 0.
1207                ENDIF
1208
1209                !Update variables
1210                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i) 
1211                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
1212                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
1213                ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
1214                zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
1215                ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
1216
1217        ENDDO
1218ENDIF
1219
1220!>LTP
1221
1222
1223
1224     ! Initialisation de zoliq (eau condensee moyenne dans la maille)
1225     DO i = 1, klon
1226        IF (rneb(i,k)>0.0) THEN
1227           zoliq(i) = zcond(i)
1228           zrho(i) = pplay(i,k) / zt(i) / RD
1229           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
1230        ENDIF
1231     ENDDO
1232!AJ<
1233     IF (.NOT. ice_thermo) THEN
1234       IF (iflag_t_glace==0) THEN
1235         DO i = 1, klon
1236            IF (rneb(i,k)>0.0) THEN
1237               zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
1238               zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
1239               zfice(i) = zfice(i)**exposant_glace_old
1240!              zfice(i) = zfice(i)**nexpo
1241         !!      zfice(i)=0.
1242            ENDIF
1243         ENDDO
1244       ELSE ! of IF (iflag_t_glace.EQ.0)
1245         CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
1246!         DO i = 1, klon
1247!            IF (rneb(i,k).GT.0.0) THEN
1248! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
1249!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
1250!                                     t_glace_max, exposant_glace)
1251!            ENDIF
1252!         ENDDO
1253       ENDIF
1254     ENDIF
1255
1256     ! Calcul de radliq (eau condensee pour le rayonnement)
1257     ! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
1258     ! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
1259     ! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
1260     ! transmise au rayonnement;
1261     ! ----------------------------------------------------------------
1262     DO i = 1, klon
1263        IF (rneb(i,k)>0.0) THEN
1264           zneb(i) = MAX(rneb(i,k), seuil_neb)
1265     !      zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) 
1266!           print*,zt(i),'fractionglace'
1267!>AJ
1268           radliq(i,k) = zoliq(i)/REAL(ninter+1)
1269        ENDIF
1270     ENDDO
1271     !
1272     DO n = 1, ninter
1273        DO i = 1, klon
1274           IF (rneb(i,k)>0.0) THEN
1275              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
1276              ! Initialization of zpluie and zice:
1277              zpluie=0
1278              zice=0
1279              IF (zneb(i)==seuil_neb) THEN
1280                 ztot = 0.0
1281              ELSE
1282                 !  quantite d'eau a eliminer: zchau (Sundqvist, 1978)
1283                 !  meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
1284                 if (ptconv(i,k)) then
1285                    zcl   =cld_lc_con
1286                    zct   =1./cld_tau_con
1287                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1288                         *fallvc(zrhol(i)) * zfice(i)
1289                 else
1290                    zcl   =cld_lc_lsc
1291                    zct   =1./cld_tau_lsc
1292                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
1293                         *fallvs(zrhol(i)) * zfice(i)
1294                 endif
1295
1296                 ! si l'heterogeneite verticale est active, on utilise
1297                 ! la fraction volumique "vraie" plutot que la fraction
1298                 ! surfacique modifiee, qui est plus grande et reduit
1299                 ! sinon l'eau in-cloud de facon artificielle
1300                 if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then
1301                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
1302                         *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl   )**2)) *(1.-zfice(i))
1303                 else
1304                    zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
1305                         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
1306                 endif
1307!AJ<
1308                 IF (.NOT. ice_thermo) THEN
1309                   ztot    = zchau + zfroi
1310                 ELSE
1311                   zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
1312                   zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
1313                   ztot    = zpluie    + zice
1314                 ENDIF
1315!>AJ
1316                 ztot    = MAX(ztot   ,0.0)
1317              ENDIF
1318              ztot    = MIN(ztot,zoliq(i))
1319!AJ<
1320     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
1321     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
1322!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
1323!      temporaires et ne doivent pas etre calcule (alors qu'elles le sont
1324!      si iflag_bergeron <> 2
1325!      A SUPPRIMER A TERME
1326              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
1327              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
1328              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
1329!>AJ
1330              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
1331           ENDIF
1332        ENDDO  ! i = 1,klon
1333     ENDDO     ! n = 1,ninter
1334
1335     ! ----------------------------------------------------------------
1336     !
1337     IF (.NOT. ice_thermo) THEN
1338       DO i = 1, klon
1339         IF (rneb(i,k)>0.0) THEN
1340           d_ql(i,k) = zoliq(i)
1341           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
1342                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1343         ENDIF
1344       ENDDO
1345     ELSE
1346!
1347!CR&JYG<
1348! On prend en compte l'effet Bergeron dans les flux de precipitation :
1349! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
1350! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
1351! et les precipitations est grossierement pris en compte en linearisant les equations
1352! et en approximant le processus de precipitation liquide par un processus a seuil.
1353! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
1354! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
1355! Le condensat precipitant solide est augmente.
1356! La vapeur d'eau est augmentee.
1357!
1358       IF ((iflag_bergeron == 2)) THEN
1359         DO i = 1, klon
1360           IF (rneb(i,k) > 0.0) THEN
1361             zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1362             zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1363           if (fl_cor_ebil > 0) then
1364             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1365             coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1366!            Calcul de DT si toute les precips liquides congelent
1367             DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1368!            On ne veut pas que T devienne superieur a la temp. de congelation.
1369!            donc que Delta > RTT-zt(i
1370             DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1371             zt(i)      = zt(i)      + DeltaT
1372!            Eau vaporisee du fait de l'augmentation de T
1373             Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1374!            on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
1375             zq(i) = zq(i) +  Deltaq
1376!            Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
1377             zcond(i)   = max( zcond(i)- Deltaq, 0. )
1378!            precip liquide qui congele ou qui s'evapore
1379             Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1380             zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1381!            bilan eau glacee
1382             zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1383           else ! if (fl_cor_ebil .GT. 0)
1384!            ancien calcul
1385             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
1386             coef1 = RLMLT*zdqs(i)/RLVTT
1387             DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
1388             zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
1389             zqprecl(i) = max( zqprecl(i) - zcp/RLMLT*(1.+coef1)*DeltaT, 0. )
1390             zcond(i)   = max( zcond(i)   - zcp/RLVTT*zdqs(i)*DeltaT, 0. )
1391             zq(i)      = zq(i)      + zcp/RLVTT*zdqs(i)*DeltaT
1392             zt(i)      = zt(i)      + DeltaT
1393           end if ! if (fl_cor_ebil .GT. 0)
1394           ENDIF  ! rneb(i,k) .GT. 0.0
1395         ENDDO
1396         DO i = 1, klon
1397           IF (rneb(i,k)>0.0) THEN
1398             d_ql(i,k) = (1-zfice(i))*zoliq(i)
1399             d_qi(i,k) = zfice(i)*zoliq(i)
1400!<LTP
1401             IF (iflag_evap_prec == 4) THEN
1402                zrflcld(i) = zrflcld(i)+zqprecl(i) &
1403                 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1404                ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1405                      *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1406                znebprecipcld(i) = rneb(i,k)
1407                zrfl(i) = zrflcld(i) + zrflclr(i)
1408                zifl(i) = ziflcld(i) + ziflclr(i)       
1409!>LTP
1410             ELSE
1411                zrfl(i) = zrfl(i)+ zqprecl(i) &
1412                 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1413                zifl(i) = zifl(i)+ zqpreci(i) &
1414                      *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
1415             
1416             ENDIF !iflag_evap_prec==4
1417
1418           ENDIF                     
1419         ENDDO
1420!!
1421       ELSE  ! iflag_bergeron
1422!>CR&JYG
1423!!
1424       DO i = 1, klon
1425         IF (rneb(i,k)>0.0) THEN
1426!CR on prend en compte la phase glace
1427!JLD inutile car on ne passe jamais ici si .not.ice_thermo
1428!           if (.not.ice_thermo) then
1429!           d_ql(i,k) = zoliq(i)
1430!           d_qi(i,k) = 0.
1431!           else
1432           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1433           d_qi(i,k) = zfice(i)*zoliq(i)
1434!           endif
1435!<LTP
1436             IF (iflag_evap_prec == 4) THEN
1437                zrflcld(i) = zrflcld(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
1438                       *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1439                ziflcld(i) = ziflcld(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
1440                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
1441                znebprecipcld(i) = rneb(i,k)
1442                zrfl(i) = zrflcld(i) + zrflclr(i)
1443                zifl(i) = ziflcld(i) + ziflclr(i)       
1444!>LTP
1445             ELSE
1446!AJ<
1447                   zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
1448                       *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1449                        zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
1450                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
1451     !      zrfl(i) = zrfl(i)+  zpluie                         &
1452     !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1453     !      zifl(i) = zifl(i)+  zice                    &
1454     !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                   
1455             ENDIF !iflag_evap_prec == 4             
1456
1457!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
1458           IF ((iflag_bergeron == 1) .AND. (zt(i) < 273.15)) THEN
1459!<LTP
1460                IF (iflag_evap_prec == 4) THEN
1461                     zsolid = zrfl(i)
1462                     ziflclr(i) = ziflclr(i) +zrflclr(i)
1463                     ziflcld(i) = ziflcld(i) +zrflcld(i)
1464                     zifl(i) = ziflclr(i)+ziflcld(i)
1465                     zrflcld(i)=0.
1466                     zrflclr(i)=0.   
1467                     zrfl(i) = zrflclr(i)+zrflcld(i)
1468!>LTP
1469                ELSE
1470                     zsolid = zrfl(i)
1471                     zifl(i) = zifl(i)+zrfl(i)
1472                     zrfl(i) = 0.
1473                 ENDIF!iflag_evap_prec==4
1474
1475           if (fl_cor_ebil > 0) then
1476              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1477                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
1478           else
1479              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1480                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
1481           end if
1482           ENDIF  ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
1483!RC   
1484
1485         ENDIF ! rneb(i,k).GT.0.0
1486       ENDDO
1487
1488       ENDIF  ! iflag_bergeron .EQ. 2
1489     ENDIF  ! .NOT. ice_thermo
1490
1491!CR: la fonte est faite au debut
1492!      IF (ice_thermo) THEN
1493!       DO i = 1, klon
1494!           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
1495!           zmelt = MIN(MAX(zmelt,0.),1.)
1496!           zrfl(i)=zrfl(i)+zmelt*zifl(i)
1497!           zifl(i)=zifl(i)*(1.-zmelt)
1498!           print*,zt(i),'octavio1'
1499!           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1500!                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
1501!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
1502!       ENDDO
1503!     ENDIF
1504
1505
1506!<LTP
1507
1508!Limitation de la fraction surfacique couverte par les précipitations lorsque l'intensité locale du flux de précipitation descend en
1509!dessous de rain_int_min
1510    IF (iflag_evap_prec==4) THEN
1511        DO i=1, klon
1512            IF (zrflclr(i) + ziflclr(i) > 0 ) THEN
1513               znebprecipclr(i) = min(znebprecipclr(i), max(zrflclr(i) &
1514                    /(znebprecipclr(i)*rain_int_min), &
1515                    ziflclr(i)/(znebprecipclr(i)*rain_int_min)))
1516            ELSE
1517                znebprecipclr(i)=0.
1518            ENDIF
1519
1520            IF (zrflcld(i) + ziflcld(i) > 0 ) THEN
1521               znebprecipcld(i) = min(znebprecipcld(i), &
1522                    max(zrflcld(i)/(znebprecipcld(i)*rain_int_min), &
1523                    ziflcld(i)/(znebprecipcld(i)*rain_int_min)))
1524            ELSE
1525                znebprecipcld(i)=0.
1526            ENDIF
1527       ENDDO
1528    ENDIf
1529
1530!>LTP
1531
1532
1533
1534
1535       
1536     IF (.NOT. ice_thermo) THEN
1537       DO i = 1, klon
1538         IF (zt(i)<RTT) THEN
1539           psfl(i,k)=zrfl(i)
1540         ELSE
1541           prfl(i,k)=zrfl(i)
1542         ENDIF
1543       ENDDO
1544     ELSE
1545     ! JAM*************************************************
1546     ! Revoir partie ci-dessous: a quoi servent psfl et prfl?
1547     ! *****************************************************
1548
1549       DO i = 1, klon
1550     !   IF (zt(i).LT.RTT) THEN
1551           psfl(i,k)=zifl(i)
1552     !   ELSE
1553           prfl(i,k)=zrfl(i)
1554     !   ENDIF
1555!>AJ
1556       ENDDO
1557     ENDIF
1558     ! ----------------------------------------------------------------
1559     ! Fin de formation des precipitations
1560     ! ----------------------------------------------------------------
1561     !
1562     ! Calculer les tendances de q et de t:
1563     !
1564     DO i = 1, klon
1565        d_q(i,k) = zq(i) - q(i,k)
1566        d_t(i,k) = zt(i) - t(i,k)
1567     ENDDO
1568     !
1569     !AA--------------- Calcul du lessivage stratiforme  -------------
1570
1571     DO i = 1,klon
1572        !
1573        if(zcond(i)>zoliq(i)+1.e-10) then
1574         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1575        else
1576         beta(i,k) = 0.
1577        endif
1578        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1579             * (paprs(i,k)-paprs(i,k+1))/RG
1580        IF (rneb(i,k)>0.0.and.zprec_cond(i)>0.) THEN
1581           !AA lessivage nucleation LMD5 dans la couche elle-meme
1582          IF (iflag_t_glace==0) THEN
1583           if (t(i,k) >= t_glace_min_old) THEN
1584              zalpha_tr = a_tr_sca(3)
1585           else
1586              zalpha_tr = a_tr_sca(4)
1587           endif
1588          ELSE ! of IF (iflag_t_glace.EQ.0)
1589           if (t(i,k) >= t_glace_min) THEN
1590              zalpha_tr = a_tr_sca(3)
1591           else
1592              zalpha_tr = a_tr_sca(4)
1593           endif
1594          ENDIF
1595           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1596           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1597           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1598           !
1599           ! nucleation avec un facteur -1 au lieu de -0.5
1600           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1601           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1602        ENDIF
1603        !
1604     ENDDO      ! boucle sur i
1605     !
1606     !AA Lessivage par impaction dans les couches en-dessous
1607     DO kk = k-1, 1, -1
1608        DO i = 1, klon
1609           IF (rneb(i,k)>0.0.and.zprec_cond(i)>0.) THEN
1610             IF (iflag_t_glace==0) THEN
1611              if (t(i,kk) >= t_glace_min_old) THEN
1612                 zalpha_tr = a_tr_sca(1)
1613              else
1614                 zalpha_tr = a_tr_sca(2)
1615              endif
1616             ELSE ! of IF (iflag_t_glace.EQ.0)
1617              if (t(i,kk) >= t_glace_min) THEN
1618                 zalpha_tr = a_tr_sca(1)
1619              else
1620                 zalpha_tr = a_tr_sca(2)
1621              endif
1622             ENDIF
1623              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1624              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
1625              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1626           ENDIF
1627        ENDDO
1628     ENDDO
1629     !
1630     !AA===============================================================
1631     !                     FIN DE LA BOUCLE VERTICALE 
1632  end DO
1633  !
1634  !AA==================================================================
1635  !
1636  ! Pluie ou neige au sol selon la temperature de la 1ere couche
1637  !
1638!CR: si la thermo de la glace est active, on calcule zifl directement
1639  IF (.NOT.ice_thermo) THEN
1640  DO i = 1, klon
1641     IF ((t(i,1)+d_t(i,1)) < RTT) THEN
1642!AJ<
1643!        snow(i) = zrfl(i)
1644        snow(i) = zrfl(i)+zifl(i)
1645!>AJ
1646        zlh_solid(i) = RLSTT-RLVTT
1647     ELSE
1648        rain(i) = zrfl(i)
1649        zlh_solid(i) = 0.
1650     ENDIF
1651  ENDDO
1652
1653  ELSE
1654     DO i = 1, klon
1655        snow(i) = zifl(i)
1656        rain(i) = zrfl(i)
1657     ENDDO
1658   
1659   ENDIF
1660  !
1661  ! For energy conservation : when snow is present, the solification
1662  ! latent heat is considered.
1663!CR: si thermo de la glace, neige deja prise en compte
1664  IF (.not.ice_thermo) THEN
1665  DO k = 1, klev
1666     DO i = 1, klon
1667        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
1668        zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
1669        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1670        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i))
1671     END DO
1672  END DO
1673  ENDIF
1674  !
1675
1676  if (ncoreczq>0) then
1677     WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1678  endif
1679
1680RETURN
1681END SUBROUTINE fisrtilp
1682END MODULE lmdz_lscp_old
Note: See TracBrowser for help on using the repository browser.