source: LMDZ6/trunk/libf/phylmd/lmdz_lscp_old.f90 @ 5408

Last change on this file since 5408 was 5285, checked in by abarral, 7 weeks ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

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