Context Navigation

thermcell_dq.F90 @ 2102

Last change on this file since 2102 was 2102, checked in by aboissinot, 6 years ago
Now detr is an argument of thermcell_dq subroutine (which is called in thermcell_main). thermcell_dq is cleaned up.
File size: 11.0 KB

Line
1	!
2	!
3	!
4	SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,detr,masse, &
5	q,dq,qa,lmin,lev_out)
6
7
8	!==============================================================================
9	! Calcul du transport verticale dans la couche limite en presence
10	! de "thermiques" explicitement representes
11	! calcul du dq/dt une fois qu'on connait les ascendances
12	!
13	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
14	! Introduction of an implicit computation of vertical advection in
15	! the environment of thermal plumes in thermcell_dq
16	! impl = 0 : explicit, 1 : implicit, -1 : old version
17	!
18	! Modif 2018/11/05 (AB alexandre.boissinot@lmd.jussieu.fr)
19	!
20	!
21	!==============================================================================
22
23	USE print_control_mod, ONLY: prt_level
24
25	IMPLICIT NONE
26
27
28	!==============================================================================
29	! Declaration
30	!==============================================================================
31
32	! inputs:
33	! -------
34
35	INTEGER ngrid, nlay
36	INTEGER impl
37	INTEGER lmin(ngrid)
38	INTEGER lev_out ! niveau pour les print
39
40	REAL ptimestep
41	REAL masse(ngrid,nlay)
42	REAL fm(ngrid,nlay+1)
43	REAL entr(ngrid,nlay)
44	REAL detr(ngrid,nlay)
45	REAL q(ngrid,nlay)
46	REAL qa(ngrid,nlay)
47
48	! outputs:
49	! --------
50
51	REAL dq(ngrid,nlay)
52
53	! local:
54	! ------
55
56	INTEGER ig, l
57	INTEGER niter, iter
58
59	REAL cfl
60	REAL qold(ngrid,nlay)
61	REAL fqa(ngrid,nlay+1)
62	REAL wqd(ngrid,nlay+1)
63	REAL zzm
64
65	CHARACTER (LEN=20) :: modname
66	CHARACTER (LEN=80) :: abort_message
67
68	!==============================================================================
69	! Initialization
70	!==============================================================================
71
72	qold = q
73
74	modname = 'thermcell_dq'
75
76	IF (impl.lt.0) THEN
77	print *, 'OLD EXPLICIT SCHEME IS USED'
78	CALL thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, &
79	& masse,q,dq,qa,lev_out)
80	RETURN
81	ENDIF
82
83	!------------------------------------------------------------------------------
84	! CFL criterion computation for advection in downdraft
85	!------------------------------------------------------------------------------
86
87	cfl = 0.
88
89	DO l=1,nlay
90	DO ig=1,ngrid
91	zzm = masse(ig,l) / ptimestep
92	cfl = max(cfl, fm(ig,l) / zzm)
93
94	IF (entr(ig,l).gt.zzm) THEN
95	print *, 'ERROR: entrainment is greater than the layer mass!'
96	print *, 'ig,l,entr', ig, l, entr(ig,l)
97	print *, '-------------------------------'
98	print , 'entrdt,mass', entr(ig,l)*ptimestep, masse(ig,l)
99	print *, '-------------------------------'
100	print *, 'fm+', fm(ig,l+1)
101	print *, 'entr,detr', entr(ig,l), detr(ig,l)
102	print *, 'fm ', fm(ig,l)
103	print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
104	print *, 'fm-', fm(ig,l-1)
105	abort_message = 'entr dt > m, 1st'
106	CALL abort_physic(modname,abort_message,1)
107	ENDIF
108	ENDDO
109	ENDDO
110
111	!------------------------------------------------------------------------------
112	! Computation of tracer concentrations in the ascending plume
113	!------------------------------------------------------------------------------
114
115	DO ig=1,ngrid
116	DO l=1,lmin(ig)
117	qa(ig,l) = q(ig,l)
118	ENDDO
119	ENDDO
120
121	DO ig=1,ngrid
122	DO l=lmin(ig)+1,nlay
123	if ((fm(ig,l+1)+detr(ig,l))ptimestep.gt.1.e-5masse(ig,l)) then
124	qa(ig,l) = (fm(ig,l) * qa(ig,l-1) + entr(ig,l) * q(ig,l)) &
125	& / (fm(ig,l+1) + detr(ig,l))
126	else
127	qa(ig,l) = q(ig,l)
128	endif
129
130	! IF (qa(ig,l).lt.0.) THEN
131	! print *, 'WARNING: qa is negative!'
132	! print *, 'qa', qa(ig,l)
133	! ENDIF
134
135	! IF (q(ig,l).lt.0.) THEN
136	! print *, 'WARNING: q is negative!'
137	! print *, 'q', q(ig,l)
138	! ENDIF
139	ENDDO
140	ENDDO
141
142	!------------------------------------------------------------------------------
143	! Plume vertical flux of q
144	!------------------------------------------------------------------------------
145
146	DO l=2,nlay-1
147	fqa(:,l) = fm(:,l) * qa(:,l-1)
148	ENDDO
149
150	fqa(:,1) = 0.
151	fqa(:,nlay) = 0.
152
153	!------------------------------------------------------------------------------
154	! Trace species evolution
155	!------------------------------------------------------------------------------
156
157	IF (impl==0) THEN
158	do l=1,nlay-1
159	q(:,l) = q(:,l) + (fqa(:,l) - fqa(:,l+1) - fm(:,l) * q(:,l) &
160	& + fm(:,l+1) * q(:,l+1)) * ptimestep / masse(:,l)
161	enddo
162	ELSE
163	do l=nlay-1,1,-1
164	q(:,l) = ( q(:,l) + ptimestep / masse(:,l) &
165	& * ( fqa(:,l) - fqa(:,l+1) + fm(:,l+1) * q(:,l+1) ) ) &
166	& / ( 1. + fm(:,l) * ptimestep / masse(:,l) )
167	ENDDO
168	ENDIF
169
170	!==============================================================================
171	! Tendencies
172	!==============================================================================
173
174	DO l=1,nlay
175	DO ig=1,ngrid
176	dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep
177	q(ig,l) = qold(ig,l)
178	ENDDO
179	ENDDO
180
181
182	RETURN
183	END
184
185
186	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
187	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
189	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
191	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
192
193
194	Subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,masse, &
195	q,dq,qa,lev_out)
196
197
198	!==============================================================================
199	!
200	! Calcul du transport verticale dans la couche limite en presence
201	! de "thermiques" explicitement representes
202	! calcul du dq/dt une fois qu'on connait les ascendances
203	!
204	!==============================================================================
205
206	USE print_control_mod, ONLY: prt_level
207
208	implicit none
209
210
211	!==============================================================================
212	! Declaration
213	!==============================================================================
214
215	integer ngrid,nlay,impl
216
217	real ptimestep
218	real masse(ngrid,nlay),fm(ngrid,nlay+1)
219	real entr(ngrid,nlay)
220	real q(ngrid,nlay)
221	real dq(ngrid,nlay)
222	integer lev_out ! niveau pour les print
223
224	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
225
226	real zzm
227
228	integer ig,l
229	real cfl
230
231	real qold(ngrid,nlay)
232	real ztimestep
233	integer niter,iter
234	CHARACTER (LEN=20) :: modname='thermcell_dq'
235	CHARACTER (LEN=80) :: abort_message
236
237
238
239	! Calcul du critere CFL pour l'advection dans la subsidence
240	cfl = 0.
241	do l=1,nlay
242	do ig=1,ngrid
243	zzm=masse(ig,l)/ptimestep
244	cfl=max(cfl,fm(ig,l)/zzm)
245	if (entr(ig,l).gt.zzm) then
246	print,'entrdt>m,2',l,entr(ig,l)*ptimestep,masse(ig,l)
247	abort_message = 'entr dt > m, 2nd'
248	CALL abort_physic (modname,abort_message,1)
249	endif
250	enddo
251	enddo
252
253	!IM 090508 print*,'CFL CFL CFL CFL ',cfl
254
255	#undef CFL
256	#ifdef CFL
257	! On subdivise le calcul en niter pas de temps.
258	niter=int(cfl)+1
259	#else
260	niter=1
261	#endif
262
263	ztimestep=ptimestep/niter
264	qold=q
265
266	do iter=1,niter
267	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
268
269	! calcul du detrainement
270	do l=1,nlay
271	do ig=1,ngrid
272	detr(ig,l)=fm(ig,l)-fm(ig,l+1)+entr(ig,l)
273	! print*,'Q2 DQ ',detr(ig,l),fm(ig,l),entr(ig,l)
274	!test
275	if (detr(ig,l).lt.0.) then
276	entr(ig,l)=entr(ig,l)-detr(ig,l)
277	detr(ig,l)=0.
278	! print*,'detr2<0!!!','ig=',ig,'l=',l,'f=',fm(ig,l),
279	! s 'f+1=',fm(ig,l+1),'e=',entr(ig,l),'d=',detr(ig,l)
280	endif
281
282	! if (fm(ig,l+1).lt.0.) then
283	! print*,'fm2<0!!!'
284	! endif
285
286	! if (entr(ig,l).lt.0.) then
287	! print*,'entr2<0!!!'
288	! endif
289	enddo
290	enddo
291
292	! calcul de la valeur dans les ascendances
293	do ig=1,ngrid
294	qa(ig,1)=q(ig,1)
295	enddo
296
297	do l=2,nlay
298	do ig=1,ngrid
299	if ((fm(ig,l+1)+detr(ig,l))ztimestep.gt.1.e-5masse(ig,l)) then
300	qa(ig,l) = (fm(ig,l)qa(ig,l-1)+entr(ig,l)q(ig,l)) &
301	& / (fm(ig,l+1)+detr(ig,l))
302	else
303	qa(ig,l)=q(ig,l)
304	endif
305
306	if (qa(ig,l).lt.0.) then
307	! print*,'qa<0!!!'
308	endif
309
310	if (q(ig,l).lt.0.) then
311	! print*,'q<0!!!'
312	endif
313	enddo
314	enddo
315
316	! Calcul du flux subsident
317
318	do l=2,nlay
319	do ig=1,ngrid
320	#undef centre
321	#ifdef centre
322	wqd(ig,l)=fm(ig,l)0.5(q(ig,l-1)+q(ig,l))
323	#else
324
325	#define plusqueun
326	#ifdef plusqueun
327	! Schema avec advection sur plus qu'une maille.
328	zzm=masse(ig,l)/ztimestep
329	if (fm(ig,l)>zzm) then
330	wqd(ig,l)=zzmq(ig,l)+(fm(ig,l)-zzm)q(ig,l+1)
331	else
332	wqd(ig,l)=fm(ig,l)*q(ig,l)
333	endif
334	#else
335	wqd(ig,l)=fm(ig,l)*q(ig,l)
336	#endif
337	#endif
338
339	if (wqd(ig,l).lt.0.) then
340	! print*,'wqd<0!!!'
341	endif
342	enddo
343	enddo
344
345	do ig=1,ngrid
346	wqd(ig,1)=0.
347	wqd(ig,nlay+1)=0.
348	enddo
349
350	! Calcul des tendances
351	do l=1,nlay
352	do ig=1,ngrid
353	q(ig,l)=q(ig,l)+(detr(ig,l)qa(ig,l)-entr(ig,l)q(ig,l) &
354	& -wqd(ig,l)+wqd(ig,l+1)) &
355	& *ztimestep/masse(ig,l)
356	! if (dq(ig,l).lt.0.) then
357	! print*,'dq<0!!!'
358	! endif
359	enddo
360	enddo
361	enddo
362
363	! Calcul des tendances
364	do l=1,nlay
365	do ig=1,ngrid
366	dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep
367	q(ig,l) = qold(ig,l)
368	enddo
369	enddo
370
371
372	return
373	end

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

Download in other formats:

Original Format