source: LMDZ5/branches/IPSLCM6.0.12/libf/phylmd/fisrtilp.F90 @ 2970

Last change on this file since 2970 was 2970, checked in by Laurent Fairhead, 7 years ago

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