Context Navigation

thermcell_dq.F90 @ 1965

Last change on this file since 1965 was 1943, checked in by idelkadi, 12 years ago

Changes concerinng the Hourdin et al. 2002 version of the thermals and
Ayotte cases

calltherm.F90

call to thermcell_2002 modified :
1) iflag_thermals changed from 1 to >= 1000 ; iflag_thermals-1000
controls sub options.
2) thermals w and fraction in output files.

hbtm.F

Singularity in the 1/heat dependency of the Monin Obukov length L when
heat=0. Since 1/L is used rather than L, it is preferable to compute
directly L. There is a dependency in 1/u* then which is treated with a
threshold value.
(+ some cleaning in the syntax).

lmdz1d.F

If nday<0, -nday is used as the total number of time steps of the
simulation.
The option with imposed wtsurf and wqsurf read in lmdz1d.def was not
active anymore.
< IF(.NOT.ok_flux_surf) THEN
changed to

IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN

before

fsens=-wtsurf*rcpd*rho(1)
flat=-wqsurf*rlvtt*rho(1)

phys_output_write.F90 et phys_local_var_mod.F90

Removing the d_u_ajsb contribution to duthe (same for dv).
Those tendencies are not computed by the model ...
< zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys - &
< d_u_ajsb(1:klon,1:klev)/pdtphys
---

zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys

! d_u_ajsb(1:klon,1:klev)/pdtphys

thermcell_dq.F90, thermcell_main.F90

Some cleaning

thermcell_old.F

Control by iflag_thermals >= 1000
wa_moy and fraca in outputs
+ cleaning

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

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

Download in other formats:

Original Format