Context Navigation

concvl.F @ 5342

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

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

Download in other formats:

Original Format