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

Last change on this file since 2865 was 2814, checked in by jyg, 8 years ago

Cosmetic changes to fisrtilp

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