Context Navigation

thermcell_dq.F90 @ 5454

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

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

Download in other formats:

Original Format