Context Navigation

conema3.F @ 1042

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

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

Download in other formats:

Original Format