Context Navigation

thermcell_main_mars.F90 @ 290

Last change on this file since 290 was 290, checked in by acolaitis, 13 years ago

--- AC 07/09/2011 ---

This is a major update for thermals and richardson layer parametrization. This update stabilizes thermals (further
studies might show that we can reduce the value of nsplit in gcm. To be tested.), and prevent downdrafts from descending into
the surface layer, which was acting as a cold feedback on the thermals. The Richardson surface layer now features
different gustiness coefficients for Richardson, heat and momentum so that u* and t* are correctly represented.

Upcoming updates will change surflayer_interpol.F90 to implement those changes in the interpolation scheme as well.

*
IMPORTANT : several parts of the vdifc code might want to use these new definitions for gustiness, u* and t*. exemple : dust devil routines
that recompute u* ? lifting routines ? TODO !

M 289 libf/phymars/thermcell_main_mars.F90
^{-----------------> Added iterations to the velocity / buoyancy / entrainment / detrainment}

computation to ensure correct convergence. Iteration number is for now set to
4, which is probably too much. This will be changed once several cases are tested.
The minimum is probably 2 iterations.

M 289 libf/phymars/vdifc.F
^{-----------------> Subgrid gustiness parametrization now uses different gustiness beta coefficients}

for heat and momentum. Comparisons with LES at Exomars landing site, matching u*
and teta* suggests values of beta=0.7 for momentum and beta=1.2 for heat, where
the formula for large scale horizontal winds in the first layer is :

U0^{2 = pu}2 + pv^{2 + (beta*wmax_th)}2

M 289 libf/phymars/vdif_cd.F
^{-----------------> LES data suggests that Richardson number distribution during daytime has a very large}

standart deviation due to highly negative Richardsons on several gridpoints. As a consequence,
the mean Richardson in daytime in the LES is much more negative than in LES. In the gcm,
parametrized subgrid turbulence prevents such cases, which can be dangerous in nearly unstable conditions.
Hence, we use beta=0.5 instead of one, so that we keep the anti-crazy-hfx function of beta and we increase the
norm of negative Richardsons in general for daytime conditions. This is set in conjonction with
beta settings for heat and momentum in vdifc.

M 289 libf/phymars/meso_inc/meso_inc_les.F
^{-----------------> HFX and USTM computations now uses the different betas for heat and momentum.}

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

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

Download in other formats:

Original Format