Context Navigation

thermcell_main_mars.F90 @ 624

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

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

Download in other formats:

Original Format