Context Navigation

thermcell_main_mars.F90 @ 619

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

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

Download in other formats:

Original Format