Context Navigation

thermcell_main_mars.F90 @ 652

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

Thermals:

minor changes

Diffusion with Yamada4:

corrected numerical oscillation in daytime first-level temperature and momentum predictions through the implementation of a sub-timestep (40 sub-steps) of one loop. (cheap in cpu time). This loop is the one that computes the tke at t+dt from the tke at t.
corrected numerical oscillations in nighttime first-level temperature and momentum predictions through the use of a minimum mixing coefficient, computed as in Holtslag and Boville (93), tuned for Mars.
yamada4 now takes into account the diffusion of tke (iflag_pbl=9), as it tends to smooth numerical oscillations. Results are very similar to those obtained without diffusion.

File size: 54.0 KB

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

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

Download in other formats:

Original Format