Context Navigation

thermcell_dq.F90 @ 1939

Last change on this file since 1939 was 142, checked in by lfita, 10 years ago

Incorporing all the work done in LMDZ_WRFmeas about:

1.- Chekcing for NaNs?
2.- Instabilities in thermcell (negative SQRT & robustness IF)

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

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

Download in other formats:

Original Format