Context Navigation

thermcell_dq.F90 @ 4133

Last change on this file since 4133 was 4133, checked in by fhourdin, 2 years ago
Corrections thermiques pour replay
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.8 KB

Line
1	subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr, &
2	& masse,q,dq,qa,lev_out)
3	USE print_control_mod, ONLY: prt_level
4	implicit none
5
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	! arguments
20	integer, intent(in) :: ngrid,nlay,impl
21	real, intent(in) :: ptimestep
22	real, intent(in), dimension(ngrid,nlay) :: masse
23	real, intent(inout), dimension(ngrid,nlay) :: entr,q
24	real, intent(in), dimension(ngrid,nlay+1) :: fm
25	real, intent(out), dimension(ngrid,nlay) :: dq,qa
26	integer, intent(in) :: lev_out ! niveau pour les print
27
28	! Local
29	real, dimension(ngrid,nlay) :: detr,qold
30	real, dimension(ngrid,nlay+1) :: wqd,fqa
31	real zzm
32	integer ig,k
33	real cfl
34
35	integer niter,iter
36	CHARACTER (LEN=20) :: modname='thermcell_dq'
37	CHARACTER (LEN=80) :: abort_message
38
39
40	! Old explicite scheme
41	if (impl<=-1) then
42	call thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, &
43	& masse,q,dq,qa,lev_out)
44	goto 1000
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,'entrdt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k)
55	abort_message = 'entr dt > m, 1st'
56	CALL abort_physic (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
117	! Trace species evolution
118	if (impl==0) then
119	do k=1,nlay-1
120	q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)q(:,k)+fm(:,k+1)q(:,k+1)) &
121	& *ptimestep/masse(:,k)
122	enddo
123	else
124	do k=nlay-1,1,-1
125	! FH debut de modif : le calcul ci dessous modifiait numériquement
126	! la concentration quand le flux de masse etait nul car on divisait
127	! puis multipliait par masse/ptimestep.
128	! q(:,k)=(masse(:,k)q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)q(:,k+1)) &
129	! & /(fm(:,k)+masse(:,k)/ptimestep)
130	q(:,k)=(q(:,k)+ptimestep/masse(:,k)(fqa(:,k)-fqa(:,k+1)+fm(:,k+1)q(:,k+1))) &
131	& /(1.+fm(:,k)*ptimestep/masse(:,k))
132	! FH fin de modif.
133	enddo
134	endif
135
136	! Tendencies
137	do k=1,nlay
138	do ig=1,ngrid
139	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
140	q(ig,k)=qold(ig,k)
141	enddo
142	enddo
143
144	1000 continue
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	USE print_control_mod, ONLY: prt_level
155	implicit none
156
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_physic (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