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

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

Further modifications to enforce physics/dynamics separation:

  • moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
  • created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
  • created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
  • reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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