Context Navigation

thermcell_main_mars.F90 @ 607

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

Minor change to internal parameters of the thermal: smoothed detrainment profiles, reduced the contribution of upward subsidence to entrainment drag, and increased pressure torque induced drag. Consequences : more precise prediction of zmax (it was overestimated by 5 to 10 percent during the late afternoon), no consequence on the thermal structure, better predictions for internal variables (much better w profile and detrainment rate profile). Tested on both high resolution and low resolution (33 level grid from me and 36 level grid from Forget)

File size: 52.9 KB

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

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

Download in other formats:

Original Format