Context Navigation

thermcell_dq.F90 @ 5440

Last change on this file since 5440 was 1985, checked in by fhourdin, 11 years ago

Modification de du schéma implicite des thermique pour que la concentration

soit inchangée en cas de flux de masse uniformément nul.

Implicit scheme for thermal plume modified to keep unchanged the concentration

when the mass flux is zero.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 8.6 KB

Line
1	subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr, &
2	& masse,q,dq,qa,lev_out)
3	implicit none
4
5	#include "iniprint.h"
6	!=======================================================================
7	!
8	! Calcul du transport verticale dans la couche limite en presence
9	! de "thermiques" explicitement representes
10	! calcul du dq/dt une fois qu'on connait les ascendances
11	!
12	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
13	! Introduction of an implicit computation of vertical advection in
14	! the environment of thermal plumes in thermcell_dq
15	! impl = 0 : explicit, 1 : implicit, -1 : old version
16	!
17	!=======================================================================
18
19	integer ngrid,nlay,impl
20
21	real ptimestep
22	real masse(ngrid,nlay),fm(ngrid,nlay+1)
23	real entr(ngrid,nlay)
24	real q(ngrid,nlay)
25	real dq(ngrid,nlay)
26	integer lev_out ! niveau pour les print
27
28	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
29
30	real zzm
31
32	integer ig,k
33	real cfl
34
35	real qold(ngrid,nlay),fqa(ngrid,nlay+1)
36	integer niter,iter
37	CHARACTER (LEN=20) :: modname='thermcell_dq'
38	CHARACTER (LEN=80) :: abort_message
39
40
41	! Old explicite scheme
42	if (impl<=-1) then
43	call thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, &
44	& masse,q,dq,qa,lev_out)
45	return
46	endif
47
48	! Calcul du critere CFL pour l'advection dans la subsidence
49	cfl = 0.
50	do k=1,nlay
51	do ig=1,ngrid
52	zzm=masse(ig,k)/ptimestep
53	cfl=max(cfl,fm(ig,k)/zzm)
54	if (entr(ig,k).gt.zzm) then
55	print,'entrdt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k)
56	abort_message = 'entr dt > m, 1st'
57	CALL abort_gcm (modname,abort_message,1)
58	endif
59	enddo
60	enddo
61
62	qold=q
63
64
65	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
66
67	! calcul du detrainement
68	do k=1,nlay
69	do ig=1,ngrid
70	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
71	! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
72	!test
73	if (detr(ig,k).lt.0.) then
74	entr(ig,k)=entr(ig,k)-detr(ig,k)
75	detr(ig,k)=0.
76	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
77	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
78	endif
79	if (fm(ig,k+1).lt.0.) then
80	! print*,'fm2<0!!!'
81	endif
82	if (entr(ig,k).lt.0.) then
83	! print*,'entr2<0!!!'
84	endif
85	enddo
86	enddo
87
88	! Computation of tracer concentrations in the ascending plume
89	do ig=1,ngrid
90	qa(ig,1)=q(ig,1)
91	enddo
92
93	do k=2,nlay
94	do ig=1,ngrid
95	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
96	& 1.e-5*masse(ig,k)) then
97	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
98	& /(fm(ig,k+1)+detr(ig,k))
99	else
100	qa(ig,k)=q(ig,k)
101	endif
102	if (qa(ig,k).lt.0.) then
103	! print*,'qa<0!!!'
104	endif
105	if (q(ig,k).lt.0.) then
106	! print*,'q<0!!!'
107	endif
108	enddo
109	enddo
110
111	! Plume vertical flux
112	do k=2,nlay-1
113	fqa(:,k)=fm(:,k)*qa(:,k-1)
114	enddo
115	fqa(:,1)=0. ; fqa(:,nlay)=0.
116
117
118	! Trace species evolution
119	if (impl==0) then
120	do k=1,nlay-1
121	q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)q(:,k)+fm(:,k+1)q(:,k+1)) &
122	& *ptimestep/masse(:,k)
123	enddo
124	else
125	do k=nlay-1,1,-1
126	! FH debut de modif : le calcul ci dessous modifiait numériquement
127	! la concentration quand le flux de masse etait nul car on divisait
128	! puis multipliait par masse/ptimestep.
129	! q(:,k)=(masse(:,k)q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)q(:,k+1)) &
130	! & /(fm(:,k)+masse(:,k)/ptimestep)
131	q(:,k)=(q(:,k)+ptimestep/masse(:,k)(fqa(:,k)-fqa(:,k+1)+fm(:,k+1)q(:,k+1))) &
132	& /(1.+fm(:,k)*ptimestep/masse(:,k))
133	! FH fin de modif.
134	enddo
135	endif
136
137	! Tendencies
138	do k=1,nlay
139	do ig=1,ngrid
140	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
141	q(ig,k)=qold(ig,k)
142	enddo
143	enddo
144
145	return
146	end
147
148	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
149	! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
150	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
151
152	subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, &
153	& masse,q,dq,qa,lev_out)
154	implicit none
155
156	#include "iniprint.h"
157	!=======================================================================
158	!
159	! Calcul du transport verticale dans la couche limite en presence
160	! de "thermiques" explicitement representes
161	! calcul du dq/dt une fois qu'on connait les ascendances
162	!
163	!=======================================================================
164
165	integer ngrid,nlay,impl
166
167	real ptimestep
168	real masse(ngrid,nlay),fm(ngrid,nlay+1)
169	real entr(ngrid,nlay)
170	real q(ngrid,nlay)
171	real dq(ngrid,nlay)
172	integer lev_out ! niveau pour les print
173
174	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
175
176	real zzm
177
178	integer ig,k
179	real cfl
180
181	real qold(ngrid,nlay)
182	real ztimestep
183	integer niter,iter
184	CHARACTER (LEN=20) :: modname='thermcell_dq'
185	CHARACTER (LEN=80) :: abort_message
186
187
188
189	! Calcul du critere CFL pour l'advection dans la subsidence
190	cfl = 0.
191	do k=1,nlay
192	do ig=1,ngrid
193	zzm=masse(ig,k)/ptimestep
194	cfl=max(cfl,fm(ig,k)/zzm)
195	if (entr(ig,k).gt.zzm) then
196	print,'entrdt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
197	abort_message = 'entr dt > m, 2nd'
198	CALL abort_gcm (modname,abort_message,1)
199	endif
200	enddo
201	enddo
202
203	!IM 090508 print*,'CFL CFL CFL CFL ',cfl
204
205	#undef CFL
206	#ifdef CFL
207	! On subdivise le calcul en niter pas de temps.
208	niter=int(cfl)+1
209	#else
210	niter=1
211	#endif
212
213	ztimestep=ptimestep/niter
214	qold=q
215
216
217	do iter=1,niter
218	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
219
220	! calcul du detrainement
221	do k=1,nlay
222	do ig=1,ngrid
223	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
224	! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
225	!test
226	if (detr(ig,k).lt.0.) then
227	entr(ig,k)=entr(ig,k)-detr(ig,k)
228	detr(ig,k)=0.
229	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
230	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
231	endif
232	if (fm(ig,k+1).lt.0.) then
233	! print*,'fm2<0!!!'
234	endif
235	if (entr(ig,k).lt.0.) then
236	! print*,'entr2<0!!!'
237	endif
238	enddo
239	enddo
240
241	! calcul de la valeur dans les ascendances
242	do ig=1,ngrid
243	qa(ig,1)=q(ig,1)
244	enddo
245
246	do k=2,nlay
247	do ig=1,ngrid
248	if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. &
249	& 1.e-5*masse(ig,k)) then
250	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
251	& /(fm(ig,k+1)+detr(ig,k))
252	else
253	qa(ig,k)=q(ig,k)
254	endif
255	if (qa(ig,k).lt.0.) then
256	! print*,'qa<0!!!'
257	endif
258	if (q(ig,k).lt.0.) then
259	! print*,'q<0!!!'
260	endif
261	enddo
262	enddo
263
264	! Calcul du flux subsident
265
266	do k=2,nlay
267	do ig=1,ngrid
268	#undef centre
269	#ifdef centre
270	wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
271	#else
272
273	#define plusqueun
274	#ifdef plusqueun
275	! Schema avec advection sur plus qu'une maille.
276	zzm=masse(ig,k)/ztimestep
277	if (fm(ig,k)>zzm) then
278	wqd(ig,k)=zzmq(ig,k)+(fm(ig,k)-zzm)q(ig,k+1)
279	else
280	wqd(ig,k)=fm(ig,k)*q(ig,k)
281	endif
282	#else
283	wqd(ig,k)=fm(ig,k)*q(ig,k)
284	#endif
285	#endif
286
287	if (wqd(ig,k).lt.0.) then
288	! print*,'wqd<0!!!'
289	endif
290	enddo
291	enddo
292	do ig=1,ngrid
293	wqd(ig,1)=0.
294	wqd(ig,nlay+1)=0.
295	enddo
296
297
298	! Calcul des tendances
299	do k=1,nlay
300	do ig=1,ngrid
301	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
302	& -wqd(ig,k)+wqd(ig,k+1)) &
303	& *ztimestep/masse(ig,k)
304	! if (dq(ig,k).lt.0.) then
305	! print*,'dq<0!!!'
306	! endif
307	enddo
308	enddo
309
310
311	enddo
312
313
314	! Calcul des tendances
315	do k=1,nlay
316	do ig=1,ngrid
317	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
318	q(ig,k)=qold(ig,k)
319	enddo
320	enddo
321
322	return
323	end

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

Download in other formats:

Original Format