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

Last change on this file since 4564 was 4535, checked in by evignon, 18 months ago

poursuite de la replay-isation de lscp en vue de la session
de reecriture de lscp_mod en juin

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