Context Navigation

thermcell_dq.F90 @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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.2 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,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,'entr dt > m ',entr(ig,k)ptimestep,masse(ig,k)
56	abort_message = ''
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	q(:,k)=(masse(:,k)q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)q(:,k+1)) &
127	& /(fm(:,k)+masse(:,k)/ptimestep)
128	enddo
129	endif
130
131	! Tendencies
132	do k=1,nlay
133	do ig=1,ngrid
134	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
135	q(ig,k)=qold(ig,k)
136	enddo
137	enddo
138
139	return
140	end
141
142	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
143	! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
144	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
145
146	subroutine thermcell_dq_o(ngrid,nlay,ptimestep,fm,entr, &
147	& masse,q,dq,qa,lev_out)
148	implicit none
149
150	#include "iniprint.h"
151	!=======================================================================
152	!
153	! Calcul du transport verticale dans la couche limite en presence
154	! de "thermiques" explicitement representes
155	! calcul du dq/dt une fois qu'on connait les ascendances
156	!
157	!=======================================================================
158
159	integer ngrid,nlay
160
161	real ptimestep
162	real masse(ngrid,nlay),fm(ngrid,nlay+1)
163	real entr(ngrid,nlay)
164	real q(ngrid,nlay)
165	real dq(ngrid,nlay)
166	integer lev_out ! niveau pour les print
167
168	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
169
170	real zzm
171
172	integer ig,k
173	real cfl
174
175	real qold(ngrid,nlay)
176	real ztimestep
177	integer niter,iter
178	CHARACTER (LEN=20) :: modname='thermcell_dq'
179	CHARACTER (LEN=80) :: abort_message
180
181
182
183	! Calcul du critere CFL pour l'advection dans la subsidence
184	cfl = 0.
185	do k=1,nlay
186	do ig=1,ngrid
187	zzm=masse(ig,k)/ptimestep
188	cfl=max(cfl,fm(ig,k)/zzm)
189	if (entr(ig,k).gt.zzm) then
190	print,'entr dt > m ',entr(ig,k)ptimestep,masse(ig,k)
191	abort_message = ''
192	CALL abort_gcm (modname,abort_message,1)
193	endif
194	enddo
195	enddo
196
197	!IM 090508 print*,'CFL CFL CFL CFL ',cfl
198
199	#undef CFL
200	#ifdef CFL
201	! On subdivise le calcul en niter pas de temps.
202	niter=int(cfl)+1
203	#else
204	niter=1
205	#endif
206
207	ztimestep=ptimestep/niter
208	qold=q
209
210
211	do iter=1,niter
212	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
213
214	! calcul du detrainement
215	do k=1,nlay
216	do ig=1,ngrid
217	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
218	! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
219	!test
220	if (detr(ig,k).lt.0.) then
221	entr(ig,k)=entr(ig,k)-detr(ig,k)
222	detr(ig,k)=0.
223	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
224	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
225	endif
226	if (fm(ig,k+1).lt.0.) then
227	! print*,'fm2<0!!!'
228	endif
229	if (entr(ig,k).lt.0.) then
230	! print*,'entr2<0!!!'
231	endif
232	enddo
233	enddo
234
235	! calcul de la valeur dans les ascendances
236	do ig=1,ngrid
237	qa(ig,1)=q(ig,1)
238	enddo
239
240	do k=2,nlay
241	do ig=1,ngrid
242	if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. &
243	& 1.e-5*masse(ig,k)) then
244	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
245	& /(fm(ig,k+1)+detr(ig,k))
246	else
247	qa(ig,k)=q(ig,k)
248	endif
249	if (qa(ig,k).lt.0.) then
250	! print*,'qa<0!!!'
251	endif
252	if (q(ig,k).lt.0.) then
253	! print*,'q<0!!!'
254	endif
255	enddo
256	enddo
257
258	! Calcul du flux subsident
259
260	do k=2,nlay
261	do ig=1,ngrid
262	#undef centre
263	#ifdef centre
264	wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
265	#else
266
267	#define plusqueun
268	#ifdef plusqueun
269	! Schema avec advection sur plus qu'une maille.
270	zzm=masse(ig,k)/ztimestep
271	if (fm(ig,k)>zzm) then
272	wqd(ig,k)=zzmq(ig,k)+(fm(ig,k)-zzm)q(ig,k+1)
273	else
274	wqd(ig,k)=fm(ig,k)*q(ig,k)
275	endif
276	#else
277	wqd(ig,k)=fm(ig,k)*q(ig,k)
278	#endif
279	#endif
280
281	if (wqd(ig,k).lt.0.) then
282	! print*,'wqd<0!!!'
283	endif
284	enddo
285	enddo
286	do ig=1,ngrid
287	wqd(ig,1)=0.
288	wqd(ig,nlay+1)=0.
289	enddo
290
291
292	! Calcul des tendances
293	do k=1,nlay
294	do ig=1,ngrid
295	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
296	& -wqd(ig,k)+wqd(ig,k+1)) &
297	& *ztimestep/masse(ig,k)
298	! if (dq(ig,k).lt.0.) then
299	! print*,'dq<0!!!'
300	! endif
301	enddo
302	enddo
303
304
305	enddo
306
307
308	! Calcul des tendances
309	do k=1,nlay
310	do ig=1,ngrid
311	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
312	q(ig,k)=qold(ig,k)
313	enddo
314	enddo
315
316	return
317	end

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

Download in other formats:

Original Format