Context Navigation

concvl.F @ 1742

Last change on this file since 1742 was 1742, checked in by idelkadi, 11 years ago

1- Inclusion des developpements de la these de Romain Pilon sur le
lessivage des aerosols :

a/ par les pluies convectives (modifs cv30_routines et cv3_routines pour

sortir les champs nécessaires au calcul off-line ; modif cvltr)

b/ par les pluies stratiformes (modifs phytrac et introduction

lsc_scav).

2- Choix entre plusieurs schemas pour les pluies stratiformes, commande
par iflag_lscav.

3- Quelques corrections dans la convection "Nouvelle Physique" pour
assurer la conservation des traceurs (cv3p1_mixing et cva_driver) (travail
de Robin Locatelli).

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 17.3 KB

Line
1	SUBROUTINE concvl (iflag_con,iflag_clos,
2	. dtime,paprs,pplay,
3	. t,q,t_wake,q_wake,s_wake,u,v,tra,ntra,
4	. ALE,ALP,work1,work2,
5	. d_t,d_q,d_u,d_v,d_tra,
6	. rain, snow, kbas, ktop, sigd,
7	. cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwdbis,
8	. Ma,mip,Vprecip,
9	. cape,cin,tvp,Tconv,iflag,
10	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
11	. qcondc,wd,pmflxr,pmflxs,
12	! RomP >>>
13	!! . da,phi,mp,dd_t,dd_q,lalim_conv,wght_th)
14	. da,phi,mp,phi2,d1a,dam,sij,clw,elij, ! RomP
15	. dd_t,dd_q,lalim_conv,wght_th, ! RomP
16	. evap, ep, epmlmMm,eplaMm, ! RomP
17	. wdtrainA,wdtrainM) ! RomP
18	! RomP <<<
19	***************************************************************
20	* *
21	* CONCVL *
22	* *
23	* *
24	* written by : Sandrine Bony-Lena, 17/05/2003, 11.16.04 *
25	* modified by : *
26	***************************************************************
27	*
28	c
29	USE dimphy
30	USE infotrac, ONLY : nbtr
31	IMPLICIT none
32	c======================================================================
33	c Auteur(s): S. Bony-Lena (LMD/CNRS) date: ???
34	c Objet: schema de convection de Emanuel (1991) interface
35	c======================================================================
36	c Arguments:
37	c dtime--input-R-pas d'integration (s)
38	c s-------input-R-la valeur "s" pour chaque couche
39	c sigs----input-R-la valeur "sigma" de chaque couche
40	c sig-----input-R-la valeur de "sigma" pour chaque niveau
41	c psolpa--input-R-la pression au sol (en Pa)
42	C pskapa--input-R-exponentiel kappa de psolpa
43	c h-------input-R-enthalpie potentielle (CpT/P*kappa)
44	c q-------input-R-vapeur d'eau (en kg/kg)
45	c
46	c work*: input et output: deux variables de travail,
47	c on peut les mettre a 0 au debut
48	c ALE-----input-R-energie disponible pour soulevement
49	c ALP-----input-R-puissance disponible pour soulevement
50	c
51	C d_h-----output-R-increment de l'enthalpie potentielle (h)
52	c d_q-----output-R-increment de la vapeur d'eau
53	c rain----output-R-la pluie (mm/s)
54	c snow----output-R-la neige (mm/s)
55	c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
56	c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
57	c dnwd0---output-R-unsaturated downdraft mass flux (kg/m**2/s)
58	c Ma------output-R-adiabatic ascent mass flux (kg/m2/s)
59	c mip-----output-R-mass flux shed by adiabatic ascent (kg/m2/s)
60	c Vprecip-output-R-vertical profile of precipitations (kg/m2/s)
61	c Tconv---output-R-environment temperature seen by convective scheme (K)
62	c Cape----output-R-CAPE (J/kg)
63	c Cin ----output-R-CIN (J/kg)
64	c Tvp-----output-R-Temperature virtuelle d'une parcelle soulevee
65	c adiabatiquement a partir du niveau 1 (K)
66	c deltapb-output-R-distance entre LCL et base de la colonne (<0 ; Pa)
67	c Ice_flag-input-L-TRUE->prise en compte de la thermodynamique de la glace
68	c dd_t-----output-R-increment de la temperature du aux descentes precipitantes
69	c dd_q-----output-R-increment de la vapeur d'eau du aux desc precip
70	c======================================================================
71	c
72	#include "dimensions.h"
73	c
74	INTEGER iflag_con,iflag_clos
75	c
76	REAL dtime, paprs(klon,klev+1),pplay(klon,klev)
77	REAL t(klon,klev),q(klon,klev),u(klon,klev),v(klon,klev)
78	REAL t_wake(klon,klev),q_wake(klon,klev)
79	Real s_wake(klon)
80	REAL tra(klon,klev,nbtr)
81	INTEGER ntra
82	REAL work1(klon,klev),work2(klon,klev),ptop2(klon)
83	REAL pmflxr(klon,klev+1),pmflxs(klon,klev+1)
84	REAL ALE(klon),ALP(klon)
85	c
86	REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
87	REAL dd_t(klon,klev),dd_q(klon,klev)
88	REAL d_tra(klon,klev,nbtr)
89	REAL rain(klon),snow(klon)
90	c
91	INTEGER kbas(klon),ktop(klon)
92	REAL em_ph(klon,klev+1),em_p(klon,klev)
93	REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev)
94
95	!! REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev) !jyg
96	REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev+1) !jyg
97
98	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
99	! RomP >>>
100	real phi2(klon,klev,klev)
101	real d1a(klon,klev),dam(klon,klev)
102	real sij(klon,klev,klev),clw(klon,klev),elij(klon,klev,klev)
103	REAL wdtrainA(klon,klev),wdtrainM(klon,klev)
104	REAL evap(klon,klev),ep(klon,klev)
105	REAL epmlmMm(klon,klev,klev),eplaMm(klon,klev)
106	! RomP <<<
107	REAL cape(klon),cin(klon),tvp(klon,klev)
108	REAL Tconv(klon,klev)
109	c
110	cCR:test: on passe lentr et alim_star des thermiques
111	INTEGER lalim_conv(klon)
112	REAL wght_th(klon,klev)
113	REAL em_sig1feed ! sigma at lower bound of feeding layer
114	REAL em_sig2feed ! sigma at upper bound of feeding layer
115	REAL em_wght(klev) ! weight density determining the feeding mixture
116	con enleve le save
117	c SAVE em_sig1feed,em_sig2feed,em_wght
118	c
119	INTEGER iflag(klon)
120	REAL rflag(klon)
121	REAL pbase(klon),bbase(klon)
122	REAL dtvpdt1(klon,klev),dtvpdq1(klon,klev)
123	REAL dplcldt(klon),dplcldr(klon)
124	REAL qcondc(klon,klev)
125	REAL wd(klon)
126	REAL Plim1(klon),Plim2(klon),asupmax(klon,klev)
127	REAL supmax0(klon),asupmaxmin(klon)
128	c
129	REAL sigd(klon)
130	REAL zx_t,zdelta,zx_qs,zcor
131	c
132	! INTEGER iflag_mix
133	! SAVE iflag_mix
134	INTEGER noff, minorig
135	INTEGER i,k,itra
136	REAL qs(klon,klev),qs_wake(klon,klev)
137	REAL cbmf(klon),plcl(klon),plfc(klon),wbeff(klon)
138	cLF SAVE cbmf
139	cIM/JYG REAL, SAVE, ALLOCATABLE :: cbmf(:)
140	ccc$OMP THREADPRIVATE(cbmf)!
141	REAL cbmflast(klon)
142	INTEGER ifrst
143	SAVE ifrst
144	DATA ifrst /0/
145	c$OMP THREADPRIVATE(ifrst)
146
147	c
148	C Variables supplementaires liees au bilan d'energie
149	c Real paire(klon)
150	cLF Real ql(klon,klev)
151	c Save paire
152	cLF Save ql
153	cLF Real t1(klon,klev),q1(klon,klev)
154	cLF Save t1,q1
155	c Data paire /1./
156	REAL, SAVE, ALLOCATABLE :: ql(:,:), q1(:,:), t1(:,:)
157	c$OMP THREADPRIVATE(ql, q1, t1)
158	c
159	C Variables liees au bilan d'energie et d'enthalpi
160	REAL ztsol(klon)
161	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
162	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
163	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
164	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
165	c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot)
166	c$OMP THREADPRIVATE(h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
167	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
168	REAL d_h_vcol_phy
169	REAL fs_bound, fq_bound
170	SAVE d_h_vcol_phy
171	c$OMP THREADPRIVATE(d_h_vcol_phy)
172	REAL zero_v(klon)
173	CHARACTER*15 ztit
174	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
175	SAVE ip_ebil
176	DATA ip_ebil/2/
177	c$OMP THREADPRIVATE(ip_ebil)
178	INTEGER if_ebil ! level for energy conserv. dignostics
179	SAVE if_ebil
180	DATA if_ebil/2/
181	c$OMP THREADPRIVATE(if_ebil)
182	c+jld ec_conser
183	REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique
184	REAL ZRCPD
185	c-jld ec_conser
186	cLF
187	INTEGER nloc
188	logical, save :: first=.true.
189	c$OMP THREADPRIVATE(first)
190	INTEGER, SAVE :: itap, igout
191	c$OMP THREADPRIVATE(itap, igout)
192	c
193	#include "YOMCST.h"
194	#include "YOMCST2.h"
195	#include "YOETHF.h"
196	#include "FCTTRE.h"
197	#include "iniprint.h"
198	c
199	if (first) then
200	c Allocate some variables LF 04/2008
201	c
202	cIM/JYG allocate(cbmf(klon))
203	allocate(ql(klon,klev))
204	allocate(t1(klon,klev))
205	allocate(q1(klon,klev))
206	itap=0
207	igout=klon/2+1/klon
208	endif
209	c Incrementer le compteur de la physique
210	itap = itap + 1
211
212	c Copy T into Tconv
213	DO k = 1,klev
214	DO i = 1,klon
215	Tconv(i,k) = T(i,k)
216	ENDDO
217	ENDDO
218	c
219	IF (if_ebil.ge.1) THEN
220	DO i=1,klon
221	ztsol(i) = t(i,1)
222	zero_v(i)=0.
223	Do k = 1,klev
224	ql(i,k) = 0.
225	ENDDO
226	END DO
227	END IF
228	c
229	cym
230	snow(:)=0
231
232	c IF (ifrst .EQ. 0) THEN
233	c ifrst = 1
234	if (first) then
235	first=.false.
236	c
237	C===========================================================================
238	C READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
239	C===========================================================================
240	C
241	if (iflag_con.eq.3) then
242	c CALL cv3_inicp()
243	CALL cv3_inip()
244	endif
245	c
246	C===========================================================================
247	C READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
248	C===========================================================================
249	C
250	cc$$$ open (56,file='supcrit.data')
251	cc$$$ read (56,*) Supcrit1, Supcrit2
252	cc$$$ close (56)
253	c
254	IF (prt_level .ge. 10)
255	& WRITE(lunout,*) 'supcrit1, supcrit2' ,supcrit1, supcrit2
256	C
257	C===========================================================================
258	C Initialisation pour les bilans d'eau et d'energie
259	C===========================================================================
260	IF (if_ebil.ge.1) d_h_vcol_phy=0.
261	c
262	DO i = 1, klon
263	cbmf(i) = 0.
264	!! plcl(i) = 0.
265	sigd(i) = 0.
266	ENDDO
267	ENDIF !(ifrst .EQ. 0)
268
269	c Initialisation a chaque pas de temps
270	plfc(:) = 0.
271	wbeff(:) = 100.
272	plcl(:) = 0.
273
274	DO k = 1, klev+1
275	DO i=1,klon
276	em_ph(i,k) = paprs(i,k) / 100.0
277	pmflxr(i,k)=0.
278	pmflxs(i,k)=0.
279	ENDDO
280	ENDDO
281	c
282	DO k = 1, klev
283	DO i=1,klon
284	em_p(i,k) = pplay(i,k) / 100.0
285	ENDDO
286	ENDDO
287	c
288	!
289	! Feeding layer
290	!
291	em_sig1feed = 1.
292	em_sig2feed = 0.97
293	c em_sig2feed = 0.8
294	! Relative Weight densities
295	do k=1,klev
296	em_wght(k)=1.
297	end do
298	cCRtest: couche alim des tehrmiques ponderee par a*
299	c DO i = 1, klon
300	c do k=1,lalim_conv(i)
301	c em_wght(k)=wght_th(i,k)
302	c print*,'em_wght=',em_wght(k),wght_th(i,k)
303	c end do
304	c END DO
305
306	if (iflag_con .eq. 4) then
307	DO k = 1, klev
308	DO i = 1, klon
309	zx_t = t(i,k)
310	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
311	zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
312	zcor=1./(1.-retv*zx_qs)
313	qs(i,k)=zx_qs*zcor
314	ENDDO
315	DO i = 1, klon
316	zx_t = t_wake(i,k)
317	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
318	zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
319	zcor=1./(1.-retv*zx_qs)
320	qs_wake(i,k)=zx_qs*zcor
321	ENDDO
322	ENDDO
323	else ! iflag_con=3 (modif de puristes qui fait la diffce pour la convergence numerique)
324	DO k = 1, klev
325	DO i = 1, klon
326	zx_t = t(i,k)
327	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
328	zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
329	zx_qs= MIN(0.5,zx_qs)
330	zcor=1./(1.-retv*zx_qs)
331	zx_qs=zx_qs*zcor
332	qs(i,k)=zx_qs
333	ENDDO
334	DO i = 1, klon
335	zx_t = t_wake(i,k)
336	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
337	zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
338	zx_qs= MIN(0.5,zx_qs)
339	zcor=1./(1.-retv*zx_qs)
340	zx_qs=zx_qs*zcor
341	qs_wake(i,k)=zx_qs
342	ENDDO
343	ENDDO
344	endif ! iflag_con
345	c
346	C------------------------------------------------------------------
347
348	C Main driver for convection:
349	C iflag_con=3 -> nvlle version de KE (JYG)
350	C iflag_con = 30 -> equivalent to convect3
351	C iflag_con = 4 -> equivalent to convect1/2
352	c
353	c
354	if (iflag_con.eq.30) then
355
356	print *, '-> cv_driver' !jyg
357	CALL cv_driver(klon,klev,klevp1,ntra,iflag_con,
358	: t,q,qs,u,v,tra,
359	$ em_p,em_ph,iflag,
360	$ d_t,d_q,d_u,d_v,d_tra,rain,
361	$ Vprecip,cbmf,work1,work2, !jyg
362	$ kbas,ktop,
363	$ dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
364	$ da,phi,mp,phi2,d1a,dam,sij,clw,elij, !RomP
365	$ evap,ep,epmlmMm,eplaMm, !RomP
366	$ wdtrainA,wdtrainM) !RomP
367	print *, 'cv_driver ->' !jyg
368	c
369	DO i = 1,klon
370	cbmf(i) = Ma(i,kbas(i))
371	ENDDO
372	c
373	else
374
375	cLF necessary for gathered fields
376	nloc=klon
377	CALL cva_driver(klon,klev,klev+1,ntra,nloc,
378	$ iflag_con,iflag_mix,iflag_clos,dtime,
379	: t,q,qs,t_wake,q_wake,qs_wake,s_wake,u,v,tra,
380	$ em_p,em_ph,
381	. ALE,ALP,
382	. em_sig1feed,em_sig2feed,em_wght,
383	. iflag,d_t,d_q,d_u,d_v,d_tra,rain,kbas,ktop,
384	$ cbmf,plcl,plfc,wbeff,work1,work2,ptop2,sigd,
385	$ Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
386	$ cape,cin,tvp,
387	$ dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
388	$ asupmaxmin,lalim_conv,
389	!AC!+!RomP
390	$ da,phi,mp,phi2,d1a,dam,sij,clw, ! RomP
391	$ elij,evap,ep,wdtrainA,wdtrainM) ! RomP
392	!AC!+!RomP
393	endif
394	C------------------------------------------------------------------
395	IF (prt_level .ge. 10)
396	. WRITE(lunout,*) ' cva_driver -> cbmf,plcl,plfc,wbeff ',
397	. cbmf(1),plcl(1),plfc(1),wbeff(1)
398
399	DO i = 1,klon
400	rain(i) = rain(i)/86400.
401	rflag(i)=iflag(i)
402	ENDDO
403
404	DO k = 1, klev
405	DO i = 1, klon
406	d_t(i,k) = dtime*d_t(i,k)
407	d_q(i,k) = dtime*d_q(i,k)
408	d_u(i,k) = dtime*d_u(i,k)
409	d_v(i,k) = dtime*d_v(i,k)
410	ENDDO
411	ENDDO
412	c
413	if (iflag_con.eq.30) then
414	DO itra = 1,ntra
415	DO k = 1, klev
416	DO i = 1, klon
417	d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
418	ENDDO
419	ENDDO
420	ENDDO
421	endif
422
423	c!AC!
424	if (iflag_con.eq.3) then
425	DO itra = 1,ntra
426	DO k = 1, klev
427	DO i = 1, klon
428	d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
429	ENDDO
430	ENDDO
431	ENDDO
432	endif
433	c!AC!
434
435	DO k = 1, klev
436	DO i = 1, klon
437	t1(i,k) = t(i,k)+ d_t(i,k)
438	q1(i,k) = q(i,k)+ d_q(i,k)
439	ENDDO
440	ENDDO
441	c !jyg
442	c--Separation neige/pluie (pour diagnostics) !jyg
443	DO k = 1, klev !jyg
444	DO i = 1, klon !jyg
445	IF (t1(i,k).LT.RTT) THEN !jyg
446	pmflxs(i,k)=Vprecip(i,k) !jyg
447	ELSE !jyg
448	pmflxr(i,k)=Vprecip(i,k) !jyg
449	ENDIF !jyg
450	ENDDO !jyg
451	ENDDO !jyg
452	c
453	cc IF (if_ebil.ge.2) THEN
454	cc ztit='after convect'
455	cc CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
456	cc e , t1,q1,ql,qs,u,v,paprs,pplay
457	cc s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
458	cc call diagphy(paire,ztit,ip_ebil
459	cc e , zero_v, zero_v, zero_v, zero_v, zero_v
460	cc e , zero_v, rain, zero_v, ztsol
461	cc e , d_h_vcol, d_qt, d_ec
462	cc s , fs_bound, fq_bound )
463	cc END IF
464	C
465	c
466	c les traceurs ne sont pas mis dans cette version de convect4:
467	if (iflag_con.eq.4) then
468	DO itra = 1,ntra
469	DO k = 1, klev
470	DO i = 1, klon
471	d_tra(i,k,itra) = 0.
472	ENDDO
473	ENDDO
474	ENDDO
475	endif
476	c print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
477
478	DO k = 1, klev
479	DO i = 1, klon
480	dtvpdt1(i,k) = 0.
481	dtvpdq1(i,k) = 0.
482	ENDDO
483	ENDDO
484	DO i = 1, klon
485	dplcldt(i) = 0.
486	dplcldr(i) = 0.
487	ENDDO
488	c
489	if(prt_level.GE.20) THEN
490	DO k=1,klev
491	! print*,'physiq apres_add_con i k it d_u d_v d_t d_q qdl0',igout
492	! .,k,itap,d_u_con(igout,k) ,d_v_con(igout,k), d_t_con(igout,k),
493	! .d_q_con(igout,k),dql0(igout,k)
494	! print*,'phys apres_add_con itap Ma cin ALE ALP wak t q undi t q'
495	! .,itap,Ma(igout,k),cin(igout),ALE(igout), ALP(igout),
496	! . t_wake(igout,k),q_wake(igout,k),t_undi(igout,k),q_undi(igout,k)
497	! print*,'phy apres_add_con itap CON rain snow EMA wk1 wk2 Vpp mip'
498	! .,itap,rain_con(igout),snow_con(igout),ema_work1(igout,k),
499	! .ema_work2(igout,k),Vprecip(igout,k), mip(igout,k)
500	! print*,'phy apres_add_con itap upwd dnwd dnwd0 cape tvp Tconv '
501	! .,itap,upwd(igout,k),dnwd(igout,k),dnwd0(igout,k),cape(igout),
502	! .tvp(igout,k),Tconv(igout,k)
503	! print*,'phy apres_add_con itap dtvpdt dtvdq dplcl dplcldr qcondc'
504	! .,itap,dtvpdt1(igout,k),dtvpdq1(igout,k),dplcldt(igout),
505	! .dplcldr(igout),qcondc(igout,k)
506	! print*,'phy apres_add_con itap wd pmflxr Kpmflxr Kp1 Kpmflxs Kp1'
507	! .,itap,wd(igout),pmflxr(igout,k),pmflxr(igout,k+1),pmflxs(igout,k)
508	! .,pmflxs(igout,k+1)
509	! print*,'phy apres_add_con itap da phi mp ftd fqd lalim wgth',
510	! .itap,da(igout,k),phi(igout,k,k),mp(igout,k),ftd(igout,k),
511	! . fqd(igout,k),lalim_conv(igout),wght_th(igout,k)
512	ENDDO
513	endif !(prt_level.EQ.20) THEN
514	c
515	RETURN
516	END
517

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

Download in other formats:

Original Format