Context Navigation

thermcell_main_mars.F90 @ 313

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

* AC 10/10/2011 *

*
This commit aims at increasing the thermals speed. Using these corrections, gcm performances in 64x48x32 using 1 tracer goes from 27.9% elapsed time in thermals to 18.76%.

*
Additional work needs to be done in tracer advection to gain speed in high tracer number configuration. (tracer advection (but not momentum nor temperature) could be decoupled from sub-timestep, as they do not act on the thermals scheme (water vapor is neglected as we use theta and not theta_v, and radiative effect of dust is not computed in the thermals.))

*
=> TOP 5 of routine contributions to gcm runtime :

Each sample counts as 0.01 seconds.

% cumulative self self total

time seconds seconds calls s/call s/call name
18.76 6.33 6.33 960 0.01 0.01 thermcell_main_mars_
17.19 12.13 5.80 svml_powf4.A
13.72 16.76 4.63 10369 0.00 0.00 filtreg_

3.94 18.09 1.33 intel_new_memset
3.73 19.35 1.26 2880 0.00 0.00 thermcell_dqupdown_

note: thermcell_main_mars_ does call quite a lot power computations (see svml_powf4.A), but this number will not increase with tracer numbers.

*
=> LOG:

M 312 libf/phymars/thermcell_main_mars.F90
^{------------------- removed (commented) computations on buoyancy which is purely diagnostic}

tuned internal convergence loop and added convergence criterion

M 312 libf/phymars/thermcell_dqupdown.F90
^{------------------- removed (commented) downdraft-related if-loops (as we do not advect tracers and momentum in downdrafts for now)}

M 312 libf/phymars/calltherm_mars.F90
^{------------------- removed (commented) diagnostic-related computations}

changed default thermals spliting and aspect ratio
corrected a bug where maximum height was not correctly computed and could result in convective adjustment used in place of thermals
when using certains sets of nsplit and r_aspect (was not happening with the baseline version, so that this correction is transparent to
users)

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

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

Download in other formats:

Original Format