source: LMDZ.3.3/trunk/libf/phylmd/fisrtilp.F @ 3661

Last change on this file since 3661 was 60, checked in by lmdz, 25 years ago

Calcul de prfl et psfl en ciel clair aussi. O. Boucher
LF
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 9.4 KB
RevLine 
[2]1      SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,
[23]2     s                   d_t, d_q, d_ql, rneb, radliq, rain, snow,
3     s                   prfl, psfl)
[2]4c
5      IMPLICIT none
6c======================================================================
7c Auteur(s): Z.X. Li (LMD/CNRS)
8c Date: le 20 mars 1995
9c Objet: condensation et precipitation stratiforme.
10c        schema de nuage
11c======================================================================
12c======================================================================
13#include "dimensions.h"
14#include "dimphy.h"
15#include "YOMCST.h"
16c
17c Arguments:
18c
19      REAL dtime ! intervalle du temps (s)
20      REAL paprs(klon,klev+1) ! pression a inter-couche
21      REAL pplay(klon,klev) ! pression au milieu de couche
22      REAL t(klon,klev) ! temperature (K)
23      REAL q(klon,klev) ! humidite specifique (kg/kg)
24      REAL d_t(klon,klev) ! incrementation de la temperature (K)
25      REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
26      REAL d_ql(klon,klev) ! incrementation de l'eau liquide
27      REAL rneb(klon,klev) ! fraction nuageuse
28      REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
29      REAL rain(klon) ! pluies (mm/s)
30      REAL snow(klon) ! neige (mm/s)
[23]31      REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
32      REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
[2]33c
34c Options du programme:
35c
36      REAL seuil_neb ! un nuage existe vraiment au-dela
37      PARAMETER (seuil_neb=0.001)
38      REAL ct ! inverse du temps pour qu'un nuage precipite
39      PARAMETER (ct=1./1800.)
40      REAL cl ! seuil de precipitation
[43]41      PARAMETER (cl=2.6e-4)
42ccc      PARAMETER (cl=2.3e-4)
43ccc      PARAMETER (cl=2.0e-4)
[2]44      INTEGER ninter ! sous-intervals pour la precipitation
45      PARAMETER (ninter=5)
46      LOGICAL evap_prec ! evaporation de la pluie
47      PARAMETER (evap_prec=.TRUE.)
48      REAL coef_eva
49      PARAMETER (coef_eva=2.0E-05)
50      LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
51      REAL ratqs ! determine la largeur de distribution de vapeur
52      PARAMETER (calcrat=.TRUE.)
53      REAL zx_min, rat_max
54      PARAMETER (zx_min=1.0, rat_max=0.01)
55      REAL zx_max, rat_min
56      PARAMETER (zx_max=0.1, rat_min=0.3)
57      REAL zx
58c
59      LOGICAL cpartiel ! condensation partielle
60      PARAMETER (cpartiel=.TRUE.)
61      REAL t_coup
62      PARAMETER (t_coup=234.0)
63c
64c Variables locales:
65c
66      INTEGER i, k, n
67      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
68      REAL zrfl(klon), zrfln(klon), zqev, zqevt
69      REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
70      REAL ztglace, zt(klon)
71      INTEGER nexpo ! exponentiel pour glace/eau
72      REAL zdz(klon),zrho(klon),ztot(klon), zrhol(klon)
73      REAL zchau(klon),zfroi(klon),zfice(klon),zneb(klon)
74c
75      LOGICAL appel1er
76      SAVE appel1er
77c
78c Fonctions en ligne:
79c
80      REAL fallv ! vitesse de chute pour crystaux de glace
81      REAL zzz
82#include "YOETHF.h"
83#include "FCTTRE.h"
[43]84      fallv (zzz) = 3.29/2.0 * ((zzz)**0.16)
85ccc      fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
86ccc      fallv (zzz) = 3.29 * ((zzz)**0.16)
[2]87c
88      DATA appel1er /.TRUE./
89c
90      IF (appel1er) THEN
91         PRINT*, 'fisrtilp, calcrat:', calcrat
92         PRINT*, 'fisrtilp, ninter:', ninter
93         PRINT*, 'fisrtilp, evap_prec:', evap_prec
94         PRINT*, 'fisrtilp, cpartiel:', cpartiel
95         IF (ABS(dtime/FLOAT(ninter)-360.0).GT.0.001) THEN
96          PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
97          PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
98          CALL abort
99         ENDIF
100         appel1er = .FALSE.
101      ENDIF
102c
103c Determiner les nuages froids par leur temperature
104c
105      ztglace = RTT - 15.0
106      nexpo = 6
107ccc      nexpo = 1
108c
109c Initialiser les sorties:
110c
[23]111      DO k = 1, klev+1
112      DO i = 1, klon
113         prfl(i,k) = 0.0
114         psfl(i,k) = 0.0
115      ENDDO
116      ENDDO
[2]117      DO k = 1, klev
118      DO i = 1, klon
119         d_t(i,k) = 0.0
120         d_q(i,k) = 0.0
121         d_ql(i,k) = 0.0
122         rneb(i,k) = 0.0
123         radliq(i,k) = 0.0
124      ENDDO
125      ENDDO
126      DO i = 1, klon
127         rain(i) = 0.0
128         snow(i) = 0.0
129      ENDDO
130c
131c Initialiser le flux de precipitation a zero
132c
133      DO i = 1, klon
134         zrfl(i) = 0.0
135         zneb(i) = seuil_neb
136      ENDDO
137c
138c Boucle verticale (du haut vers le bas)
139c
140      DO 9999 k = klev, 1, -1
141c
142      DO i = 1, klon
143         zt(i)=t(i,k)
144         zq(i)=q(i,k)
145      ENDDO
146c
147c Calculer l'evaporation de la precipitation
148c
149      IF (evap_prec) THEN
150      DO i = 1, klon
151      IF (zrfl(i) .GT.0.) THEN
152         IF (thermcep) THEN
153           zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
154           zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
155           zqs(i)=MIN(0.5,zqs(i))
156           zcor=1./(1.-RETV*zqs(i))
157           zqs(i)=zqs(i)*zcor
158         ELSE
159           IF (zt(i) .LT. t_coup) THEN
160              zqs(i) = qsats(zt(i)) / pplay(i,k)
161           ELSE
162              zqs(i) = qsatl(zt(i)) / pplay(i,k)
163           ENDIF
164         ENDIF
165         zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
166         zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
167     .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
168         zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
169     .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
170         zqev = MIN (zqev, zqevt)
171         zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
172     .                            /RG/dtime
173         zq(i) = zq(i) - (zrfln(i)-zrfl(i))
174     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
175         zt(i) = zt(i) + (zrfln(i)-zrfl(i))
176     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
177     .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
178         zrfl(i) = zrfln(i)
179      ENDIF
180      ENDDO
181      ENDIF
182c
183c Calculer Qs et L/Cp*dQs/dT:
184c
185      IF (thermcep) THEN
186         DO i = 1, klon
187           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
188           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
189           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
190           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
191           zqs(i) = MIN(0.5,zqs(i))
192           zcor = 1./(1.-RETV*zqs(i))
193           zqs(i) = zqs(i)*zcor
194           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
195         ENDDO
196      ELSE
197         DO i = 1, klon
198            IF (zt(i).LT.t_coup) THEN
199               zqs(i) = qsats(zt(i))/pplay(i,k)
200               zdqs(i) = dqsats(zt(i),zqs(i))
201            ELSE
202               zqs(i) = qsatl(zt(i))/pplay(i,k)
203               zdqs(i) = dqsatl(zt(i),zqs(i))
204            ENDIF
205         ENDDO
206      ENDIF
207c
208c Determiner la condensation partielle et calculer la quantite
209c de l'eau condensee:
210c
211      IF (cpartiel) THEN
212         DO i = 1, klon
213c
214            zx = pplay(i,k)/paprs(i,1)
215            zx = (zx_max-zx)/(zx_max-zx_min)
216            zx = MIN(MAX(zx,0.0),1.0)
217            zx = zx * zx * zx
218            ratqs = zx * (rat_max-rat_min) + rat_min
[43]219            IF (.NOT.calcrat) ratqs=0.05
220ccc            IF (.NOT.calcrat) ratqs=0.2
[2]221c
222            zdelq = ratqs * zq(i)
223            rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
224            zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
225            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
226            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
227            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
228            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
229         ENDDO
230      ELSE
231         DO i = 1, klon
232            IF (zq(i).GT.zqs(i)) THEN
233               rneb(i,k) = 1.0
234            ELSE
235               rneb(i,k) = 0.0
236            ENDIF
237            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
238         ENDDO
239      ENDIF
240c
241      DO i = 1, klon
242         zq(i) = zq(i) - zcond(i)
243         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
244      ENDDO
245c
246c Partager l'eau condensee en precipitation et eau liquide nuageuse
247c
248      DO i = 1, klon
249      IF (rneb(i,k).GT.0.0) THEN
250         zoliq(i) = zcond(i)
251         zrho(i) = pplay(i,k) / zt(i) / RD
252         zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
253         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
254         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
255         zfice(i) = zfice(i)**nexpo
256         zneb(i) = MAX(rneb(i,k), seuil_neb)
257         radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
258      ENDIF
259      ENDDO
260c
261      DO n = 1, ninter
262      DO i = 1, klon
263      IF (rneb(i,k).GT.0.0) THEN
264         zchau(i) = ct*dtime/FLOAT(ninter) * zoliq(i)
265     .          * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) *(1.-zfice(i))
266         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
267         zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
268     .              *fallv(zrhol(i)) * zfice(i)
269         ztot(i) = zchau(i) + zfroi(i)
270         IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
271         ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
272         zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
273         radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
274      ENDIF
275      ENDDO
276      ENDDO
277c
278      DO i = 1, klon
279      IF (rneb(i,k).GT.0.0) THEN
280         d_ql(i,k) = zoliq(i)
281         zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
282     .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
283      ENDIF
[60]284      IF (zt(i).LT.RTT) THEN
285        psfl(i,k)=zrfl(i)
286      ELSE
287        prfl(i,k)=zrfl(i)
288      ENDIF
[2]289      ENDDO
290c
291c Calculer les tendances de q et de t:
292c
293      DO i = 1, klon
294         d_q(i,k) = zq(i) - q(i,k)
295         d_t(i,k) = zt(i) - t(i,k)
296      ENDDO
297c
298 9999 CONTINUE
299c
300c Pluie ou neige au sol selon la temperature de la 1ere couche
301c
302      DO i = 1, klon
303      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
304         snow(i) = zrfl(i)
305      ELSE
306         rain(i) = zrfl(i)
307      ENDIF
308      ENDDO
309c
310      RETURN
311      END
Note: See TracBrowser for help on using the repository browser.