source: LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp.F @ 436

Last change on this file since 436 was 410, checked in by lmdzadmin, 22 years ago

correction bug IM
IM/LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 16.0 KB
Line 
1c $Header$
2c
3      SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs,
4     s                   d_t, d_q, d_ql, rneb, radliq, rain, snow,
5     s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
6     s                   frac_impa, frac_nucl,
7     s                   prfl, psfl, rhcl)
8
9c
10      IMPLICIT none
11c======================================================================
12c Auteur(s): Z.X. Li (LMD/CNRS)
13c Date: le 20 mars 1995
14c Objet: condensation et precipitation stratiforme.
15c        schema de nuage
16c======================================================================
17c======================================================================
18#include "dimensions.h"
19#include "dimphy.h"
20#include "YOMCST.h"
21#include "tracstoke.h"
22#include "fisrtilp.h"
23c
24c Arguments:
25c
26      REAL dtime ! intervalle du temps (s)
27      REAL paprs(klon,klev+1) ! pression a inter-couche
28      REAL pplay(klon,klev) ! pression au milieu de couche
29      REAL t(klon,klev) ! temperature (K)
30      REAL q(klon,klev) ! humidite specifique (kg/kg)
31      REAL d_t(klon,klev) ! incrementation de la temperature (K)
32      REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
33      REAL d_ql(klon,klev) ! incrementation de l'eau liquide
34      REAL rneb(klon,klev) ! fraction nuageuse
35      REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
36      REAL rhcl(klon,klev) ! humidite relative en ciel clair
37      REAL rain(klon) ! pluies (mm/s)
38      REAL snow(klon) ! neige (mm/s)
39      REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
40      REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
41cAA
42c Coeffients de fraction lessivee : pour OFF-LINE
43c
44      REAL pfrac_nucl(klon,klev)
45      REAL pfrac_1nucl(klon,klev)
46      REAL pfrac_impa(klon,klev)
47c
48c Fraction d'aerosols lessivee par impaction et par nucleation
49c POur ON-LINE
50c
51      REAL frac_impa(klon,klev)
52      REAL frac_nucl(klon,klev)
53      real zct(klon),zcl(klon)
54cAA
55c
56c Options du programme:
57c
58      REAL seuil_neb ! un nuage existe vraiment au-dela
59      PARAMETER (seuil_neb=0.001)
60
61      INTEGER ninter ! sous-intervals pour la precipitation
62      PARAMETER (ninter=5)
63      LOGICAL evap_prec ! evaporation de la pluie
64      PARAMETER (evap_prec=.TRUE.)
65      REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
66      logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
67
68      real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
69      real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
70      real erf
71c
72      LOGICAL cpartiel ! condensation partielle
73      PARAMETER (cpartiel=.TRUE.)
74      REAL t_coup
75      PARAMETER (t_coup=234.0)
76c
77c Variables locales:
78c
79      INTEGER i, k, n, kk
80      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
81      REAL zrfl(klon), zrfln(klon), zqev, zqevt
82      REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
83      REAL ztglace, zt(klon)
84      INTEGER nexpo ! exponentiel pour glace/eau
85      REAL zdz(klon),zrho(klon),ztot(klon), zrhol(klon)
86      REAL zchau(klon),zfroi(klon),zfice(klon),zneb(klon)
87c
88      LOGICAL appel1er
89      SAVE appel1er
90c
91c---------------------------------------------------------------
92c
93cAA Variables traceurs:
94cAA  Provisoire !!! Parametres alpha du lessivage
95cAA  A priori on a 4 scavenging # possibles
96c
97      REAL a_tr_sca(4)
98      save a_tr_sca
99c
100c Variables intermediaires
101c
102      REAL zalpha_tr
103      REAL zfrac_lessi
104      REAL zprec_cond(klon)
105cAA
106      REAL zmair, zcpair, zcpeau
107C     Pour la conversion eau-neige
108      REAL zlh_solid(klon), zm_solid
109cIM
110      INTEGER klevm1
111c---------------------------------------------------------------
112c
113c Fonctions en ligne:
114c
115      REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
116      REAL zzz
117#include "YOETHF.h"
118#include "FCTTRE.h"
119      fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
120      fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
121c
122      DATA appel1er /.TRUE./
123
124      IF (appel1er) THEN
125c
126         PRINT*, 'fisrtilp, ninter:', ninter
127         PRINT*, 'fisrtilp, evap_prec:', evap_prec
128         PRINT*, 'fisrtilp, cpartiel:', cpartiel
129         IF (ABS(dtime/FLOAT(ninter)-360.0).GT.0.001) THEN
130          PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
131          PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
132c         CALL abort
133         ENDIF
134         appel1er = .FALSE.
135c
136cAA initialiation provisoire
137       a_tr_sca(1) = -0.5
138       a_tr_sca(2) = -0.5
139       a_tr_sca(3) = -0.5
140       a_tr_sca(4) = -0.5
141c
142cAA Initialisation a 1 des coefs des fractions lessivees
143c
144      DO k = 1, klev
145       DO i = 1, klon
146          pfrac_nucl(i,k)=1.
147          pfrac_1nucl(i,k)=1.
148          pfrac_impa(i,k)=1.
149       ENDDO
150      ENDDO
151
152      ENDIF          !  test sur appel1er
153c
154cMAf Initialisation a 0 de zoliq
155       DO i = 1, klon
156          zoliq(i)=0.
157       ENDDO
158c Determiner les nuages froids par leur temperature
159c  nexpo regle la raideur de la transition eau liquide / eau glace.
160c
161      ztglace = RTT - 15.0
162      nexpo = 6
163ccc      nexpo = 1
164c
165c Initialiser les sorties:
166c
167      DO k = 1, klev+1
168      DO i = 1, klon
169         prfl(i,k) = 0.0
170         psfl(i,k) = 0.0
171      ENDDO
172      ENDDO
173
174      DO k = 1, klev
175      DO i = 1, klon
176         d_t(i,k) = 0.0
177         d_q(i,k) = 0.0
178         d_ql(i,k) = 0.0
179         rneb(i,k) = 0.0
180         radliq(i,k) = 0.0
181         frac_nucl(i,k) = 1.
182         frac_impa(i,k) = 1.
183      ENDDO
184      ENDDO
185      DO i = 1, klon
186         rain(i) = 0.0
187         snow(i) = 0.0
188      ENDDO
189c
190c Initialiser le flux de precipitation a zero
191c
192      DO i = 1, klon
193         zrfl(i) = 0.0
194         zneb(i) = seuil_neb
195      ENDDO
196c
197c
198cAA Pour plus de securite
199
200      zalpha_tr   = 0.
201      zfrac_lessi = 0.
202
203cAA----------------------------------------------------------
204c
205c Boucle verticale (du haut vers le bas)
206c
207cIM : klevm1
208      klevm1=klev-1
209      DO 9999 k = klev, 1, -1
210c
211cAA----------------------------------------------------------
212c
213      DO i = 1, klon
214         zt(i)=t(i,k)
215         zq(i)=q(i,k)
216      ENDDO
217c
218c Calculer la varition de temp. de l'air du a la chaleur sensible
219C transporter par la pluie.
220C Il resterait a rajouter cet effet de la chaleur sensible sur les
221C flux de surface, du a la diff. de temp. entre le 1er niveau et la
222C surface.
223C
224      DO i = 1, klon
225cIM
226       IF(k.LE.klevm1) THEN         
227        zmair=(paprs(i,k)-paprs(i,k+1))/RG
228        zcpair=RCPD*(1.0+RVTMP2*zq(i))
229        zcpeau=RCPD*RVTMP2
230        zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau
231     $      + zmair*zcpair*zt(i) )
232     $      / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
233CC        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
234       ENDIF
235      ENDDO
236c
237c
238c Calculer l'evaporation de la precipitation
239c
240
241
242      IF (evap_prec) THEN
243      DO i = 1, klon
244      IF (zrfl(i) .GT.0.) THEN
245         IF (thermcep) THEN
246           zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
247           zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
248           zqs(i)=MIN(0.5,zqs(i))
249           zcor=1./(1.-RETV*zqs(i))
250           zqs(i)=zqs(i)*zcor
251         ELSE
252           IF (zt(i) .LT. t_coup) THEN
253              zqs(i) = qsats(zt(i)) / pplay(i,k)
254           ELSE
255              zqs(i) = qsatl(zt(i)) / pplay(i,k)
256           ENDIF
257         ENDIF
258         zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
259         zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
260     .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
261         zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
262     .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
263         zqev = MIN (zqev, zqevt)
264         zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
265     .                            /RG/dtime
266
267c pour la glace, on réévapore toute la précip dans la couche du dessous
268c la glace venant de la couche du dessus est simplement dans la couche
269c du dessous.
270
271         IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
272
273         zq(i) = zq(i) - (zrfln(i)-zrfl(i))
274     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
275         zt(i) = zt(i) + (zrfln(i)-zrfl(i))
276     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
277     .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
278         zrfl(i) = zrfln(i)
279      ENDIF
280      ENDDO
281      ENDIF
282c
283c Calculer Qs et L/Cp*dQs/dT:
284c
285      IF (thermcep) THEN
286         DO i = 1, klon
287           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
288           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
289           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
290           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
291           zqs(i) = MIN(0.5,zqs(i))
292           zcor = 1./(1.-RETV*zqs(i))
293           zqs(i) = zqs(i)*zcor
294           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
295         ENDDO
296      ELSE
297         DO i = 1, klon
298            IF (zt(i).LT.t_coup) THEN
299               zqs(i) = qsats(zt(i))/pplay(i,k)
300               zdqs(i) = dqsats(zt(i),zqs(i))
301            ELSE
302               zqs(i) = qsatl(zt(i))/pplay(i,k)
303               zdqs(i) = dqsatl(zt(i),zqs(i))
304            ENDIF
305         ENDDO
306      ENDIF
307c
308c Determiner la condensation partielle et calculer la quantite
309c de l'eau condensee:
310c
311      IF (cpartiel) THEN
312
313c        print*,'Dans partiel k=',k
314c
315c   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
316c   nuageuse a partir des PDF de Sandrine Bony.
317c   rneb  : fraction nuageuse
318c   zqn   : eau totale dans le nuage
319c   zcond : eau condensee moyenne dans la maille.
320c           on prend en compte le réchauffement qui diminue la partie condensee
321c
322c   Version avec les raqts
323
324         if (iflag_pdf.eq.0) then
325
326           do i=1,klon
327            zdelq = min(ratqs(i,k),0.99) * zq(i)
328            rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
329            zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
330           enddo
331
332         else
333c
334c   Version avec les nouvelles PDFs.
335           do i=1,klon
336              if(zq(i).lt.1.e-15) then
337                print*,'ZQ(',i,',',k,')=',zq(i)
338                zq(i)=1.e-15
339              endif
340           enddo
341           do i=1,klon
342            zpdf_sig(i)=ratqs(i,k)*zq(i)
343            zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
344            zpdf_delta(i)=log(zq(i)/zqs(i))
345            zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
346            zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
347            zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
348            zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
349            zpdf_e1(i)=1.-erf(zpdf_e1(i))
350            zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
351            zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
352            zpdf_e2(i)=1.-erf(zpdf_e2(i))
353            if (zpdf_e1(i).lt.1.e-10) then
354               rneb(i,k)=0.
355               zqn(i)=zqs(i)
356            else
357               rneb(i,k)=0.5*zpdf_e1(i)
358               zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
359            endif
360           
361           enddo
362
363        endif ! iflag_pdf
364
365         do i=1,klon
366            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
367            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
368            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
369c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
370c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
371c  la convection.
372c  ATTENTION !!! Il va falloir verifier tout ca.
373            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
374c           print*,'ZDQS ',zdqs(i)
375c--Olivier
376            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
377            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
378            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
379c--fin
380           ENDDO
381      ELSE
382         DO i = 1, klon
383            IF (zq(i).GT.zqs(i)) THEN
384               rneb(i,k) = 1.0
385            ELSE
386               rneb(i,k) = 0.0
387            ENDIF
388            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
389         ENDDO
390      ENDIF
391c
392      DO i = 1, klon
393         zq(i) = zq(i) - zcond(i)
394c         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
395         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
396      ENDDO
397c
398c Partager l'eau condensee en precipitation et eau liquide nuageuse
399c
400      DO i = 1, klon
401      IF (rneb(i,k).GT.0.0) THEN
402         zoliq(i) = zcond(i)
403         zrho(i) = pplay(i,k) / zt(i) / RD
404         zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
405         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
406         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
407         zfice(i) = zfice(i)**nexpo
408         zneb(i) = MAX(rneb(i,k), seuil_neb)
409         radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
410      ENDIF
411      ENDDO
412c
413      DO n = 1, ninter
414      DO i = 1, klon
415      IF (rneb(i,k).GT.0.0) THEN
416         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
417
418         if (ptconv(i,k)) then
419            zcl(i)=cld_lc_con
420            zct(i)=1./cld_tau_con
421         else
422            zcl(i)=cld_lc_lsc
423            zct(i)=1./cld_tau_lsc
424         endif
425c  quantité d'eau à élminier.
426         zchau(i) = zct(i)*dtime/FLOAT(ninter) * zoliq(i)
427     .         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl(i))**2)) *(1.-zfice(i))
428c  meme chose pour la glace.
429         if (ptconv(i,k)) then
430            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
431     .              *fallvc(zrhol(i)) * zfice(i)
432         else
433            zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
434     .              *fallvs(zrhol(i)) * zfice(i)
435         endif
436         ztot(i) = zchau(i) + zfroi(i)
437         IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
438         ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
439         zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
440         radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
441      ENDIF
442      ENDDO
443      ENDDO
444c
445      DO i = 1, klon
446      IF (rneb(i,k).GT.0.0) THEN
447         d_ql(i,k) = zoliq(i)
448         zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
449     .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
450      ENDIF
451      IF (zt(i).LT.RTT) THEN
452        psfl(i,k)=zrfl(i)
453      ELSE
454        prfl(i,k)=zrfl(i)
455      ENDIF
456      ENDDO
457c
458c Calculer les tendances de q et de t:
459c
460      DO i = 1, klon
461         d_q(i,k) = zq(i) - q(i,k)
462         d_t(i,k) = zt(i) - t(i,k)
463      ENDDO
464c
465cAA--------------- Calcul du lessivage stratiforme  -------------
466
467      DO i = 1,klon
468c
469         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
470     .                * (paprs(i,k)-paprs(i,k+1))/RG
471         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
472cAA lessivage nucleation LMD5 dans la couche elle-meme
473            if (t(i,k) .GE. ztglace) THEN
474               zalpha_tr = a_tr_sca(3)
475            else
476               zalpha_tr = a_tr_sca(4)
477            endif
478            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
479            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
480            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
481c
482c nucleation avec un facteur -1 au lieu de -0.5
483            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
484            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
485         ENDIF
486c
487      ENDDO      ! boucle sur i
488c
489cAA Lessivage par impaction dans les couches en-dessous
490      DO kk = k-1, 1, -1
491        DO i = 1, klon
492          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
493            if (t(i,kk) .GE. ztglace) THEN
494              zalpha_tr = a_tr_sca(1)
495            else
496              zalpha_tr = a_tr_sca(2)
497            endif
498            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
499            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
500            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
501          ENDIF
502        ENDDO
503      ENDDO
504c
505cAA----------------------------------------------------------
506c                     FIN DE BOUCLE SUR K   
507 9999 CONTINUE
508c
509cAA-----------------------------------------------------------
510c
511c Pluie ou neige au sol selon la temperature de la 1ere couche
512c
513      DO i = 1, klon
514      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
515         snow(i) = zrfl(i)
516         zlh_solid(i) = RLSTT-RLVTT
517      ELSE
518         rain(i) = zrfl(i)
519         zlh_solid(i) = 0.
520      ENDIF
521      ENDDO
522C
523C For energy conservation : when snow is present, the solification
524c latent heat is considered.
525      DO k = 1, klev
526        DO i = 1, klon
527          zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
528          zmair=(paprs(i,k)-paprs(i,k+1))/RG
529          zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
530          d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
531        END DO
532      END DO
533c
534      RETURN
535      END
Note: See TracBrowser for help on using the repository browser.