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

Last change on this file since 379 was 373, checked in by lmdzadmin, 22 years ago

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