source: LMDZ5/trunk/libf/phylmd/fisrtilp_tr.F @ 1969

Last change on this file since 1969 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

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