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

Last change on this file since 1990 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

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