Context Navigation

fisrtilp.F @ 877

Last change on this file since 877 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

Line
1	SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,
2	s d_t, d_q, d_ql, rneb, radliq, rain, snow,
3	s prfl, psfl)
4	c
5	IMPLICIT none
6	c======================================================================
7	c Auteur(s): Z.X. Li (LMD/CNRS)
8	c Date: le 20 mars 1995
9	c Objet: condensation et precipitation stratiforme.
10	c schema de nuage
11	c======================================================================
12	c======================================================================
13	#include "dimensions.h"
14	#include "dimphy.h"
15	#include "YOMCST.h"
16	c
17	c Arguments:
18	c
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)
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)
33	c
34	c Options du programme:
35	c
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
41	PARAMETER (cl=2.6e-4)
42	ccc PARAMETER (cl=2.3e-4)
43	ccc PARAMETER (cl=2.0e-4)
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
58	c
59	LOGICAL cpartiel ! condensation partielle
60	PARAMETER (cpartiel=.TRUE.)
61	REAL t_coup
62	PARAMETER (t_coup=234.0)
63	c
64	c Variables locales:
65	c
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)
74	c
75	LOGICAL appel1er
76	SAVE appel1er
77	c
78	c Fonctions en ligne:
79	c
80	REAL fallv ! vitesse de chute pour crystaux de glace
81	REAL zzz
82	#include "YOETHF.h"
83	#include "FCTTRE.h"
84	fallv (zzz) = 3.29/2.0 * ((zzz)**0.16)
85	ccc fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
86	ccc fallv (zzz) = 3.29 * ((zzz)**0.16)
87	c
88	DATA appel1er /.TRUE./
89	c
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
102	c
103	c Determiner les nuages froids par leur temperature
104	c
105	ztglace = RTT - 15.0
106	nexpo = 6
107	ccc nexpo = 1
108	c
109	c Initialiser les sorties:
110	c
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
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
130	c
131	c Initialiser le flux de precipitation a zero
132	c
133	DO i = 1, klon
134	zrfl(i) = 0.0
135	zneb(i) = seuil_neb
136	ENDDO
137	c
138	c Boucle verticale (du haut vers le bas)
139	c
140	DO 9999 k = klev, 1, -1
141	c
142	DO i = 1, klon
143	zt(i)=t(i,k)
144	zq(i)=q(i,k)
145	ENDDO
146	c
147	c Calculer l'evaporation de la precipitation
148	c
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
182	c
183	c Calculer Qs et L/Cp*dQs/dT:
184	c
185	IF (thermcep) THEN
186	DO i = 1, klon
187	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
188	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
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
207	c
208	c Determiner la condensation partielle et calculer la quantite
209	c de l'eau condensee:
210	c
211	IF (cpartiel) THEN
212	DO i = 1, klon
213	c
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
219	IF (.NOT.calcrat) ratqs=0.05
220	ccc IF (.NOT.calcrat) ratqs=0.2
221	c
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
240	c
241	DO i = 1, klon
242	zq(i) = zq(i) - zcond(i)
243	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
244	ENDDO
245	c
246	c Partager l'eau condensee en precipitation et eau liquide nuageuse
247	c
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
260	c
261	DO n = 1, ninter
262	DO i = 1, klon
263	IF (rneb(i,k).GT.0.0) THEN
264	zchau(i) = ctdtime/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
277	c
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
284	IF (zt(i).LT.RTT) THEN
285	psfl(i,k)=zrfl(i)
286	ELSE
287	prfl(i,k)=zrfl(i)
288	ENDIF
289	ENDDO
290	c
291	c Calculer les tendances de q et de t:
292	c
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
297	c
298	9999 CONTINUE
299	c
300	c Pluie ou neige au sol selon la temperature de la 1ere couche
301	c
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
309	c
310	RETURN
311	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format