source: LMDZ.3.3/tags/IPSL-CM4_v1/libf/phylmd/fisrtilp.F @ 402

Last change on this file since 402 was 393, checked in by lmdzadmin, 22 years ago

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