Context Navigation

thermcell_main_mars.F90 @ 342

Last change on this file since 342 was 337, checked in by acolaitis, 14 years ago

10/10/2011 == AC

*
This commit aims at increasing the thermals speed, especially for large tracer number configurations. The idea behind this commit is to advect non-active conserved variables outside of the sub-timestep of the thermals. Because these variables are not used in thermals computation, we can decouple them:

momentum: can be decoupled because we assume a constant ratio between horizontal velocity in alimentation layer and maximum vertical velocity in the thermals.

tracer: can be decoupled because we do not take condensation of any tracer into account and hence do not liberate latent heat nor form clouds in the thermals.

temperature: cannot be decoupled (of course)
*

D 336 libf/phymars/thermcell_dqupdown.F90
^{---------------- Deleted and replaced by a simpler version. Notes about downdraft advection are still available from revision 336 of SVN in thermcell_dqupdown.}

A 0 libf/phymars/thermcell_dqup.F90
^{---------------- New upward advection for tracers and momentum in thermals. Several changes are done compared to the old approach:}

Updraft quantities are not longer computed by making hypothesis on the amount of advected air.
- In general, the formalism for updraft computation is much simpler and clearer.
Tracer tendancies are no longer computed using the conservation equation. Instead, we use the divergence of an approximated turbulent flux of the concerned quantity, where downdraft are also neglected.

M 336 libf/phymars/thermcell_main_mars.F90
^{---------------- The Main does not call anymore thermcell_dqupdown, which it was doing 2+tracer_number times per subtimestep (140 times per physical step for a 2 tracer config)

M 336 libf/phymars/calltherm_mars.F90}---------------- Entrainment, detrainment and mass-fluxes are recomputed in the sub-timestep loop. Their final value after iterations is used by the new advection routine to compute tracer and momentum fluxes.

* Results

Conservation of tracers has been assessed over 1 yr in 1D and found to be comparable to that obtained with the simple convective adjustment. (it actually seems to be better by a factor of 10%!)
GCM speed-up is of about 20% compared to the old thermal configuration, for a 2 tracer case.
Advection of sharp tracer profiles has been successfully observed, similar to the old method.

File size: 49.0 KB

Line
1	!
2	!
3	SUBROUTINE thermcell_main_mars(ptimestep &
4	& ,pplay,pplev,pphi,zlev,zlay &
5	& ,pu,pv,pt,pq,pq2 &
6	& ,pduadj,pdvadj,pdtadj,pdqadj,pdq2adj &
7	& ,fm,entr,detr,lmax,zmax &
8	& ,r_aspect &
9	& ,zw2,fraca &
10	& ,zpopsk,ztla,heatFlux,heatFlux_down &
11	& ,buoyancyOut, buoyancyEst)
12
13	IMPLICIT NONE
14
15	!=======================================================================
16	! Mars version of thermcell_main. Author : A Colaitis
17	!=======================================================================
18
19	#include "dimensions.h"
20	#include "dimphys.h"
21	#include "comcstfi.h"
22
23	! ============== INPUTS ==============
24
25	REAL, INTENT(IN) :: ptimestep,r_aspect
26	REAL, INTENT(IN) :: pt(ngridmx,nlayermx)
27	REAL, INTENT(IN) :: pu(ngridmx,nlayermx)
28	REAL, INTENT(IN) :: pv(ngridmx,nlayermx)
29	REAL, INTENT(IN) :: pq(ngridmx,nlayermx,nqmx)
30	REAL, INTENT(IN) :: pq2(ngridmx,nlayermx)
31	REAL, INTENT(IN) :: pplay(ngridmx,nlayermx)
32	REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1)
33	REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
34	REAL, INTENT(IN) :: zlay(ngridmx,nlayermx)
35	REAL, INTENT(IN) :: zlev(ngridmx,nlayermx+1)
36
37	! ============== OUTPUTS ==============
38
39	REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx)
40	REAL :: pduadj(ngridmx,nlayermx)
41	REAL :: pdvadj(ngridmx,nlayermx)
42	REAL :: pdqadj(ngridmx,nlayermx,nqmx)
43	! REAL, INTENT(OUT) :: pdq2adj(ngridmx,nlayermx)
44	REAL :: pdq2adj(ngridmx,nlayermx)
45	REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1)
46
47	! Diagnostics
48	REAL, INTENT(OUT) :: heatFlux(ngridmx,nlayermx) ! interface heatflux
49	REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
50	! REAL, INTENT(OUT) :: buoyancyOut(ngridmx,nlayermx) ! interlayer buoyancy term
51	! REAL, INTENT(OUT) :: buoyancyEst(ngridmx,nlayermx) ! interlayer estimated buoyancy term
52
53	! dummy variables when output not needed :
54
55	! REAL :: heatFlux(ngridmx,nlayermx) ! interface heatflux
56	! REAL :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
57	REAL :: buoyancyOut(ngridmx,nlayermx) ! interlayer buoyancy term
58	REAL :: buoyancyEst(ngridmx,nlayermx) ! interlayer estimated buoyancy term
59
60
61	! ============== LOCAL ================
62
63	INTEGER ig,k,l,ll,iq
64	INTEGER lmax(ngridmx),lmin(ngridmx),lalim(ngridmx)
65	INTEGER lmix(ngridmx)
66	INTEGER lmix_bis(ngridmx)
67	REAL linter(ngridmx)
68	REAL zmix(ngridmx)
69	REAL zmax(ngridmx)
70	REAL ztva(ngridmx,nlayermx),zw_est(ngridmx,nlayermx+1),ztva_est(ngridmx,nlayermx)
71	REAL zmax_sec(ngridmx)
72	REAL zh(ngridmx,nlayermx)
73	REAL zdthladj(ngridmx,nlayermx)
74	REAL zdthladj_down(ngridmx,nlayermx)
75	REAL ztvd(ngridmx,nlayermx)
76	REAL ztv(ngridmx,nlayermx)
77	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx),zo(ngridmx,nlayermx)
78	REAL zva(ngridmx,nlayermx)
79	REAL zua(ngridmx,nlayermx)
80
81	REAL zta(ngridmx,nlayermx)
82	REAL fraca(ngridmx,nlayermx+1)
83	REAL q2(ngridmx,nlayermx)
84	REAL rho(ngridmx,nlayermx),rhobarz(ngridmx,nlayermx),masse(ngridmx,nlayermx)
85	REAL zpopsk(ngridmx,nlayermx)
86
87	REAL wmax(ngridmx)
88	REAL wmax_sec(ngridmx)
89	REAL fm(ngridmx,nlayermx+1),entr(ngridmx,nlayermx),detr(ngridmx,nlayermx)
90
91	REAL fm_down(ngridmx,nlayermx+1)
92
93	REAL ztla(ngridmx,nlayermx)
94
95	REAL f_star(ngridmx,nlayermx+1),entr_star(ngridmx,nlayermx)
96	REAL detr_star(ngridmx,nlayermx)
97	REAL alim_star_tot(ngridmx)
98	REAL alim_star(ngridmx,nlayermx)
99	REAL alim_star_clos(ngridmx,nlayermx)
100	REAL f(ngridmx)
101
102	REAL teta_th_int(ngridmx,nlayermx)
103	REAL teta_env_int(ngridmx,nlayermx)
104	REAL teta_down_int(ngridmx,nlayermx)
105
106	CHARACTER (LEN=20) :: modname
107	CHARACTER (LEN=80) :: abort_message
108
109	! ============= PLUME VARIABLES ============
110
111	REAL w_est(ngridmx,nlayermx+1)
112	REAL wa_moy(ngridmx,nlayermx+1)
113	REAL wmaxa(ngridmx)
114	REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
115	LOGICAL active(ngridmx),activetmp(ngridmx)
116	REAL a1,b1,ae,be,ad,bd
117	INTEGER tic
118
119	! ==========================================
120
121	! ============= HEIGHT VARIABLES ===========
122
123	REAL num(ngridmx)
124	REAL denom(ngridmx)
125	REAL zlevinter(ngridmx)
126
127	! =========================================
128
129	! ============= DRY VARIABLES =============
130
131	REAL zw2_dry(ngridmx,nlayermx+1)
132	REAL f_star_dry(ngridmx,nlayermx+1)
133	REAL ztva_dry(ngridmx,nlayermx+1)
134	REAL wmaxa_dry(ngridmx)
135	REAL wa_moy_dry(ngridmx,nlayermx+1)
136	REAL linter_dry(ngridmx),zlevinter_dry(ngridmx)
137	INTEGER lmix_dry(ngridmx),lmax_dry(ngridmx)
138
139	! =========================================
140
141	! ============= CLOSURE VARIABLES =========
142
143	REAL zdenom(ngridmx)
144	REAL alim_star2(ngridmx)
145	REAL alim_star_tot_clos(ngridmx)
146	INTEGER llmax
147
148	! =========================================
149
150	! ============= FLUX2 VARIABLES ===========
151
152	INTEGER ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
153	INTEGER ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
154	REAL zfm
155	REAL f_old,ddd0,eee0,ddd,eee,zzz
156	REAL fomass_max,alphamax
157
158	! =========================================
159
160	! ============= DTETA VARIABLES ===========
161
162	! rien : on prends la divergence du flux turbulent
163
164	! =========================================
165
166	! ============= DV2 VARIABLES =============
167	! not used for now
168
169	REAL qa(ngridmx,nlayermx),detr_dv2(ngridmx,nlayermx),zf,zf2
170	REAL wvd(ngridmx,nlayermx+1),wud(ngridmx,nlayermx+1)
171	REAL gamma0(ngridmx,nlayermx+1),gamma(ngridmx,nlayermx+1)
172	REAL ue(ngridmx,nlayermx),ve(ngridmx,nlayermx)
173	LOGICAL ltherm(ngridmx,nlayermx)
174	REAL dua(ngridmx,nlayermx),dva(ngridmx,nlayermx)
175	INTEGER iter
176	INTEGER nlarga0
177
178	! =========================================
179
180	!-----------------------------------------------------------------------
181	! initialisation:
182	! ---------------
183
184	entr(:,:)=0.
185	detr(:,:)=0.
186	fm(:,:)=0.
187	! zu(:,:)=pu(:,:)
188	! zv(:,:)=pv(:,:)
189	ztv(:,:)=pt(:,:)/zpopsk(:,:)
190
191	!------------------------------------------------------------------------
192	! --------------------
193	!
194	!
195	! + + + + + + + + + + +
196	!
197	!
198	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
199	! wh,wt,wo ...
200	!
201	! + + + + + + + + + + + zh,zu,zv,zo,rho
202	!
203	!
204	! -------------------- zlev(1)
205	! \\\\\\\\\\\\\\\\\\\\
206	!
207	!
208
209	!-----------------------------------------------------------------------
210	! Calcul des altitudes des couches
211	!-----------------------------------------------------------------------
212
213	! do l=2,nlayermx
214	! zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/g
215	! enddo
216	! zlev(:,1)=0.
217	! zlev(:,nlayermx+1)=(2.*pphi(:,nlayermx)-pphi(:,nlayermx-1))/g
218
219	! zlay(:,:)=pphi(:,:)/g
220	!-----------------------------------------------------------------------
221	! Calcul des densites
222	!-----------------------------------------------------------------------
223
224	rho(:,:)=pplay(:,:)/(r*pt(:,:))
225
226	rhobarz(:,1)=rho(:,1)
227
228	do l=2,nlayermx
229	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
230	enddo
231
232	!calcul de la masse
233	do l=1,nlayermx
234	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
235	enddo
236
237
238	!------------------------------------------------------------------
239	!
240	! /\|\
241	! -------- \| F_k+1 -------
242	! ----> D_k
243	! /\|\ <---- E_k , A_k
244	! -------- \| F_k ---------
245	! ----> D_k-1
246	! <---- E_k-1 , A_k-1
247	!
248	!
249	! ---------------------------
250	!
251	! ----- F_lmax+1=0 ---------- \
252	! lmax (zmax) \|
253	! --------------------------- \|
254	! \|
255	! --------------------------- \|
256	! \|
257	! --------------------------- \|
258	! \|
259	! --------------------------- \|
260	! \|
261	! --------------------------- \|
262	! \| E
263	! --------------------------- \| D
264	! \|
265	! --------------------------- \|
266	! \|
267	! --------------------------- \ \|
268	! lalim \| \|
269	! --------------------------- \| \|
270	! \| \|
271	! --------------------------- \| \|
272	! \| A \|
273	! --------------------------- \| \|
274	! \| \|
275	! --------------------------- \| \|
276	! lmin (=1 pour le moment) \| \|
277	! ----- F_lmin=0 ------------ / /
278	!
279	! ---------------------------
280	! //////////////////////////
281	!
282
283	!=============================================================================
284	! Calculs initiaux ne faisant pas intervenir les changements de phase
285	!=============================================================================
286
287	!------------------------------------------------------------------
288	! 1. alim_star est le profil vertical de l'alimentation a la base du
289	! panache thermique, calcule a partir de la flotabilite de l'air sec
290	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
291	!------------------------------------------------------------------
292	!
293	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
294	lmin=1
295
296	!-----------------------------------------------------------------------------
297	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
298	! panache sec conservatif (e=d=0) alimente selon alim_star
299	! Il s'agit d'un calcul de type CAPE
300	! zmax_sec est utilise pour determiner la geometrie du thermique.
301	!------------------------------------------------------------------------------
302	!---------------------------------------------------------------------------------
303	!calcul du melange et des variables dans le thermique
304	!--------------------------------------------------------------------------------
305
306	! ===========================================================================
307	! ===================== PLUME ===============================================
308	! ===========================================================================
309
310	! Initialisations des variables reeles
311	ztva(:,:)=ztv(:,:)
312	ztva_est(:,:)=ztva(:,:)
313	ztla(:,:)=0.
314	zdz=0.
315	zbuoy(:,:)=0.
316	w_est(:,:)=0.
317	f_star(:,:)=0.
318	wa_moy(:,:)=0.
319	linter(:)=1.
320
321	! a1=2.5 ; b1=0.0015 ; ae=0.045 ; be = 0.6 ! svn baseline
322
323	a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57 !improved fits
324	ad = 0.0005114 ; bd = -0.662
325
326	! Initialisation des variables entieres
327	lmix(:)=1
328	lmix_bis(:)=2
329	wmaxa(:)=0.
330	lalim(:)=1
331
332	!-------------------------------------------------------------------------
333	! On ne considere comme actif que les colonnes dont les deux premieres
334	! couches sont instables.
335	!-------------------------------------------------------------------------
336	active(:)=ztv(:,1)>ztv(:,2)
337	do ig=1,ngridmx
338	if (ztv(ig,1)>=(ztv(ig,2))) then
339	alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.) &
340	& *sqrt(zlev(ig,2))
341	! & /sqrt(zlev(ig,2))
342	! & *zlev(ig,2)
343	lalim(ig)=2
344	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1)
345	endif
346	enddo
347
348	do l=2,nlayermx-1
349	! do l=2,4
350	do ig=1,ngridmx
351	if (ztv(ig,l)>(ztv(ig,l+1)+0.) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then
352	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
353	& *sqrt(zlev(ig,l+1))
354	! & *zlev(ig,2)
355	lalim(ig)=l+1
356	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
357	endif
358	enddo
359	enddo
360	do l=1,nlayermx
361	do ig=1,ngridmx
362	if (alim_star_tot(ig) > 1.e-10 ) then
363	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
364	endif
365	enddo
366	enddo
367
368	alim_star_tot(:)=1.
369	! if(alim_star(1,1) .ne. 0.) then
370	! print*, 'lalim star' ,lalim(1)
371	! endif
372
373	!------------------------------------------------------------------------------
374	! Calcul dans la premiere couche
375	! On decide dans cette version que le thermique n'est actif que si la premiere
376	! couche est instable.
377	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
378	! dans une couche l>1
379	!------------------------------------------------------------------------------
380
381	do ig=1,ngridmx
382	! Le panache va prendre au debut les caracteristiques de l'air contenu
383	! dans cette couche.
384	if (active(ig)) then
385	ztla(ig,1)=ztv(ig,1)
386	!cr: attention, prise en compte de f(1)=1 => AC : what ? f(1) =0. ! (d'ou f(2)=a(1)
387	! dans un panache conservatif
388	f_star(ig,1)=0.
389
390	f_star(ig,2)=alim_star(ig,1)
391
392	zw2(ig,2)=2.g(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
393	& *(zlev(ig,2)-zlev(ig,1)) &
394	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
395	w_est(ig,2)=zw2(ig,2)
396
397	endif
398	enddo
399
400
401	!==============================================================================
402	!boucle de calcul de la vitesse verticale dans le thermique
403	!==============================================================================
404	do l=2,nlayermx-1
405	!==============================================================================
406
407
408	! On decide si le thermique est encore actif ou non
409	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
410	do ig=1,ngridmx
411	active(ig)=active(ig) &
412	& .and. zw2(ig,l)>1.e-10 &
413	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
414	enddo
415
416	!---------------------------------------------------------------------------
417	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
418	! sans tenir compte du detrainement et de l'entrainement dans cette
419	! couche
420	! C'est a dire qu'on suppose
421	! ztla(l)=ztla(l-1)
422	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
423	! avant) a l'alimentation pour avoir un calcul plus propre
424	!---------------------------------------------------------------------------
425
426	do ig=1,ngridmx
427	if(active(ig)) then
428
429	! if(l .lt. lalim(ig)) then
430	! ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
431	! & alim_star(ig,l)*ztv(ig,l)) &
432	! & /(f_star(ig,l)+alim_star(ig,l))
433	! else
434	ztva_est(ig,l)=ztla(ig,l-1)
435	! endif
436
437	zdz=zlev(ig,l+1)-zlev(ig,l)
438	zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
439
440	if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
441	w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.zdza1zbuoy(ig,l)-2.zdzw_est(ig,l)b1 &
442	& -2.zdzw_est(ig,l)ae(a1zbuoy(ig,l)/w_est(ig,l)-b1)*be)
443	else
444	w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.zdza1zbuoy(ig,l)-2.zdzw_est(ig,l)b1)
445	endif
446	if (w_est(ig,l+1).lt.0.) then
447	w_est(ig,l+1)=zw2(ig,l)
448	endif
449	endif
450	enddo
451
452	!-------------------------------------------------
453	!calcul des taux d'entrainement et de detrainement
454	!-------------------------------------------------
455
456	do ig=1,ngridmx
457	if (active(ig)) then
458
459	zw2m=w_est(ig,l+1)
460
461	if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
462	entr_star(ig,l)=f_star(ig,l)zdz &
463	& MAX(0.,ae(a1(zbuoy(ig,l)/zw2m)-b1)**be)
464	else
465	entr_star(ig,l)=0.
466	endif
467
468	if(zbuoy(ig,l) .gt. 0.) then
469	if(l .lt. lalim(ig)) then
470	detr_star(ig,l)=0.
471	else
472
473	! detr_star(ig,l) = f_star(ig,l)zdz &
474	! & 0.0105((zbuoy(ig,l)/zw2m)/0.048)*(1./1.7)
475	! detr_star(ig,l) = f_star(ig,l)zdz &
476	! & 0.0085((zbuoy(ig,l)/zw2m)/0.05)*(1./1.55)
477
478	! last baseline from direct les
479	! detr_star(ig,l) = f_star(ig,l)zdz &
480	! & 0.065(2.5(zbuoy(ig,l)/zw2m))**0.75
481
482	! new param from continuity eq with a fit on dfdz
483	detr_star(ig,l) = f_star(ig,l)zdz &
484	& ad
485
486	! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) !svn baseline
487	! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0008)
488
489	! & 0.014((zbuoy(ig,l)/zw2m)/0.05)*(1./1.35)
490	! detr_star(ig,l) = f_star(ig,l)zdz &
491	! & ((zbuoy(ig,l)/zw2m)/2.222)! + 0.0002)
492
493	endif
494	else
495	detr_star(ig,l)=f_star(ig,l)zdz &
496	& bd*zbuoy(ig,l)/zw2m
497
498	! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) !svn baseline
499
500	! & 5.(-afactzbetalphazbuoy(ig,l)/zw2m)
501	! & 5.(-afact*zbuoy(ig,l)/zw2m)
502
503	! last baseline from direct les
504	! & 0.065(-2.5(zbuoy(ig,l)/zw2m))**0.75
505
506	! new param from continuity eq with a fit on dfdz
507
508
509	endif
510
511	! En dessous de lalim, on prend le max de alim_star et entr_star pour
512	! alim_star et 0 sinon
513
514	if (l.lt.lalim(ig)) then
515	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
516	entr_star(ig,l)=0.
517	endif
518
519	! Calcul du flux montant normalise
520
521	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
522	& -detr_star(ig,l)
523
524	endif
525	enddo
526
527
528	!----------------------------------------------------------------------------
529	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
530	!---------------------------------------------------------------------------
531
532	DO tic=0,6 ! internal convergence loop
533	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
534	do ig=1,ngridmx
535	if (activetmp(ig)) then
536
537	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
538	& (alim_star(ig,l)+entr_star(ig,l))*ztv(ig,l)) &
539	& /(f_star(ig,l+1)+detr_star(ig,l))
540
541	endif
542	enddo
543
544	activetmp(:)=activetmp(:).and.(abs(ztla(:,l)-ztva(:,l)).gt.0.01)
545
546	do ig=1,ngridmx
547	if (activetmp(ig)) then
548	ztva(ig,l) = ztla(ig,l)
549	zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
550
551	if (((a1*zbuoy(ig,l)/zw2(ig,l)-b1) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then
552	zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.zdza1*zbuoy(ig,l)- &
553	& 2.zdzzw2(ig,l)b1-2.zdzzw2(ig,l)ae(a1zbuoy(ig,l)/zw2(ig,l)-b1)**be)
554	else
555	zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.zdza1zbuoy(ig,l)-2.zdzzw2(ig,l)b1)
556	endif
557	endif
558	enddo
559
560	! ================ RECOMPUTE ENTR, DETR, and F FROM NEW W2 ===================
561
562
563	do ig=1,ngridmx
564	if (activetmp(ig)) then
565
566	zw2m=zw2(ig,l+1)
567	if(zw2m .gt. 0) then
568	if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
569	entr_star(ig,l)=f_star(ig,l)zdz &
570	& MAX(0.,ae(a1(zbuoy(ig,l)/zw2m)-b1)**be)
571	else
572	entr_star(ig,l)=0.
573	endif
574
575	if(zbuoy(ig,l) .gt. 0.) then
576	if(l .lt. lalim(ig)) then
577	detr_star(ig,l)=0.
578	else
579	detr_star(ig,l) = f_star(ig,l)zdz &
580	& ad
581	endif
582	else
583	detr_star(ig,l)=f_star(ig,l)zdz &
584	& bd*zbuoy(ig,l)/zw2m
585	endif
586	else
587	entr_star(ig,l)=0.
588	detr_star(ig,l)=0.
589	endif
590
591	! En dessous de lalim, on prend le max de alim_star et entr_star pour
592	! alim_star et 0 sinon
593
594	if (l.lt.lalim(ig)) then
595	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
596	entr_star(ig,l)=0.
597	endif
598
599	! Calcul du flux montant normalise
600
601	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
602	& -detr_star(ig,l)
603
604	endif
605	enddo
606
607	ENDDO ! of tic
608
609	!---------------------------------------------------------------------------
610	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
611	!---------------------------------------------------------------------------
612
613	do ig=1,ngridmx
614	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
615	print*,'On tombe sur le cas particulier de thermcell_plume'
616	zw2(ig,l+1)=0.
617	linter(ig)=l+1
618	endif
619
620	if (zw2(ig,l+1).lt.0.) then
621	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
622	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
623	zw2(ig,l+1)=0.
624	endif
625	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
626
627	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
628	! lmix est le niveau de la couche ou w (wa_moy) est maximum
629	!on rajoute le calcul de lmix_bis
630	lmix(ig)=l+1
631	wmaxa(ig)=wa_moy(ig,l+1)
632	endif
633	enddo
634
635	!=========================================================================
636	! FIN DE LA BOUCLE VERTICALE
637	enddo
638	!=========================================================================
639
640	!on recalcule alim_star_tot
641	do ig=1,ngridmx
642	alim_star_tot(ig)=0.
643	enddo
644	do ig=1,ngridmx
645	do l=1,lalim(ig)-1
646	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
647	enddo
648	enddo
649
650	do l=1,nlayermx
651	do ig=1,ngridmx
652	if (alim_star_tot(ig) > 1.e-10 ) then
653	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
654	endif
655	enddo
656	enddo
657
658	! ===========================================================================
659	! ================= FIN PLUME ===============================================
660	! ===========================================================================
661
662
663	! ===========================================================================
664	! ================= HEIGHT ==================================================
665	! ===========================================================================
666
667	! Attention, w2 est transforme en sa racine carree dans cette routine
668
669	!-------------------------------------------------------------------------------
670	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
671	!-------------------------------------------------------------------------------
672
673	!calcul de la hauteur max du thermique
674	do ig=1,ngridmx
675	lmax(ig)=lalim(ig)
676	enddo
677	do ig=1,ngridmx
678	do l=nlayermx,lalim(ig)+1,-1
679	if (zw2(ig,l).le.1.e-10) then
680	lmax(ig)=l-1
681	endif
682	enddo
683	enddo
684
685	! On traite le cas particulier qu'il faudrait éviter ou le thermique
686	! atteind le haut du modele ...
687	do ig=1,ngridmx
688	if ( zw2(ig,nlayermx) > 1.e-10 ) then
689	print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
690	lmax(ig)=nlayermx
691	endif
692	enddo
693
694	! pas de thermique si couche 1 stable
695	! do ig=1,ngridmx
696	! if (lmin(ig).gt.1) then
697	! lmax(ig)=1
698	! lmin(ig)=1
699	! lalim(ig)=1
700	! endif
701	! enddo
702	!
703	! Determination de zw2 max
704	do ig=1,ngridmx
705	wmax(ig)=0.
706	enddo
707
708	do l=1,nlayermx
709	do ig=1,ngridmx
710	if (l.le.lmax(ig)) then
711	if (zw2(ig,l).lt.0.)then
712	print*,'pb2 zw2<0'
713	endif
714	zw2(ig,l)=sqrt(zw2(ig,l))
715	wmax(ig)=max(wmax(ig),zw2(ig,l))
716	else
717	zw2(ig,l)=0.
718	endif
719	enddo
720	enddo
721	! Longueur caracteristique correspondant a la hauteur des thermiques.
722	do ig=1,ngridmx
723	zmax(ig)=0.
724	zlevinter(ig)=zlev(ig,1)
725	enddo
726
727	num(:)=0.
728	denom(:)=0.
729	do ig=1,ngridmx
730	do l=1,nlayermx
731	num(ig)=num(ig)+zw2(ig,l)zlev(ig,l)(zlev(ig,l+1)-zlev(ig,l))
732	denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
733	enddo
734	enddo
735	do ig=1,ngridmx
736	if (denom(ig).gt.1.e-10) then
737	zmax(ig)=2.*num(ig)/denom(ig)
738	endif
739	enddo
740
741	! def de zmix continu (profil parabolique des vitesses)
742	do ig=1,ngridmx
743	if (lmix(ig).gt.1) then
744	! test
745	if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) &
746	& *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) &
747	& -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) &
748	& *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10) &
749	& then
750	!
751	zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) &
752	& ((zlev(ig,lmix(ig))zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)*zlev(ig,lmix(ig)+1))) &
753	& -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) &
754	& ((zlev(ig,lmix(ig)-1)zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))*zlev(ig,lmix(ig))))) &
755	& /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) &
756	& *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) &
757	& -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) &
758	& *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
759	else
760	zmix(ig)=zlev(ig,lmix(ig))
761	print*,'pb zmix'
762	endif
763	else
764	zmix(ig)=0.
765	endif
766	!test
767	if ((zmax(ig)-zmix(ig)).le.0.) then
768	zmix(ig)=0.9*zmax(ig)
769	endif
770	enddo
771	!
772	! calcul du nouveau lmix correspondant
773	do ig=1,ngridmx
774	do l=1,nlayermx
775	if (zmix(ig).ge.zlev(ig,l).and. &
776	& zmix(ig).lt.zlev(ig,l+1)) then
777	lmix(ig)=l
778	endif
779	enddo
780	enddo
781
782
783	! Attention, w2 est transforme en sa racine carree dans cette routine
784
785	! ===========================================================================
786	! ================= FIN HEIGHT ==============================================
787	! ===========================================================================
788
789	! Choix de la fonction d'alimentation utilisee pour la fermeture.
790
791	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
792
793	! ===========================================================================
794	! ============= CLOSURE =====================================================
795	! ===========================================================================
796
797	!-------------------------------------------------------------------------------
798	! Fermeture,determination de f
799	!-------------------------------------------------------------------------------
800	! Appel avec la version seche
801
802	alim_star2(:)=0.
803	alim_star_tot_clos(:)=0.
804	f(:)=0.
805
806	! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
807	llmax=1
808	do ig=1,ngridmx
809	if (lalim(ig)>llmax) llmax=lalim(ig)
810	enddo
811
812
813	! Calcul des integrales sur la verticale de alim_star et de
814	! alim_star^2/(rho dz)
815	do k=1,llmax-1
816	do ig=1,ngridmx
817	if (k<lalim(ig)) then
818	alim_star2(ig)=alim_star2(ig)+alim_star_clos(ig,k)*alim_star_clos(ig,k) &
819	& /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
820	alim_star_tot_clos(ig)=alim_star_tot_clos(ig)+alim_star_clos(ig,k)
821	endif
822	enddo
823	enddo
824
825	! WARNING : MARS MODIF : we have added 2. : ratio of wmax/vmoy
826	! True ratio is 3.5 but wetake into account the vmoy is the one alimentating
827	! the thermal, so there are vs=0 into the vmoy... the true vmoy is lower. (a la louche)
828	! And r_aspect has been changed from 2 to 1.5 from observations
829	do ig=1,ngridmx
830	if (alim_star2(ig)>1.e-10) then
831	! f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/ &
832	! & (max(500.,zmax_sec(ig))r_aspectalim_star2(ig))
833	f(ig)=wmax(ig)*alim_star_tot_clos(ig)/ &
834	& (max(500.,zmax(ig))r_aspectalim_star2(ig))
835
836	endif
837	enddo
838
839	! ===========================================================================
840	! ============= FIN CLOSURE =================================================
841	! ===========================================================================
842
843
844	! ===========================================================================
845	! ============= FLUX2 =======================================================
846	! ===========================================================================
847
848	!-------------------------------------------------------------------------------
849	!deduction des flux
850	!-------------------------------------------------------------------------------
851
852	fomass_max=0.8
853	alphamax=0.5
854
855	ncorecfm1=0
856	ncorecfm2=0
857	ncorecfm3=0
858	ncorecfm4=0
859	ncorecfm5=0
860	ncorecfm6=0
861	ncorecfm7=0
862	ncorecfm8=0
863	ncorecalpha=0
864
865	!-------------------------------------------------------------------------
866	! Multiplication par le flux de masse issu de la femreture
867	!-------------------------------------------------------------------------
868
869	do l=1,nlayermx
870	entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
871	detr(:,l)=f(:)*detr_star(:,l)
872	enddo
873
874	do l=1,nlayermx
875	do ig=1,ngridmx
876	if (l.lt.lmax(ig)) then
877	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
878	elseif(l.eq.lmax(ig)) then
879	fm(ig,l+1)=0.
880	detr(ig,l)=fm(ig,l)+entr(ig,l)
881	else
882	fm(ig,l+1)=0.
883	endif
884	enddo
885	enddo
886
887	! Test provisoire : pour comprendre pourquoi on corrige plein de fois
888	! le cas fm6, on commence par regarder une premiere fois avant les
889	! autres corrections.
890
891	! do l=1,nlayermx
892	! do ig=1,ngridmx
893	! if (detr(ig,l).gt.fm(ig,l)) then
894	! ncorecfm8=ncorecfm8+1
895	! endif
896	! enddo
897	! enddo
898
899	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
900	! FH Version en cours de test;
901	! par rapport a thermcell_flux, on fait une grande boucle sur "l"
902	! et on modifie le flux avec tous les contrï¿½les appliques d'affilee
903	! pour la meme couche
904	! Momentanement, on duplique le calcule du flux pour pouvoir comparer
905	! les flux avant et apres modif
906	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
907
908	do l=1,nlayermx
909
910	do ig=1,ngridmx
911	if (l.lt.lmax(ig)) then
912	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
913	elseif(l.eq.lmax(ig)) then
914	fm(ig,l+1)=0.
915	detr(ig,l)=fm(ig,l)+entr(ig,l)
916	else
917	fm(ig,l+1)=0.
918	endif
919	enddo
920
921
922	!-------------------------------------------------------------------------
923	! Verification de la positivite des flux de masse
924	!-------------------------------------------------------------------------
925
926	do ig=1,ngridmx
927
928	if (fm(ig,l+1).lt.0.) then
929	if((l+1) .eq. lmax(ig)) then
930	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
931	fm(ig,l+1)=0.
932	entr(ig,l+1)=0.
933	ncorecfm2=ncorecfm2+1
934	else
935	print*,'fm(l+1)<0 : ig, l+1,lmax :',ig,l+1,lmax(ig),fm(ig,l+1)
936	ncorecfm1=ncorecfm1+1
937	fm(ig,l+1)=fm(ig,l)
938	detr(ig,l)=entr(ig,l)
939	endif
940	endif
941
942	enddo
943
944	! Les "optimisations" du flux sont desactivees : moins de bidouilles
945	! je considere que celles ci ne sont pas justifiees ou trop delicates
946	! pour MARS, d'apres des observations LES.
947	!-------------------------------------------------------------------------
948	!Test sur fraca croissant
949	!-------------------------------------------------------------------------
950	! if (iflag_thermals_optflux==0) then
951	! do ig=1,ngridmx
952	! if (l.ge.lalim(ig).and.l.le.lmax(ig) &
953	! & .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
954	!! zzz est le flux en l+1 a frac constant
955	! zzz=fm(ig,l)rhobarz(ig,l+1)zw2(ig,l+1) &
956	! & /(rhobarz(ig,l)*zw2(ig,l))
957	! if (fm(ig,l+1).gt.zzz) then
958	! detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
959	! fm(ig,l+1)=zzz
960	! ncorecfm4=ncorecfm4+1
961	! endif
962	! endif
963	! enddo
964	! endif
965	!
966	!-------------------------------------------------------------------------
967	!test sur flux de masse croissant
968	!-------------------------------------------------------------------------
969	! if (iflag_thermals_optflux==0) then
970	! do ig=1,ngridmx
971	! if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
972	! f_old=fm(ig,l+1)
973	! fm(ig,l+1)=fm(ig,l)
974	! detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
975	! ncorecfm5=ncorecfm5+1
976	! endif
977	! enddo
978	! endif
979	!
980	!-------------------------------------------------------------------------
981	!detr ne peut pas etre superieur a fm
982	!-------------------------------------------------------------------------
983
984	do ig=1,ngridmx
985	if (detr(ig,l).gt.fm(ig,l)) then
986	ncorecfm6=ncorecfm6+1
987	detr(ig,l)=fm(ig,l)
988	entr(ig,l)=fm(ig,l+1)
989
990	! Dans le cas ou on est au dessus de la couche d'alimentation et que le
991	! detrainement est plus fort que le flux de masse, on stope le thermique.
992	! endif
993
994	if(l.gt.lmax(ig)) then
995	! if(l.gt.lalim(ig)) then
996	detr(ig,l)=0.
997	fm(ig,l+1)=0.
998	entr(ig,l)=0.
999	endif
1000
1001	endif
1002
1003	enddo
1004
1005	!-------------------------------------------------------------------------
1006	!fm ne peut pas etre negatif
1007	!-------------------------------------------------------------------------
1008
1009	do ig=1,ngridmx
1010	if (fm(ig,l+1).lt.0.) then
1011	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
1012	fm(ig,l+1)=0.
1013	ncorecfm2=ncorecfm2+1
1014	endif
1015	enddo
1016
1017	!-----------------------------------------------------------------------
1018	!la fraction couverte ne peut pas etre superieure a 1
1019	!-----------------------------------------------------------------------
1020	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1021	! FH Partie a revisiter.
1022	! Il semble qu'etaient codees ici deux optiques dans le cas
1023	! F/ (rho *w) > 1
1024	! soit limiter la hauteur du thermique en considerant que c'est
1025	! la derniere chouche, soit limiter F a rho w.
1026	! Dans le second cas, il faut en fait limiter a un peu moins
1027	! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
1028	! dans thermcell_main et qu'il semble de toutes facons deraisonable
1029	! d'avoir des fractions de 1..
1030	! Ci dessous, et dans l'etat actuel, le premier des deux if est
1031	! sans doute inutile.
1032	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1033
1034	do ig=1,ngridmx
1035	if (zw2(ig,l+1).gt.1.e-10) then
1036	zfm=rhobarz(ig,l+1)zw2(ig,l+1)alphamax
1037	if ( fm(ig,l+1) .gt. zfm) then
1038	f_old=fm(ig,l+1)
1039	fm(ig,l+1)=zfm
1040	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
1041	ncorecalpha=ncorecalpha+1
1042	endif
1043	endif
1044
1045	enddo
1046
1047	! Fin de la grande boucle sur les niveaux verticaux
1048	enddo
1049
1050	!-----------------------------------------------------------------------
1051	! On fait en sorte que la quantite totale d'air entraine dans le
1052	! panache ne soit pas trop grande comparee a la masse de la maille
1053	!-----------------------------------------------------------------------
1054
1055	do l=1,nlayermx-1
1056	do ig=1,ngridmx
1057	eee0=entr(ig,l)
1058	ddd0=detr(ig,l)
1059	eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
1060	ddd=detr(ig,l)-eee
1061	if (eee.gt.0.) then
1062	ncorecfm3=ncorecfm3+1
1063	entr(ig,l)=entr(ig,l)-eee
1064	if ( ddd.gt.0.) then
1065	! l'entrainement est trop fort mais l'exces peut etre compense par une
1066	! diminution du detrainement)
1067	detr(ig,l)=ddd
1068	else
1069	! l'entrainement est trop fort mais l'exces doit etre compense en partie
1070	! par un entrainement plus fort dans la couche superieure
1071	if(l.eq.lmax(ig)) then
1072	detr(ig,l)=fm(ig,l)+entr(ig,l)
1073	else
1074	entr(ig,l+1)=entr(ig,l+1)-ddd
1075	detr(ig,l)=0.
1076	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
1077	detr(ig,l)=0.
1078	endif
1079	endif
1080	endif
1081	enddo
1082	enddo
1083	!
1084	! ddd=detr(ig)-entre
1085	!on s assure que tout s annule bien en zmax
1086	do ig=1,ngridmx
1087	fm(ig,lmax(ig)+1)=0.
1088	entr(ig,lmax(ig))=0.
1089	detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
1090	enddo
1091
1092	!-----------------------------------------------------------------------
1093	! Impression du nombre de bidouilles qui ont ete necessaires
1094	!-----------------------------------------------------------------------
1095
1096	!IM 090508 beg
1097	if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngridmx/4. ) then
1098	print*,'thermcell warning : large number of corrections'
1099	print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
1100	& ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
1101	& ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
1102	& ncorecfm6,'x fm6', &
1103	& ncorecfm7,'x fm7', &
1104	& ncorecfm8,'x fm8', &
1105	& ncorecalpha,'x alpha'
1106	endif
1107
1108	! ===========================================================================
1109	! ============= FIN FLUX2 ===================================================
1110	! ===========================================================================
1111
1112
1113	! ===========================================================================
1114	! ============= TRANSPORT ===================================================
1115	! ===========================================================================
1116
1117	!------------------------------------------------------------------
1118	! calcul du transport vertical
1119	!------------------------------------------------------------------
1120
1121	! ------------------------------------------------------------------
1122	! Transport de teta dans l'updraft : (remplace thermcell_dq)
1123	! ------------------------------------------------------------------
1124
1125	zdthladj(:,:)=0.
1126
1127	do ig=1,ngridmx
1128	if(lmax(ig) .gt. 1) then
1129	do k=1,lmax(ig)
1130	zdthladj(ig,k)=(1./masse(ig,k))(fm(ig,k+1)ztv(ig,k+1)- &
1131	& fm(ig,k)ztv(ig,k)+fm(ig,k)ztva(ig,k)-fm(ig,k+1)*ztva(ig,k+1))
1132	if (ztv(ig,k) + ptimestep*zdthladj(ig,k) .le. 0.) then
1133	print,'Teta<0 in thermcell_dTeta up: qenv .. dq : ', ztv(ig,k),ptimestepzdthladj(ig,k)
1134	if(ztv(ig,k) .gt. 0.) then
1135	zdthladj(ig,k)=0.
1136	endif
1137	endif
1138	enddo
1139	endif
1140	enddo
1141
1142	! ------------------------------------------------------------------
1143	! Prescription des proprietes du downdraft
1144	! ------------------------------------------------------------------
1145
1146	ztvd(:,:)=ztv(:,:)
1147	fm_down(:,:)=0.
1148	do ig=1,ngridmx
1149	if (lmax(ig) .gt. 1) then
1150	do l=1,lmax(ig)
1151	if(zlay(ig,l) .le. 0.7*zmax(ig)) then
1152	fm_down(ig,l) =-1.9*fm(ig,l)
1153	endif
1154
1155	if(zlay(ig,l) .le. 0.06*zmax(ig)) then
1156	ztvd(ig,l)=ztv(ig,l)max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)(ztva(ig,l)/ztv(ig,l) - 1.)))
1157	elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
1158	ztvd(ig,l)=ztv(ig,l)max(0.,1.-0.25(ztva(ig,l)/ztv(ig,l) - 1.))
1159	elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
1160	ztvd(ig,l)=ztv(ig,l)max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)(ztva(ig,l)/ztv(ig,l) - 1.)))
1161	else
1162	ztvd(ig,l)=ztv(ig,l)
1163	endif
1164
1165	! if(zlay(ig,l) .le. 0.06*zmax(ig)) then
1166	! ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig))/19.231 + 0.9938)
1167	! elseif(zlay(ig,l) .le. 0.6*zmax(ig)) then
1168	! ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig)-0.075)/187.931 + 0.9982)
1169	! else
1170	!! elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
1171	! ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig) -0.60)/(-1333) + 1.00025)
1172	!! else
1173	!! ztvd(ig,l)=ztv(ig,l)
1174	! endif
1175
1176	enddo
1177	endif
1178	enddo
1179
1180	! ------------------------------------------------------------------
1181	! Transport de teta dans le downdraft : (remplace thermcell_dq)
1182	! ------------------------------------------------------------------
1183
1184	zdthladj_down(:,:)=0.
1185
1186	do ig=1,ngridmx
1187	if(lmax(ig) .gt. 1) then
1188	! No downdraft in the very-near surface layer, we begin at k=3
1189	do k=2,lmax(ig)
1190	zdthladj_down(ig,k)=(1./masse(ig,k))(fm_down(ig,k+1)ztv(ig,k+1)- &
1191	& fm_down(ig,k)ztv(ig,k)+fm_down(ig,k)ztvd(ig,k)-fm_down(ig,k+1)*ztvd(ig,k+1))
1192	if (ztv(ig,k) + ptimestep*zdthladj_down(ig,k) .le. 0.) then
1193	print,'q<0 in thermcell_dTeta down: qenv .. dq : ', ztv(ig,k),ptimestepzdthladj_down(ig,k)
1194	if(ztv(ig,k) .gt. 0.) then
1195	zdthladj(ig,k)=0.
1196	endif
1197	endif
1198	enddo
1199	endif
1200	enddo
1201
1202	!------------------------------------------------------------------
1203	! Calcul de la fraction de l'ascendance
1204	!------------------------------------------------------------------
1205	do ig=1,ngridmx
1206	fraca(ig,1)=0.
1207	fraca(ig,nlayermx+1)=0.
1208	enddo
1209	do l=2,nlayermx
1210	do ig=1,ngridmx
1211	if (zw2(ig,l).gt.1.e-10) then
1212	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
1213	else
1214	fraca(ig,l)=0.
1215	endif
1216	enddo
1217	enddo
1218
1219
1220
1221	! ===========================================================================
1222	! ============= DV2 =========================================================
1223	! ===========================================================================
1224	! ROUTINE OVERIDE : ne prends pas en compte le downdraft
1225	! de plus, le gradient de pression horizontal semble tout deregler... A VOIR
1226
1227	if (0 .eq. 1) then
1228
1229	!------------------------------------------------------------------
1230	! calcul du transport vertical du moment horizontal
1231	!------------------------------------------------------------------
1232
1233	! Calcul du transport de V tenant compte d'echange par gradient
1234	! de pression horizontal avec l'environnement
1235
1236	! calcul du detrainement
1237	!---------------------------
1238
1239	nlarga0=0.
1240
1241	do k=1,nlayermx
1242	do ig=1,ngridmx
1243	detr_dv2(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
1244	enddo
1245	enddo
1246
1247	! calcul de la valeur dans les ascendances
1248	do ig=1,ngridmx
1249	zua(ig,1)=zu(ig,1)
1250	zva(ig,1)=zv(ig,1)
1251	ue(ig,1)=zu(ig,1)
1252	ve(ig,1)=zv(ig,1)
1253	enddo
1254
1255	gamma(1:ngridmx,1)=0.
1256	do k=2,nlayermx
1257	do ig=1,ngridmx
1258	ltherm(ig,k)=(fm(ig,k+1)+detr_dv2(ig,k))ptimestep > 1.e-5masse(ig,k)
1259	if(ltherm(ig,k).and.zmax(ig)>0.) then
1260	gamma0(ig,k)=masse(ig,k) &
1261	& sqrt( 0.5(fraca(ig,k+1)+fraca(ig,k)) ) &
1262	& *0.5/zmax(ig) &
1263	& *1.
1264	else
1265	gamma0(ig,k)=0.
1266	endif
1267	if (ltherm(ig,k).and.zmax(ig)<=0.) nlarga0=nlarga0+1
1268	enddo
1269	enddo
1270
1271	gamma(:,:)=0.
1272
1273	do k=2,nlayermx
1274
1275	do ig=1,ngridmx
1276
1277	if (ltherm(ig,k)) then
1278	dua(ig,k)=zua(ig,k-1)-zu(ig,k-1)
1279	dva(ig,k)=zva(ig,k-1)-zv(ig,k-1)
1280	else
1281	zua(ig,k)=zu(ig,k)
1282	zva(ig,k)=zv(ig,k)
1283	ue(ig,k)=zu(ig,k)
1284	ve(ig,k)=zv(ig,k)
1285	endif
1286	enddo
1287
1288
1289	! Debut des iterations
1290	!----------------------
1291
1292	! AC WARNING : SALE !
1293
1294	do iter=1,5
1295	do ig=1,ngridmx
1296	! Pour memoire : calcul prenant en compte la fraction reelle
1297	! zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
1298	! zf2=1./(1.-zf)
1299	! Calcul avec fraction infiniement petite
1300	zf=0.
1301	zf2=1.
1302
1303	! la première fois on multiplie le coefficient de freinage
1304	! par le module du vent dans la couche en dessous.
1305	! Mais pourquoi donc ???
1306	if (ltherm(ig,k)) then
1307	! On choisit une relaxation lineaire.
1308	! gamma(ig,k)=gamma0(ig,k)
1309	! On choisit une relaxation quadratique.
1310	gamma(ig,k)=gamma0(ig,k)sqrt(dua(ig,k)2+dva(ig,k)*2)
1311	zua(ig,k)=(fm(ig,k)*zua(ig,k-1) &
1312	& +(zf2entr(ig,k)+gamma(ig,k))zu(ig,k)) &
1313	& /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)zfzf2 &
1314	& +gamma(ig,k))
1315	zva(ig,k)=(fm(ig,k)*zva(ig,k-1) &
1316	& +(zf2entr(ig,k)+gamma(ig,k))zv(ig,k)) &
1317	& /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)zfzf2 &
1318	& +gamma(ig,k))
1319
1320	! print*,' OUTPUT DV2 !!!!!!!!!!!!!!!!!!!!',k,zua(ig,k),zva(ig,k),zu(ig,k),zv(ig,k),dua(ig,k),dva(ig,k)
1321	dua(ig,k)=zua(ig,k)-zu(ig,k)
1322	dva(ig,k)=zva(ig,k)-zv(ig,k)
1323	ue(ig,k)=(zu(ig,k)-zfzua(ig,k))zf2
1324	ve(ig,k)=(zv(ig,k)-zfzva(ig,k))zf2
1325	endif
1326	enddo
1327	! Fin des iterations
1328	!--------------------
1329	enddo
1330
1331	enddo ! k=2,nlayermx
1332
1333	! Calcul du flux vertical de moment dans l'environnement.
1334	!---------------------------------------------------------
1335	do k=2,nlayermx
1336	do ig=1,ngridmx
1337	wud(ig,k)=fm(ig,k)*ue(ig,k)
1338	wvd(ig,k)=fm(ig,k)*ve(ig,k)
1339	enddo
1340	enddo
1341	do ig=1,ngridmx
1342	wud(ig,1)=0.
1343	wud(ig,nlayermx+1)=0.
1344	wvd(ig,1)=0.
1345	wvd(ig,nlayermx+1)=0.
1346	enddo
1347
1348	! calcul des tendances.
1349	!-----------------------
1350	do k=1,nlayermx
1351	do ig=1,ngridmx
1352	pduadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zua(ig,k) &
1353	& -(entr(ig,k)+gamma(ig,k))*ue(ig,k) &
1354	& -wud(ig,k)+wud(ig,k+1)) &
1355	& /masse(ig,k)
1356	pdvadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zva(ig,k) &
1357	& -(entr(ig,k)+gamma(ig,k))*ve(ig,k) &
1358	& -wvd(ig,k)+wvd(ig,k+1)) &
1359	& /masse(ig,k)
1360	enddo
1361	enddo
1362
1363
1364	! Sorties eventuelles.
1365	!----------------------
1366
1367	! if (nlarga0>0) then
1368	! print*,'WARNING !!!!!! DANS THERMCELL_DV2 '
1369	! print*,nlarga0,' points pour lesquels laraga=0. dans un thermique'
1370	! print*,'Il faudrait decortiquer ces points'
1371	! endif
1372
1373	! ===========================================================================
1374	! ============= FIN DV2 =====================================================
1375	! ===========================================================================
1376
1377	else
1378
1379	! modname='momentum'
1380	! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, &
1381	! & masse,zu,pduadj,ztvd,fm_down,ztv,modname,lmax)
1382	!
1383	! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, &
1384	! & masse,zv,pdvadj,ztvd,fm_down,ztv,modname,lmax)
1385
1386	endif
1387
1388	!------------------------------------------------------------------
1389	! incrementation dt
1390	!------------------------------------------------------------------
1391
1392	do l=1,nlayermx
1393	do ig=1,ngridmx
1394	pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)
1395	enddo
1396	enddo
1397
1398	!------------------------------------------------------------------
1399	! calcul du transport vertical de traceurs
1400	!------------------------------------------------------------------
1401
1402	! if (nqmx .ne. 0.) then
1403	! modname='tracer'
1404	! DO iq=1,nqmx
1405	! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, &
1406	! & masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax)
1407	!
1408	! ENDDO
1409	! endif
1410
1411	!------------------------------------------------------------------
1412	! calcul du transport vertical de la tke
1413	!------------------------------------------------------------------
1414
1415	! modname='tke'
1416	! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, &
1417	! & masse,pq2,pdq2adj,ztvd,fm_down,ztv,modname,lmax)
1418
1419	! ===========================================================================
1420	! ============= FIN TRANSPORT ===============================================
1421	! ===========================================================================
1422
1423
1424	!------------------------------------------------------------------
1425	! Calculs de diagnostiques pour les sorties
1426	!------------------------------------------------------------------
1427	! DIAGNOSTIQUE
1428	! We compute interface values for teta env and th. The last interface
1429	! value does not matter, as the mass flux is 0 there.
1430
1431
1432	do l=1,nlayermx-1
1433	do ig=1,ngridmx
1434	teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l))
1435	teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l))
1436	teta_env_int(ig,l)=0.5*(ztv(ig,l+1)+ztv(ig,l))
1437	enddo
1438	enddo
1439	do ig=1,ngridmx
1440	teta_th_int(ig,nlayermx)=teta_th_int(ig,nlayermx-1)
1441	teta_env_int(ig,nlayermx)=teta_env_int(ig,nlayermx-1)
1442	teta_down_int(ig,nlayermx)=teta_down_int(ig,nlayermx-1)
1443	enddo
1444	do l=1,nlayermx
1445	do ig=1,ngridmx
1446	heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
1447	! buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1448	! buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1449	heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l)
1450	enddo
1451	enddo
1452
1453	return
1454	end

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

Download in other formats:

Original Format