Context Navigation

thermcell.F @ 3722

Last change on this file since 3722 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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

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

Download in other formats:

Original Format