Context Navigation

thermcell_main_mars.F90 @ 496

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

M 489 libf/phymars/thermcell_main_mars.F90
^{--------------------> New parameters for thermals, following latest version of the relevant article}

M 489 libf/phymars/physiq.F
^{--------------------> Minor modifications}

D 489 libf/phymars/surflayer_interpol.F
A 0 libf/phymars/pbl_parameters.F
^{--------------------> Replaced surflayer_interpol.F by a cleaner and richer version called pbl_parameters.F}

This routine is the base of what will later be implemented in the MCD and as a tool to
compute surface properties from output fields (so that it will ultimately disappear from libf
and might become an utility)

M 489 libf/phymars/vdif_cd.F
^{--------------------> Minor modification to the Richardson computation}

M 489 libf/phymars/vdifc.F
^{--------------------> Minor modification to the aerodynamic conductances computation, replacing wmax by hfmax, which}

is a better proxy for convective activity. Might be replaced by w_star later.

File size: 47.3 KB

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

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

Download in other formats:

Original Format