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

Last change on this file since 741 was 559, checked in by lmdzadmin, 20 years ago

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