Context Navigation

thermcell.F @ 4104

Last change on this file since 4104 was 987, checked in by Laurent Fairhead, 16 years ago
Du nettoyage sur le parallelisme, inclusion de nouvelles interfaces pour OPA9 et ORCHIDEE YM LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 39.0 KB

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

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

Download in other formats:

Original Format