Context Navigation

concvl.F90 @ 2699

Last change on this file since 2699 was 2007, checked in by fhourdin, 11 years ago

Modification pour garantir la conservation des espèces traces.

Modification for tracer conservation

Jean-Yves Grandpeix

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: 17.4 KB

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

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

Download in other formats:

Original Format