Context Navigation

conema3.F @ 2570

Last change on this file since 2570 was 634, checked in by Laurent Fairhead, 20 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 12.5 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	cym REAL em_t(klev)
89	REAL,ALLOCATABLE,SAVE :: em_t(:)
90	cym REAL em_q(klev)
91	REAL,ALLOCATABLE,SAVE :: em_q(:)
92	cym REAL em_qs(klev)
93	REAL,ALLOCATABLE,SAVE :: em_qs(:)
94	cym REAL em_u(klev), em_v(klev), em_tra(klev,ntrac)
95	REAL,ALLOCATABLE,SAVE :: em_u(:),em_v(:),em_tra(:,:)
96	cym REAL em_ph(klev+1), em_p(klev)
97	REAL,ALLOCATABLE,SAVE ::em_ph(:),em_p(:)
98	cym REAL em_work1(klev), em_work2(klev)
99	REAL,ALLOCATABLE,SAVE ::em_work1(:),em_work2(:)
100	cym REAL em_precip, em_d_t(klev), em_d_q(klev)
101	REAL,SAVE :: em_precip
102	REAL,ALLOCATABLE,SAVE :: em_d_t(:),em_d_q(:)
103	cym REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,ntrac)
104	REAL,ALLOCATABLE,SAVE ::em_d_u(:),em_d_v(:),em_d_tra(:,:)
105	cym REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
106	REAL,ALLOCATABLE,SAVE :: em_upwd(:),em_dnwd(:),em_dnwdbis(:)
107	REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
108	REAL em_dplcldt, em_dplcldr
109	cym SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
110	cym SAVE em_u,em_v, em_tra
111	cym SAVE em_d_u,em_d_v, em_d_tra
112	cym SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
113
114	INTEGER em_bas, em_top
115	SAVE em_bas, em_top
116
117	REAL em_wd
118	REAL em_qcond(klev)
119	REAL em_qcondc(klev)
120	c
121	REAL zx_t, zx_qs, zdelta, zcor
122	INTEGER iflag
123	REAL sigsum
124	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
125	c VARIABLES A SORTIR
126	cccccccccccccccccccccccccccccccccccccccccccccccccc
127
128	cym REAL emmip(klev) !variation de flux ascnon dilue i et i+1
129	REAL,ALLOCATABLE,SAVE ::emmip(:)
130	cym SAVE emmip
131	cym real emMke(klev)
132	REAL,ALLOCATABLE,SAVE ::emMke(:)
133	cym save emMke
134	real top
135	real bas
136	cym real emMa(klev)
137	REAL,ALLOCATABLE,SAVE ::emMa(:)
138	cym save emMa
139	real Ma(klon,klev)
140	real Ment(klev,klev)
141	real Qent(klev,klev)
142	real TPS(klev),TLS(klev)
143	real SIJ(klev,klev)
144	real em_CAPE, em_TVP(klev)
145	real em_pbase, em_bbase
146	integer iw,j,k,ix,iy
147
148	c -- sb: pour schema nuages:
149
150	integer iflagcon
151	integer em_ifc(klev)
152
153	real em_pradj
154	real em_cldf(klev), em_cldq(klev)
155	real em_ftadj(klev), em_fradj(klev)
156
157	integer ifc(klon,klev)
158	real pradj(klon)
159	real cldf(klon,klev), cldq(klon,klev)
160	real ftadj(klon,klev), fqadj(klon,klev)
161
162	c sb --
163
164	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
165	c
166	#include "YOMCST.h"
167	#include "YOETHF.h"
168	#include "FCTTRE.h"
169
170	if (first) then
171
172	allocate(em_t(klev))
173	allocate(em_q(klev))
174	allocate(em_qs(klev))
175	allocate(em_u(klev), em_v(klev), em_tra(klev,ntrac))
176	allocate(em_ph(klev+1), em_p(klev))
177	allocate(em_work1(klev), em_work2(klev))
178	allocate(em_d_t(klev), em_d_q(klev))
179	allocate(em_d_u(klev), em_d_v(klev), em_d_tra(klev,ntrac))
180	allocate(em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev))
181	allocate(emmip(klev))
182	allocate(emMke(klev))
183	allocate(emMa(klev))
184
185	first=.false.
186	endif
187
188	qcond_incld(:,:) = 0.
189	c
190	c$$$ print*,'debut conema'
191
192	DO 999 i = 1, klon
193	DO l = 1, klev+1
194	em_ph(l) = paprs(i,l) / 100.0
195	ENDDO
196	c
197	DO l = 1, klev
198	em_p(l) = pplay(i,l) / 100.0
199	em_t(l) = t(i,l)
200	em_q(l) = q(i,l)
201	em_u(l) = u(i,l)
202	em_v(l) = v(i,l)
203	do itra = 1, ntra
204	em_tra(l,itra) = tra(i,l,itra)
205	enddo
206	c$$$ print*,'em_t',em_t
207	c$$$ print*,'em_q',em_q
208	c$$$ print*,'em_qs',em_qs
209	c$$$ print*,'em_u',em_u
210	c$$$ print*,'em_v',em_v
211	c$$$ print*,'em_tra',em_tra
212	c$$$ print*,'em_p',em_p
213
214
215	c
216	zx_t = em_t(l)
217	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
218	zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(l)/100.0
219	zx_qs=MIN(0.5,zx_qs)
220	c$$$ print*,'zx_qs',zx_qs
221	zcor=1./(1.-retv*zx_qs)
222	zx_qs=zx_qs*zcor
223	em_qs(l) = zx_qs
224	c$$$ print*,'em_qs',em_qs
225	c
226	em_work1(l) = work1(i,l)
227	em_work2(l) = work2(i,l)
228	emMke(l)=0
229	c emMa(l)=0
230	c Ma(i,l)=0
231
232	em_dtvpdt1(l) = 0.
233	em_dtvpdq1(l) = 0.
234	dtvpdt1(i,l) = 0.
235	dtvpdq1(i,l) = 0.
236	ENDDO
237	c
238	em_dplcldt = 0.
239	em_dplcldr = 0.
240	rain(i) = 0.0
241	snow(i) = 0.0
242	kbas(i) = 1
243	ktop(i) = 1
244	c ajout SB:
245	bas = 1
246	top = 1
247
248
249	c sb3d write(*,1792) (em_work1(m),m=1,klev)
250	1792 format('sig avant convect ',/,10(1X,E13.5))
251	c
252	c sb d write(*,1793) (em_work2(m),m=1,klev)
253	1793 format('w avant convect ',/,10(1X,E13.5))
254
255	c$$$ print*,'avant convect'
256	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
257	c
258
259	c print*,'avant convect i=',i
260	CALL convect3(dtime,epmax,ok_adj_ema,
261	. em_t, em_q, em_qs,em_u ,em_v ,
262	. em_tra, em_p, em_ph,
263	. klev, klev+1, klev-1,ntra, dtime, iflag,
264	. em_d_t, em_d_q,em_d_u,em_d_v,
265	. em_d_tra, em_precip,
266	. em_bas, em_top,em_upwd, em_dnwd, em_dnwdbis,
267	. em_work1, em_work2,emmip,emMke,emMa,Ment,
268	. Qent,TPS,TLS,SIJ,em_CAPE,em_TVP,em_pbase,em_bbase,
269	. em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr, ! sbl
270	. em_d_t2,em_d_q2,em_d_u2,em_d_v2,em_wd,em_qcond,em_qcondc)!sbl
271	c print*,'apres convect '
272	c
273	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
274	c
275	c -- sb: Appel schema statistique de nuages couple a la convection
276	c (Bony et Emanuel 2001):
277
278	c -- creer cvthermo.h qui contiendra les cstes thermo de LMDZ:
279
280	iflagcon = 3
281	c CALL cv_thermo(iflagcon)
282
283	c -- appel schema de nuages:
284
285	c CALL CLOUDS_SUB_LS(klev,em_q,em_qs,em_t
286	c i ,em_p,em_ph,dtime,em_qcondc
287	c o ,em_cldf,em_cldq,em_pradj,em_ftadj,em_fradj,em_ifc)
288
289	do k = 1, klev
290	cldf(i,k) = em_cldf(k) ! cloud fraction (0-1)
291	cldq(i,k) = em_cldq(k) ! in-cloud water content (kg/kg)
292	ftadj(i,k) = em_ftadj(k) ! (dT/dt)_{LS adj} (K/s)
293	fqadj(i,k) = em_fradj(k) ! (dq/dt)_{LS adj} (kg/kg/s)
294	ifc(i,k) = em_ifc(k) ! flag convergence clouds_gno (1 ou 2)
295	enddo
296	pradj(i) = em_pradj ! precip from LS supersat adj (mm/day)
297
298	c sb --
299	c
300	c SB:
301	if (iflag.ne.1 .and. iflag.ne.4) then
302	em_CAPE = 0.
303	do l = 1, klev
304	em_upwd(l) = 0.
305	em_dnwd(l) = 0.
306	em_dnwdbis(l) = 0.
307	emMa(l) = 0.
308	em_TVP(l) = 0.
309	enddo
310	endif
311	c fin SB
312	c
313	c If sig has been set to zero, then set Ma to zero
314	c
315	sigsum = 0.
316	do k = 1,klev
317	sigsum = sigsum + em_work1(k)
318	enddo
319	if (sigsum .eq. 0.0) then
320	do k = 1,klev
321	emMa(k) = 0.
322	enddo
323	endif
324	c
325	c sb3d print*,'i, iflag=',i,iflag
326	c
327	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
328	c
329	c SORTIE DES ICB ET INB
330	c en fait inb et icb correspondent au niveau ou se trouve
331	c le nuage,le numero d'interface
332	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
333
334	c modif SB:
335	if (iflag.EQ.1 .or. iflag.EQ.4) then
336	top=em_top
337	bas=em_bas
338	kbas(i) = em_bas
339	ktop(i) = em_top
340	endif
341
342	pbase(i) = em_pbase
343	bbase(i) = em_bbase
344	rain(i) = em_precip/ 86400.0
345	snow(i) = 0.0
346	cape(i) = em_CAPE
347	wd(i) = em_wd
348	rflag(i) = float(iflag)
349	c SB kbas(i) = em_bas
350	c SB ktop(i) = em_top
351	dplcldt(i) = em_dplcldt
352	dplcldr(i) = em_dplcldr
353	DO l = 1, klev
354	d_t2(i,l) = dtime * em_d_t2(l)
355	d_q2(i,l) = dtime * em_d_q2(l)
356	d_u2(i,l) = dtime * em_d_u2(l)
357	d_v2(i,l) = dtime * em_d_v2(l)
358
359	d_t(i,l) = dtime * em_d_t(l)
360	d_q(i,l) = dtime * em_d_q(l)
361	d_u(i,l) = dtime * em_d_u(l)
362	d_v(i,l) = dtime * em_d_v(l)
363	do itra = 1, ntra
364	d_tra(i,l,itra) = dtime * em_d_tra(l,itra)
365	enddo
366	upwd(i,l) = em_upwd(l)
367	dnwd(i,l) = em_dnwd(l)
368	dnwdbis(i,l) = em_dnwdbis(l)
369	work1(i,l) = em_work1(l)
370	work2(i,l) = em_work2(l)
371	Ma(i,l)=emMa(l)
372	tvp(i,l)=em_TVP(l)
373	dtvpdt1(i,l) = em_dtvpdt1(l)
374	dtvpdq1(i,l) = em_dtvpdq1(l)
375
376	if (iflag_clw.eq.0) then
377	qcond_incld(i,l) = em_qcondc(l)
378	else if (iflag_clw.eq.1) then
379	qcond_incld(i,l) = em_qcond(l)
380	endif
381	ENDDO
382	999 CONTINUE
383
384	c On calcule une eau liquide diagnostique en fonction de la
385	c precip.
386	if ( iflag_clw.eq.2 ) then
387	do l=1,klev
388	do i=1,klon
389	if (ktop(i)-kbas(i).gt.0.and.
390	s l.ge.kbas(i).and.l.le.ktop(i)) then
391	qcond_incld(i,l)=rain(i)*8.e4
392	c s *(pplay(i,l )-paprs(i,ktop(i)+1))
393	s /(pplay(i,kbas(i))-pplay(i,ktop(i)))
394	c s **2
395	else
396	qcond_incld(i,l)=0.
397	endif
398	enddo
399	print*,'l=',l,', qcond_incld=',qcond_incld(1,l)
400	enddo
401	endif
402
403
404	RETURN
405	END
406

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

Download in other formats:

Original Format