Context Navigation

thermcell.F90 @ 5875

Last change on this file since 5875 was 1992, checked in by lguez, 12 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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: 32.6 KB

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

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

Download in other formats:

Original Format