source: LMDZ4/trunk/libf/phylmd/fisrtilp.F @ 792

Last change on this file since 792 was 766, checked in by Laurent Fairhead, 17 years ago

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