source: LMDZ5/trunk/libf/phylmd/fisrtilp.F90 @ 2490

Last change on this file since 2490 was 2466, checked in by fhourdin, 9 years ago

Correction de la fonte des precipitations glacees. Jean-Yves Grandpeix.

  • 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: 38.5 KB
Line 
1!
2! $Id: fisrtilp.F90 2466 2016-03-14 16:56:37Z jyg $
3!
4!
5SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
6     d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow,          &
7     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
8     frac_impa, frac_nucl, beta,                        &
9     prfl, psfl, rhcl, zqta, fraca,                     &
10     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
11     iflag_ice_thermo)
12
13  !
14  USE dimphy
15  USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14)
16  USE print_control_mod, ONLY: prt_level, lunout
17  IMPLICIT none
18  !======================================================================
19  ! Auteur(s): Z.X. Li (LMD/CNRS)
20  ! Date: le 20 mars 1995
21  ! Objet: condensation et precipitation stratiforme.
22  !        schema de nuage
23  !======================================================================
24  !======================================================================
25  include "YOMCST.h"
26  include "fisrtilp.h"
27  include "nuage.h" ! JBM (3/14)
28
29  !
30  ! Arguments:
31  !
32  REAL dtime ! intervalle du temps (s)
33  REAL paprs(klon,klev+1) ! pression a inter-couche
34  REAL pplay(klon,klev) ! pression au milieu de couche
35  REAL t(klon,klev) ! temperature (K)
36  REAL q(klon,klev) ! humidite specifique (kg/kg)
37  REAL d_t(klon,klev) ! incrementation de la temperature (K)
38  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
39  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
40  REAL d_qi(klon,klev) ! incrementation de l'eau glace
41  REAL rneb(klon,klev) ! fraction nuageuse
42  REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
43  REAL rhcl(klon,klev) ! humidite relative en ciel clair
44  REAL rain(klon) ! pluies (mm/s)
45  REAL snow(klon) ! neige (mm/s)
46  REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
47  REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
48  REAL ztv(klon,klev)
49  REAL zqta(klon,klev),fraca(klon,klev)
50  REAL sigma1(klon,klev),sigma2(klon,klev)
51  REAL qltot(klon,klev),ctot(klon,klev)
52  REAL zpspsk(klon,klev),ztla(klon,klev)
53  REAL zthl(klon,klev)
54  REAL ztfondue, qsl, qsi
55
56  logical lognormale(klon)
57  logical ice_thermo
58
59  !AA
60  ! Coeffients de fraction lessivee : pour OFF-LINE
61  !
62  REAL pfrac_nucl(klon,klev)
63  REAL pfrac_1nucl(klon,klev)
64  REAL pfrac_impa(klon,klev)
65  !
66  ! Fraction d'aerosols lessivee par impaction et par nucleation
67  ! POur ON-LINE
68  !
69  REAL frac_impa(klon,klev)
70  REAL frac_nucl(klon,klev)
71  real zct      ,zcl
72  !AA
73  !
74  ! Options du programme:
75  !
76  REAL seuil_neb ! un nuage existe vraiment au-dela
77  PARAMETER (seuil_neb=0.001)
78
79  INTEGER ninter ! sous-intervals pour la precipitation
80  INTEGER ncoreczq 
81  INTEGER iflag_cld_th
82  INTEGER iflag_ice_thermo
83  PARAMETER (ninter=5)
84  LOGICAL evap_prec ! evaporation de la pluie
85  PARAMETER (evap_prec=.TRUE.)
86  REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
87  logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
88
89  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
90  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
91  real erf   
92  REAL qcloud(klon)
93  !
94  LOGICAL cpartiel ! condensation partielle
95  PARAMETER (cpartiel=.TRUE.)
96  REAL t_coup
97  PARAMETER (t_coup=234.0)
98  !
99  ! Variables locales:
100  !
101  INTEGER i, k, n, kk
102  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 
103  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
104  LOGICAL convergence(klon)
105  REAL DDT0
106  PARAMETER (DDT0=.01)
107  INTEGER n_i(klon), iter
108  REAL cste
109 
110  REAL zrfl(klon), zrfln(klon), zqev, zqevt
111  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
112  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
113  REAL zoliqp(klon), zoliqi(klon)
114  REAL zt(klon)
115! JBM (3/14) nexpo is replaced by exposant_glace
116! REAL nexpo ! exponentiel pour glace/eau
117! INTEGER, PARAMETER :: nexpo=6
118  INTEGER exposant_glace_old
119  REAL t_glace_min_old
120  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
121  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
122  REAL zmelt, zpluie, zice, zcondold
123  PARAMETER (ztfondue=278.15)
124  REAL dzfice(klon)
125  REAL zsolid
126!!!!
127!  Variables pour Bergeron
128  REAL zcp, coef1, DeltaT
129  REAL zqpreci(klon), zqprecl(klon)
130  !
131  LOGICAL appel1er
132  SAVE appel1er
133  !$OMP THREADPRIVATE(appel1er)
134  !
135  !---------------------------------------------------------------
136  !
137  !AA Variables traceurs:
138  !AA  Provisoire !!! Parametres alpha du lessivage
139  !AA  A priori on a 4 scavenging # possibles
140  !
141  REAL a_tr_sca(4)
142  save a_tr_sca
143  !$OMP THREADPRIVATE(a_tr_sca)
144  !
145  ! Variables intermediaires
146  !
147  REAL zalpha_tr
148  REAL zfrac_lessi
149  REAL zprec_cond(klon)
150  !AA
151! RomP >>> 15 nov 2012
152  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
153! RomP <<<
154  REAL zmair, zcpair, zcpeau
155  !     Pour la conversion eau-neige
156  REAL zlh_solid(klon), zm_solid
157  !IM
158  !ym      INTEGER klevm1
159  !---------------------------------------------------------------
160  !
161  ! Fonctions en ligne:
162  !
163  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
164  REAL zzz
165  include "YOETHF.h"
166  include "FCTTRE.h"
167  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
168  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
169  !
170  DATA appel1er /.TRUE./
171  !ym
172!CR: pour iflag_ice_thermo=2, on active que la convection
173!  ice_thermo = iflag_ice_thermo .GE. 1
174   ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
175  zdelq=0.0
176
177  if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
178  IF (appel1er) THEN
179     !
180     WRITE(lunout,*) 'fisrtilp, ninter:', ninter
181     WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
182     WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
183     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
184        WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
185        WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
186        !         CALL abort
187     ENDIF
188     appel1er = .FALSE.
189     !
190     !AA initialiation provisoire
191     a_tr_sca(1) = -0.5
192     a_tr_sca(2) = -0.5
193     a_tr_sca(3) = -0.5
194     a_tr_sca(4) = -0.5
195     !
196     !AA Initialisation a 1 des coefs des fractions lessivees
197     !
198     !cdir collapse
199     DO k = 1, klev
200        DO i = 1, klon
201           pfrac_nucl(i,k)=1.
202           pfrac_1nucl(i,k)=1.
203           pfrac_impa(i,k)=1.
204           beta(i,k)=0.  !RomP initialisation
205        ENDDO
206     ENDDO
207
208  ENDIF          !  test sur appel1er
209  !
210  !MAf Initialisation a 0 de zoliq
211  !      DO i = 1, klon
212  !         zoliq(i)=0.
213  !      ENDDO
214  ! Determiner les nuages froids par leur temperature
215  !  nexpo regle la raideur de la transition eau liquide / eau glace.
216  !
217!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
218  IF (iflag_t_glace.EQ.0) THEN
219!   ztglace = RTT - 15.0
220    t_glace_min_old = RTT - 15.0
221    !AJ<
222    IF (ice_thermo) THEN
223!     nexpo = 2
224      exposant_glace_old = 2
225    ELSE
226!     nexpo = 6
227      exposant_glace_old = 6
228    ENDIF
229   
230  ENDIF
231 
232!!  RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
233!!  RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
234!>AJ
235  !cc      nexpo = 1
236  !
237  ! Initialiser les sorties:
238  !
239  !cdir collapse
240  DO k = 1, klev+1
241     DO i = 1, klon
242        prfl(i,k) = 0.0
243        psfl(i,k) = 0.0
244     ENDDO
245  ENDDO
246
247  !cdir collapse
248  DO k = 1, klev
249     DO i = 1, klon
250        d_t(i,k) = 0.0
251        d_q(i,k) = 0.0
252        d_ql(i,k) = 0.0
253        d_qi(i,k) = 0.0
254        rneb(i,k) = 0.0
255        radliq(i,k) = 0.0
256        frac_nucl(i,k) = 1.
257        frac_impa(i,k) = 1.
258     ENDDO
259  ENDDO
260  DO i = 1, klon
261     rain(i) = 0.0
262     snow(i) = 0.0
263     zoliq(i)=0.
264     !     ENDDO
265     !
266     ! Initialiser le flux de precipitation a zero
267     !
268     !     DO i = 1, klon
269     zrfl(i) = 0.0
270     zifl(i) = 0.0
271     zneb(i) = seuil_neb
272  ENDDO
273  !
274  !
275  !AA Pour plus de securite
276
277  zalpha_tr   = 0.
278  zfrac_lessi = 0.
279
280  !AA----------------------------------------------------------
281  !
282  ncoreczq=0
283  ! Boucle verticale (du haut vers le bas)
284  !
285  !IM : klevm1
286  !ym      klevm1=klev-1
287  DO k = klev, 1, -1
288     !
289     !AA----------------------------------------------------------
290     !
291     DO i = 1, klon
292        zt(i)=t(i,k)
293        zq(i)=q(i,k)
294     ENDDO
295     !
296     ! Calculer la varition de temp. de l'air du a la chaleur sensible
297     ! transporter par la pluie.
298     ! Il resterait a rajouter cet effet de la chaleur sensible sur les
299     ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
300     ! surface.
301     !
302     IF(k.LE.klevm1) THEN         
303        DO i = 1, klon
304           !IM
305           zmair=(paprs(i,k)-paprs(i,k+1))/RG
306           zcpair=RCPD*(1.0+RVTMP2*zq(i))
307           zcpeau=RCPD*RVTMP2
308           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
309                + zmair*zcpair*zt(i) ) &
310                / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
311           !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
312        ENDDO
313     ENDIF
314     !
315     !
316     ! Calculer l'evaporation de la precipitation
317     !
318
319
320     ! Calculer l'evaporation de la precipitation
321     !
322
323
324     IF (evap_prec) THEN
325        DO i = 1, klon
326!AJ<
327!!           IF (zrfl(i) .GT.0.) THEN
328           IF (zrfl(i)+zifl(i).GT.0.) THEN
329!>AJ
330              IF (thermcep) THEN
331                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
332                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
333                 zqs(i)=MIN(0.5,zqs(i))
334                 zcor=1./(1.-RETV*zqs(i))
335                 zqs(i)=zqs(i)*zcor
336              ELSE
337                 IF (zt(i) .LT. t_coup) THEN
338                    zqs(i) = qsats(zt(i)) / pplay(i,k)
339                 ELSE
340                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
341                 ENDIF
342              ENDIF
343           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
344        ENDDO
345!AJ<
346       IF (.NOT. ice_thermo) THEN
347        DO i = 1, klon
348!AJ<
349!!           IF (zrfl(i) .GT.0.) THEN
350           IF (zrfl(i)+zifl(i).GT.0.) THEN
351!>AJ
352                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
353                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
354                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
355                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
356                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
357                zqev = MIN (zqev, zqevt)
358                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
359                     /RG/dtime
360       
361                ! pour la glace, on ré-évapore toute la précip dans la
362                ! couche du dessous
363                ! la glace venant de la couche du dessus est simplement
364                ! dans la couche du dessous.
365       
366                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
367       
368                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
369                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
370                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
371                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
372                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
373                zrfl(i) = zrfln(i)
374                zifl(i) = 0.
375           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
376        ENDDO
377!
378       ELSE ! (.NOT. ice_thermo)
379!
380        DO i = 1, klon
381!AJ<
382!!           IF (zrfl(i) .GT.0.) THEN
383           IF (zrfl(i)+zifl(i).GT.0.) THEN
384!>AJ
385     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
386     ! Modification de l'évaporation avec la glace
387     ! Différentiation entre précipitation liquide et solide
388     ! On suppose que coef_evai=2*coef_eva
389     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
390     
391         zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
392     !    zqev0 = MAX (0.0, zqs(i)-zq(i) )
393
394     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
395     ! On différencie qsat pour l'eau et la glace
396     ! Si zdelta=1. --> glace
397     ! Si zdelta=0. --> eau liquide
398     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
399         
400         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
401         qsl= MIN(0.5,qsl)
402         zcor= 1./(1.-RETV*qsl)
403         qsl= qsl*zcor
404         
405         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
406              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
407         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
408              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
409
410         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
411         qsi= MIN(0.5,qsi)
412         zcor= 1./(1.-RETV*qsi)
413         qsi= qsi*zcor
414
415         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
416              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
417         zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
418              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
419
420             
421     !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
422     ! Vérification sur l'évaporation
423     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
424     
425         IF (zqevt+zqevti.GT.zqev0) THEN
426                  zqev=zqev0*zqevt/(zqevt+zqevti)
427                  zqevi=zqev0*zqevti/(zqevt+zqevti)
428             
429         ELSE
430             IF (zqevt+zqevti.GT.0.) THEN
431                  zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
432                  zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
433             ELSE
434             zqev=0.
435             zqevi=0.
436             ENDIF
437         ENDIF
438     
439         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
440                                 /RG/dtime)
441         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
442                                 /RG/dtime)
443         
444     ! Pour la glace, on révapore toute la précip dans la couche du dessous
445     ! la glace venant de la couche du dessus est simplement dans la couche
446     ! du dessous.
447     
448     !    IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
449!         print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
450         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
451                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
452         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
453                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
454                  * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
455                  + (zifln(i)-zifl(i)) &
456                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
457                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
458     
459         zrfl(i) = zrfln(i)
460         zifl(i) = zifln(i)
461
462!CR ATTENTION: deplacement de la fonte de la glace
463!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
464!!!        zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2  !!!!!!!!! jyg
465!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
466           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
467           zmelt = MIN(MAX(zmelt,0.),1.)
468           zrfl(i)=zrfl(i)+zmelt*zifl(i)
469           zifl(i)=zifl(i)*(1.-zmelt)
470!           print*,zt(i),'octavio1'
471           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
472                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
473!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
474!fin CR
475
476
477
478           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
479        ENDDO
480
481      ENDIF ! (.NOT. ice_thermo)
482     
483     ENDIF ! (evap_prec)
484     !
485     ! Calculer Qs et L/Cp*dQs/dT:
486     !
487     IF (thermcep) THEN
488        DO i = 1, klon
489           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
490           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
491           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
492           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
493           zqs(i) = MIN(0.5,zqs(i))
494           zcor = 1./(1.-RETV*zqs(i))
495           zqs(i) = zqs(i)*zcor
496           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
497        ENDDO
498     ELSE
499        DO i = 1, klon
500           IF (zt(i).LT.t_coup) THEN
501              zqs(i) = qsats(zt(i))/pplay(i,k)
502              zdqs(i) = dqsats(zt(i),zqs(i))
503           ELSE
504              zqs(i) = qsatl(zt(i))/pplay(i,k)
505              zdqs(i) = dqsatl(zt(i),zqs(i))
506           ENDIF
507        ENDDO
508     ENDIF
509     !
510     ! Determiner la condensation partielle et calculer la quantite
511     ! de l'eau condensee:
512     !
513!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
514!       if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
515!         write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
516!        " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
517!         stop
518!       endif
519
520     IF (cpartiel) THEN
521
522        !        print*,'Dans partiel k=',k
523        !
524        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
525        !   nuageuse a partir des PDF de Sandrine Bony.
526        !   rneb  : fraction nuageuse
527        !   zqn   : eau totale dans le nuage
528        !   zcond : eau condensee moyenne dans la maille.
529        !  on prend en compte le réchauffement qui diminue la partie
530        ! condensee
531        !
532        !   Version avec les raqts
533
534        if (iflag_pdf.eq.0) then
535
536           do i=1,klon
537              zdelq = min(ratqs(i,k),0.99) * zq(i)
538              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
539              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
540           enddo
541
542        else
543           !
544           !   Version avec les nouvelles PDFs.
545           do i=1,klon
546              if(zq(i).lt.1.e-15) then
547                 ncoreczq=ncoreczq+1
548                 zq(i)=1.e-15
549              endif
550           enddo
551
552           if (iflag_cld_th>=5) then
553
554              call cloudth(klon,klev,k,ztv, &
555                   zq,zqta,fraca, &
556                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
557                   ratqs,zqs,t)
558
559              do i=1,klon
560                 rneb(i,k)=ctot(i,k)
561                 zqn(i)=qcloud(i)
562              enddo
563
564           endif
565
566           if (iflag_cld_th <= 4) then
567              lognormale = .true.
568           elseif (iflag_cld_th >= 6) then
569              ! lognormale en l'absence des thermiques
570              lognormale = fraca(:,k) < 1e-10
571           else
572              ! Dans le cas iflag_cld_th=5, on prend systématiquement la
573              ! bi-gaussienne
574              lognormale = .false.
575           end if
576
577!CR: variation de qsat avec T en présence de glace ou non
578!initialisations
579           do i=1,klon
580              DT(i) = 0.
581              n_i(i)=0
582              Tbef(i)=zt(i)
583              qlbef(i)=0.
584           enddo
585
586
587!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
588           if (iflag_fisrtilp_qsat.ge.0) then
589             do iter=1,iflag_fisrtilp_qsat+1
590               
591               do i=1,klon
592!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
593                 convergence(i)=abs(DT(i)).gt.DDT0
594                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
595                 Tbef(i)=Tbef(i)+DT(i)
596                 if (.not.ice_thermo) then   
597                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
598                 else
599                    if (iflag_t_glace.eq.0) then
600                    zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
601                    else if (iflag_t_glace.eq.1) then
602                    zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
603                    endif
604                 endif
605                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
606                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
607                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
608                 zqs(i) = MIN(0.5,zqs(i))
609                 zcor = 1./(1.-RETV*zqs(i))
610                 zqs(i) = zqs(i)*zcor
611                 zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
612                 zpdf_sig(i)=ratqs(i,k)*zq(i)
613                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
614                 zpdf_delta(i)=log(zq(i)/zqs(i))
615                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
616                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
617                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
618                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
619                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
620                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
621                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
622                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
623             
624                 if (zpdf_e1(i).lt.1.e-10) then
625                    rneb(i,k)=0.
626                    zqn(i)=zqs(i)
627                 else
628                    rneb(i,k)=0.5*zpdf_e1(i)
629                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
630                 endif
631
632                 endif !convergence
633               enddo ! boucle en i
634
635                 if (.not. ice_thermo) then
636
637                 do i=1,klon
638                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
639
640                 qlbef(i)=max(0.,zqn(i)-zqs(i))
641                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
642                 denom=1.+rneb(i,k)*zdqs(i)
643                 DT(i)=num/denom
644                 n_i(i)=n_i(i)+1
645                 endif
646                 enddo
647
648                 else
649
650!calcul de la fraction de glace
651!CR: on utilise la nouvelle fonction de JBM pour l ancien calcul
652!                 zfice(i) = icefrac_lsc(Tbef(i), t_glace_min, &
653!                                     t_glace_max, exposant_glace)
654!                 zfice(i) = 1.0 - (Tbef(i)-ztglace) / (RTT-ztglace)
655!                 zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
656!                 zfice(i) = zfice(i)**nexpo
657                 if (iflag_t_glace.eq.1) then
658                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
659                 endif
660
661                 do i=1,klon
662                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
663                 
664                 if (iflag_t_glace.eq.0) then
665                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
666                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
667                    zfice(i) = zfice(i)**exposant_glace_old
668                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - RTT)
669                 endif
670                 
671                 if (iflag_t_glace.eq.1) then
672                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) / (t_glace_min - t_glace_max)
673                 endif
674               
675                 if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
676                    dzfice(i)=0.
677                 endif
678
679                 if (zfice(i).lt.1) then
680                    cste=RLVTT
681                 else
682                    cste=RLSTT
683                 endif
684
685                 qlbef(i)=max(0.,zqn(i)-zqs(i))
686                 num=-Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
687                 denom=1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
688                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
689                 DT(i)=num/denom
690                 n_i(i)=n_i(i)+1
691
692                 endif ! fin convergence
693                 enddo ! fin boucle i
694
695                 endif !ice_thermo
696
697!                 endif
698!               enddo
699 
700         
701             enddo
702           endif
703
704
705        endif ! iflag_pdf
706
707
708!        if (iflag_fisrtilp_qsat.eq.-1) then
709!CR: ATTENTION option fausse mais a existe: pour la re-activer, prendre iflag_fisrtilp_qsat=0 et activer les lignes suivantes:
710       IF (1.eq.0) THEN
711       DO i=1,klon
712           IF (rneb(i,k) .LE. 0.0) THEN
713              zqn(i) = 0.0
714              rneb(i,k) = 0.0
715              zcond(i) = 0.0
716              rhcl(i,k)=zq(i)/zqs(i)
717           ELSE IF (rneb(i,k) .GE. 1.0) THEN
718              zqn(i) = zq(i)
719              rneb(i,k) = 1.0                 
720              zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
721              rhcl(i,k)=1.0
722           ELSE
723              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
724              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
725           ENDIF
726        ENDDO
727        ENDIF
728
729!        ELSE
730
731        DO i=1,klon
732           IF (rneb(i,k) .LE. 0.0) THEN
733              zqn(i) = 0.0
734              rneb(i,k) = 0.0
735              zcond(i) = 0.0
736              rhcl(i,k)=zq(i)/zqs(i)
737           ELSE IF (rneb(i,k) .GE. 1.0) THEN
738              zqn(i) = zq(i)
739              rneb(i,k) = 1.0
740              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
741              rhcl(i,k)=1.0
742           ELSE
743              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
744              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
745           ENDIF
746        ENDDO
747
748
749!        ENDIF
750
751        !         do i=1,klon
752        !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
753        !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
754        !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
755        !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
756        !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
757        !c  la convection.
758        !c  ATTENTION !!! Il va falloir verifier tout ca.
759        !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
760        !c           print*,'ZDQS ',zdqs(i)
761        !c--Olivier
762        !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
763        !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
764        !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
765        !c--fin
766        !           ENDDO
767     ELSE
768        DO i = 1, klon
769           IF (zq(i).GT.zqs(i)) THEN
770              rneb(i,k) = 1.0
771           ELSE
772              rneb(i,k) = 0.0
773           ENDIF
774           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
775        ENDDO
776     ENDIF
777     !
778     DO i = 1, klon
779        zq(i) = zq(i) - zcond(i)
780        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
781     ENDDO
782!AJ<
783     IF (.NOT. ice_thermo) THEN
784        if (iflag_fisrtilp_qsat.lt.1) then
785           DO i = 1, klon
786             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
787           ENDDO
788        else if (iflag_fisrtilp_qsat.gt.0) then
789           DO i= 1, klon
790             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
791           ENDDO
792        endif
793     ELSE
794         if (iflag_t_glace.eq.1) then
795            CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
796         endif
797         if (iflag_fisrtilp_qsat.lt.1) then
798           DO i = 1, klon
799! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
800!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
801!                                     t_glace_max, exposant_glace)
802              if (iflag_t_glace.eq.0) then
803                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
804                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
805                    zfice(i) = zfice(i)**exposant_glace_old
806              endif
807              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
808                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
809           ENDDO
810         else
811           DO i=1, klon
812! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
813!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
814!                                     t_glace_max, exposant_glace)
815              if (iflag_t_glace.eq.0) then
816                    zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old)
817                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
818                    zfice(i) = zfice(i)**exposant_glace_old
819              endif
820              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
821                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
822           ENDDO
823         endif
824!         print*,zt(i),zrfl(i),zifl(i),'temp1'
825       ENDIF
826!>AJ
827     !
828     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
829     !
830     DO i = 1, klon
831        IF (rneb(i,k).GT.0.0) THEN
832           zoliq(i) = zcond(i)
833           zrho(i) = pplay(i,k) / zt(i) / RD
834           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
835        ENDIF
836     ENDDO
837!AJ<
838     IF (.NOT. ice_thermo) THEN
839       IF (iflag_t_glace.EQ.0) THEN
840         DO i = 1, klon
841            IF (rneb(i,k).GT.0.0) THEN
842               zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
843               zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
844               zfice(i) = zfice(i)**exposant_glace_old
845!              zfice(i) = zfice(i)**nexpo
846         !!      zfice(i)=0.
847            ENDIF
848         ENDDO
849       ELSE ! of IF (iflag_t_glace.EQ.0)
850         CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
851!         DO i = 1, klon
852!            IF (rneb(i,k).GT.0.0) THEN
853! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
854!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
855!                                     t_glace_max, exposant_glace)
856!            ENDIF
857!         ENDDO
858       ENDIF
859     ENDIF
860     DO i = 1, klon
861        IF (rneb(i,k).GT.0.0) THEN
862           zneb(i) = MAX(rneb(i,k), seuil_neb)
863     !      zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) 
864!           print*,zt(i),'fractionglace'
865!>AJ
866           radliq(i,k) = zoliq(i)/REAL(ninter+1)
867        ENDIF
868     ENDDO
869     !
870     DO n = 1, ninter
871        DO i = 1, klon
872           IF (rneb(i,k).GT.0.0) THEN
873              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
874              ! Initialization of zpluie and zice:
875              zpluie=0
876              zice=0
877              IF (zneb(i).EQ.seuil_neb) THEN
878                 ztot = 0.0
879              ELSE
880                 !  quantite d'eau a eliminer: zchau
881                 !  meme chose pour la glace: zfroi
882                 if (ptconv(i,k)) then
883                    zcl   =cld_lc_con
884                    zct   =1./cld_tau_con
885                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
886                         *fallvc(zrhol(i)) * zfice(i)
887                 else
888                    zcl   =cld_lc_lsc
889                    zct   =1./cld_tau_lsc
890                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
891                         *fallvs(zrhol(i)) * zfice(i)
892                 endif
893                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
894                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
895!AJ<
896                 IF (.NOT. ice_thermo) THEN
897                   ztot    = zchau + zfroi
898                 ELSE
899                   zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
900                   zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
901                   ztot    = zpluie    + zice
902                 ENDIF
903!>AJ
904                 ztot    = MAX(ztot   ,0.0)
905              ENDIF
906              ztot    = MIN(ztot,zoliq(i))
907!AJ<
908     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
909     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
910              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
911              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
912              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
913!>AJ
914              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
915           ENDIF
916        ENDDO  ! i = 1,klon
917     ENDDO     ! n = 1,ninter
918     !
919     IF (.NOT. ice_thermo) THEN
920       DO i = 1, klon
921         IF (rneb(i,k).GT.0.0) THEN
922           d_ql(i,k) = zoliq(i)
923           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
924                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
925         ENDIF
926       ENDDO
927     ELSE
928!
929!CR&JYG<
930! On prend en compte l'effet Bergeron dans les flux de precipitation :
931! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui
932! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat
933! et les precipitations est grossierement pris en compte en linearisant les equations
934! et en approximant le processus de precipitation liquide par un processus a seuil.
935! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération.
936! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T).
937! Le condensat precipitant solide est augmente.
938! La vapeur d'eau est augmentee.
939!
940       IF ((iflag_bergeron .EQ. 2)) THEN
941         DO i = 1, klon
942           IF (rneb(i,k) .GT. 0.0) THEN
943             zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
944             zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
945             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
946             coef1 = RLMLT*zdqs(i)/RLVTT
947             DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.)
948             zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT
949             zqprecl(i) = max( zqprecl(i) - zcp/RLMLT*(1.+coef1)*DeltaT, 0. )
950             zcond(i)   = max( zcond(i)   - zcp/RLVTT*zdqs(i)*DeltaT, 0. )
951             zq(i)      = zq(i)      + zcp/RLVTT*zdqs(i)*DeltaT
952             zt(i)      = zt(i)      + DeltaT
953           ENDIF  ! rneb(i,k) .GT. 0.0
954         ENDDO
955         DO i = 1, klon
956           IF (rneb(i,k).GT.0.0) THEN
957             d_ql(i,k) = (1-zfice(i))*zoliq(i)
958             d_qi(i,k) = zfice(i)*zoliq(i)
959             zrfl(i) = zrfl(i)+ zqprecl(i) &
960                 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
961             zifl(i) = zifl(i)+ zqpreci(i) &
962                      *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
963           ENDIF                     
964         ENDDO
965!!
966       ELSE  ! iflag_bergeron
967!>CR&JYG
968!!
969       DO i = 1, klon
970         IF (rneb(i,k).GT.0.0) THEN
971!CR on prend en compte la phase glace
972           if (.not.ice_thermo) then
973           d_ql(i,k) = zoliq(i)
974           d_qi(i,k) = 0.
975           else
976           d_ql(i,k) = (1-zfice(i))*zoliq(i)
977           d_qi(i,k) = zfice(i)*zoliq(i)
978           endif
979!AJ<
980           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
981               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
982           zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
983                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
984     !      zrfl(i) = zrfl(i)+  zpluie                         &
985     !          *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
986     !      zifl(i) = zifl(i)+  zice                    &
987     !               *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                   
988
989!CR : on prend en compte l'effet Bergeron dans les flux de precipitation
990           IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN
991              zsolid = zrfl(i)
992              zifl(i) = zifl(i)+zrfl(i)
993              zrfl(i) = 0.
994              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
995                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
996           ENDIF  ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
997!RC   
998
999         ENDIF ! rneb(i,k).GT.0.0
1000       ENDDO
1001
1002       ENDIF  ! iflag_bergeron .EQ. 2
1003     ENDIF  ! .NOT. ice_thermo
1004
1005!CR: la fonte est faite au debut
1006!      IF (ice_thermo) THEN
1007!       DO i = 1, klon
1008!           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
1009!           zmelt = MIN(MAX(zmelt,0.),1.)
1010!           zrfl(i)=zrfl(i)+zmelt*zifl(i)
1011!           zifl(i)=zifl(i)*(1.-zmelt)
1012!           print*,zt(i),'octavio1'
1013!           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
1014!                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
1015!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
1016!       ENDDO
1017!     ENDIF
1018
1019       
1020     IF (.NOT. ice_thermo) THEN
1021       DO i = 1, klon
1022         IF (zt(i).LT.RTT) THEN
1023           psfl(i,k)=zrfl(i)
1024         ELSE
1025           prfl(i,k)=zrfl(i)
1026         ENDIF
1027       ENDDO
1028     ELSE
1029     ! JAM*************************************************
1030     ! Revoir partie ci-dessous: à quoi servent psfl et prfl?
1031     ! *****************************************************
1032
1033       DO i = 1, klon
1034     !   IF (zt(i).LT.RTT) THEN
1035           psfl(i,k)=zifl(i)
1036     !   ELSE
1037           prfl(i,k)=zrfl(i)
1038     !   ENDIF
1039!>AJ
1040       ENDDO
1041     ENDIF
1042     !
1043     !
1044     ! Calculer les tendances de q et de t:
1045     !
1046     DO i = 1, klon
1047        d_q(i,k) = zq(i) - q(i,k)
1048        d_t(i,k) = zt(i) - t(i,k)
1049     ENDDO
1050     !
1051     !AA--------------- Calcul du lessivage stratiforme  -------------
1052
1053     DO i = 1,klon
1054        !
1055        if(zcond(i).gt.zoliq(i)+1.e-10) then
1056         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1057        else
1058         beta(i,k) = 0.
1059        endif
1060        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
1061             * (paprs(i,k)-paprs(i,k+1))/RG
1062        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1063           !AA lessivage nucleation LMD5 dans la couche elle-meme
1064          IF (iflag_t_glace.EQ.0) THEN
1065           if (t(i,k) .GE. t_glace_min_old) THEN
1066              zalpha_tr = a_tr_sca(3)
1067           else
1068              zalpha_tr = a_tr_sca(4)
1069           endif
1070          ELSE ! of IF (iflag_t_glace.EQ.0)
1071           if (t(i,k) .GE. t_glace_min) THEN
1072              zalpha_tr = a_tr_sca(3)
1073           else
1074              zalpha_tr = a_tr_sca(4)
1075           endif
1076          ENDIF
1077           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1078           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1079           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1080           !
1081           ! nucleation avec un facteur -1 au lieu de -0.5
1082           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1083           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
1084        ENDIF
1085        !
1086     ENDDO      ! boucle sur i
1087     !
1088     !AA Lessivage par impaction dans les couches en-dessous
1089     DO kk = k-1, 1, -1
1090        DO i = 1, klon
1091           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
1092             IF (iflag_t_glace.EQ.0) THEN
1093              if (t(i,kk) .GE. t_glace_min_old) THEN
1094                 zalpha_tr = a_tr_sca(1)
1095              else
1096                 zalpha_tr = a_tr_sca(2)
1097              endif
1098             ELSE ! of IF (iflag_t_glace.EQ.0)
1099              if (t(i,kk) .GE. t_glace_min) THEN
1100                 zalpha_tr = a_tr_sca(1)
1101              else
1102                 zalpha_tr = a_tr_sca(2)
1103              endif
1104             ENDIF
1105              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1106              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
1107              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1108           ENDIF
1109        ENDDO
1110     ENDDO
1111     !
1112     !AA----------------------------------------------------------
1113     !                     FIN DE BOUCLE SUR K   
1114  end DO
1115  !
1116  !AA-----------------------------------------------------------
1117  !
1118  ! Pluie ou neige au sol selon la temperature de la 1ere couche
1119  !
1120!CR: si la thermo de la glace est active, on calcule zifl directement
1121  IF (.NOT.ice_thermo) THEN
1122  DO i = 1, klon
1123     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
1124!AJ<
1125!        snow(i) = zrfl(i)
1126        snow(i) = zrfl(i)+zifl(i)
1127!>AJ
1128        zlh_solid(i) = RLSTT-RLVTT
1129     ELSE
1130        rain(i) = zrfl(i)
1131        zlh_solid(i) = 0.
1132     ENDIF
1133  ENDDO
1134
1135  ELSE
1136     DO i = 1, klon
1137        snow(i) = zifl(i)
1138        rain(i) = zrfl(i)
1139     ENDDO
1140   
1141   ENDIF
1142  !
1143  ! For energy conservation : when snow is present, the solification
1144  ! latent heat is considered.
1145!CR: si thermo de la glace, neige deja prise en compte
1146  IF (.not.ice_thermo) THEN
1147  DO k = 1, klev
1148     DO i = 1, klon
1149        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
1150        zmair=(paprs(i,k)-paprs(i,k+1))/RG
1151        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
1152        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
1153     END DO
1154  END DO
1155  ENDIF
1156  !
1157
1158  if (ncoreczq>0) then
1159     WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
1160  endif
1161
1162END SUBROUTINE fisrtilp
Note: See TracBrowser for help on using the repository browser.