Context Navigation

thermcell_main_mars.F90 @ 659

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

Thermals: minor changes to some parameters, according to latest high domain size simulations. This really does not change the result when using yamada4 for diffusion anyway, but is consistent with the most precise LES observations. On the conceptual level, one of the important change is that the mass-flux compensation for updrafts and downdrafts is now a subsidence and no more an uplift, as in terrestrial modeling.

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

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

Download in other formats:

Original Format