source: LMDZ4/trunk/libf/phylmd/fisrtilp_tr.F @ 928

Last change on this file since 928 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: 12.8 KB
RevLine 
[524]1!
2! $Header$
3!
4c
5      SUBROUTINE fisrtilp_tr(dtime,paprs,pplay,t,q,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,
10     s                   RHcl) ! relative humidity in clear sky (needed for aer optical properties; aeropt.F)
11
12c
[766]13      USE dimphy
[524]14      IMPLICIT none
15c======================================================================
16c Auteur(s): Z.X. Li (LMD/CNRS)
17c Date: le 20 mars 1995
18c Objet: condensation et precipitation stratiforme.
19c        schema de nuage
20c======================================================================
21c======================================================================
[766]22cym#include "dimensions.h"
23cym#include "dimphy.h"
[524]24#include "YOMCST.h"
25#include "tracstoke.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 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)
43     
44Cjq   For aerosol opt properties needed (see aeropt.F)
45      REAL RHcl(klon,klev)
46     
47cAA
48c Coeffients de fraction lessivee : pour OFF-LINE
49c
50      REAL pfrac_nucl(klon,klev)
51      REAL pfrac_1nucl(klon,klev)
52      REAL pfrac_impa(klon,klev)
53c
54c Fraction d'aerosols lessivee par impaction et par nucleation
55c POur ON-LINE
56c
57      REAL frac_impa(klon,klev)
58      REAL frac_nucl(klon,klev)
59cAA
60c
61c Options du programme:
62c
63      REAL seuil_neb ! un nuage existe vraiment au-dela
64      PARAMETER (seuil_neb=0.001)
65      REAL ct ! inverse du temps pour qu'un nuage precipite
66      PARAMETER (ct=1./1800.)
67      REAL cl ! seuil de precipitation
68      PARAMETER (cl=2.6e-4)
69ccc      PARAMETER (cl=2.3e-4)
70ccc      PARAMETER (cl=2.0e-4)
71      INTEGER ninter ! sous-intervals pour la precipitation
72      PARAMETER (ninter=5)
73      LOGICAL evap_prec ! evaporation de la pluie
74      PARAMETER (evap_prec=.TRUE.)
75      REAL coef_eva
76      PARAMETER (coef_eva=2.0E-05)
77      LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
78      REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
79      PARAMETER (calcrat=.TRUE.)
80      REAL zx_min, rat_max
81      PARAMETER (zx_min=1.0, rat_max=0.01)
82      REAL zx_max, rat_min
83      PARAMETER (zx_max=0.1, rat_min=0.3)
84      REAL zx
85c
86      LOGICAL cpartiel ! condensation partielle
87      PARAMETER (cpartiel=.TRUE.)
88      REAL t_coup
89      PARAMETER (t_coup=234.0)
90c
91c Variables locales:
92c
93      INTEGER i, k, n, kk
94      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
95      REAL zrfl(klon), zrfln(klon), zqev, zqevt
96      REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
97      REAL ztglace, zt(klon)
98      INTEGER nexpo ! exponentiel pour glace/eau
99      REAL zdz(klon),zrho(klon),ztot(klon), zrhol(klon)
100      REAL zchau(klon),zfroi(klon),zfice(klon),zneb(klon)
101c
102      LOGICAL appel1er
103      SAVE appel1er
[766]104c$OMP THREADPRIVATE(appel1er)
[524]105c
106c---------------------------------------------------------------
107c
108cAA Variables traceurs:
109cAA  Provisoire !!! Parametres alpha du lessivage
110cAA  A priori on a 4 scavenging # possibles
111c
112      REAL a_tr_sca(4)
113      save a_tr_sca
[766]114c$OMP THREADPRIVATE(a_tr_sca)
[524]115c
116c Variables intermediaires
117c
118      REAL zalpha_tr
119      REAL zfrac_lessi
120      REAL zprec_cond(klon)
121cAA
122c---------------------------------------------------------------
123c
124c Fonctions en ligne:
125c
126      REAL fallv ! vitesse de chute pour crystaux de glace
127      REAL zzz
128#include "YOETHF.h"
129#include "FCTTRE.h"
130      fallv (zzz) = 3.29/2.0 * ((zzz)**0.16)
131ccc      fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
132ccc      fallv (zzz) = 3.29 * ((zzz)**0.16)
133c
134      DATA appel1er /.TRUE./
135c
136      IF (appel1er) THEN
137c
138         PRINT*, 'fisrtilp, calcrat:', calcrat
139         PRINT*, 'fisrtilp, ninter:', ninter
140         PRINT*, 'fisrtilp, evap_prec:', evap_prec
141         PRINT*, 'fisrtilp, cpartiel:', cpartiel
142         IF (ABS(dtime/FLOAT(ninter)-360.0).GT.0.001) THEN
143          PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
144          PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
145          CALL abort
146         ENDIF
147         appel1er = .FALSE.
148c
149cAA initialiation provisoire
150       a_tr_sca(1) = -0.5
151       a_tr_sca(2) = -0.5
152       a_tr_sca(3) = -0.5
153       a_tr_sca(4) = -0.5
154c
155cAA Initialisation a 1 des coefs des fractions lessivees
156c
157      DO k = 1, klev
158       DO i = 1, klon
159          pfrac_nucl(i,k)=1.
160          pfrac_1nucl(i,k)=1.
161          pfrac_impa(i,k)=1.
162       ENDDO
163      ENDDO
164
165      ENDIF          !  test sur appel1er
166c
167cMAf Initialisation a 0 de zoliq
168       DO i = 1, klon
169          zoliq(i)=0.
170       ENDDO
171c Determiner les nuages froids par leur temperature
172c
173      ztglace = RTT - 15.0
174      nexpo = 6
175ccc      nexpo = 1
176c
177c Initialiser les sorties:
178c
179      DO k = 1, klev+1
180      DO i = 1, klon
181         prfl(i,k) = 0.0
182         psfl(i,k) = 0.0
183      ENDDO
184      ENDDO
185
186      DO k = 1, klev
187      DO i = 1, klon
188         d_t(i,k) = 0.0
189         d_q(i,k) = 0.0
190         d_ql(i,k) = 0.0
191         rneb(i,k) = 0.0
192         radliq(i,k) = 0.0
193         frac_nucl(i,k) = 1.
194         frac_impa(i,k) = 1.
195      ENDDO
196      ENDDO
197      DO i = 1, klon
198         rain(i) = 0.0
199         snow(i) = 0.0
200      ENDDO
201c
202c Initialiser le flux de precipitation a zero
203c
204      DO i = 1, klon
205         zrfl(i) = 0.0
206         zneb(i) = seuil_neb
207      ENDDO
208c
209c
210cAA Pour plus de securite
211
212      zalpha_tr   = 0.
213      zfrac_lessi = 0.
214
215cAA----------------------------------------------------------
216c
217c Boucle verticale (du haut vers le bas)
218c
219      DO 9999 k = klev, 1, -1
220c
221cAA----------------------------------------------------------
222c
223      DO i = 1, klon
224         zt(i)=t(i,k)
225         zq(i)=q(i,k)
226      ENDDO
227c
228c Calculer l'evaporation de la precipitation
229c
230      IF (evap_prec) THEN
231      DO i = 1, klon
232      IF (zrfl(i) .GT.0.) THEN
233         IF (thermcep) THEN
234           zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
235           zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
236           zqs(i)=MIN(0.5,zqs(i))
237           zcor=1./(1.-RETV*zqs(i))
238           zqs(i)=zqs(i)*zcor
239         ELSE
240           IF (zt(i) .LT. t_coup) THEN
241              zqs(i) = qsats(zt(i)) / pplay(i,k)
242           ELSE
243              zqs(i) = qsatl(zt(i)) / pplay(i,k)
244           ENDIF
245         ENDIF
246         zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
247         zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
248     .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
249         zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
250     .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
251         zqev = MIN (zqev, zqevt)
252         zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
253     .                            /RG/dtime
254         zq(i) = zq(i) - (zrfln(i)-zrfl(i))
255     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
256         zt(i) = zt(i) + (zrfln(i)-zrfl(i))
257     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
258     .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
259         zrfl(i) = zrfln(i)
260      ENDIF
261      ENDDO
262      ENDIF
263c
264c Calculer Qs et L/Cp*dQs/dT:
265c
266      IF (thermcep) THEN
267         DO i = 1, klon
268           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
269           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
270           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
271           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
272           zqs(i) = MIN(0.5,zqs(i))
273           zcor = 1./(1.-RETV*zqs(i))
274           zqs(i) = zqs(i)*zcor
275           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
276         ENDDO
277      ELSE
278         DO i = 1, klon
279            IF (zt(i).LT.t_coup) THEN
280               zqs(i) = qsats(zt(i))/pplay(i,k)
281               zdqs(i) = dqsats(zt(i),zqs(i))
282            ELSE
283               zqs(i) = qsatl(zt(i))/pplay(i,k)
284               zdqs(i) = dqsatl(zt(i),zqs(i))
285            ENDIF
286         ENDDO
287      ENDIF
288c
289c Determiner la condensation partielle et calculer la quantite
290c de l'eau condensee:
291c
292      IF (cpartiel) THEN
293         DO i = 1, klon
294c
295            zdelq = ratqs(i,k) * zq(i)
296            rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
297            zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
298            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
299            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
300            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
301            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
302           
303c--Olivier
304            RHcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
305            IF (rneb(i,k) .LE. 0.0) RHcl(i,k)=zq(i)/zqs(i)
306            IF (rneb(i,k) .GE. 1.0) RHcl(i,k)=1.0
307c--fin
308           
309         ENDDO
310      ELSE
311         DO i = 1, klon
312            IF (zq(i).GT.zqs(i)) THEN
313               rneb(i,k) = 1.0
314            ELSE
315               rneb(i,k) = 0.0
316            ENDIF
317            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
318         ENDDO
319      ENDIF
320c
321      DO i = 1, klon
322         zq(i) = zq(i) - zcond(i)
323         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
324      ENDDO
325c
326c Partager l'eau condensee en precipitation et eau liquide nuageuse
327c
328      DO i = 1, klon
329      IF (rneb(i,k).GT.0.0) THEN
330         zoliq(i) = zcond(i)
331         zrho(i) = pplay(i,k) / zt(i) / RD
332         zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
333         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
334         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
335         zfice(i) = zfice(i)**nexpo
336         zneb(i) = MAX(rneb(i,k), seuil_neb)
337         radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
338      ENDIF
339      ENDDO
340c
341      DO n = 1, ninter
342      DO i = 1, klon
343      IF (rneb(i,k).GT.0.0) THEN
344         zchau(i) = ct*dtime/FLOAT(ninter) * zoliq(i)
345     .          * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) *(1.-zfice(i))
346         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
347         zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
348     .              *fallv(zrhol(i)) * zfice(i)
349         ztot(i) = zchau(i) + zfroi(i)
350         IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
351         ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
352         zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
353         radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
354      ENDIF
355      ENDDO
356      ENDDO
357c
358      DO i = 1, klon
359      IF (rneb(i,k).GT.0.0) THEN
360         d_ql(i,k) = zoliq(i)
361         zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
362     .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
363      ENDIF
364      IF (zt(i).LT.RTT) THEN
365        psfl(i,k)=zrfl(i)
366      ELSE
367        prfl(i,k)=zrfl(i)
368      ENDIF
369      ENDDO
370c
371c Calculer les tendances de q et de t:
372c
373      DO i = 1, klon
374         d_q(i,k) = zq(i) - q(i,k)
375         d_t(i,k) = zt(i) - t(i,k)
376      ENDDO
377c
378cAA--------------- Calcul du lessivage stratiforme  -------------
379
380      DO i = 1,klon
381c
382         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
383     .                * (paprs(i,k)-paprs(i,k+1))/RG
384         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
385cAA lessivage nucleation LMD5 dans la couche elle-meme
386            if (t(i,k) .GE. ztglace) THEN
387               zalpha_tr = a_tr_sca(3)
388            else
389               zalpha_tr = a_tr_sca(4)
390            endif
391            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
392            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
393            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
394c
395c nucleation avec un facteur -1 au lieu de -0.5
396            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
397            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
398         ENDIF
399c
400      ENDDO      ! boucle sur i
401c
402cAA Lessivage par impaction dans les couches en-dessous
403      DO kk = k-1, 1, -1
404        DO i = 1, klon
405          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
406            if (t(i,kk) .GE. ztglace) THEN
407              zalpha_tr = a_tr_sca(1)
408            else
409              zalpha_tr = a_tr_sca(2)
410            endif
411            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
412            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
413            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
414          ENDIF
415        ENDDO
416      ENDDO
417c
418cAA----------------------------------------------------------
419c                     FIN DE BOUCLE SUR K   
420 9999 CONTINUE
421c
422cAA-----------------------------------------------------------
423c
424c Pluie ou neige au sol selon la temperature de la 1ere couche
425c
426      DO i = 1, klon
427      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
428         snow(i) = zrfl(i)
429      ELSE
430         rain(i) = zrfl(i)
431      ENDIF
432      ENDDO
433c
434      RETURN
435      END
Note: See TracBrowser for help on using the repository browser.