Context Navigation

conema3.F @ 5169

Last change on this file since 5169 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.0 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE conema3 (dtime,paprs,pplay,t,q,u,v,tra,ntra,
5	. work1,work2,d_t,d_q,d_u,d_v,d_tra,
6	. rain, snow, kbas, ktop,
7	. upwd,dnwd,dnwdbis,bas,top,Ma,cape,tvp,rflag,
8	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
9	. qcond_incld)
10
11	USE dimphy
12	USE infotrac, ONLY : nbtr
13	IMPLICIT none
14	c======================================================================
15	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
16	c Objet: schema de convection de Emanuel (1991) interface
17	c Mai 1998: Interface modifiee pour implementation dans LMDZ
18	c======================================================================
19	c Arguments:
20	c dtime---input-R-pas d'integration (s)
21	c paprs---input-R-pression inter-couches (Pa)
22	c pplay---input-R-pression au milieu des couches (Pa)
23	c t-------input-R-temperature (K)
24	c q-------input-R-humidite specifique (kg/kg)
25	c u-------input-R-vitesse du vent zonal (m/s)
26	c v-------input-R-vitesse duvent meridien (m/s)
27	c tra-----input-R-tableau de rapport de melange des traceurs
28	c work*: input et output: deux variables de travail,
29	c on peut les mettre a 0 au debut
30	c
31	C d_t-----output-R-increment de la temperature
32	c d_q-----output-R-increment de la vapeur d'eau
33	c d_u-----output-R-increment de la vitesse zonale
34	c d_v-----output-R-increment de la vitesse meridienne
35	c d_tra---output-R-increment du contenu en traceurs
36	c rain----output-R-la pluie (mm/s)
37	c snow----output-R-la neige (mm/s)
38	c kbas----output-R-bas du nuage (integer)
39	c ktop----output-R-haut du nuage (integer)
40	c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
41	c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
42	c dnwdbis-output-R-unsaturated downdraft mass flux (kg/m**2/s)
43	c bas-----output-R-bas du nuage (real)
44	c top-----output-R-haut du nuage (real)
45	c Ma------output-R-flux ascendant non dilue (kg/m**2/s)
46	c cape----output-R-CAPE
47	c tvp-----output-R-virtual temperature of the lifted parcel
48	c rflag---output-R-flag sur le fonctionnement de convect
49	c pbase---output-R-pression a la base du nuage (Pa)
50	c bbase---output-R-buoyancy a la base du nuage (K)
51	c dtvpdt1-output-R-derivative of parcel virtual temp wrt T1
52	c dtvpdq1-output-R-derivative of parcel virtual temp wrt Q1
53	c dplcldt-output-R-derivative of the PCP pressure wrt T1
54	c dplcldr-output-R-derivative of the PCP pressure wrt Q1
55	c======================================================================
56	c
57	#include "dimensions.h"
58	#include "conema3.h"
59	INTEGER i, l,m,itra
60	INTEGER ntra ! if no tracer transport
61	! is needed, set ntra = 1 (or 0)
62	REAL dtime
63	c
64	REAL d_t2(klon,klev), d_q2(klon,klev) ! sbl
65	REAL d_u2(klon,klev), d_v2(klon,klev) ! sbl
66	REAL em_d_t2(klev), em_d_q2(klev) ! sbl
67	REAL em_d_u2(klev), em_d_v2(klev) ! sbl
68	c
69	REAL paprs(klon,klev+1), pplay(klon,klev)
70	REAL t(klon,klev), q(klon,klev), d_t(klon,klev), d_q(klon,klev)
71	REAL u(klon,klev), v(klon,klev), tra(klon,klev,ntra)
72	REAL d_u(klon,klev), d_v(klon,klev), d_tra(klon,klev,ntra)
73	REAL work1(klon,klev), work2(klon,klev)
74	REAL upwd(klon,klev), dnwd(klon,klev), dnwdbis(klon,klev)
75	REAL rain(klon)
76	REAL snow(klon)
77	REAL cape(klon), tvp(klon,klev), rflag(klon)
78	REAL pbase(klon), bbase(klon)
79	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
80	REAL dplcldt(klon), dplcldr(klon)
81	INTEGER kbas(klon), ktop(klon)
82
83	REAL wd(klon)
84	REAL qcond_incld(klon,klev)
85	c
86	LOGICAL,SAVE :: first=.true.
87	c$OMP THREADPRIVATE(first)
88
89	cym REAL em_t(klev)
90	REAL,ALLOCATABLE,SAVE :: em_t(:)
91	c$OMP THREADPRIVATE(em_t)
92	cym REAL em_q(klev)
93	REAL,ALLOCATABLE,SAVE :: em_q(:)
94	c$OMP THREADPRIVATE(em_q)
95	cym REAL em_qs(klev)
96	REAL,ALLOCATABLE,SAVE :: em_qs(:)
97	c$OMP THREADPRIVATE(em_qs)
98	cym REAL em_u(klev), em_v(klev), em_tra(klev,nbtr)
99	REAL,ALLOCATABLE,SAVE :: em_u(:),em_v(:),em_tra(:,:)
100	c$OMP THREADPRIVATE(em_u,em_v,em_tra)
101	cym REAL em_ph(klev+1), em_p(klev)
102	REAL,ALLOCATABLE,SAVE ::em_ph(:),em_p(:)
103	c$OMP THREADPRIVATE(em_ph,em_p)
104	cym REAL em_work1(klev), em_work2(klev)
105	REAL,ALLOCATABLE,SAVE ::em_work1(:),em_work2(:)
106	c$OMP THREADPRIVATE(em_work1,em_work2)
107	cym REAL em_precip, em_d_t(klev), em_d_q(klev)
108	REAL,SAVE :: em_precip
109	c$OMP THREADPRIVATE(em_precip)
110	REAL,ALLOCATABLE,SAVE :: em_d_t(:),em_d_q(:)
111	c$OMP THREADPRIVATE(em_d_t,em_d_q)
112	cym REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr)
113	REAL,ALLOCATABLE,SAVE ::em_d_u(:),em_d_v(:),em_d_tra(:,:)
114	c$OMP THREADPRIVATE(em_d_u,em_d_v,em_d_tra)
115	cym REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
116	REAL,ALLOCATABLE,SAVE :: em_upwd(:),em_dnwd(:),em_dnwdbis(:)
117	c$OMP THREADPRIVATE(em_upwd,em_dnwd,em_dnwdbis)
118	REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
119	REAL em_dplcldt, em_dplcldr
120	cym SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
121	cym SAVE em_u,em_v, em_tra
122	cym SAVE em_d_u,em_d_v, em_d_tra
123	cym SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
124
125	INTEGER em_bas, em_top
126	SAVE em_bas, em_top
127	c$OMP THREADPRIVATE(em_bas,em_top)
128	REAL em_wd
129	REAL em_qcond(klev)
130	REAL em_qcondc(klev)
131	c
132	REAL zx_t, zx_qs, zdelta, zcor
133	INTEGER iflag
134	REAL sigsum
135	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
136	c VARIABLES A SORTIR
137	cccccccccccccccccccccccccccccccccccccccccccccccccc
138
139	cym REAL emmip(klev) !variation de flux ascnon dilue i et i+1
140	REAL,ALLOCATABLE,SAVE ::emmip(:)
141	c$OMP THREADPRIVATE(emmip)
142	cym SAVE emmip
143	cym real emMke(klev)
144	REAL,ALLOCATABLE,SAVE ::emMke(:)
145	c$OMP THREADPRIVATE(emMke)
146	cym save emMke
147	real top
148	real bas
149	cym real emMa(klev)
150	REAL,ALLOCATABLE,SAVE ::emMa(:)
151	c$OMP THREADPRIVATE(emMa)
152	cym save emMa
153	real Ma(klon,klev)
154	real Ment(klev,klev)
155	real Qent(klev,klev)
156	real TPS(klev),TLS(klev)
157	real SIJ(klev,klev)
158	real em_CAPE, em_TVP(klev)
159	real em_pbase, em_bbase
160	integer iw,j,k,ix,iy
161
162	c -- sb: pour schema nuages:
163
164	integer iflagcon
165	integer em_ifc(klev)
166
167	real em_pradj
168	real em_cldf(klev), em_cldq(klev)
169	real em_ftadj(klev), em_fradj(klev)
170
171	integer ifc(klon,klev)
172	real pradj(klon)
173	real cldf(klon,klev), cldq(klon,klev)
174	real ftadj(klon,klev), fqadj(klon,klev)
175
176	c sb --
177
178	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
179	c
180	#include "YOMCST.h"
181	#include "YOETHF.h"
182	#include "FCTTRE.h"
183
184	if (first) then
185
186	allocate(em_t(klev))
187	allocate(em_q(klev))
188	allocate(em_qs(klev))
189	allocate(em_u(klev), em_v(klev), em_tra(klev,nbtr))
190	allocate(em_ph(klev+1), em_p(klev))
191	allocate(em_work1(klev), em_work2(klev))
192	allocate(em_d_t(klev), em_d_q(klev))
193	allocate(em_d_u(klev), em_d_v(klev), em_d_tra(klev,nbtr))
194	allocate(em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev))
195	allocate(emmip(klev))
196	allocate(emMke(klev))
197	allocate(emMa(klev))
198
199	first=.false.
200	endif
201
202	qcond_incld(:,:) = 0.
203	c
204	c@$$ print*,'debut conema'
205
206	DO 999 i = 1, klon
207	DO l = 1, klev+1
208	em_ph(l) = paprs(i,l) / 100.0
209	ENDDO
210	c
211	DO l = 1, klev
212	em_p(l) = pplay(i,l) / 100.0
213	em_t(l) = t(i,l)
214	em_q(l) = q(i,l)
215	em_u(l) = u(i,l)
216	em_v(l) = v(i,l)
217	do itra = 1, ntra
218	em_tra(l,itra) = tra(i,l,itra)
219	enddo
220	c@$$ print*,'em_t',em_t
221	c@$$ print*,'em_q',em_q
222	c@$$ print*,'em_qs',em_qs
223	c@$$ print*,'em_u',em_u
224	c@$$ print*,'em_v',em_v
225	c@$$ print*,'em_tra',em_tra
226	c@$$ print*,'em_p',em_p
227
228
229	c
230	zx_t = em_t(l)
231	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
232	zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(l)/100.0
233	zx_qs=MIN(0.5,zx_qs)
234	c@$$ print*,'zx_qs',zx_qs
235	zcor=1./(1.-retv*zx_qs)
236	zx_qs=zx_qs*zcor
237	em_qs(l) = zx_qs
238	c@$$ print*,'em_qs',em_qs
239	c
240	em_work1(l) = work1(i,l)
241	em_work2(l) = work2(i,l)
242	emMke(l)=0
243	c emMa(l)=0
244	c Ma(i,l)=0
245
246	em_dtvpdt1(l) = 0.
247	em_dtvpdq1(l) = 0.
248	dtvpdt1(i,l) = 0.
249	dtvpdq1(i,l) = 0.
250	ENDDO
251	c
252	em_dplcldt = 0.
253	em_dplcldr = 0.
254	rain(i) = 0.0
255	snow(i) = 0.0
256	kbas(i) = 1
257	ktop(i) = 1
258	c ajout SB:
259	bas = 1
260	top = 1
261
262
263	c sb3d write(*,1792) (em_work1(m),m=1,klev)
264	1792 format('sig avant convect ',/,10(1X,E13.5))
265	c
266	c sb d write(*,1793) (em_work2(m),m=1,klev)
267	1793 format('w avant convect ',/,10(1X,E13.5))
268
269	c@$$ print*,'avant convect'
270	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
271	c
272
273	c print*,'avant convect i=',i
274	CALL convect3(dtime,epmax,ok_adj_ema,
275	. em_t, em_q, em_qs,em_u ,em_v ,
276	. em_tra, em_p, em_ph,
277	. klev, klev+1, klev-1,ntra, dtime, iflag,
278	. em_d_t, em_d_q,em_d_u,em_d_v,
279	. em_d_tra, em_precip,
280	. em_bas, em_top,em_upwd, em_dnwd, em_dnwdbis,
281	. em_work1, em_work2,emmip,emMke,emMa,Ment,
282	. Qent,TPS,TLS,SIJ,em_CAPE,em_TVP,em_pbase,em_bbase,
283	. em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr, ! sbl
284	. em_d_t2,em_d_q2,em_d_u2,em_d_v2,em_wd,em_qcond,em_qcondc)!sbl
285	c print*,'apres convect '
286	c
287	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
288	c
289	c -- sb: Appel schema statistique de nuages couple a la convection
290	c (Bony et Emanuel 2001):
291
292	c -- creer cvthermo.h qui contiendra les cstes thermo de LMDZ:
293
294	iflagcon = 3
295	c CALL cv_thermo(iflagcon)
296
297	c -- appel schema de nuages:
298
299	c CALL CLOUDS_SUB_LS(klev,em_q,em_qs,em_t
300	c i ,em_p,em_ph,dtime,em_qcondc
301	c o ,em_cldf,em_cldq,em_pradj,em_ftadj,em_fradj,em_ifc)
302
303	do k = 1, klev
304	cldf(i,k) = em_cldf(k) ! cloud fraction (0-1)
305	cldq(i,k) = em_cldq(k) ! in-cloud water content (kg/kg)
306	ftadj(i,k) = em_ftadj(k) ! (dT/dt)_{LS adj} (K/s)
307	fqadj(i,k) = em_fradj(k) ! (dq/dt)_{LS adj} (kg/kg/s)
308	ifc(i,k) = em_ifc(k) ! flag convergence clouds_gno (1 ou 2)
309	enddo
310	pradj(i) = em_pradj ! precip from LS supersat adj (mm/day)
311
312	c sb --
313	c
314	c SB:
315	if (iflag.ne.1 .and. iflag.ne.4) then
316	em_CAPE = 0.
317	do l = 1, klev
318	em_upwd(l) = 0.
319	em_dnwd(l) = 0.
320	em_dnwdbis(l) = 0.
321	emMa(l) = 0.
322	em_TVP(l) = 0.
323	enddo
324	endif
325	c fin SB
326	c
327	c If sig has been set to zero, then set Ma to zero
328	c
329	sigsum = 0.
330	do k = 1,klev
331	sigsum = sigsum + em_work1(k)
332	enddo
333	if (sigsum .eq. 0.0) then
334	do k = 1,klev
335	emMa(k) = 0.
336	enddo
337	endif
338	c
339	c sb3d print*,'i, iflag=',i,iflag
340	c
341	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
342	c
343	c SORTIE DES ICB ET INB
344	c en fait inb et icb correspondent au niveau ou se trouve
345	c le nuage,le numero d'interface
346	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
347
348	c modif SB:
349	if (iflag.EQ.1 .or. iflag.EQ.4) then
350	top=em_top
351	bas=em_bas
352	kbas(i) = em_bas
353	ktop(i) = em_top
354	endif
355
356	pbase(i) = em_pbase
357	bbase(i) = em_bbase
358	rain(i) = em_precip/ 86400.0
359	snow(i) = 0.0
360	cape(i) = em_CAPE
361	wd(i) = em_wd
362	rflag(i) = float(iflag)
363	c SB kbas(i) = em_bas
364	c SB ktop(i) = em_top
365	dplcldt(i) = em_dplcldt
366	dplcldr(i) = em_dplcldr
367	DO l = 1, klev
368	d_t2(i,l) = dtime * em_d_t2(l)
369	d_q2(i,l) = dtime * em_d_q2(l)
370	d_u2(i,l) = dtime * em_d_u2(l)
371	d_v2(i,l) = dtime * em_d_v2(l)
372
373	d_t(i,l) = dtime * em_d_t(l)
374	d_q(i,l) = dtime * em_d_q(l)
375	d_u(i,l) = dtime * em_d_u(l)
376	d_v(i,l) = dtime * em_d_v(l)
377	do itra = 1, ntra
378	d_tra(i,l,itra) = dtime * em_d_tra(l,itra)
379	enddo
380	upwd(i,l) = em_upwd(l)
381	dnwd(i,l) = em_dnwd(l)
382	dnwdbis(i,l) = em_dnwdbis(l)
383	work1(i,l) = em_work1(l)
384	work2(i,l) = em_work2(l)
385	Ma(i,l)=emMa(l)
386	tvp(i,l)=em_TVP(l)
387	dtvpdt1(i,l) = em_dtvpdt1(l)
388	dtvpdq1(i,l) = em_dtvpdq1(l)
389
390	if (iflag_clw.eq.0) then
391	qcond_incld(i,l) = em_qcondc(l)
392	else if (iflag_clw.eq.1) then
393	qcond_incld(i,l) = em_qcond(l)
394	endif
395	ENDDO
396	999 CONTINUE
397
398	c On calcule une eau liquide diagnostique en fonction de la
399	c precip.
400	if ( iflag_clw.eq.2 ) then
401	do l=1,klev
402	do i=1,klon
403	if (ktop(i)-kbas(i).gt.0.and.
404	s l.ge.kbas(i).and.l.le.ktop(i)) then
405	qcond_incld(i,l)=rain(i)*8.e4
406	c s *(pplay(i,l )-paprs(i,ktop(i)+1))
407	s /(pplay(i,kbas(i))-pplay(i,ktop(i)))
408	c s **2
409	else
410	qcond_incld(i,l)=0.
411	endif
412	enddo
413	print*,'l=',l,', qcond_incld=',qcond_incld(1,l)
414	enddo
415	endif
416
417
418	RETURN
419	END
420

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

Download in other formats:

Original Format