Context Navigation

thermcell_main_mars.F90 @ 1009

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

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

Download in other formats:

Original Format