Context Navigation

fisrtilp.F90 @ 1905

Last change on this file since 1905 was 1901, checked in by idelkadi, 11 years ago
Prise en compte de la variation de l'humidite a saturation (qsat) avec la variation de temperature associee a la condensation dans le calcul de l'eau nuageuse dans fisrtilp. Introduction d'une boucle de convergence avec un nombre maximum d'iterations fixe a iflag_fisrtilp_qsat. Si iflag_fisrtilp_qsat=0: calcul sans variation de qsat avec T Si iflag_fisrtilp_qsat=-1: calcul bogue
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 30.8 KB

Line
1	!
2	! $Id: fisrtilp.F90 1901 2013-11-19 09:29:06Z acozic $
3	!
4	!
5	SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
6	d_t, d_q, d_ql, rneb, radliq, rain, snow, &
7	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
8	frac_impa, frac_nucl, beta, &
9	prfl, psfl, rhcl, zqta, fraca, &
10	ztv, zpspsk, ztla, zthl, iflag_cldcon, &
11	iflag_ice_thermo)
12
13	!
14	USE dimphy
15	IMPLICIT none
16	!======================================================================
17	! Auteur(s): Z.X. Li (LMD/CNRS)
18	! Date: le 20 mars 1995
19	! Objet: condensation et precipitation stratiforme.
20	! schema de nuage
21	!======================================================================
22	!======================================================================
23	!ym include "dimensions.h"
24	!ym include "dimphy.h"
25	include "YOMCST.h"
26	include "tracstoke.h"
27	include "fisrtilp.h"
28	include "iniprint.h"
29
30	!
31	! Arguments:
32	!
33	REAL dtime ! intervalle du temps (s)
34	REAL paprs(klon,klev+1) ! pression a inter-couche
35	REAL pplay(klon,klev) ! pression au milieu de couche
36	REAL t(klon,klev) ! temperature (K)
37	REAL q(klon,klev) ! humidite specifique (kg/kg)
38	REAL d_t(klon,klev) ! incrementation de la temperature (K)
39	REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
40	REAL d_ql(klon,klev) ! incrementation de l'eau liquide
41	REAL rneb(klon,klev) ! fraction nuageuse
42	REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
43	REAL rhcl(klon,klev) ! humidite relative en ciel clair
44	REAL rain(klon) ! pluies (mm/s)
45	REAL snow(klon) ! neige (mm/s)
46	REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
47	REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
48	REAL ztv(klon,klev)
49	REAL zqta(klon,klev),fraca(klon,klev)
50	REAL sigma1(klon,klev),sigma2(klon,klev)
51	REAL qltot(klon,klev),ctot(klon,klev)
52	REAL zpspsk(klon,klev),ztla(klon,klev)
53	REAL zthl(klon,klev)
54	REAL ztfondue, qsl, qsi
55
56	logical lognormale(klon)
57	logical ice_thermo
58
59	!AA
60	! Coeffients de fraction lessivee : pour OFF-LINE
61	!
62	REAL pfrac_nucl(klon,klev)
63	REAL pfrac_1nucl(klon,klev)
64	REAL pfrac_impa(klon,klev)
65	!
66	! Fraction d'aerosols lessivee par impaction et par nucleation
67	! POur ON-LINE
68	!
69	REAL frac_impa(klon,klev)
70	REAL frac_nucl(klon,klev)
71	real zct ,zcl
72	!AA
73	!
74	! Options du programme:
75	!
76	REAL seuil_neb ! un nuage existe vraiment au-dela
77	PARAMETER (seuil_neb=0.001)
78
79	INTEGER ninter ! sous-intervals pour la precipitation
80	INTEGER ncoreczq
81	INTEGER iflag_cldcon
82	INTEGER iflag_ice_thermo
83	PARAMETER (ninter=5)
84	LOGICAL evap_prec ! evaporation de la pluie
85	PARAMETER (evap_prec=.TRUE.)
86	REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
87	logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
88
89	real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
90	real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
91	real erf
92	REAL qcloud(klon)
93	!
94	LOGICAL cpartiel ! condensation partielle
95	PARAMETER (cpartiel=.TRUE.)
96	REAL t_coup
97	PARAMETER (t_coup=234.0)
98	!
99	! Variables locales:
100	!
101	INTEGER i, k, n, kk
102	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
103	REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
104	LOGICAL convergence(klon)
105	REAL DDT0
106	PARAMETER (DDT0=.01)
107	INTEGER n_i, iter
108
109	REAL zrfl(klon), zrfln(klon), zqev, zqevt
110	REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
111	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
112	REAL zoliqp(klon), zoliqi(klon)
113	REAL ztglace, zt(klon)
114	INTEGER nexpo ! exponentiel pour glace/eau
115	REAL zdz(klon),zrho(klon),ztot , zrhol(klon)
116	REAL zchau ,zfroi ,zfice(klon),zneb(klon)
117	REAL zmelt, zpluie, zice, zcondold
118	PARAMETER (ztfondue=278.15)
119	!
120	LOGICAL appel1er
121	SAVE appel1er
122	!$OMP THREADPRIVATE(appel1er)
123	!
124	!---------------------------------------------------------------
125	!
126	!AA Variables traceurs:
127	!AA Provisoire !!! Parametres alpha du lessivage
128	!AA A priori on a 4 scavenging # possibles
129	!
130	REAL a_tr_sca(4)
131	save a_tr_sca
132	!$OMP THREADPRIVATE(a_tr_sca)
133	!
134	! Variables intermediaires
135	!
136	REAL zalpha_tr
137	REAL zfrac_lessi
138	REAL zprec_cond(klon)
139	!AA
140	! RomP >>> 15 nov 2012
141	REAL beta(klon,klev) ! taux de conversion de l'eau cond
142	! RomP <<<
143	REAL zmair, zcpair, zcpeau
144	! Pour la conversion eau-neige
145	REAL zlh_solid(klon), zm_solid
146	!IM
147	!ym INTEGER klevm1
148	!---------------------------------------------------------------
149	!
150	! Fonctions en ligne:
151	!
152	REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
153	REAL zzz
154	include "YOETHF.h"
155	include "FCTTRE.h"
156	fallvc (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_con
157	fallvs (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_lsc
158	!
159	DATA appel1er /.TRUE./
160	!ym
161	ice_thermo = iflag_ice_thermo .GE. 1
162	zdelq=0.0
163
164	if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
165	IF (appel1er) THEN
166	!
167	WRITE(lunout,*) 'fisrtilp, ninter:', ninter
168	WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
169	WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
170	IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
171	WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
172	WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
173	! CALL abort
174	ENDIF
175	appel1er = .FALSE.
176	!
177	!AA initialiation provisoire
178	a_tr_sca(1) = -0.5
179	a_tr_sca(2) = -0.5
180	a_tr_sca(3) = -0.5
181	a_tr_sca(4) = -0.5
182	!
183	!AA Initialisation a 1 des coefs des fractions lessivees
184	!
185	!cdir collapse
186	DO k = 1, klev
187	DO i = 1, klon
188	pfrac_nucl(i,k)=1.
189	pfrac_1nucl(i,k)=1.
190	pfrac_impa(i,k)=1.
191	beta(i,k)=0. !RomP initialisation
192	ENDDO
193	ENDDO
194
195	ENDIF ! test sur appel1er
196	!
197	!MAf Initialisation a 0 de zoliq
198	! DO i = 1, klon
199	! zoliq(i)=0.
200	! ENDDO
201	! Determiner les nuages froids par leur temperature
202	! nexpo regle la raideur de la transition eau liquide / eau glace.
203	!
204	ztglace = RTT - 15.0
205	!AJ<
206	IF (ice_thermo) THEN
207	nexpo = 2
208	ELSE
209	nexpo = 6
210	ENDIF
211	!! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
212	!! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
213	!>AJ
214	!cc nexpo = 1
215	!
216	! Initialiser les sorties:
217	!
218	!cdir collapse
219	DO k = 1, klev+1
220	DO i = 1, klon
221	prfl(i,k) = 0.0
222	psfl(i,k) = 0.0
223	ENDDO
224	ENDDO
225
226	!cdir collapse
227	DO k = 1, klev
228	DO i = 1, klon
229	d_t(i,k) = 0.0
230	d_q(i,k) = 0.0
231	d_ql(i,k) = 0.0
232	rneb(i,k) = 0.0
233	radliq(i,k) = 0.0
234	frac_nucl(i,k) = 1.
235	frac_impa(i,k) = 1.
236	ENDDO
237	ENDDO
238	DO i = 1, klon
239	rain(i) = 0.0
240	snow(i) = 0.0
241	zoliq(i)=0.
242	! ENDDO
243	!
244	! Initialiser le flux de precipitation a zero
245	!
246	! DO i = 1, klon
247	zrfl(i) = 0.0
248	zifl(i) = 0.0
249	zneb(i) = seuil_neb
250	ENDDO
251	!
252	!
253	!AA Pour plus de securite
254
255	zalpha_tr = 0.
256	zfrac_lessi = 0.
257
258	!AA----------------------------------------------------------
259	!
260	ncoreczq=0
261	! Boucle verticale (du haut vers le bas)
262	!
263	!IM : klevm1
264	!ym klevm1=klev-1
265	DO k = klev, 1, -1
266	!
267	!AA----------------------------------------------------------
268	!
269	DO i = 1, klon
270	zt(i)=t(i,k)
271	zq(i)=q(i,k)
272	ENDDO
273	!
274	! Calculer la varition de temp. de l'air du a la chaleur sensible
275	! transporter par la pluie.
276	! Il resterait a rajouter cet effet de la chaleur sensible sur les
277	! flux de surface, du a la diff. de temp. entre le 1er niveau et la
278	! surface.
279	!
280	IF(k.LE.klevm1) THEN
281	DO i = 1, klon
282	!IM
283	zmair=(paprs(i,k)-paprs(i,k+1))/RG
284	zcpair=RCPD(1.0+RVTMP2zq(i))
285	zcpeau=RCPD*RVTMP2
286	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zrfl(i)dtime*zcpeau &
287	+ zmairzcpairzt(i) ) &
288	/ (zmairzcpair + zrfl(i)dtime*zcpeau)
289	! C WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
290	ENDDO
291	ENDIF
292	!
293	!
294	! Calculer l'evaporation de la precipitation
295	!
296
297
298	! Calculer l'evaporation de la precipitation
299	!
300
301
302	IF (evap_prec) THEN
303	DO i = 1, klon
304	!AJ<
305	!! IF (zrfl(i) .GT.0.) THEN
306	IF (zrfl(i)+zifl(i).GT.0.) THEN
307	!>AJ
308	IF (thermcep) THEN
309	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
310	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
311	zqs(i)=MIN(0.5,zqs(i))
312	zcor=1./(1.-RETV*zqs(i))
313	zqs(i)=zqs(i)*zcor
314	ELSE
315	IF (zt(i) .LT. t_coup) THEN
316	zqs(i) = qsats(zt(i)) / pplay(i,k)
317	ELSE
318	zqs(i) = qsatl(zt(i)) / pplay(i,k)
319	ENDIF
320	ENDIF
321	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
322	ENDDO
323	!AJ<
324	IF (.NOT. ice_thermo) THEN
325	DO i = 1, klon
326	!AJ<
327	!! IF (zrfl(i) .GT.0.) THEN
328	IF (zrfl(i)+zifl(i).GT.0.) THEN
329	!>AJ
330	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
331	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
332	* (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
333	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
334	* RG*dtime/(paprs(i,k)-paprs(i,k+1))
335	zqev = MIN (zqev, zqevt)
336	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
337	/RG/dtime
338
339	! pour la glace, on ré-évapore toute la précip dans la
340	! couche du dessous
341	! la glace venant de la couche du dessus est simplement
342	! dans la couche du dessous.
343
344	IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
345
346	zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
347	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
348	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
349	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
350	* RLVTT/RCPD/(1.0+RVTMP2*zq(i))
351	zrfl(i) = zrfln(i)
352	zifl(i) = 0.
353	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
354	ENDDO
355	!
356	ELSE ! (.NOT. ice_thermo)
357	!
358	DO i = 1, klon
359	!AJ<
360	!! IF (zrfl(i) .GT.0.) THEN
361	IF (zrfl(i)+zifl(i).GT.0.) THEN
362	!>AJ
363	!JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
364	! Modification de l'évaporation avec la glace
365	! Différentiation entre précipitation liquide et solide
366	! On suppose que coef_evai=2*coef_eva
367	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
368
369	zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
370	! zqev0 = MAX (0.0, zqs(i)-zq(i) )
371
372	!JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
373	! On différencie qsat pour l'eau et la glace
374	! Si zdelta=1. --> glace
375	! Si zdelta=0. --> eau liquide
376	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
377
378	qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
379	qsl= MIN(0.5,qsl)
380	zcor= 1./(1.-RETV*qsl)
381	qsl= qsl*zcor
382
383	zqevt = 1.coef_eva(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
384	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
385	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
386	RGdtime/(paprs(i,k)-paprs(i,k+1))
387
388	qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
389	qsi= MIN(0.5,qsi)
390	zcor= 1./(1.-RETV*qsi)
391	qsi= qsi*zcor
392
393	zqevti = 1.coef_eva(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
394	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
395	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
396	RGdtime/(paprs(i,k)-paprs(i,k+1))
397
398
399	!JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
400	! Vérification sur l'évaporation
401	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
402
403	IF (zqevt+zqevti.GT.zqev0) THEN
404	zqev=zqev0*zqevt/(zqevt+zqevti)
405	zqevi=zqev0*zqevti/(zqevt+zqevti)
406
407	ELSE
408	IF (zqevt+zqevti.GT.0.) THEN
409	zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
410	zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
411	ELSE
412	zqev=0.
413	zqevi=0.
414	ENDIF
415	ENDIF
416
417	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
418	/RG/dtime)
419	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
420	/RG/dtime)
421
422	! Pour la glace, on révapore toute la précip dans la couche du dessous
423	! la glace venant de la couche du dessus est simplement dans la couche
424	! du dessous.
425
426	! IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
427	! print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
428	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
429	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
430	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
431	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
432	* RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
433	+ (zifln(i)-zifl(i)) &
434	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
435	* RLSTT/RCPD/(1.0+RVTMP2*zq(i))
436
437	zrfl(i) = zrfln(i)
438	zifl(i) = zifln(i)
439
440	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
441	ENDDO
442
443	ENDIF ! (.NOT. ice_thermo)
444
445	ENDIF ! (evap_prec)
446	!
447	! Calculer Qs et L/Cp*dQs/dT:
448	!
449	IF (thermcep) THEN
450	DO i = 1, klon
451	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
452	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
453	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
454	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
455	zqs(i) = MIN(0.5,zqs(i))
456	zcor = 1./(1.-RETV*zqs(i))
457	zqs(i) = zqs(i)*zcor
458	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
459	ENDDO
460	ELSE
461	DO i = 1, klon
462	IF (zt(i).LT.t_coup) THEN
463	zqs(i) = qsats(zt(i))/pplay(i,k)
464	zdqs(i) = dqsats(zt(i),zqs(i))
465	ELSE
466	zqs(i) = qsatl(zt(i))/pplay(i,k)
467	zdqs(i) = dqsatl(zt(i),zqs(i))
468	ENDIF
469	ENDDO
470	ENDIF
471	!
472	! Determiner la condensation partielle et calculer la quantite
473	! de l'eau condensee:
474	!
475	!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
476	if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
477	write(,) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
478	" but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
479	stop
480	endif
481
482	IF (cpartiel) THEN
483
484	! print*,'Dans partiel k=',k
485	!
486	! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
487	! nuageuse a partir des PDF de Sandrine Bony.
488	! rneb : fraction nuageuse
489	! zqn : eau totale dans le nuage
490	! zcond : eau condensee moyenne dans la maille.
491	! on prend en compte le réchauffement qui diminue la partie
492	! condensee
493	!
494	! Version avec les raqts
495
496	if (iflag_pdf.eq.0) then
497
498	do i=1,klon
499	zdelq = min(ratqs(i,k),0.99) * zq(i)
500	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
501	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
502	enddo
503
504	else
505	!
506	! Version avec les nouvelles PDFs.
507	do i=1,klon
508	if(zq(i).lt.1.e-15) then
509	ncoreczq=ncoreczq+1
510	zq(i)=1.e-15
511	endif
512	enddo
513
514	if (iflag_cldcon>=5) then
515
516	call cloudth(klon,klev,k,ztv, &
517	zq,zqta,fraca, &
518	qcloud,ctot,zpspsk,paprs,ztla,zthl, &
519	ratqs,zqs,t)
520
521	do i=1,klon
522	rneb(i,k)=ctot(i,k)
523	zqn(i)=qcloud(i)
524	enddo
525
526	endif
527
528	if (iflag_cldcon <= 4) then
529	lognormale = .true.
530	elseif (iflag_cldcon >= 6) then
531	! lognormale en l'absence des thermiques
532	lognormale = fraca(:,k) < 1e-10
533	else
534	! Dans le cas iflag_cldcon=5, on prend systématiquement la
535	! bi-gaussienne
536	lognormale = .false.
537	end if
538
539	do i=1,klon
540	Tbef(i)=zt(i)
541	qlbef(i)=0.
542	if (lognormale(i)) then
543	zpdf_sig(i)=ratqs(i,k)*zq(i)
544	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
545	zpdf_delta(i)=log(zq(i)/zqs(i))
546	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
547	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
548	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
549	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
550	zpdf_e1(i)=1.-erf(zpdf_e1(i))
551	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
552	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
553	zpdf_e2(i)=1.-erf(zpdf_e2(i))
554
555	if (zpdf_e1(i).lt.1.e-10) then
556	rneb(i,k)=0.
557	zqn(i)=zqs(i)
558	else
559	rneb(i,k)=0.5*zpdf_e1(i)
560	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
561	endif
562
563	qlbef(i)=max(0.,zqn(i)-zqs(i))
564	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
565	denom=1.+rneb(i,k)*zdqs(i)
566	DT(i)=num/denom
567	endif
568	enddo
569
570	n_i=1
571	! do while ((abs(DT(i)).gt.DDT0).and.((zqn(i)-zqs(i)).gt.0.))
572	! iflag_fisrtilp_qsat=0: qsat ne varie pas avec T
573	! iflag_fisrtilp_qsat > 1 : nombre d iterations max pour converger sur le calcul de qsat(T)
574	! do while (n_i.le.iflag_fisrtilp_qsat)
575	if (iflag_fisrtilp_qsat.ge.1) then
576	do iter=1,iflag_fisrtilp_qsat
577	do i=1,klon
578	convergence(i)=abs(DT(i)).gt.DDT0
579	if (convergence(i).and.lognormale(i)) then
580	Tbef(i)=Tbef(i)+DT(i)
581	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(i)))
582	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
583	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
584	zqs(i)= R2ES * FOEEW(Tbef(i),zdelta)/pplay(i,k)
585	zqs(i)=MIN(0.5,zqs(i))
586	zcor=1./(1.-retv*zqs(i))
587	zqs(i)=zqs(i)*zcor
588	zdqs(i)=FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
589
590	zpdf_sig(i)=ratqs(i,k)*zq(i)
591	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
592	zpdf_delta(i)=log(zq(i)/zqs(i))
593	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
594	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
595	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
596	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
597	zpdf_e1(i)=1.-erf(zpdf_e1(i))
598	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
599	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
600	zpdf_e2(i)=1.-erf(zpdf_e2(i))
601
602	if (zpdf_e1(i).lt.1.e-10) then
603	rneb(i,k)=0.
604	zqn(i)=zqs(i)
605	else
606	rneb(i,k)=0.5*zpdf_e1(i)
607	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
608	endif
609
610	qlbef(i)=max(0.,zqn(i)-zqs(i))
611	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
612	denom=1.+rneb(i,k)*zdqs(i)
613	DT(i)=num/denom
614	n_i=n_i+1
615	endif
616	enddo
617	enddo
618
619	endif
620
621
622	endif ! iflag_pdf
623
624	if (iflag_fisrtilp_qsat.eq.-1) then
625	!CR: ATTENTION option fausse mais a existe
626	DO i=1,klon
627	IF (rneb(i,k) .LE. 0.0) THEN
628	zqn(i) = 0.0
629	rneb(i,k) = 0.0
630	zcond(i) = 0.0
631	rhcl(i,k)=zq(i)/zqs(i)
632	ELSE IF (rneb(i,k) .GE. 1.0) THEN
633	zqn(i) = zq(i)
634	rneb(i,k) = 1.0
635	zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
636	rhcl(i,k)=1.0
637	ELSE
638	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
639	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
640	ENDIF
641	ENDDO
642
643	ELSE
644
645	DO i=1,klon
646	IF (rneb(i,k) .LE. 0.0) THEN
647	zqn(i) = 0.0
648	rneb(i,k) = 0.0
649	zcond(i) = 0.0
650	rhcl(i,k)=zq(i)/zqs(i)
651	ELSE IF (rneb(i,k) .GE. 1.0) THEN
652	zqn(i) = zq(i)
653	rneb(i,k) = 1.0
654	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
655	rhcl(i,k)=1.0
656	ELSE
657	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
658	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
659	ENDIF
660	ENDDO
661
662
663	ENDIF
664
665	! do i=1,klon
666	! IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
667	! IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
668	! rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
669	!c zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
670	!c On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
671	!c la convection.
672	!c ATTENTION !!! Il va falloir verifier tout ca.
673	! zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
674	!c print*,'ZDQS ',zdqs(i)
675	!c--Olivier
676	! rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
677	! IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
678	! IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
679	!c--fin
680	! ENDDO
681	ELSE
682	DO i = 1, klon
683	IF (zq(i).GT.zqs(i)) THEN
684	rneb(i,k) = 1.0
685	ELSE
686	rneb(i,k) = 0.0
687	ENDIF
688	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
689	ENDDO
690	ENDIF
691	!
692	DO i = 1, klon
693	zq(i) = zq(i) - zcond(i)
694	! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
695	ENDDO
696	!AJ<
697	IF (.NOT. ice_thermo) THEN
698	if (iflag_fisrtilp_qsat.lt.1) then
699	DO i = 1, klon
700	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
701	ENDDO
702	else if (iflag_fisrtilp_qsat.gt.0) then
703	DO i= 1, klon
704	if (lognormale(i)) then
705	zt(i)=Tbef(i)
706	else
707	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
708	endif
709	ENDDO
710	endif
711	ELSE
712	if (iflag_fisrtilp_qsat.lt.1) then
713	DO i = 1, klon
714	zfice(i) = 1.0 - (zt(i)-ztglace) / (273.15-ztglace)
715	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
716	zfice(i) = zfice(i)**nexpo
717	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
718	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*zq(i))
719	ENDDO
720	else
721	DO i=1, klon
722	zfice(i) = 1.0 - (zt(i)-ztglace) / (273.15-ztglace)
723	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
724	zfice(i) = zfice(i)**nexpo
725	!CR: ATTENTION zt different de Tbef: à corriger
726	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
727	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
728	ENDDO
729	endif
730	! print*,zt(i),zrfl(i),zifl(i),'temp1'
731	ENDIF
732	!>AJ
733	!
734	! Partager l'eau condensee en precipitation et eau liquide nuageuse
735	!
736	DO i = 1, klon
737	IF (rneb(i,k).GT.0.0) THEN
738	zoliq(i) = zcond(i)
739	zrho(i) = pplay(i,k) / zt(i) / RD
740	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
741	ENDIF
742	ENDDO
743	!AJ<
744	IF (.NOT. ice_thermo) THEN
745	DO i = 1, klon
746	IF (rneb(i,k).GT.0.0) THEN
747	zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
748	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
749	zfice(i) = zfice(i)**nexpo
750	!! zfice(i)=0.
751	ENDIF
752	ENDDO
753	ENDIF
754	DO i = 1, klon
755	IF (rneb(i,k).GT.0.0) THEN
756	zneb(i) = MAX(rneb(i,k), seuil_neb)
757	! zt(i) = zt(i)+zcond(i)zfice(i)RLMLT/RCPD/(1.0+RVTMP2*zq(i))
758	! print*,zt(i),'fractionglace'
759	!>AJ
760	radliq(i,k) = zoliq(i)/REAL(ninter+1)
761	ENDIF
762	ENDDO
763	!
764	DO n = 1, ninter
765	DO i = 1, klon
766	IF (rneb(i,k).GT.0.0) THEN
767	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
768	! Initialization of zpluie and zice:
769	zpluie=0
770	zice=0
771	IF (zneb(i).EQ.seuil_neb) THEN
772	ztot = 0.0
773	ELSE
774	! quantite d'eau a eliminer: zchau
775	! meme chose pour la glace: zfroi
776	if (ptconv(i,k)) then
777	zcl =cld_lc_con
778	zct =1./cld_tau_con
779	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
780	fallvc(zrhol(i)) zfice(i)
781	else
782	zcl =cld_lc_lsc
783	zct =1./cld_tau_lsc
784	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
785	fallvs(zrhol(i)) zfice(i)
786	endif
787	zchau = zct dtime/REAL(ninter) zoliq(i) &
788	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )2)) (1.-zfice(i))
789	!AJ<
790	IF (.NOT. ice_thermo) THEN
791	ztot = zchau + zfroi
792	ELSE
793	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
794	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
795	ztot = zpluie + zice
796	ENDIF
797	!>AJ
798	ztot = MAX(ztot ,0.0)
799	ENDIF
800	ztot = MIN(ztot,zoliq(i))
801	!AJ<
802	! zoliqp = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
803	! zoliqi = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
804	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
805	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
806	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
807	!>AJ
808	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
809	ENDIF
810	ENDDO
811	ENDDO
812	!
813	IF (.NOT. ice_thermo) THEN
814	DO i = 1, klon
815	IF (rneb(i,k).GT.0.0) THEN
816	d_ql(i,k) = zoliq(i)
817	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
818	* (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
819	ENDIF
820	ENDDO
821	ELSE
822	DO i = 1, klon
823	IF (rneb(i,k).GT.0.0) THEN
824	d_ql(i,k) = zoliq(i)
825	!AJ<
826	zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
827	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
828	zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
829	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
830	! zrfl(i) = zrfl(i)+ zpluie &
831	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
832	! zifl(i) = zifl(i)+ zice &
833	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
834
835	ENDIF
836	ENDDO
837	ENDIF
838
839	IF (ice_thermo) THEN
840	DO i = 1, klon
841	zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
842	zmelt = MIN(MAX(zmelt,0.),1.)
843	zrfl(i)=zrfl(i)+zmelt*zifl(i)
844	zifl(i)=zifl(i)*(1.-zmelt)
845	! print*,zt(i),'octavio1'
846	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
847	RLMLT/RCPD/(1.0+RVTMP2zq(i))
848	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
849	ENDDO
850	ENDIF
851
852
853	IF (.NOT. ice_thermo) THEN
854	DO i = 1, klon
855	IF (zt(i).LT.RTT) THEN
856	psfl(i,k)=zrfl(i)
857	ELSE
858	prfl(i,k)=zrfl(i)
859	ENDIF
860	ENDDO
861	ELSE
862	! JAM*************************************************
863	! Revoir partie ci-dessous: à quoi servent psfl et prfl?
864	! *****************************************************
865
866	DO i = 1, klon
867	! IF (zt(i).LT.RTT) THEN
868	psfl(i,k)=zifl(i)
869	! ELSE
870	prfl(i,k)=zrfl(i)
871	! ENDIF
872	!>AJ
873	ENDDO
874	ENDIF
875	!
876	!
877	! Calculer les tendances de q et de t:
878	!
879	DO i = 1, klon
880	d_q(i,k) = zq(i) - q(i,k)
881	d_t(i,k) = zt(i) - t(i,k)
882	ENDDO
883	!
884	!AA--------------- Calcul du lessivage stratiforme -------------
885
886	DO i = 1,klon
887	!
888	if(zcond(i).gt.zoliq(i)+1.e-10) then
889	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
890	else
891	beta(i,k) = 0.
892	endif
893	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
894	* (paprs(i,k)-paprs(i,k+1))/RG
895	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
896	!AA lessivage nucleation LMD5 dans la couche elle-meme
897	if (t(i,k) .GE. ztglace) THEN
898	zalpha_tr = a_tr_sca(3)
899	else
900	zalpha_tr = a_tr_sca(4)
901	endif
902	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
903	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
904	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
905	!
906	! nucleation avec un facteur -1 au lieu de -0.5
907	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
908	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
909	ENDIF
910	!
911	ENDDO ! boucle sur i
912	!
913	!AA Lessivage par impaction dans les couches en-dessous
914	DO kk = k-1, 1, -1
915	DO i = 1, klon
916	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
917	if (t(i,kk) .GE. ztglace) THEN
918	zalpha_tr = a_tr_sca(1)
919	else
920	zalpha_tr = a_tr_sca(2)
921	endif
922	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
923	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
924	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
925	ENDIF
926	ENDDO
927	ENDDO
928	!
929	!AA----------------------------------------------------------
930	! FIN DE BOUCLE SUR K
931	end DO
932	!
933	!AA-----------------------------------------------------------
934	!
935	! Pluie ou neige au sol selon la temperature de la 1ere couche
936	!
937	DO i = 1, klon
938	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
939	!AJ<
940	!! snow(i) = zrfl(i)
941	snow(i) = zrfl(i)+zifl(i)
942	!>AJ
943	zlh_solid(i) = RLSTT-RLVTT
944	ELSE
945	rain(i) = zrfl(i)
946	zlh_solid(i) = 0.
947	ENDIF
948	ENDDO
949	!
950	! For energy conservation : when snow is present, the solification
951	! latent heat is considered.
952	DO k = 1, klev
953	DO i = 1, klon
954	zcpair=RCPD(1.0+RVTMP2(q(i,k)+d_q(i,k)))
955	zmair=(paprs(i,k)-paprs(i,k+1))/RG
956	zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
957	d_t(i,k) = d_t(i,k) + zlh_solid(i) zm_solid / (zcpairzmair)
958	END DO
959	END DO
960	!
961
962	if (ncoreczq>0) then
963	WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
964	endif
965
966	END SUBROUTINE fisrtilp

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

Download in other formats:

Original Format