source: LMDZ6/trunk/libf/phylmd/fisrtilp.F90 @ 4009

Last change on this file since 4009 was 3992, checked in by musat, 3 years ago

Correction de bug
IM

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