Context Navigation

thermcell_dq.F90 @ 103

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

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

Download in other formats:

Original Format