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

Last change on this file since 2343 was 2343, checked in by Ehouarn Millour, 9 years ago

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