source: LMDZ4/branches/LMDZ4_par_0/libf/phylmd/fisrtilp_tr.F @ 5456

Last change on this file since 5456 was 634, checked in by Laurent Fairhead, 20 years ago

Modifications faites à la physique pour la rendre parallele YM
Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester
les modifs pleinement avant leurs inclusions dans le tronc principal
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.7 KB
Line 
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
13      USE dimphy
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======================================================================
22cym#include "dimensions.h"
23cym#include "dimphy.h"
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
104c
105c---------------------------------------------------------------
106c
107cAA Variables traceurs:
108cAA  Provisoire !!! Parametres alpha du lessivage
109cAA  A priori on a 4 scavenging # possibles
110c
111      REAL a_tr_sca(4)
112      save a_tr_sca
113c
114c Variables intermediaires
115c
116      REAL zalpha_tr
117      REAL zfrac_lessi
118      REAL zprec_cond(klon)
119cAA
120c---------------------------------------------------------------
121c
122c Fonctions en ligne:
123c
124      REAL fallv ! vitesse de chute pour crystaux de glace
125      REAL zzz
126#include "YOETHF.h"
127#include "FCTTRE.h"
128      fallv (zzz) = 3.29/2.0 * ((zzz)**0.16)
129ccc      fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
130ccc      fallv (zzz) = 3.29 * ((zzz)**0.16)
131c
132      DATA appel1er /.TRUE./
133c
134      IF (appel1er) THEN
135c
136         PRINT*, 'fisrtilp, calcrat:', calcrat
137         PRINT*, 'fisrtilp, ninter:', ninter
138         PRINT*, 'fisrtilp, evap_prec:', evap_prec
139         PRINT*, 'fisrtilp, cpartiel:', cpartiel
140         IF (ABS(dtime/FLOAT(ninter)-360.0).GT.0.001) THEN
141          PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
142          PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
143          CALL abort
144         ENDIF
145         appel1er = .FALSE.
146c
147cAA initialiation provisoire
148       a_tr_sca(1) = -0.5
149       a_tr_sca(2) = -0.5
150       a_tr_sca(3) = -0.5
151       a_tr_sca(4) = -0.5
152c
153cAA Initialisation a 1 des coefs des fractions lessivees
154c
155      DO k = 1, klev
156       DO i = 1, klon
157          pfrac_nucl(i,k)=1.
158          pfrac_1nucl(i,k)=1.
159          pfrac_impa(i,k)=1.
160       ENDDO
161      ENDDO
162
163      ENDIF          !  test sur appel1er
164c
165cMAf Initialisation a 0 de zoliq
166       DO i = 1, klon
167          zoliq(i)=0.
168       ENDDO
169c Determiner les nuages froids par leur temperature
170c
171      ztglace = RTT - 15.0
172      nexpo = 6
173ccc      nexpo = 1
174c
175c Initialiser les sorties:
176c
177      DO k = 1, klev+1
178      DO i = 1, klon
179         prfl(i,k) = 0.0
180         psfl(i,k) = 0.0
181      ENDDO
182      ENDDO
183
184      DO k = 1, klev
185      DO i = 1, klon
186         d_t(i,k) = 0.0
187         d_q(i,k) = 0.0
188         d_ql(i,k) = 0.0
189         rneb(i,k) = 0.0
190         radliq(i,k) = 0.0
191         frac_nucl(i,k) = 1.
192         frac_impa(i,k) = 1.
193      ENDDO
194      ENDDO
195      DO i = 1, klon
196         rain(i) = 0.0
197         snow(i) = 0.0
198      ENDDO
199c
200c Initialiser le flux de precipitation a zero
201c
202      DO i = 1, klon
203         zrfl(i) = 0.0
204         zneb(i) = seuil_neb
205      ENDDO
206c
207c
208cAA Pour plus de securite
209
210      zalpha_tr   = 0.
211      zfrac_lessi = 0.
212
213cAA----------------------------------------------------------
214c
215c Boucle verticale (du haut vers le bas)
216c
217      DO 9999 k = klev, 1, -1
218c
219cAA----------------------------------------------------------
220c
221      DO i = 1, klon
222         zt(i)=t(i,k)
223         zq(i)=q(i,k)
224      ENDDO
225c
226c Calculer l'evaporation de la precipitation
227c
228      IF (evap_prec) THEN
229      DO i = 1, klon
230      IF (zrfl(i) .GT.0.) THEN
231         IF (thermcep) THEN
232           zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
233           zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
234           zqs(i)=MIN(0.5,zqs(i))
235           zcor=1./(1.-RETV*zqs(i))
236           zqs(i)=zqs(i)*zcor
237         ELSE
238           IF (zt(i) .LT. t_coup) THEN
239              zqs(i) = qsats(zt(i)) / pplay(i,k)
240           ELSE
241              zqs(i) = qsatl(zt(i)) / pplay(i,k)
242           ENDIF
243         ENDIF
244         zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
245         zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
246     .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
247         zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
248     .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
249         zqev = MIN (zqev, zqevt)
250         zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
251     .                            /RG/dtime
252         zq(i) = zq(i) - (zrfln(i)-zrfl(i))
253     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
254         zt(i) = zt(i) + (zrfln(i)-zrfl(i))
255     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
256     .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
257         zrfl(i) = zrfln(i)
258      ENDIF
259      ENDDO
260      ENDIF
261c
262c Calculer Qs et L/Cp*dQs/dT:
263c
264      IF (thermcep) THEN
265         DO i = 1, klon
266           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
267           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
268           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
269           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
270           zqs(i) = MIN(0.5,zqs(i))
271           zcor = 1./(1.-RETV*zqs(i))
272           zqs(i) = zqs(i)*zcor
273           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
274         ENDDO
275      ELSE
276         DO i = 1, klon
277            IF (zt(i).LT.t_coup) THEN
278               zqs(i) = qsats(zt(i))/pplay(i,k)
279               zdqs(i) = dqsats(zt(i),zqs(i))
280            ELSE
281               zqs(i) = qsatl(zt(i))/pplay(i,k)
282               zdqs(i) = dqsatl(zt(i),zqs(i))
283            ENDIF
284         ENDDO
285      ENDIF
286c
287c Determiner la condensation partielle et calculer la quantite
288c de l'eau condensee:
289c
290      IF (cpartiel) THEN
291         DO i = 1, klon
292c
293            zdelq = ratqs(i,k) * zq(i)
294            rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
295            zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
296            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
297            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
298            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
299            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
300           
301c--Olivier
302            RHcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
303            IF (rneb(i,k) .LE. 0.0) RHcl(i,k)=zq(i)/zqs(i)
304            IF (rneb(i,k) .GE. 1.0) RHcl(i,k)=1.0
305c--fin
306           
307         ENDDO
308      ELSE
309         DO i = 1, klon
310            IF (zq(i).GT.zqs(i)) THEN
311               rneb(i,k) = 1.0
312            ELSE
313               rneb(i,k) = 0.0
314            ENDIF
315            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
316         ENDDO
317      ENDIF
318c
319      DO i = 1, klon
320         zq(i) = zq(i) - zcond(i)
321         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
322      ENDDO
323c
324c Partager l'eau condensee en precipitation et eau liquide nuageuse
325c
326      DO i = 1, klon
327      IF (rneb(i,k).GT.0.0) THEN
328         zoliq(i) = zcond(i)
329         zrho(i) = pplay(i,k) / zt(i) / RD
330         zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
331         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
332         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
333         zfice(i) = zfice(i)**nexpo
334         zneb(i) = MAX(rneb(i,k), seuil_neb)
335         radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
336      ENDIF
337      ENDDO
338c
339      DO n = 1, ninter
340      DO i = 1, klon
341      IF (rneb(i,k).GT.0.0) THEN
342         zchau(i) = ct*dtime/FLOAT(ninter) * zoliq(i)
343     .          * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) *(1.-zfice(i))
344         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
345         zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
346     .              *fallv(zrhol(i)) * zfice(i)
347         ztot(i) = zchau(i) + zfroi(i)
348         IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
349         ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
350         zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
351         radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
352      ENDIF
353      ENDDO
354      ENDDO
355c
356      DO i = 1, klon
357      IF (rneb(i,k).GT.0.0) THEN
358         d_ql(i,k) = zoliq(i)
359         zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
360     .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
361      ENDIF
362      IF (zt(i).LT.RTT) THEN
363        psfl(i,k)=zrfl(i)
364      ELSE
365        prfl(i,k)=zrfl(i)
366      ENDIF
367      ENDDO
368c
369c Calculer les tendances de q et de t:
370c
371      DO i = 1, klon
372         d_q(i,k) = zq(i) - q(i,k)
373         d_t(i,k) = zt(i) - t(i,k)
374      ENDDO
375c
376cAA--------------- Calcul du lessivage stratiforme  -------------
377
378      DO i = 1,klon
379c
380         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
381     .                * (paprs(i,k)-paprs(i,k+1))/RG
382         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
383cAA lessivage nucleation LMD5 dans la couche elle-meme
384            if (t(i,k) .GE. ztglace) THEN
385               zalpha_tr = a_tr_sca(3)
386            else
387               zalpha_tr = a_tr_sca(4)
388            endif
389            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
390            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
391            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
392c
393c nucleation avec un facteur -1 au lieu de -0.5
394            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
395            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
396         ENDIF
397c
398      ENDDO      ! boucle sur i
399c
400cAA Lessivage par impaction dans les couches en-dessous
401      DO kk = k-1, 1, -1
402        DO i = 1, klon
403          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
404            if (t(i,kk) .GE. ztglace) THEN
405              zalpha_tr = a_tr_sca(1)
406            else
407              zalpha_tr = a_tr_sca(2)
408            endif
409            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
410            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
411            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
412          ENDIF
413        ENDDO
414      ENDDO
415c
416cAA----------------------------------------------------------
417c                     FIN DE BOUCLE SUR K   
418 9999 CONTINUE
419c
420cAA-----------------------------------------------------------
421c
422c Pluie ou neige au sol selon la temperature de la 1ere couche
423c
424      DO i = 1, klon
425      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
426         snow(i) = zrfl(i)
427      ELSE
428         rain(i) = zrfl(i)
429      ENDIF
430      ENDDO
431c
432      RETURN
433      END
Note: See TracBrowser for help on using the repository browser.