Context Navigation

thermcell_main_mars.F90 @ 161

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

================================================
======== IMPLEMENTATION OF THERMALS ============
================================================

Author: A. Colaitis (2011-06-16)

The main goal of this revision is to start including the thermals into the model
for development purposes. Users should not use the thermals yet, as
several major configuration changes still need to be done.

This version includes :

updraft and downdraft parametrizations
velocity in the thermal, including drag
plume height analysis
closure equation
updraft transport of heat, tracers and momentum
downdraft transport of heat

This model should not be used without upcoming developments, namely :

downdraft transport of tracers and momentum
updraft & downdraft transport of q2 (tke)
revision of vdif_kc to compute q2 for non-stratified cases

Thermals could also include in a later revision :

momentum loss during transport (horizontal drag)

Compilation of the thermals has been successfully tested on ifort, gfortran and pgf90

================================================
================================================

M libf/phymars/callkeys.h
M libf/phymars/inifis.F

Added new control flags to call the thermals :

calltherm (false by default) <- to call thermals
outptherm (false by default) <- to output thermal-related diagnostics (for dev purposes)

================================================

M libf/phymars/vdifc.F
^{------> added a temporary output for thermal-related diagnostics}

M libf/phymars/testphys1d.F
^{------> added treatment for a initialization from a profile of neutral gas (ar)}

-> will be transformed in a decaying tracer for thermal diagnostics

M libf/phymars/physiq.F
^{------> added a section to call the thermals}

-> changed the call to convadj
-> added thermal-related outputs for diagnostics

M libf/phymars/convadj.F
^{------> takes now into account the height of thermals to execute convective adjustment}

=> note : convective adjustment needs to be activated when using thermals, in case of a

second instable layer above the thermals

================================================

A libf/phymars/calltherm_interface.F90
^{------> Interface between physiq.F and the thermals}

A libf/phymars/calltherm_mars.F90
^{------> Routine running the sub-timestep of the thermals}

A libf/phymars/thermcell_main_mars.F90
^{------> Main thermals routine specific to Martian physics}

A libf/phymars/thermcell_dqupdown.F90
^{------> Thermals subroutine computing transport of quantities by updrafts and downdrafts}

A libf/phymars/thermcell.F90
^{------> Module including parameters from the Earth to Mars importation. Will disappear in future dev}

================================================
================================================

File size: 49.9 KB

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

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

Download in other formats:

Original Format