Context Navigation

thermcell.F @ 1972

Last change on this file since 1972 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 39.2 KB

Line
1	!
2	! $Id: thermcell.F 1907 2013-11-26 13:10:46Z fhourdin $
3	!
4	SUBROUTINE calcul_sec(ngrid,nlay,ptimestep
5	s ,pplay,pplev,pphi,zlev
6	s ,pu,pv,pt,po
7	s ,zmax,wmax,zw2,lmix
8	c s ,pu_therm,pv_therm
9	s ,r_aspect,l_mix,w2di,tho)
10
11	USE dimphy
12	IMPLICIT NONE
13
14	c=======================================================================
15	c
16	c Calcul du transport verticale dans la couche limite en presence
17	c de "thermiques" explicitement representes
18	c
19	c Réécriture à partir d'un listing papier à Habas, le 14/02/00
20	c
21	c le thermique est supposé homogène et dissipé par mélange avec
22	c son environnement. la longueur l_mix contrôle l'efficacité du
23	c mélange
24	c
25	c Le calcul du transport des différentes espèces se fait en prenant
26	c en compte:
27	c 1. un flux de masse montant
28	c 2. un flux de masse descendant
29	c 3. un entrainement
30	c 4. un detrainement
31	c
32	c=======================================================================
33
34	c-----------------------------------------------------------------------
35	c declarations:
36	c -------------
37
38	#include "dimensions.h"
39	cccc#include "dimphy.h"
40	#include "YOMCST.h"
41
42	c arguments:
43	c ----------
44
45	INTEGER ngrid,nlay,w2di
46	REAL tho
47	real ptimestep,l_mix,r_aspect
48	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
49	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
50	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
51	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
52	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
53	real pphi(ngrid,nlay)
54
55	integer idetr
56	save idetr
57	data idetr/3/
58	c$OMP THREADPRIVATE(idetr)
59	c local:
60	c ------
61
62	INTEGER ig,k,l,lmaxa(klon),lmix(klon)
63	real zsortie1d(klon)
64	c CR: on remplace lmax(klon,klev+1)
65	INTEGER lmax(klon),lmin(klon),lentr(klon)
66	real linter(klon)
67	real zmix(klon), fracazmix(klon)
68	c RC
69	real zmax(klon),zw,zw2(klon,klev+1),ztva(klon,klev)
70
71	real zlev(klon,klev+1),zlay(klon,klev)
72	REAL zh(klon,klev),zdhadj(klon,klev)
73	REAL ztv(klon,klev)
74	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
75	REAL wh(klon,klev+1)
76	real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
77	real zla(klon,klev+1)
78	real zwa(klon,klev+1)
79	real zld(klon,klev+1)
80	! real zwd(klon,klev+1)
81	real zsortie(klon,klev)
82	real zva(klon,klev)
83	real zua(klon,klev)
84	real zoa(klon,klev)
85
86	real zha(klon,klev)
87	real wa_moy(klon,klev+1)
88	real fraca(klon,klev+1)
89	real fracc(klon,klev+1)
90	real zf,zf2
91	real thetath2(klon,klev),wth2(klon,klev)
92	! common/comtherm/thetath2,wth2
93
94	real count_time
95	! integer isplit,nsplit
96	integer isplit,nsplit,ialt
97	parameter (nsplit=10)
98	data isplit/0/
99	save isplit
100	c$OMP THREADPRIVATE(isplit)
101
102	logical sorties
103	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
104	real zpspsk(klon,klev)
105
106	c real wmax(klon,klev),wmaxa(klon)
107	real wmax(klon),wmaxa(klon)
108	real wa(klon,klev,klev+1)
109	real wd(klon,klev+1)
110	real larg_part(klon,klev,klev+1)
111	real fracd(klon,klev+1)
112	real xxx(klon,klev+1)
113	real larg_cons(klon,klev+1)
114	real larg_detr(klon,klev+1)
115	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
116	real pu_therm(klon,klev),pv_therm(klon,klev)
117	real fm(klon,klev+1),entr(klon,klev)
118	real fmc(klon,klev+1)
119
120	cCR:nouvelles variables
121	real f_star(klon,klev+1),entr_star(klon,klev)
122	real entr_star_tot(klon),entr_star2(klon)
123	real zalim(klon)
124	integer lalim(klon)
125	real norme(klon)
126	real f(klon), f0(klon)
127	real zlevinter(klon)
128	logical therm
129	logical first
130	data first /.false./
131	save first
132	c$OMP THREADPRIVATE(first)
133	cRC
134
135	character*2 str2
136	character*10 str10
137
138	character (len=20) :: modname='calcul_sec'
139	character (len=80) :: abort_message
140
141
142	! LOGICAL vtest(klon),down
143
144	EXTERNAL SCOPY
145
146	integer ncorrec
147	save ncorrec
148	data ncorrec/0/
149	c$OMP THREADPRIVATE(ncorrec)
150
151	c
152	c-----------------------------------------------------------------------
153	c initialisation:
154	c ---------------
155	c
156	sorties=.true.
157	IF(ngrid.NE.klon) THEN
158	PRINT*
159	PRINT*,'STOP dans convadj'
160	PRINT*,'ngrid =',ngrid
161	PRINT*,'klon =',klon
162	ENDIF
163	c
164	c-----------------------------------------------------------------------
165	c incrementation eventuelle de tendances precedentes:
166	c ---------------------------------------------------
167
168	c print*,'0 OK convect8'
169
170	DO 1010 l=1,nlay
171	DO 1015 ig=1,ngrid
172	zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
173	zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
174	zu(ig,l)=pu(ig,l)
175	zv(ig,l)=pv(ig,l)
176	zo(ig,l)=po(ig,l)
177	ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
178	1015 CONTINUE
179	1010 CONTINUE
180
181	c print*,'1 OK convect8'
182	c --------------------
183	c
184	c
185	c + + + + + + + + + + +
186	c
187	c
188	c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
189	c wh,wt,wo ...
190	c
191	c + + + + + + + + + + + zh,zu,zv,zo,rho
192	c
193	c
194	c -------------------- zlev(1)
195	c \\\\\\\\\\\\\\\\\\\\
196	c
197	c
198
199	c-----------------------------------------------------------------------
200	c Calcul des altitudes des couches
201	c-----------------------------------------------------------------------
202
203	do l=2,nlay
204	do ig=1,ngrid
205	zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
206	enddo
207	enddo
208	do ig=1,ngrid
209	zlev(ig,1)=0.
210	zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
211	enddo
212	do l=1,nlay
213	do ig=1,ngrid
214	zlay(ig,l)=pphi(ig,l)/RG
215	enddo
216	enddo
217
218	c print*,'2 OK convect8'
219	c-----------------------------------------------------------------------
220	c Calcul des densites
221	c-----------------------------------------------------------------------
222
223	do l=1,nlay
224	do ig=1,ngrid
225	rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
226	enddo
227	enddo
228
229	do l=2,nlay
230	do ig=1,ngrid
231	rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
232	enddo
233	enddo
234
235	do k=1,nlay
236	do l=1,nlay+1
237	do ig=1,ngrid
238	wa(ig,k,l)=0.
239	enddo
240	enddo
241	enddo
242
243	c print*,'3 OK convect8'
244	c------------------------------------------------------------------
245	c Calcul de w2, quarre de w a partir de la cape
246	c a partir de w2, on calcule wa, vitesse de l'ascendance
247	c
248	c ATTENTION: Dans cette version, pour cause d'economie de memoire,
249	c w2 est stoke dans wa
250	c
251	c ATTENTION: dans convect8, on n'utilise le calcule des wa
252	c independants par couches que pour calculer l'entrainement
253	c a la base et la hauteur max de l'ascendance.
254	c
255	c Indicages:
256	c l'ascendance provenant du niveau k traverse l'interface l avec
257	c une vitesse wa(k,l).
258	c
259	c --------------------
260	c
261	c + + + + + + + + + +
262	c
263	c wa(k,l) ---- -------------------- l
264	c /\
265	c /\|\|\ + + + + + + + + + +
266	c \|\|
267	c \|\| --------------------
268	c \|\|
269	c \|\| + + + + + + + + + +
270	c \|\|
271	c \|\| --------------------
272	c \|\|__
273	c \|___ + + + + + + + + + + k
274	c
275	c --------------------
276	c
277	c
278	c
279	c------------------------------------------------------------------
280
281	cCR: ponderation entrainement des couches instables
282	cdef des entr_star tels que entr=f*entr_star
283	do l=1,klev
284	do ig=1,ngrid
285	entr_star(ig,l)=0.
286	enddo
287	enddo
288	c determination de la longueur de la couche d entrainement
289	do ig=1,ngrid
290	lentr(ig)=1
291	enddo
292
293	con ne considere que les premieres couches instables
294	therm=.false.
295	do k=nlay-2,1,-1
296	do ig=1,ngrid
297	if (ztv(ig,k).gt.ztv(ig,k+1).and.
298	s ztv(ig,k+1).le.ztv(ig,k+2)) then
299	lentr(ig)=k+1
300	therm=.true.
301	endif
302	enddo
303	enddo
304	climitation de la valeur du lentr
305	c do ig=1,ngrid
306	c lentr(ig)=min(5,lentr(ig))
307	c enddo
308	c determination du lmin: couche d ou provient le thermique
309	do ig=1,ngrid
310	lmin(ig)=1
311	enddo
312	do ig=1,ngrid
313	do l=nlay,2,-1
314	if (ztv(ig,l-1).gt.ztv(ig,l)) then
315	lmin(ig)=l-1
316	endif
317	enddo
318	enddo
319	cinitialisations
320	do ig=1,ngrid
321	zalim(ig)=0.
322	norme(ig)=0.
323	lalim(ig)=1
324	enddo
325	do k=1,klev-1
326	do ig=1,ngrid
327	zalim(ig)=zalim(ig)+zlev(ig,k)*MAX(0.,(ztv(ig,k)-ztv(ig,k+1))
328	s /(zlev(ig,k+1)-zlev(ig,k)))
329	c s *(zlev(ig,k+1)-zlev(ig,k))
330	norme(ig)=norme(ig)+MAX(0.,(ztv(ig,k)-ztv(ig,k+1))
331	s /(zlev(ig,k+1)-zlev(ig,k)))
332	c s *(zlev(ig,k+1)-zlev(ig,k))
333	enddo
334	enddo
335	do ig=1,ngrid
336	if (norme(ig).gt.1.e-10) then
337	zalim(ig)=max(10.*zalim(ig)/norme(ig),zlev(ig,2))
338	c zalim(ig)=min(zalim(ig),zlev(ig,lentr(ig)))
339	endif
340	enddo
341	cdétermination du lalim correspondant
342	do k=1,klev-1
343	do ig=1,ngrid
344	if ((zalim(ig).gt.zlev(ig,k)).and.(zalim(ig).le.zlev(ig,k+1)))
345	s then
346	lalim(ig)=k
347	endif
348	enddo
349	enddo
350	c
351	c definition de l'entrainement des couches
352	do l=1,klev-1
353	do ig=1,ngrid
354	if (ztv(ig,l).gt.ztv(ig,l+1).and.
355	s l.ge.lmin(ig).and.l.lt.lentr(ig)) then
356	entr_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)
357	c s *(zlev(ig,l+1)-zlev(ig,l))
358	s *sqrt(zlev(ig,l+1))
359	cautre def
360	c entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
361	c s /zlev(ig,lentr(ig)+2)))**(3./2.)
362	endif
363	enddo
364	enddo
365	cnouveau test
366	c if (therm) then
367	do l=1,klev-1
368	do ig=1,ngrid
369	if (ztv(ig,l).gt.ztv(ig,l+1).and.
370	s l.ge.lmin(ig).and.l.le.lalim(ig)
371	s .and.zalim(ig).gt.1.e-10) then
372	c if (l.le.lentr(ig)) then
373	c entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
374	c s /zalim(ig)))**(3./2.)
375	c write(10,*)zlev(ig,l),entr_star(ig,l)
376	endif
377	enddo
378	enddo
379	c endif
380	c pas de thermique si couche 1 stable
381	do ig=1,ngrid
382	if (lmin(ig).gt.5) then
383	do l=1,klev
384	entr_star(ig,l)=0.
385	enddo
386	endif
387	enddo
388	c calcul de l entrainement total
389	do ig=1,ngrid
390	entr_star_tot(ig)=0.
391	enddo
392	do ig=1,ngrid
393	do k=1,klev
394	entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
395	enddo
396	enddo
397	c Calcul entrainement normalise
398	do ig=1,ngrid
399	if (entr_star_tot(ig).gt.1.e-10) then
400	c do l=1,lentr(ig)
401	do l=1,klev
402	cdef possibles pour entr_star: zdthetadz, dthetadz, zdtheta
403	entr_star(ig,l)=entr_star(ig,l)/entr_star_tot(ig)
404	enddo
405	endif
406	enddo
407	c
408	c print*,'fin calcul entr_star'
409	do k=1,klev
410	do ig=1,ngrid
411	ztva(ig,k)=ztv(ig,k)
412	enddo
413	enddo
414	cRC
415	c print*,'7 OK convect8'
416	do k=1,klev+1
417	do ig=1,ngrid
418	zw2(ig,k)=0.
419	fmc(ig,k)=0.
420	cCR
421	f_star(ig,k)=0.
422	cRC
423	larg_cons(ig,k)=0.
424	larg_detr(ig,k)=0.
425	wa_moy(ig,k)=0.
426	enddo
427	enddo
428
429	c print*,'8 OK convect8'
430	do ig=1,ngrid
431	linter(ig)=1.
432	lmaxa(ig)=1
433	lmix(ig)=1
434	wmaxa(ig)=0.
435	enddo
436
437	cCR:
438	do l=1,nlay-2
439	do ig=1,ngrid
440	if (ztv(ig,l).gt.ztv(ig,l+1)
441	s .and.entr_star(ig,l).gt.1.e-10
442	s .and.zw2(ig,l).lt.1e-10) then
443	f_star(ig,l+1)=entr_star(ig,l)
444	ctest:calcul de dteta
445	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
446	s *(zlev(ig,l+1)-zlev(ig,l))
447	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
448	larg_detr(ig,l)=0.
449	else if ((zw2(ig,l).ge.1e-10).and.
450	s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
451	f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
452	ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
453	s *ztv(ig,l))/f_star(ig,l+1)
454	zw2(ig,l+1)=zw2(ig,l)(f_star(ig,l)/f_star(ig,l+1))*2+
455	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
456	s *(zlev(ig,l+1)-zlev(ig,l))
457	endif
458	c determination de zmax continu par interpolation lineaire
459	if (zw2(ig,l+1).lt.0.) then
460	ctest
461	if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
462	c print*,'pb linter'
463	endif
464	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
465	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
466	zw2(ig,l+1)=0.
467	lmaxa(ig)=l
468	else
469	if (zw2(ig,l+1).lt.0.) then
470	c print*,'pb1 zw2<0'
471	endif
472	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
473	endif
474	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
475	c lmix est le niveau de la couche ou w (wa_moy) est maximum
476	lmix(ig)=l+1
477	wmaxa(ig)=wa_moy(ig,l+1)
478	endif
479	enddo
480	enddo
481	c print*,'fin calcul zw2'
482	c
483	c Calcul de la couche correspondant a la hauteur du thermique
484	do ig=1,ngrid
485	lmax(ig)=lentr(ig)
486	c lmax(ig)=lalim(ig)
487	enddo
488	do ig=1,ngrid
489	do l=nlay,lentr(ig)+1,-1
490	c do l=nlay,lalim(ig)+1,-1
491	if (zw2(ig,l).le.1.e-10) then
492	lmax(ig)=l-1
493	endif
494	enddo
495	enddo
496	c pas de thermique si couche 1 stable
497	do ig=1,ngrid
498	if (lmin(ig).gt.5) then
499	lmax(ig)=1
500	lmin(ig)=1
501	lentr(ig)=1
502	lalim(ig)=1
503	endif
504	enddo
505	c
506	c Determination de zw2 max
507	do ig=1,ngrid
508	wmax(ig)=0.
509	enddo
510
511	do l=1,nlay
512	do ig=1,ngrid
513	if (l.le.lmax(ig)) then
514	if (zw2(ig,l).lt.0.)then
515	c print*,'pb2 zw2<0'
516	endif
517	zw2(ig,l)=sqrt(zw2(ig,l))
518	wmax(ig)=max(wmax(ig),zw2(ig,l))
519	else
520	zw2(ig,l)=0.
521	endif
522	enddo
523	enddo
524
525	c Longueur caracteristique correspondant a la hauteur des thermiques.
526	do ig=1,ngrid
527	zmax(ig)=0.
528	zlevinter(ig)=zlev(ig,1)
529	enddo
530	do ig=1,ngrid
531	c calcul de zlevinter
532	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
533	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
534	s -zlev(ig,lmax(ig)))
535	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
536	enddo
537	do ig=1,ngrid
538	c write(8,*)zmax(ig),lmax(ig),lentr(ig),lmin(ig)
539	enddo
540	con stope après les calculs de zmax et wmax
541	RETURN
542
543	c print*,'avant fermeture'
544	c Fermeture,determination de f
545	cAttention! entrainement normalisé ou pas?
546	do ig=1,ngrid
547	entr_star2(ig)=0.
548	enddo
549	do ig=1,ngrid
550	if (entr_star_tot(ig).LT.1.e-10) then
551	f(ig)=0.
552	else
553	do k=lmin(ig),lentr(ig)
554	c do k=lmin(ig),lalim(ig)
555	entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
556	s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
557	enddo
558	c Nouvelle fermeture
559	f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
560	s *entr_star2(ig))
561	c s *entr_star_tot(ig)
562	ctest
563	c if (first) then
564	f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
565	s *wmax(ig))
566	c endif
567	endif
568	f0(ig)=f(ig)
569	c first=.true.
570	enddo
571	c print*,'apres fermeture'
572	con stoppe après la fermeture
573	RETURN
574	c Calcul de l'entrainement
575	do k=1,klev
576	do ig=1,ngrid
577	entr(ig,k)=f(ig)*entr_star(ig,k)
578	enddo
579	enddo
580	con stoppe après le calcul de entr
581	c RETURN
582	cCR:test pour entrainer moins que la masse
583	c do ig=1,ngrid
584	c do l=1,lentr(ig)
585	c if ((entr(ig,l)ptimestep).gt.(0.9masse(ig,l))) then
586	c entr(ig,l+1)=entr(ig,l+1)+entr(ig,l)
587	c s -0.9*masse(ig,l)/ptimestep
588	c entr(ig,l)=0.9*masse(ig,l)/ptimestep
589	c endif
590	c enddo
591	c enddo
592	cCR: fin test
593	c Calcul des flux
594	do ig=1,ngrid
595	do l=1,lmax(ig)-1
596	fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
597	enddo
598	enddo
599
600	cRC
601
602
603	c print*,'9 OK convect8'
604	c print*,'WA1 ',wa_moy
605
606	c determination de l'indice du debut de la mixed layer ou w decroit
607
608	c calcul de la largeur de chaque ascendance dans le cas conservatif.
609	c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
610	c d'une couche est égale à la hauteur de la couche alimentante.
611	c La vitesse maximale dans l'ascendance est aussi prise comme estimation
612	c de la vitesse d'entrainement horizontal dans la couche alimentante.
613
614	do l=2,nlay
615	do ig=1,ngrid
616	if (l.le.lmaxa(ig)) then
617	zw=max(wa_moy(ig,l),1.e-10)
618	larg_cons(ig,l)=zmax(ig)*r_aspect
619	s fmc(ig,l)/(rhobarz(ig,l)zw)
620	endif
621	enddo
622	enddo
623
624	do l=2,nlay
625	do ig=1,ngrid
626	if (l.le.lmaxa(ig)) then
627	c if (idetr.eq.0) then
628	c cette option est finalement en dur.
629	if ((l_mix*zlev(ig,l)).lt.0.)then
630	c print,'pb l_mixzlev<0'
631	endif
632	cCR: test: nouvelle def de lambda
633	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
634	if (zw2(ig,l).gt.1.e-10) then
635	larg_detr(ig,l)=sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
636	else
637	larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
638	endif
639	cRC
640	c else if (idetr.eq.1) then
641	c larg_detr(ig,l)=larg_cons(ig,l)
642	c s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
643	c else if (idetr.eq.2) then
644	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
645	c s *sqrt(wa_moy(ig,l))
646	c else if (idetr.eq.4) then
647	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
648	c s *wa_moy(ig,l)
649	c endif
650	endif
651	enddo
652	enddo
653
654	c print*,'10 OK convect8'
655	c print*,'WA2 ',wa_moy
656	c calcul de la fraction de la maille concernée par l'ascendance en tenant
657	c compte de l'epluchage du thermique.
658	c
659	cCR def de zmix continu (profil parabolique des vitesses)
660	do ig=1,ngrid
661	if (lmix(ig).gt.1.) then
662	c test
663	if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
664	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
665	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
666	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
667	s then
668	c
669	zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
670	s ((zlev(ig,lmix(ig)))2-(zlev(ig,lmix(ig)+1))*2)
671	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
672	s ((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))
673	s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
674	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
675	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
676	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
677	else
678	zmix(ig)=zlev(ig,lmix(ig))
679	c print*,'pb zmix'
680	endif
681	else
682	zmix(ig)=0.
683	endif
684	ctest
685	if ((zmax(ig)-zmix(ig)).lt.0.) then
686	zmix(ig)=0.99*zmax(ig)
687	c print*,'pb zmix>zmax'
688	endif
689	enddo
690	c
691	c calcul du nouveau lmix correspondant
692	do ig=1,ngrid
693	do l=1,klev
694	if (zmix(ig).ge.zlev(ig,l).and.
695	s zmix(ig).lt.zlev(ig,l+1)) then
696	lmix(ig)=l
697	endif
698	enddo
699	enddo
700	c
701	do l=2,nlay
702	do ig=1,ngrid
703	if(larg_cons(ig,l).gt.1.) then
704	c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
705	fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
706	s /(r_aspect*zmax(ig))
707	c test
708	fraca(ig,l)=max(fraca(ig,l),0.)
709	fraca(ig,l)=min(fraca(ig,l),0.5)
710	fracd(ig,l)=1.-fraca(ig,l)
711	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
712	else
713	c wa_moy(ig,l)=0.
714	fraca(ig,l)=0.
715	fracc(ig,l)=0.
716	fracd(ig,l)=1.
717	endif
718	enddo
719	enddo
720	cCR: calcul de fracazmix
721	do ig=1,ngrid
722	fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
723	s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
724	s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
725	s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
726	enddo
727	c
728	do l=2,nlay
729	do ig=1,ngrid
730	if(larg_cons(ig,l).gt.1.) then
731	if (l.gt.lmix(ig)) then
732	ctest
733	if (zmax(ig)-zmix(ig).lt.1.e-10) then
734	c print*,'pb xxx'
735	xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
736	else
737	xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
738	endif
739	if (idetr.eq.0) then
740	fraca(ig,l)=fracazmix(ig)
741	else if (idetr.eq.1) then
742	fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
743	else if (idetr.eq.2) then
744	fraca(ig,l)=fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
745	else
746	fraca(ig,l)=fracazmix(ig)xxx(ig,l)*2
747	endif
748	c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
749	fraca(ig,l)=max(fraca(ig,l),0.)
750	fraca(ig,l)=min(fraca(ig,l),0.5)
751	fracd(ig,l)=1.-fraca(ig,l)
752	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
753	endif
754	endif
755	enddo
756	enddo
757
758	c print*,'fin calcul fraca'
759	c print*,'11 OK convect8'
760	c print*,'Ea3 ',wa_moy
761	c------------------------------------------------------------------
762	c Calcul de fracd, wd
763	c somme wa - wd = 0
764	c------------------------------------------------------------------
765
766
767	do ig=1,ngrid
768	fm(ig,1)=0.
769	fm(ig,nlay+1)=0.
770	enddo
771
772	do l=2,nlay
773	do ig=1,ngrid
774	fm(ig,l)=fraca(ig,l)wa_moy(ig,l)rhobarz(ig,l)
775	cCR:test
776	if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
777	s .and.l.gt.lmix(ig)) then
778	fm(ig,l)=fm(ig,l-1)
779	c write(1,*)'ajustement fm, l',l
780	endif
781	c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
782	cRC
783	enddo
784	do ig=1,ngrid
785	if(fracd(ig,l).lt.0.1) then
786	abort_message = 'fracd trop petit'
787	CALL abort_gcm (modname,abort_message,1)
788
789	else
790	c vitesse descendante "diagnostique"
791	wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
792	endif
793	enddo
794	enddo
795
796	do l=1,nlay
797	do ig=1,ngrid
798	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
799	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
800	enddo
801	enddo
802
803	c print*,'12 OK convect8'
804	c print*,'WA4 ',wa_moy
805	cc------------------------------------------------------------------
806	c calcul du transport vertical
807	c------------------------------------------------------------------
808
809	go to 4444
810	c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
811	do l=2,nlay-1
812	do ig=1,ngrid
813	if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
814	s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
815	c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
816	c s ,fm(ig,l+1)*ptimestep
817	c s ,' M=',masse(ig,l),masse(ig,l+1)
818	endif
819	enddo
820	enddo
821
822	do l=1,nlay
823	do ig=1,ngrid
824	if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
825	c print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
826	c s ,entr(ig,l)*ptimestep
827	c s ,' M=',masse(ig,l)
828	endif
829	enddo
830	enddo
831
832	do l=1,nlay
833	do ig=1,ngrid
834	if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
835	c print*,'WARN!!! fm exagere ig=',ig,' l=',l
836	c s ,' FM=',fm(ig,l)
837	endif
838	if(.not.masse(ig,l).ge.1.e-10
839	s .or..not.masse(ig,l).le.1.e4) then
840	c print*,'WARN!!! masse exagere ig=',ig,' l=',l
841	c s ,' M=',masse(ig,l)
842	c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
843	c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
844	c print*,'zlev(ig,l+1),zlev(ig,l)'
845	c s ,zlev(ig,l+1),zlev(ig,l)
846	c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
847	c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
848	endif
849	if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
850	c print*,'WARN!!! entr exagere ig=',ig,' l=',l
851	c s ,' E=',entr(ig,l)
852	endif
853	enddo
854	enddo
855
856	4444 continue
857
858	cCR:redefinition du entr
859	do l=1,nlay
860	do ig=1,ngrid
861	detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
862	if (detr(ig,l).lt.0.) then
863	c entr(ig,l)=entr(ig,l)-detr(ig,l)
864	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
865	detr(ig,l)=0.
866	c print*,'WARNING !!! detrainement negatif ',ig,l
867	endif
868	enddo
869	enddo
870	cRC
871	if (w2di.eq.1) then
872	fm0=fm0+ptimestep*(fm-fm0)/tho
873	entr0=entr0+ptimestep*(entr-entr0)/tho
874	else
875	fm0=fm
876	entr0=entr
877	endif
878
879	if (1.eq.1) then
880	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
881	. ,zh,zdhadj,zha)
882	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
883	. ,zo,pdoadj,zoa)
884	else
885	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
886	. ,zh,zdhadj,zha)
887	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
888	. ,zo,pdoadj,zoa)
889	endif
890
891	if (1.eq.0) then
892	call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
893	. ,fraca,zmax
894	. ,zu,zv,pduadj,pdvadj,zua,zva)
895	else
896	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
897	. ,zu,pduadj,zua)
898	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
899	. ,zv,pdvadj,zva)
900	endif
901
902	do l=1,nlay
903	do ig=1,ngrid
904	zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
905	zf2=zf/(1.-zf)
906	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
907	wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
908	enddo
909	enddo
910
911
912
913	c print*,'13 OK convect8'
914	c print*,'WA5 ',wa_moy
915	do l=1,nlay
916	do ig=1,ngrid
917	pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
918	enddo
919	enddo
920
921
922	c do l=1,nlay
923	c do ig=1,ngrid
924	c if(abs(pdtadj(ig,l))*86400..gt.500.) then
925	c print*,'WARN!!! ig=',ig,' l=',l
926	c s ,' pdtadj=',pdtadj(ig,l)
927	c endif
928	c if(abs(pdoadj(ig,l))*86400..gt.1.) then
929	c print*,'WARN!!! ig=',ig,' l=',l
930	c s ,' pdoadj=',pdoadj(ig,l)
931	c endif
932	c enddo
933	c enddo
934
935	c print*,'14 OK convect8'
936	c------------------------------------------------------------------
937	c Calculs pour les sorties
938	c------------------------------------------------------------------
939
940	if(sorties) then
941	do l=1,nlay
942	do ig=1,ngrid
943	zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
944	zld(ig,l)=fracd(ig,l)*zmax(ig)
945	if(1.-fracd(ig,l).gt.1.e-10)
946	s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
947	enddo
948	enddo
949
950	cdeja fait
951	c do l=1,nlay
952	c do ig=1,ngrid
953	c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
954	c if (detr(ig,l).lt.0.) then
955	c entr(ig,l)=entr(ig,l)-detr(ig,l)
956	c detr(ig,l)=0.
957	c print*,'WARNING !!! detrainement negatif ',ig,l
958	c endif
959	c enddo
960	c enddo
961
962	c print*,'15 OK convect8'
963
964	isplit=isplit+1
965
966
967	c #define und
968	goto 123
969	#ifdef und
970	CALL writeg1d(1,nlay,wd,'wd ','wd ')
971	CALL writeg1d(1,nlay,zwa,'wa ','wa ')
972	CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')
973	CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')
974	CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')
975	CALL writeg1d(1,nlay,zla,'la ','la ')
976	CALL writeg1d(1,nlay,zld,'ld ','ld ')
977	CALL writeg1d(1,nlay,pt,'pt ','pt ')
978	CALL writeg1d(1,nlay,zh,'zh ','zh ')
979	CALL writeg1d(1,nlay,zha,'zha ','zha ')
980	CALL writeg1d(1,nlay,zu,'zu ','zu ')
981	CALL writeg1d(1,nlay,zv,'zv ','zv ')
982	CALL writeg1d(1,nlay,zo,'zo ','zo ')
983	CALL writeg1d(1,nlay,wh,'wh ','wh ')
984	CALL writeg1d(1,nlay,wu,'wu ','wu ')
985	CALL writeg1d(1,nlay,wv,'wv ','wv ')
986	CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')
987	CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')
988	CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')
989	CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')
990	CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')
991	CALL writeg1d(1,nlay,entr ,'entr ','entr ')
992	CALL writeg1d(1,nlay,detr ,'detr ','detr ')
993	CALL writeg1d(1,nlay,fm ,'fm ','fm ')
994
995	CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')
996	CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')
997	CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')
998
999	c recalcul des flux en diagnostique...
1000	c print*,'PAS DE TEMPS ',ptimestep
1001	call dt2F(pplev,pplay,pt,pdtadj,wh)
1002	CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')
1003	#endif
1004	123 continue
1005	#define troisD
1006	#ifdef troisD
1007	c if (sorties) then
1008	c print*,'Debut des wrgradsfi'
1009
1010	c print*,'16 OK convect8'
1011	call wrgradsfi(1,nlay,wd,'wd ','wd ')
1012	call wrgradsfi(1,nlay,zwa,'wa ','wa ')
1013	call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')
1014	call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')
1015	call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')
1016	call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')
1017	c print*,'WA6 ',wa_moy
1018	call wrgradsfi(1,nlay,zla,'la ','la ')
1019	call wrgradsfi(1,nlay,zld,'ld ','ld ')
1020	call wrgradsfi(1,nlay,pt,'pt ','pt ')
1021	call wrgradsfi(1,nlay,zh,'zh ','zh ')
1022	call wrgradsfi(1,nlay,zha,'zha ','zha ')
1023	call wrgradsfi(1,nlay,zua,'zua ','zua ')
1024	call wrgradsfi(1,nlay,zva,'zva ','zva ')
1025	call wrgradsfi(1,nlay,zu,'zu ','zu ')
1026	call wrgradsfi(1,nlay,zv,'zv ','zv ')
1027	call wrgradsfi(1,nlay,zo,'zo ','zo ')
1028	call wrgradsfi(1,nlay,wh,'wh ','wh ')
1029	call wrgradsfi(1,nlay,wu,'wu ','wu ')
1030	call wrgradsfi(1,nlay,wv,'wv ','wv ')
1031	call wrgradsfi(1,nlay,wo,'wo ','wo ')
1032	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
1033	call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')
1034	call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')
1035	call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')
1036	call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')
1037	call wrgradsfi(1,nlay,entr,'entr ','entr ')
1038	call wrgradsfi(1,nlay,detr,'detr ','detr ')
1039	call wrgradsfi(1,nlay,fm,'fm ','fm ')
1040	call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')
1041	call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')
1042	call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')
1043	call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')
1044
1045	call wrgradsfi(1,nlay,zo,'zo ','zo ')
1046	call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')
1047	call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')
1048
1049	cCR:nouveaux diagnostiques
1050	call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ')
1051	call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ')
1052	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
1053	call wrgradsfi(1,1,zmix,'zmix ','zmix ')
1054	zsortie1d(:)=lmax(:)
1055	call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ')
1056	call wrgradsfi(1,1,wmax,'wmax ','wmax ')
1057	zsortie1d(:)=lmix(:)
1058	call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ')
1059	zsortie1d(:)=lentr(:)
1060	call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ')
1061
1062	c print*,'17 OK convect8'
1063
1064	do k=1,klev/10
1065	write(str2,'(i2.2)') k
1066	str10='wa'//str2
1067	do l=1,nlay
1068	do ig=1,ngrid
1069	zsortie(ig,l)=wa(ig,k,l)
1070	enddo
1071	enddo
1072	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
1073	do l=1,nlay
1074	do ig=1,ngrid
1075	zsortie(ig,l)=larg_part(ig,k,l)
1076	enddo
1077	enddo
1078	str10='la'//str2
1079	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
1080	enddo
1081
1082
1083	c print*,'18 OK convect8'
1084	c endif
1085	c print*,'Fin des wrgradsfi'
1086	#endif
1087
1088	endif
1089
1090	c if(wa_moy(1,4).gt.1.e-10) stop
1091
1092	c print*,'19 OK convect8'
1093	return
1094	end
1095
1096	SUBROUTINE fermeture_seche(ngrid,nlay
1097	s ,pplay,pplev,pphi,zlev,rhobarz,f0,zpspsk
1098	s ,alim_star,zh,zo,lentr,lmin,nu_min,nu_max,r_aspect
1099	s ,zmax,wmax)
1100
1101	USE dimphy
1102	IMPLICIT NONE
1103
1104	#include "dimensions.h"
1105	cccc#include "dimphy.h"
1106	#include "YOMCST.h"
1107
1108	INTEGER ngrid,nlay
1109	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
1110	real pphi(ngrid,nlay)
1111	real zlev(klon,klev+1)
1112	real alim_star(klon,klev)
1113	real f0(klon)
1114	integer lentr(klon)
1115	integer lmin(klon)
1116	real zmax(klon)
1117	real wmax(klon)
1118	real nu_min
1119	real nu_max
1120	real r_aspect
1121	real rhobarz(klon,klev+1)
1122	REAL zh(klon,klev)
1123	real zo(klon,klev)
1124	real zpspsk(klon,klev)
1125
1126	integer ig,l
1127
1128	real f_star(klon,klev+1)
1129	real detr_star(klon,klev)
1130	real entr_star(klon,klev)
1131	real zw2(klon,klev+1)
1132	real linter(klon)
1133	integer lmix(klon)
1134	integer lmax(klon)
1135	real zlevinter(klon)
1136	real wa_moy(klon,klev+1)
1137	real wmaxa(klon)
1138	REAL ztv(klon,klev)
1139	REAL ztva(klon,klev)
1140	real nu(klon,klev)
1141	! real zmax0_sec(klon)
1142	! save zmax0_sec
1143	REAL, SAVE, ALLOCATABLE :: zmax0_sec(:)
1144	c$OMP THREADPRIVATE(zmax0_sec)
1145	logical, save :: first = .true.
1146	c$OMP THREADPRIVATE(first)
1147
1148	if (first) then
1149	allocate(zmax0_sec(klon))
1150	first=.false.
1151	endif
1152
1153	do l=1,nlay
1154	do ig=1,ngrid
1155	ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
1156	ztv(ig,l)=ztv(ig,l)(1.+RETVzo(ig,l))
1157	enddo
1158	enddo
1159	do l=1,nlay-2
1160	do ig=1,ngrid
1161	if (ztv(ig,l).gt.ztv(ig,l+1)
1162	s .and.alim_star(ig,l).gt.1.e-10
1163	s .and.zw2(ig,l).lt.1e-10) then
1164	f_star(ig,l+1)=alim_star(ig,l)
1165	ctest:calcul de dteta
1166	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
1167	s *(zlev(ig,l+1)-zlev(ig,l))
1168	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
1169	else if ((zw2(ig,l).ge.1e-10).and.
1170	s (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
1171	cestimation du detrainement a partir de la geometrie du pas precedent
1172	ctests sur la definition du detr
1173	nu(ig,l)=(nu_min+nu_max)/2.
1174	s *(1.-(nu_max-nu_min)/(nu_max+nu_min)
1175	s *tanh((((ztva(ig,l-1)-ztv(ig,l))/ztv(ig,l))/0.0005)))
1176
1177	detr_star(ig,l)=rhobarz(ig,l)
1178	s *sqrt(zw2(ig,l))
1179	s /(r_aspectzmax0_sec(ig))
1180	c s /(r_aspectzmax0(ig))
1181	s (sqrt(nu(ig,l)*zlev(ig,l+1)
1182	s /sqrt(zw2(ig,l)))
1183	s -sqrt(nu(ig,l)*zlev(ig,l)
1184	s /sqrt(zw2(ig,l))))
1185	detr_star(ig,l)=detr_star(ig,l)/f0(ig)
1186	if ((detr_star(ig,l)).gt.f_star(ig,l)) then
1187	detr_star(ig,l)=f_star(ig,l)
1188	endif
1189	entr_star(ig,l)=0.9*detr_star(ig,l)
1190	if ((l.lt.lentr(ig))) then
1191	entr_star(ig,l)=0.
1192	c detr_star(ig,l)=0.
1193	endif
1194	c print*,'ok detr_star'
1195	cprise en compte du detrainement dans le calcul du flux
1196	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
1197	s -detr_star(ig,l)
1198	ctest sur le signe de f_star
1199	if ((f_star(ig,l+1)+detr_star(ig,l)).gt.1.e-10) then
1200	cAM on melange Tl et qt du thermique
1201	ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+(entr_star(ig,l)
1202	s +alim_star(ig,l))
1203	s *ztv(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
1204	zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)
1205	s /(f_star(ig,l+1)+detr_star(ig,l)))**2+
1206	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1207	s *(zlev(ig,l+1)-zlev(ig,l))
1208	endif
1209	endif
1210	c
1211	if (zw2(ig,l+1).lt.0.) then
1212	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
1213	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1214	zw2(ig,l+1)=0.
1215	c print*,'linter=',linter(ig)
1216	else
1217	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1218	endif
1219	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1220	c lmix est le niveau de la couche ou w (wa_moy) est maximum
1221	lmix(ig)=l+1
1222	wmaxa(ig)=wa_moy(ig,l+1)
1223	endif
1224	enddo
1225	enddo
1226	c print*,'fin calcul zw2'
1227	c
1228	c Calcul de la couche correspondant a la hauteur du thermique
1229	do ig=1,ngrid
1230	lmax(ig)=lentr(ig)
1231	enddo
1232	do ig=1,ngrid
1233	do l=nlay,lentr(ig)+1,-1
1234	if (zw2(ig,l).le.1.e-10) then
1235	lmax(ig)=l-1
1236	endif
1237	enddo
1238	enddo
1239	c pas de thermique si couche 1 stable
1240	do ig=1,ngrid
1241	if (lmin(ig).gt.1) then
1242	lmax(ig)=1
1243	lmin(ig)=1
1244	lentr(ig)=1
1245	endif
1246	enddo
1247	c
1248	c Determination de zw2 max
1249	do ig=1,ngrid
1250	wmax(ig)=0.
1251	enddo
1252
1253	do l=1,nlay
1254	do ig=1,ngrid
1255	if (l.le.lmax(ig)) then
1256	if (zw2(ig,l).lt.0.)then
1257	c print*,'pb2 zw2<0'
1258	endif
1259	zw2(ig,l)=sqrt(zw2(ig,l))
1260	wmax(ig)=max(wmax(ig),zw2(ig,l))
1261	else
1262	zw2(ig,l)=0.
1263	endif
1264	enddo
1265	enddo
1266
1267	c Longueur caracteristique correspondant a la hauteur des thermiques.
1268	do ig=1,ngrid
1269	zmax(ig)=0.
1270	zlevinter(ig)=zlev(ig,1)
1271	enddo
1272	do ig=1,ngrid
1273	c calcul de zlevinter
1274	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
1275	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
1276	s -zlev(ig,lmax(ig)))
1277	cpour le cas ou on prend tjs lmin=1
1278	c zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
1279	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
1280	zmax0_sec(ig)=zmax(ig)
1281	enddo
1282
1283	RETURN
1284	END

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

Download in other formats:

Original Format