Context Navigation

concvl.F @ 923

Last change on this file since 923 was 879, checked in by Laurent Fairhead, 17 years ago
Suite de la bascule vers une physique avec thermiques, nouvelle convection, poche froide ... LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 12.3 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE concvl (iflag_con,iflag_clos,
5	. dtime,paprs,pplay,
6	. t,q,t_wake,q_wake,u,v,tra,ntra,
7	. ALE,ALP,work1,work2,
8	. d_t,d_q,d_u,d_v,d_tra,
9	. rain, snow, kbas, ktop, sigd,
10	. upwd,dnwd,dnwdbis,Ma,mip,Vprecip,
11	. cape,cin,tvp,Tconv,iflag,
12	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
13	. qcondc,wd,pmflxr,pmflxs,
14	. da,phi,mp,dd_t,dd_q,lalim_conv,wght_th)
15	***************************************************************
16	* *
17	* CONCVL *
18	* *
19	* *
20	* written by : Sandrine Bony-Lena, 17/05/2003, 11.16.04 *
21	* modified by : *
22	***************************************************************
23	*
24	c
25	c USE dimphy
26	IMPLICIT none
27	c======================================================================
28	c Auteur(s): S. Bony-Lena (LMD/CNRS) date: ???
29	c Objet: schema de convection de Emanuel (1991) interface
30	c======================================================================
31	c Arguments:
32	c dtime--input-R-pas d'integration (s)
33	c s-------input-R-la valeur "s" pour chaque couche
34	c sigs----input-R-la valeur "sigma" de chaque couche
35	c sig-----input-R-la valeur de "sigma" pour chaque niveau
36	c psolpa--input-R-la pression au sol (en Pa)
37	C pskapa--input-R-exponentiel kappa de psolpa
38	c h-------input-R-enthalpie potentielle (CpT/P*kappa)
39	c q-------input-R-vapeur d'eau (en kg/kg)
40	c
41	c work*: input et output: deux variables de travail,
42	c on peut les mettre a 0 au debut
43	c ALE-----input-R-energie disponible pour soulevement
44	c ALP-----input-R-puissance disponible pour soulevement
45	c
46	C d_h-----output-R-increment de l'enthalpie potentielle (h)
47	c d_q-----output-R-increment de la vapeur d'eau
48	c rain----output-R-la pluie (mm/s)
49	c snow----output-R-la neige (mm/s)
50	c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
51	c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
52	c dnwd0---output-R-unsaturated downdraft mass flux (kg/m**2/s)
53	c Ma------output-R-adiabatic ascent mass flux (kg/m2/s)
54	c mip-----output-R-mass flux shed by adiabatic ascent (kg/m2/s)
55	c Vprecip-output-R-vertical profile of precipitations (kg/m2/s)
56	c Tconv---output-R-environment temperature seen by convective scheme (K)
57	c Cape----output-R-CAPE (J/kg)
58	c Cin ----output-R-CIN (J/kg)
59	c Tvp-----output-R-Temperature virtuelle d'une parcelle soulevee
60	c adiabatiquement a partir du niveau 1 (K)
61	c deltapb-output-R-distance entre LCL et base de la colonne (<0 ; Pa)
62	c Ice_flag-input-L-TRUE->prise en compte de la thermodynamique de la glace
63	c dd_t-----output-R-increment de la temperature du aux descentes precipitantes
64	c dd_q-----output-R-increment de la vapeur d'eau du aux desc precip
65	c======================================================================
66	c
67	#include "dimensions.h"
68	#include "dimphy.h"
69	c
70	integer NTRAC
71	PARAMETER (NTRAC=nqmx-2)
72	c
73	INTEGER iflag_con,iflag_clos
74	c
75	REAL dtime, paprs(klon,klev+1),pplay(klon,klev)
76	REAL t(klon,klev),q(klon,klev),u(klon,klev),v(klon,klev)
77	REAL t_wake(klon,klev),q_wake(klon,klev)
78	REAL tra(klon,klev,ntrac)
79	INTEGER ntra
80	REAL work1(klon,klev),work2(klon,klev),ptop2(klon)
81	REAL pmflxr(klon,klev+1),pmflxs(klon,klev+1)
82	REAL ALE(klon),ALP(klon)
83	c
84	REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
85	REAL dd_t(klon,klev),dd_q(klon,klev)
86	REAL d_tra(klon,klev,ntrac)
87	REAL rain(klon),snow(klon)
88	c
89	INTEGER kbas(klon),ktop(klon)
90	REAL em_ph(klon,klev+1),em_p(klon,klev)
91	REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev)
92	REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev)
93	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
94	REAL cape(klon),cin(klon),tvp(klon,klev)
95	REAL Tconv(klon,klev)
96	c
97	cCR:test: on passe lentr et alim_star des thermiques
98	INTEGER lalim_conv(klon)
99	REAL wght_th(klon,klev)
100	REAL em_sig1feed ! sigma at lower bound of feeding layer
101	REAL em_sig2feed ! sigma at upper bound of feeding layer
102	REAL em_wght(klev) ! weight density determining the feeding mixture
103	con enleve le save
104	c SAVE em_sig1feed,em_sig2feed,em_wght
105	c
106	INTEGER iflag(klon)
107	REAL rflag(klon)
108	REAL pbase(klon),bbase(klon)
109	REAL dtvpdt1(klon,klev),dtvpdq1(klon,klev)
110	REAL dplcldt(klon),dplcldr(klon)
111	REAL qcondc(klon,klev)
112	REAL wd(klon)
113	REAL Plim1(klon),Plim2(klon),asupmax(klon,klev)
114	REAL supmax0(klon),asupmaxmin(klon)
115	c
116	REAL sigd(klon)
117	REAL zx_t,zdelta,zx_qs,zcor
118	c
119	INTEGER iflag_mix
120	SAVE iflag_mix
121	INTEGER noff, minorig
122	INTEGER i,k,itra
123	REAL qs(klon,klev),qs_wake(klon,klev)
124	REAL cbmf(klon)
125	SAVE 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	Real ql(klon,klev)
136	c Save paire
137	Save ql
138	Real t1(klon,klev),q1(klon,klev)
139	Save t1,q1
140	c Data paire /1./
141	c
142	C Variables liees au bilan d'energie et d'enthalpi
143	REAL ztsol(klon)
144	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
145	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
146	SAVE 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	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
149	REAL d_h_vcol_phy
150	REAL fs_bound, fq_bound
151	SAVE d_h_vcol_phy
152	REAL zero_v(klon)
153	CHARACTER*15 ztit
154	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
155	SAVE ip_ebil
156	DATA ip_ebil/2/
157	INTEGER if_ebil ! level for energy conserv. dignostics
158	SAVE if_ebil
159	DATA if_ebil/2/
160	c+jld ec_conser
161	REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique
162	REAL ZRCPD
163	c-jld ec_conser
164	c
165	#include "YOMCST.h"
166	#include "YOMCST2.h"
167	#include "YOETHF.h"
168	#include "FCTTRE.h"
169	c
170
171	c Copy T into Tconv
172	DO k = 1,klev
173	DO i = 1,klon
174	Tconv(i,k) = T(i,k)
175	ENDDO
176	ENDDO
177	c
178	IF (if_ebil.ge.1) THEN
179	DO i=1,klon
180	ztsol(i) = t(i,1)
181	zero_v(i)=0.
182	Do k = 1,klev
183	ql(i,k) = 0.
184	ENDDO
185	END DO
186	END IF
187	c
188	cym
189	snow(:)=0
190
191	IF (ifrst .EQ. 0) THEN
192	ifrst = 1
193	c
194	C===========================================================================
195	C READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
196	C===========================================================================
197	C
198	if (iflag_con.eq.3) then
199	c CALL cv3_inicp(iflag_clos,iflag_mix)
200	CALL cv3_inip(iflag_mix)
201	endif
202	c
203	C===========================================================================
204	C READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
205	C===========================================================================
206	C
207	open (56,file='supcrit.data')
208	read (56,*) Supcrit1, Supcrit2
209	close (56)
210	c
211	print*, 'supcrit1, supcrit2' ,supcrit1, supcrit2
212	C
213	C===========================================================================
214	C Initialisation pour les bilans d'eau et d'energie
215	C===========================================================================
216	IF (if_ebil.ge.1) d_h_vcol_phy=0.
217	c
218	DO i = 1, klon
219	cbmf(i) = 0.
220	ENDDO
221	ENDIF !(ifrst .EQ. 0)
222
223	DO k = 1, klev+1
224	DO i=1,klon
225	em_ph(i,k) = paprs(i,k) / 100.0
226	pmflxs(i,k)=0.
227	ENDDO
228	ENDDO
229	c
230	DO k = 1, klev
231	DO i=1,klon
232	em_p(i,k) = pplay(i,k) / 100.0
233	ENDDO
234	ENDDO
235	c
236	!
237	! Feeding layer
238	!
239	em_sig1feed = 1.
240	em_sig2feed = 0.97
241	c em_sig2feed = 0.8
242	! Relative Weight densities
243	do k=1,klev
244	em_wght(k)=1.
245	end do
246	cCRtest: couche alim des tehrmiques ponderee par a*
247	c DO i = 1, klon
248	c do k=1,lalim_conv(i)
249	c em_wght(k)=wght_th(i,k)
250	c print*,'em_wght=',em_wght(k),wght_th(i,k)
251	c end do
252	c END DO
253
254	if (iflag_con .eq. 4) then
255	DO k = 1, klev
256	DO i = 1, klon
257	zx_t = t(i,k)
258	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
259	zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
260	zcor=1./(1.-retv*zx_qs)
261	qs(i,k)=zx_qs*zcor
262	ENDDO
263	DO i = 1, klon
264	zx_t = t_wake(i,k)
265	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
266	zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
267	zcor=1./(1.-retv*zx_qs)
268	qs_wake(i,k)=zx_qs*zcor
269	ENDDO
270	ENDDO
271	else ! iflag_con=3 (modif de puristes qui fait la diffce pour la convergence numerique)
272	DO k = 1, klev
273	DO i = 1, klon
274	zx_t = t(i,k)
275	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
276	zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
277	zx_qs= MIN(0.5,zx_qs)
278	zcor=1./(1.-retv*zx_qs)
279	zx_qs=zx_qs*zcor
280	qs(i,k)=zx_qs
281	ENDDO
282	DO i = 1, klon
283	zx_t = t_wake(i,k)
284	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
285	zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
286	zx_qs= MIN(0.5,zx_qs)
287	zcor=1./(1.-retv*zx_qs)
288	zx_qs=zx_qs*zcor
289	qs_wake(i,k)=zx_qs
290	ENDDO
291	ENDDO
292	endif ! iflag_con
293	c
294	C------------------------------------------------------------------
295
296	C Main driver for convection:
297	C iflag_con=3 -> nvlle version de KE (JYG)
298	C iflag_con = 30 -> equivalent to convect3
299	C iflag_con = 4 -> equivalent to convect1/2
300	c
301	c
302	if (iflag_con.eq.30) then
303
304	CALL cv_driver(klon,klev,klev+1,ntra,iflag_con,
305	: t,q,qs,u,v,tra,
306	$ em_p,em_ph,iflag,
307	$ d_t,d_q,d_u,d_v,d_tra,rain,
308	$ pmflxr,cbmf,work1,work2,
309	$ kbas,ktop,
310	$ dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
311	$ da,phi,mp)
312
313	else
314
315	CALL cva_driver(klon,klev,klev+1,ntra,
316	$ iflag_con,iflag_mix,iflag_clos,dtime,
317	: t,q,qs,t_wake,q_wake,qs_wake,u,v,tra,
318	$ em_p,em_ph,
319	. ALE,ALP,
320	. em_sig1feed,em_sig2feed,em_wght,
321	. iflag,d_t,d_q,d_u,d_v,d_tra,rain,kbas,ktop,
322	$ cbmf,work1,work2,ptop2,sigd,
323	$ Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
324	$ cape,cin,tvp,
325	$ dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
326	$ asupmaxmin,lalim_conv)
327	endif
328	C------------------------------------------------------------------
329
330	DO i = 1,klon
331	rain(i) = rain(i)/86400.
332	rflag(i)=iflag(i)
333	ENDDO
334
335	DO k = 1, klev
336	DO i = 1, klon
337	d_t(i,k) = dtime*d_t(i,k)
338	d_q(i,k) = dtime*d_q(i,k)
339	d_u(i,k) = dtime*d_u(i,k)
340	d_v(i,k) = dtime*d_v(i,k)
341	ENDDO
342	ENDDO
343	c
344	if (iflag_con.eq.30) then
345	DO itra = 1,ntra
346	DO k = 1, klev
347	DO i = 1, klon
348	d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
349	ENDDO
350	ENDDO
351	ENDDO
352	endif
353
354	DO k = 1, klev
355	DO i = 1, klon
356	t1(i,k) = t(i,k)+ d_t(i,k)
357	q1(i,k) = q(i,k)+ d_q(i,k)
358	ENDDO
359	ENDDO
360	c
361	cc IF (if_ebil.ge.2) THEN
362	cc ztit='after convect'
363	cc CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
364	cc e , t1,q1,ql,qs,u,v,paprs,pplay
365	cc s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
366	cc call diagphy(paire,ztit,ip_ebil
367	cc e , zero_v, zero_v, zero_v, zero_v, zero_v
368	cc e , zero_v, rain, zero_v, ztsol
369	cc e , d_h_vcol, d_qt, d_ec
370	cc s , fs_bound, fq_bound )
371	cc END IF
372	C
373	c
374	c les traceurs ne sont pas mis dans cette version de convect4:
375	if (iflag_con.eq.4) then
376	DO itra = 1,ntra
377	DO k = 1, klev
378	DO i = 1, klon
379	d_tra(i,k,itra) = 0.
380	ENDDO
381	ENDDO
382	ENDDO
383	endif
384	print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
385
386	RETURN
387	END
388

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

Download in other formats:

Original Format