Context Navigation

thermcell_plume.F90 @ 2352

Last change on this file since 2352 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 41.6 KB

Line
1	!
2	! $Id: thermcell_plume.F90 2311 2015-06-25 07:45:24Z emillour $
3	!
4	SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
5	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
6	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
7	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
8	& ,lev_out,lunout1,igout)
9	! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
10	!--------------------------------------------------------------------------
11	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12	!--------------------------------------------------------------------------
13	USE IOIPSL, ONLY : getin
14
15	USE print_control_mod, ONLY: prt_level
16	IMPLICIT NONE
17
18	#include "YOMCST.h"
19	#include "YOETHF.h"
20	#include "FCTTRE.h"
21	#include "thermcell.h"
22
23	INTEGER itap
24	INTEGER lunout1,igout
25	INTEGER ngrid,klev
26	REAL ptimestep
27	REAL ztv(ngrid,klev)
28	REAL zthl(ngrid,klev)
29	REAL po(ngrid,klev)
30	REAL zl(ngrid,klev)
31	REAL rhobarz(ngrid,klev)
32	REAL zlev(ngrid,klev+1)
33	REAL pplev(ngrid,klev+1)
34	REAL pphi(ngrid,klev)
35	REAL zpspsk(ngrid,klev)
36	REAL alim_star(ngrid,klev)
37	REAL f0(ngrid)
38	INTEGER lalim(ngrid)
39	integer lev_out ! niveau pour les print
40	integer nbpb
41
42	real alim_star_tot(ngrid)
43
44	REAL ztva(ngrid,klev)
45	REAL ztla(ngrid,klev)
46	REAL zqla(ngrid,klev)
47	REAL zqta(ngrid,klev)
48	REAL zha(ngrid,klev)
49
50	REAL detr_star(ngrid,klev)
51	REAL coefc
52	REAL entr_star(ngrid,klev)
53	REAL detr(ngrid,klev)
54	REAL entr(ngrid,klev)
55
56	REAL csc(ngrid,klev)
57
58	REAL zw2(ngrid,klev+1)
59	REAL w_est(ngrid,klev+1)
60	REAL f_star(ngrid,klev+1)
61	REAL wa_moy(ngrid,klev+1)
62
63	REAL ztva_est(ngrid,klev)
64	REAL ztv_est(ngrid,klev)
65	REAL zqla_est(ngrid,klev)
66	REAL zqsatth(ngrid,klev)
67	REAL zta_est(ngrid,klev)
68	REAL ztemp(ngrid),zqsat(ngrid)
69	REAL zdw2,zdw2bis
70	REAL zw2modif
71	REAL zw2fact,zw2factbis
72	REAL zeps(ngrid,klev)
73
74	REAL linter(ngrid)
75	INTEGER lmix(ngrid)
76	INTEGER lmix_bis(ngrid)
77	REAL wmaxa(ngrid)
78
79	INTEGER ig,l,k,lt,it,lm
80
81	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
82	real zbuoyjam(ngrid,klev),zdqtjam(ngrid,klev)
83	real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
84	real ztv1,ztv2,factinv,zinv,zlmel
85	real zlmelup,zlmeldwn,zlt,zltdwn,zltup
86	real atv1,atv2,btv1,btv2
87	real ztv_est1,ztv_est2
88	real zcor,zdelta,zcvm5,qlbef
89	real zbetalpha, coefzlmel
90	real eps
91	REAL REPS,RLvCp,DDT0
92	PARAMETER (DDT0=.01)
93	logical Zsat
94	LOGICAL active(ngrid),activetmp(ngrid)
95	REAL fact_gamma,fact_gamma2,fact_epsilon2
96
97	REAL, SAVE :: fact_epsilon, fact_epsilon_omp=0.002
98	REAL, SAVE :: betalpha, betalpha_omp=0.9
99	REAL, SAVE :: afact, afact_omp=2./3.
100	REAL, SAVE :: fact_shell, fact_shell_omp=1.
101	REAL,SAVE :: detr_min,detr_min_omp=1.e-5
102	REAL,SAVE :: entr_min,entr_min_omp=1.e-5
103	REAL,SAVE :: detr_q_coef,detr_q_coef_omp=0.012
104	REAL,SAVE :: detr_q_power,detr_q_power_omp=0.5
105	REAL,SAVE :: mix0,mix0_omp=0.
106
107	LOGICAL, SAVE :: first=.true.
108
109	REAL c2(ngrid,klev)
110
111	if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
112	Zsat=.false.
113	! Initialisation
114
115	RLvCp = RLVTT/RCPD
116	IF (first) THEN
117	!$OMP MASTER
118	! FH : if ok_sync=.true. , the time axis is written at each time step
119	! in the output files. Only at the end in the opposite case
120	CALL getin('thermals_fact_epsilon',fact_epsilon_omp)
121	CALL getin('thermals_betalpha',betalpha_omp)
122	CALL getin('thermals_afact',afact_omp)
123	CALL getin('thermals_fact_shell',fact_shell_omp)
124	CALL getin('thermals_detr_min',detr_min_omp)
125	CALL getin('thermals_entr_min',entr_min_omp)
126	CALL getin('thermals_detr_q_coef',detr_q_coef_omp)
127	CALL getin('thermals_detr_q_power',detr_q_power_omp)
128	CALL getin('thermals_mix0',mix0_omp)
129	! CALL getin('thermals_X',X_omp)
130	! X=X_omp
131	!$OMP END MASTER
132	!$OMP BARRIER
133	fact_epsilon=fact_epsilon_omp
134	betalpha=betalpha_omp
135	afact=afact_omp
136	fact_shell=fact_shell_omp
137	detr_min=detr_min_omp
138	entr_min=entr_min_omp
139	detr_q_coef=detr_q_coef_omp
140	detr_q_power=detr_q_power_omp
141	mix0=mix0_omp
142	first=.false.
143	ENDIF
144
145	zbetalpha=betalpha/(1.+betalpha)
146
147
148	! Initialisations des variables r?elles
149	if (1==1) then
150	ztva(:,:)=ztv(:,:)
151	ztva_est(:,:)=ztva(:,:)
152	ztv_est(:,:)=ztv(:,:)
153	ztla(:,:)=zthl(:,:)
154	zqta(:,:)=po(:,:)
155	zqla(:,:)=0.
156	zha(:,:) = ztva(:,:)
157	else
158	ztva(:,:)=0.
159	ztv_est(:,:)=0.
160	ztva_est(:,:)=0.
161	ztla(:,:)=0.
162	zqta(:,:)=0.
163	zha(:,:) =0.
164	endif
165
166	zqla_est(:,:)=0.
167	zqsatth(:,:)=0.
168	zqla(:,:)=0.
169	detr_star(:,:)=0.
170	entr_star(:,:)=0.
171	alim_star(:,:)=0.
172	alim_star_tot(:)=0.
173	csc(:,:)=0.
174	detr(:,:)=0.
175	entr(:,:)=0.
176	zw2(:,:)=0.
177	zbuoy(:,:)=0.
178	zbuoyjam(:,:)=0.
179	gamma(:,:)=0.
180	zeps(:,:)=0.
181	w_est(:,:)=0.
182	f_star(:,:)=0.
183	wa_moy(:,:)=0.
184	linter(:)=1.
185	! linter(:)=1.
186	! Initialisation des variables entieres
187	lmix(:)=1
188	lmix_bis(:)=2
189	wmaxa(:)=0.
190	lalim(:)=1
191
192
193	!-------------------------------------------------------------------------
194	! On ne considere comme actif que les colonnes dont les deux premieres
195	! couches sont instables.
196	!-------------------------------------------------------------------------
197	active(:)=ztv(:,1)>ztv(:,2)
198
199	!-------------------------------------------------------------------------
200	! Definition de l'alimentation a l'origine dans thermcell_init
201	!-------------------------------------------------------------------------
202	do l=1,klev-1
203	do ig=1,ngrid
204	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
205	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
206	& *sqrt(zlev(ig,l+1))
207	lalim(ig)=l+1
208	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
209	! print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
210	endif
211	enddo
212	enddo
213	do l=1,klev
214	do ig=1,ngrid
215	if (alim_star_tot(ig) > 1.e-10 ) then
216	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
217	endif
218	enddo
219	enddo
220	alim_star_tot(:)=1.
221
222
223
224
225	!------------------------------------------------------------------------------
226	! Calcul dans la premiere couche
227	! On decide dans cette version que le thermique n'est actif que si la premiere
228	! couche est instable.
229	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
230	! dans une couche l>1
231	!------------------------------------------------------------------------------
232	do ig=1,ngrid
233	! Le panache va prendre au debut les caracteristiques de l'air contenu
234	! dans cette couche.
235	if (active(ig)) then
236	ztla(ig,1)=zthl(ig,1)
237	zqta(ig,1)=po(ig,1)
238	zqla(ig,1)=zl(ig,1)
239	!cr: attention, prise en compte de f*(1)=1
240	f_star(ig,2)=alim_star(ig,1)
241	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
242	& *(zlev(ig,2)-zlev(ig,1)) &
243	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
244	w_est(ig,2)=zw2(ig,2)
245	endif
246	enddo
247	!
248
249	!==============================================================================
250	!boucle de calcul de la vitesse verticale dans le thermique
251	!==============================================================================
252	do l=2,klev-1
253	!==============================================================================
254
255
256	! On decide si le thermique est encore actif ou non
257	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
258	do ig=1,ngrid
259	active(ig)=active(ig) &
260	& .and. zw2(ig,l)>1.e-10 &
261	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
262	enddo
263
264
265
266	!---------------------------------------------------------------------------
267	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
268	! sans tenir compte du detrainement et de l'entrainement dans cette
269	! couche
270	! C'est a dire qu'on suppose
271	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
272	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
273	! avant) a l'alimentation pour avoir un calcul plus propre
274	!---------------------------------------------------------------------------
275
276	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
277	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
278	do ig=1,ngrid
279	! print*,'active',active(ig),ig,l
280	if(active(ig)) then
281	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
282	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
283	zta_est(ig,l)=ztva_est(ig,l)
284	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
285	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
286	& -zqla_est(ig,l))-zqla_est(ig,l))
287
288
289	!Modif AJAM
290
291	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
292	zdz=zlev(ig,l+1)-zlev(ig,l)
293	lmel=fact_thermals_ed_dz*zlev(ig,l)
294	! lmel=0.09*zlev(ig,l)
295	zlmel=zlev(ig,l)+lmel
296	zlmelup=zlmel+(zdz/2)
297	zlmeldwn=zlmel-(zdz/2)
298
299	lt=l+1
300	zlt=zlev(ig,lt)
301	zdz3=zlev(ig,lt+1)-zlt
302	zltdwn=zlt-zdz3/2
303	zltup=zlt+zdz3/2
304
305	!=========================================================================
306	! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
307	!=========================================================================
308
309	!--------------------------------------------------
310	if (iflag_thermals_ed.lt.8) then
311	!--------------------------------------------------
312	!AJ052014: J'ai remplac?? la boucle do par un do while
313	! afin de faire moins de calcul dans la boucle
314	!--------------------------------------------------
315	do while (zlmelup.gt.zltup)
316	lt=lt+1
317	zlt=zlev(ig,lt)
318	zdz3=zlev(ig,lt+1)-zlt
319	zltdwn=zlt-zdz3/2
320	zltup=zlt+zdz3/2
321	enddo
322	!--------------------------------------------------
323	!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
324	! on cherche o?? se trouve l'altitude d'inversion
325	! en calculant ztv1 (interpolation de la valeur de
326	! theta au niveau lt en utilisant les niveaux lt-1 et
327	! lt-2) et ztv2 (interpolation avec les niveaux lt+1
328	! et lt+2). Si theta r??ellement calcul??e au niveau lt
329	! comprise entre ztv1 et ztv2, alors il y a inversion
330	! et on calcule son altitude zinv en supposant que ztv(lt)
331	! est une combinaison lineaire de ztv1 et ztv2.
332	! Ensuite, on calcule la flottabilite en comparant
333	! la temperature de la couche l a celle de l'air situe
334	! l+lmel plus haut, ce qui necessite de savoir quel fraction
335	! de cet air est au-dessus ou en-dessous de l'inversion
336	!--------------------------------------------------
337	atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2))
338	btv1=(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
339	& /(zlev(ig,lt-1)-zlev(ig,lt-2))
340	atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1))
341	btv2=(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
342	& /(zlev(ig,lt+2)-zlev(ig,lt+1))
343
344	ztv1=atv1*zlt+btv1
345	ztv2=atv2*zlt+btv2
346
347	if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
348
349	!--------------------------------------------------
350	!AJ052014: D??calage de zinv qui est entre le haut
351	! et le bas de la couche lt
352	!--------------------------------------------------
353	factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
354	zinv=zltdwn+zdz3*factinv
355
356
357	if (zlmeldwn.ge.zinv) then
358	ztv_est(ig,l)=atv2*zlmel+btv2
359	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
360	& +(1.-fact_shell)*zbuoy(ig,l)
361	elseif (zlmelup.ge.zinv) then
362	ztv_est2=atv20.5(zlmelup+zinv)+btv2
363	ztv_est1=atv10.5(zinv+zlmeldwn)+btv1
364	ztv_est(ig,l)=((zlmelup-zinv)/zdz)ztv_est2+((zinv-zlmeldwn)/zdz)ztv_est1
365
366	zbuoyjam(ig,l)=fact_shellRG(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- &
367	& ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
368	& ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
369
370	else
371	ztv_est(ig,l)=atv1*zlmel+btv1
372	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
373	& +(1.-fact_shell)*zbuoy(ig,l)
374	endif
375
376	else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
377
378	if (zlmeldwn.gt.zltdwn) then
379	zbuoyjam(ig,l)=fact_shellRG((ztva_est(ig,l)- &
380	& ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l)
381	else
382	zbuoyjam(ig,l)=fact_shellRG(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
383	& ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
384	& ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
385
386	endif
387
388	! zbuoyjam(ig,l)=fact_shellRG(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
389	! & ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
390	! & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
391	! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
392	! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
393	! & po(ig,lt-1))/po(ig,lt-1))
394	endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
395
396	else ! if (iflag_thermals_ed.lt.8) then
397	lt=l+1
398	zlt=zlev(ig,lt)
399	zdz2=zlev(ig,lt)-zlev(ig,l)
400
401	do while (lmel.gt.zdz2)
402	lt=lt+1
403	zlt=zlev(ig,lt)
404	zdz2=zlt-zlev(ig,l)
405	enddo
406	zdz3=zlev(ig,lt+1)-zlt
407	zltdwn=zlev(ig,lt)-zdz3/2
408	zlmelup=zlmel+(zdz/2)
409	coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
410	zbuoyjam(ig,l)=1.RG(coefzlmel*(ztva_est(ig,l)- &
411	& ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
412	& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
413	endif ! if (iflag_thermals_ed.lt.8) then
414
415	!------------------------------------------------
416	!AJAM:nouveau calcul de w?
417	!------------------------------------------------
418	zdz=zlev(ig,l+1)-zlev(ig,l)
419	zdzbis=zlev(ig,l)-zlev(ig,l-1)
420	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
421
422	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
423	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
424	zdw2=afact*zbuoy(ig,l)/fact_epsilon
425	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
426	! zdw2bis=0.5*(zdw2+zdw2bis)
427	lm=Max(1,l-2)
428	! zdw2=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l) &
429	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
430	! zdw2bis=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l-1) &
431	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
432	! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
433	! w_est(ig,l+1)=(zdz/zdzbis)Max(0.0001,exp(-zw2fact) &
434	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
435	! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
436	! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))zdw2+w_est(ig,l)exp(-zw2fact))
437
438	!--------------------------------------------------
439	!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
440	!--------------------------------------------------
441	if (iflag_thermals_ed==8) then
442	! Ancienne version
443	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
444	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
445	! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
446
447	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
448
449	! Nouvelle version Arnaud
450	else
451	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
452	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
453	! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
454
455	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
456
457	! w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))(exp(-zw2fact) &
458	! & (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
459	! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
460
461
462
463	! w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2biszw2fact)exp(-zw2fact))
464
465	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(zw2(ig,l)+zdw2zw2fact)*exp(-zw2fact)+ &
466	! & (zdzbis-zdz)/zdzbis(zw2(ig,l-1)+zdw2biszw2factbis)*exp(-zw2factbis))
467
468	! w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
469
470	endif
471
472
473	if (iflag_thermals_ed<6) then
474	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
475	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
476	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
477	fact_epsilon=0.0002/(zalpha+0.1)
478	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
479	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
480	zdw2=afact*zbuoy(ig,l)/fact_epsilon
481	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
482	! w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2zw2fact)exp(-zw2fact))
483
484	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
485	! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
486	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
487
488	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
489
490
491	endif
492	!--------------------------------------------------
493	!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
494	!on fait max(0.0001,.....)
495	!--------------------------------------------------
496
497	! if (w_est(ig,l+1).lt.0.) then
498	! w_est(ig,l+1)=zw2(ig,l)
499	! w_est(ig,l+1)=0.0001
500	! endif
501
502	endif
503	enddo
504
505
506	!-------------------------------------------------
507	!calcul des taux d'entrainement et de detrainement
508	!-------------------------------------------------
509
510	do ig=1,ngrid
511	if (active(ig)) then
512
513	! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
514	zw2m=w_est(ig,l+1)
515	! zw2m=zw2(ig,l)
516	zdz=zlev(ig,l+1)-zlev(ig,l)
517	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
518	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
519	zbuoybis=zbuoy(ig,l)
520	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
521	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
522
523
524	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
525	! & afact*zbuoybis/zw2m - fact_epsilon )
526
527	! entr_star(ig,l)=MAX(0.,f_star(ig,l)zdzzbetalpha* &
528	! & afact*zbuoybis/zw2m - fact_epsilon )
529
530
531
532	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
533
534	!=========================================================================
535	! 4. Calcul de l'entrainement et du detrainement
536	!=========================================================================
537
538	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
539	! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
540	! entrbis=entr_star(ig,l)
541
542	if (iflag_thermals_ed.lt.6) then
543	fact_epsilon=0.0002/(zalpha+0.1)
544	endif
545
546
547
548	detr_star(ig,l)=f_star(ig,l)*zdz &
549	& ( mix0 0.1 / (zalpha+0.001) &
550	& + MAX(detr_min, -afactzbetalphazbuoyjam(ig,l)/zw2m &
551	& + detr_q_coef(zdqt(ig,l)/zw2m)*detr_q_power))
552
553	! detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
554	! & ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
555
556	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
557
558	entr_star(ig,l)=f_star(ig,l)zdz ( &
559	& mix0 * 0.1 / (zalpha+0.001) &
560	& + zbetalpha*MAX(entr_min, &
561	& afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
562
563
564	! entr_star(ig,l)=f_star(ig,l)zdz ( &
565	! & mix0 * 0.1 / (zalpha+0.001) &
566	! & + MAX(entr_min, &
567	! & zbetalphaafactzbuoyjam(ig,l)/zw2m - fact_epsilon + &
568	! & detr_q_coef(zdqt(ig,l)/zw2m)*detr_q_power))
569
570
571	! entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
572	! & ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
573
574	! entr_star(ig,l)=Max(0.,f_star(ig,l)zdzzbetalpha* &
575	! & afact*zbuoy(ig,l)/zw2m &
576	! & - 1.*fact_epsilon)
577
578
579	! En dessous de lalim, on prend le max de alim_star et entr_star pour
580	! alim_star et 0 sinon
581	if (l.lt.lalim(ig)) then
582	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
583	entr_star(ig,l)=0.
584	endif
585	! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
586	! alim_star(ig,l)=entrbis
587	! endif
588
589	! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
590	! Calcul du flux montant normalise
591	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
592	& -detr_star(ig,l)
593
594	endif
595	enddo
596
597
598	!============================================================================
599	! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
600	!===========================================================================
601
602	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
603	do ig=1,ngrid
604	if (activetmp(ig)) then
605	Zsat=.false.
606	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
607	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
608	& /(f_star(ig,l+1)+detr_star(ig,l))
609	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
610	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
611	& /(f_star(ig,l+1)+detr_star(ig,l))
612
613	endif
614	enddo
615
616	ztemp(:)=zpspsk(:,l)*ztla(:,l)
617	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
618	do ig=1,ngrid
619	if (activetmp(ig)) then
620	! on ecrit de maniere conservative (sat ou non)
621	! T = Tl +Lv/Cp ql
622	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
623	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
624	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
625	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
626	zha(ig,l) = ztva(ig,l)
627	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
628	& -zqla(ig,l))-zqla(ig,l))
629	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
630	zdz=zlev(ig,l+1)-zlev(ig,l)
631	zdzbis=zlev(ig,l)-zlev(ig,l-1)
632	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
633	!!!!!!! fact_epsilon=0.002
634	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
635	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
636	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
637	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
638	! zdw2=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l) &
639	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
640	! lm=Max(1,l-2)
641	! zdw2bis=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l-1) &
642	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
643	! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
644	! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(zw2(ig,l)+zdw2zw2fact)*exp(-zw2fact)+ &
645	! & (zdzbis-zdz)/zdzbis(zw2(ig,l-1)+zdw2biszw2factbis)*exp(-zw2factbis))
646	! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2zw2fact)exp(-zw2fact))
647	! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
648	! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
649	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
650	if (iflag_thermals_ed==8) then
651	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
652	else
653	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
654	endif
655	! zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))(exp(-zw2fact) &
656	! & (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
657	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
658
659
660	if (iflag_thermals_ed.lt.6) then
661	zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
662	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
663	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
664	fact_epsilon=0.0002/(zalpha+0.1)**1
665	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
666	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
667	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
668	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
669
670	! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
671	! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
672	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
673	! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2zw2fact)exp(-zw2fact))
674	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
675
676	endif
677
678
679	endif
680	enddo
681
682	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
683	!
684	!===========================================================================
685	! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
686	!===========================================================================
687
688	nbpb=0
689	do ig=1,ngrid
690	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
691	! stop'On tombe sur le cas particulier de thermcell_dry'
692	! print*,'On tombe sur le cas particulier de thermcell_plume'
693	nbpb=nbpb+1
694	zw2(ig,l+1)=0.
695	linter(ig)=l+1
696	endif
697
698	if (zw2(ig,l+1).lt.0.) then
699	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
700	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
701	zw2(ig,l+1)=0.
702	!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
703	elseif (f_star(ig,l+1).lt.0.) then
704	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
705	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
706	zw2(ig,l+1)=0.
707	!fin CR:04/05/12
708	endif
709
710	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
711
712	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
713	! lmix est le niveau de la couche ou w (wa_moy) est maximum
714	!on rajoute le calcul de lmix_bis
715	if (zqla(ig,l).lt.1.e-10) then
716	lmix_bis(ig)=l+1
717	endif
718	lmix(ig)=l+1
719	wmaxa(ig)=wa_moy(ig,l+1)
720	endif
721	enddo
722
723	if (nbpb>0) then
724	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
725	endif
726
727	!=========================================================================
728	! FIN DE LA BOUCLE VERTICALE
729	enddo
730	!=========================================================================
731
732	!on recalcule alim_star_tot
733	do ig=1,ngrid
734	alim_star_tot(ig)=0.
735	enddo
736	do ig=1,ngrid
737	do l=1,lalim(ig)-1
738	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
739	enddo
740	enddo
741
742
743	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
744
745	#undef wrgrads_thermcell
746	#ifdef wrgrads_thermcell
747	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
748	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
749	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
750	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
751	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
752	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
753	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
754	#endif
755
756
757	return
758	end
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
795	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
796	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
797	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
798	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
799	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
800	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
801	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
802	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
803	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
804	& ,lev_out,lunout1,igout)
805	!& ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
806
807	!--------------------------------------------------------------------------
808	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
809	! Version conforme a l'article de Rio et al. 2010.
810	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
811	!--------------------------------------------------------------------------
812
813	USE print_control_mod, ONLY: prt_level
814	IMPLICIT NONE
815
816	#include "YOMCST.h"
817	#include "YOETHF.h"
818	#include "FCTTRE.h"
819	#include "thermcell.h"
820
821	INTEGER itap
822	INTEGER lunout1,igout
823	INTEGER ngrid,klev
824	REAL ptimestep
825	REAL ztv(ngrid,klev)
826	REAL zthl(ngrid,klev)
827	REAL po(ngrid,klev)
828	REAL zl(ngrid,klev)
829	REAL rhobarz(ngrid,klev)
830	REAL zlev(ngrid,klev+1)
831	REAL pplev(ngrid,klev+1)
832	REAL pphi(ngrid,klev)
833	REAL zpspsk(ngrid,klev)
834	REAL alim_star(ngrid,klev)
835	REAL f0(ngrid)
836	INTEGER lalim(ngrid)
837	integer lev_out ! niveau pour les print
838	integer nbpb
839
840	real alim_star_tot(ngrid)
841
842	REAL ztva(ngrid,klev)
843	REAL ztla(ngrid,klev)
844	REAL zqla(ngrid,klev)
845	REAL zqta(ngrid,klev)
846	REAL zha(ngrid,klev)
847
848	REAL detr_star(ngrid,klev)
849	REAL coefc
850	REAL entr_star(ngrid,klev)
851	REAL detr(ngrid,klev)
852	REAL entr(ngrid,klev)
853
854	REAL csc(ngrid,klev)
855
856	REAL zw2(ngrid,klev+1)
857	REAL w_est(ngrid,klev+1)
858	REAL f_star(ngrid,klev+1)
859	REAL wa_moy(ngrid,klev+1)
860
861	REAL ztva_est(ngrid,klev)
862	REAL zqla_est(ngrid,klev)
863	REAL zqsatth(ngrid,klev)
864	REAL zta_est(ngrid,klev)
865	REAL zbuoyjam(ngrid,klev)
866	REAL ztemp(ngrid),zqsat(ngrid)
867	REAL zdw2
868	REAL zw2modif
869	REAL zw2fact
870	REAL zeps(ngrid,klev)
871
872	REAL linter(ngrid)
873	INTEGER lmix(ngrid)
874	INTEGER lmix_bis(ngrid)
875	REAL wmaxa(ngrid)
876
877	INTEGER ig,l,k
878
879	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
880	real zbuoybis
881	real zcor,zdelta,zcvm5,qlbef,zdz2
882	real betalpha,zbetalpha
883	real eps, afact
884	REAL REPS,RLvCp,DDT0
885	PARAMETER (DDT0=.01)
886	logical Zsat
887	LOGICAL active(ngrid),activetmp(ngrid)
888	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
889	REAL c2(ngrid,klev)
890	Zsat=.false.
891	! Initialisation
892
893	RLvCp = RLVTT/RCPD
894	fact_epsilon=0.002
895	betalpha=0.9
896	afact=2./3.
897
898	zbetalpha=betalpha/(1.+betalpha)
899
900
901	! Initialisations des variables reeles
902	if (1==1) then
903	ztva(:,:)=ztv(:,:)
904	ztva_est(:,:)=ztva(:,:)
905	ztla(:,:)=zthl(:,:)
906	zqta(:,:)=po(:,:)
907	zha(:,:) = ztva(:,:)
908	else
909	ztva(:,:)=0.
910	ztva_est(:,:)=0.
911	ztla(:,:)=0.
912	zqta(:,:)=0.
913	zha(:,:) =0.
914	endif
915
916	zqla_est(:,:)=0.
917	zqsatth(:,:)=0.
918	zqla(:,:)=0.
919	detr_star(:,:)=0.
920	entr_star(:,:)=0.
921	alim_star(:,:)=0.
922	alim_star_tot(:)=0.
923	csc(:,:)=0.
924	detr(:,:)=0.
925	entr(:,:)=0.
926	zw2(:,:)=0.
927	zbuoy(:,:)=0.
928	zbuoyjam(:,:)=0.
929	gamma(:,:)=0.
930	zeps(:,:)=0.
931	w_est(:,:)=0.
932	f_star(:,:)=0.
933	wa_moy(:,:)=0.
934	linter(:)=1.
935	! linter(:)=1.
936	! Initialisation des variables entieres
937	lmix(:)=1
938	lmix_bis(:)=2
939	wmaxa(:)=0.
940	lalim(:)=1
941
942
943	!-------------------------------------------------------------------------
944	! On ne considere comme actif que les colonnes dont les deux premieres
945	! couches sont instables.
946	!-------------------------------------------------------------------------
947	active(:)=ztv(:,1)>ztv(:,2)
948
949	!-------------------------------------------------------------------------
950	! Definition de l'alimentation a l'origine dans thermcell_init
951	!-------------------------------------------------------------------------
952	do l=1,klev-1
953	do ig=1,ngrid
954	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
955	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
956	& *sqrt(zlev(ig,l+1))
957	lalim(ig)=l+1
958	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
959	endif
960	enddo
961	enddo
962	do l=1,klev
963	do ig=1,ngrid
964	if (alim_star_tot(ig) > 1.e-10 ) then
965	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
966	endif
967	enddo
968	enddo
969	alim_star_tot(:)=1.
970
971
972
973	!------------------------------------------------------------------------------
974	! Calcul dans la premiere couche
975	! On decide dans cette version que le thermique n'est actif que si la premiere
976	! couche est instable.
977	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
978	! dans une couche l>1
979	!------------------------------------------------------------------------------
980	do ig=1,ngrid
981	! Le panache va prendre au debut les caracteristiques de l'air contenu
982	! dans cette couche.
983	if (active(ig)) then
984	ztla(ig,1)=zthl(ig,1)
985	zqta(ig,1)=po(ig,1)
986	zqla(ig,1)=zl(ig,1)
987	!cr: attention, prise en compte de f*(1)=1
988	f_star(ig,2)=alim_star(ig,1)
989	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
990	& *(zlev(ig,2)-zlev(ig,1)) &
991	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
992	w_est(ig,2)=zw2(ig,2)
993	endif
994	enddo
995	!
996
997	!==============================================================================
998	!boucle de calcul de la vitesse verticale dans le thermique
999	!==============================================================================
1000	do l=2,klev-1
1001	!==============================================================================
1002
1003
1004	! On decide si le thermique est encore actif ou non
1005	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
1006	do ig=1,ngrid
1007	active(ig)=active(ig) &
1008	& .and. zw2(ig,l)>1.e-10 &
1009	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
1010	enddo
1011
1012
1013
1014	!---------------------------------------------------------------------------
1015	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
1016	! sans tenir compte du detrainement et de l'entrainement dans cette
1017	! couche
1018	! C'est a dire qu'on suppose
1019	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
1020	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
1021	! avant) a l'alimentation pour avoir un calcul plus propre
1022	!---------------------------------------------------------------------------
1023
1024	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
1025	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
1026
1027	do ig=1,ngrid
1028	! print*,'active',active(ig),ig,l
1029	if(active(ig)) then
1030	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
1031	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
1032	zta_est(ig,l)=ztva_est(ig,l)
1033	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
1034	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
1035	& -zqla_est(ig,l))-zqla_est(ig,l))
1036
1037	!------------------------------------------------
1038	!AJAM:nouveau calcul de w?
1039	!------------------------------------------------
1040	zdz=zlev(ig,l+1)-zlev(ig,l)
1041	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1042
1043	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
1044	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
1045	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
1046
1047
1048	if (w_est(ig,l+1).lt.0.) then
1049	w_est(ig,l+1)=zw2(ig,l)
1050	endif
1051	endif
1052	enddo
1053
1054
1055	!-------------------------------------------------
1056	!calcul des taux d'entrainement et de detrainement
1057	!-------------------------------------------------
1058
1059	do ig=1,ngrid
1060	if (active(ig)) then
1061
1062	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
1063	zw2m=w_est(ig,l+1)
1064	zdz=zlev(ig,l+1)-zlev(ig,l)
1065	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1066	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
1067	zbuoybis=zbuoy(ig,l)
1068	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
1069	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
1070
1071
1072	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
1073	& afact*zbuoybis/zw2m - fact_epsilon )
1074
1075
1076	detr_star(ig,l)=f_star(ig,l)*zdz &
1077	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
1078	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
1079
1080	! En dessous de lalim, on prend le max de alim_star et entr_star pour
1081	! alim_star et 0 sinon
1082	if (l.lt.lalim(ig)) then
1083	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
1084	entr_star(ig,l)=0.
1085	endif
1086
1087	! Calcul du flux montant normalise
1088	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
1089	& -detr_star(ig,l)
1090
1091	endif
1092	enddo
1093
1094
1095	!----------------------------------------------------------------------------
1096	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
1097	!---------------------------------------------------------------------------
1098	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
1099	do ig=1,ngrid
1100	if (activetmp(ig)) then
1101	Zsat=.false.
1102	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
1103	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
1104	& /(f_star(ig,l+1)+detr_star(ig,l))
1105	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
1106	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
1107	& /(f_star(ig,l+1)+detr_star(ig,l))
1108
1109	endif
1110	enddo
1111
1112	ztemp(:)=zpspsk(:,l)*ztla(:,l)
1113	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
1114
1115	do ig=1,ngrid
1116	if (activetmp(ig)) then
1117	! on ecrit de maniere conservative (sat ou non)
1118	! T = Tl +Lv/Cp ql
1119	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
1120	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
1121	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1122	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1123	zha(ig,l) = ztva(ig,l)
1124	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
1125	& -zqla(ig,l))-zqla(ig,l))
1126	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1127	zdz=zlev(ig,l+1)-zlev(ig,l)
1128	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
1129
1130	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
1131	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
1132	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
1133	endif
1134	enddo
1135
1136	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
1137	!
1138	!---------------------------------------------------------------------------
1139	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
1140	!---------------------------------------------------------------------------
1141
1142	nbpb=0
1143	do ig=1,ngrid
1144	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1145	! stop'On tombe sur le cas particulier de thermcell_dry'
1146	! print*,'On tombe sur le cas particulier de thermcell_plume'
1147	nbpb=nbpb+1
1148	zw2(ig,l+1)=0.
1149	linter(ig)=l+1
1150	endif
1151
1152	if (zw2(ig,l+1).lt.0.) then
1153	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
1154	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1155	zw2(ig,l+1)=0.
1156	elseif (f_star(ig,l+1).lt.0.) then
1157	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
1158	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
1159	! print*,"linter plume", linter(ig)
1160	zw2(ig,l+1)=0.
1161	endif
1162
1163	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1164
1165	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1166	! lmix est le niveau de la couche ou w (wa_moy) est maximum
1167	!on rajoute le calcul de lmix_bis
1168	if (zqla(ig,l).lt.1.e-10) then
1169	lmix_bis(ig)=l+1
1170	endif
1171	lmix(ig)=l+1
1172	wmaxa(ig)=wa_moy(ig,l+1)
1173	endif
1174	enddo
1175
1176	if (nbpb>0) then
1177	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
1178	endif
1179
1180	!=========================================================================
1181	! FIN DE LA BOUCLE VERTICALE
1182	enddo
1183	!=========================================================================
1184
1185	!on recalcule alim_star_tot
1186	do ig=1,ngrid
1187	alim_star_tot(ig)=0.
1188	enddo
1189	do ig=1,ngrid
1190	do l=1,lalim(ig)-1
1191	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
1192	enddo
1193	enddo
1194
1195
1196	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
1197
1198	#undef wrgrads_thermcell
1199	#ifdef wrgrads_thermcell
1200	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
1201	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
1202	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
1203	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
1204	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
1205	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
1206	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
1207	#endif
1208
1209
1210	return
1211	end

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

Download in other formats:

Original Format