Context Navigation

thermcell_main.F90 @ 1294

Last change on this file since 1294 was 1294, checked in by Laurent Fairhead, 14 years ago
Modifications pour la nouvelle version des thermiques (2009/2010) CR et FH
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 26.5 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep &
5	& ,pplay,pplev,pphi,debut &
6	& ,pu,pv,pt,po &
7	& ,pduadj,pdvadj,pdtadj,pdoadj &
8	& ,fm0,entr0,detr0,zqta,zqla,lmax &
9	& ,ratqscth,ratqsdiff,zqsatth &
10	& ,r_aspect,l_mix,tau_thermals &
11	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
12	& ,zmax0, f0,zw2,fraca)
13
14	USE dimphy
15	USE comgeomphy , ONLY:rlond,rlatd
16	IMPLICIT NONE
17
18	!=======================================================================
19	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
20	! Version du 09.02.07
21	! Calcul du transport vertical dans la couche limite en presence
22	! de "thermiques" explicitement representes avec processus nuageux
23	!
24	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
25	!
26	! le thermique est suppose homogene et dissipe par melange avec
27	! son environnement. la longueur l_mix controle l'efficacite du
28	! melange
29	!
30	! Le calcul du transport des differentes especes se fait en prenant
31	! en compte:
32	! 1. un flux de masse montant
33	! 2. un flux de masse descendant
34	! 3. un entrainement
35	! 4. un detrainement
36	!
37	!=======================================================================
38
39	!-----------------------------------------------------------------------
40	! declarations:
41	! -------------
42
43	#include "dimensions.h"
44	#include "YOMCST.h"
45	#include "YOETHF.h"
46	#include "FCTTRE.h"
47	#include "iniprint.h"
48
49	! arguments:
50	! ----------
51
52	!IM 140508
53	INTEGER itap
54
55	INTEGER ngrid,nlay,w2di
56	real tau_thermals
57	real ptimestep,l_mix,r_aspect
58	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
59	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
60	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
61	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
62	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
63	real pphi(ngrid,nlay)
64
65	! local:
66	! ------
67
68	integer icount
69	data icount/0/
70	save icount
71	!$OMP THREADPRIVATE(icount)
72
73	integer,save :: igout=1
74	!$OMP THREADPRIVATE(igout)
75	integer,save :: lunout1=6
76	!$OMP THREADPRIVATE(lunout1)
77	integer,save :: lev_out=10
78	!$OMP THREADPRIVATE(lev_out)
79
80	INTEGER ig,k,l,ll
81	real zsortie1d(klon)
82	INTEGER lmax(klon),lmin(klon),lalim(klon)
83	INTEGER lmix(klon)
84	INTEGER lmix_bis(klon)
85	real linter(klon)
86	real zmix(klon)
87	real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
88	! real fraca(klon,klev)
89
90	real zmax_sec(klon)
91	!on garde le zmax du pas de temps precedent
92	real zmax0(klon)
93	!FH/IM save zmax0
94
95	real lambda
96
97	real zlev(klon,klev+1),zlay(klon,klev)
98	real deltaz(klon,klev)
99	REAL zh(klon,klev)
100	real zthl(klon,klev),zdthladj(klon,klev)
101	REAL ztv(klon,klev)
102	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
103	real zl(klon,klev)
104	real zsortie(klon,klev)
105	real zva(klon,klev)
106	real zua(klon,klev)
107	real zoa(klon,klev)
108
109	real zta(klon,klev)
110	real zha(klon,klev)
111	real fraca(klon,klev+1)
112	real zf,zf2
113	real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
114	real q2(klon,klev)
115	! FH probleme de dimensionnement avec l'allocation dynamique
116	! common/comtherm/thetath2,wth2
117	real wq(klon,klev)
118	real wthl(klon,klev)
119	real wthv(klon,klev)
120
121	real ratqscth(klon,klev)
122	real var
123	real vardiff
124	real ratqsdiff(klon,klev)
125
126	logical sorties
127	real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
128	real zpspsk(klon,klev)
129
130	real wmax(klon)
131	real wmax_sec(klon)
132	real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
133	real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
134
135	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
136	!niveau de condensation
137	integer nivcon(klon)
138	real zcon(klon)
139	REAL CHI
140	real zcon2(klon)
141	real pcon(klon)
142	real zqsat(klon,klev)
143	real zqsatth(klon,klev)
144
145	real f_star(klon,klev+1),entr_star(klon,klev)
146	real detr_star(klon,klev)
147	real alim_star_tot(klon)
148	real alim_star(klon,klev)
149	real alim_star_clos(klon,klev)
150	real f(klon), f0(klon)
151	!FH/IM save f0
152	real zlevinter(klon)
153	logical debut
154	real seuil
155	real csc(klon,klev)
156
157	!
158	!nouvelles variables pour la convection
159	real Ale_bl(klon)
160	real Alp_bl(klon)
161	real alp_int(klon)
162	real ale_int(klon)
163	integer n_int(klon)
164	real fm_tot(klon)
165	real wght_th(klon,klev)
166	integer lalim_conv(klon)
167	!v1d logical therm
168	!v1d save therm
169
170	character*2 str2
171	character*10 str10
172
173	EXTERNAL SCOPY
174	!
175
176	!-----------------------------------------------------------------------
177	! initialisation:
178	! ---------------
179	!
180
181	seuil=0.25
182
183	if (debut) then
184	fm0=0.
185	entr0=0.
186	detr0=0.
187
188
189	! #define wrgrads_thermcell
190	#undef wrgrads_thermcell
191	#ifdef wrgrads_thermcell
192	! Initialisation des sorties grads pour les thermiques.
193	! Pour l'instant en 1D sur le point igout.
194	! Utilise par thermcell_out3d.h
195	str10='therm'
196	call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
197	& rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., &
198	& ptimestep,str10,'therm ')
199	#endif
200
201
202
203	endif
204
205	fm=0. ; entr=0. ; detr=0.
206
207	icount=icount+1
208
209	!IM 090508 beg
210	!print*,'====================================================================='
211	!print*,'====================================================================='
212	!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
213	!print*,'====================================================================='
214	!print*,'====================================================================='
215	!IM 090508 end
216
217	if (prt_level.ge.1) print*,'thermcell_main V4'
218
219	sorties=.true.
220	IF(ngrid.NE.klon) THEN
221	PRINT*
222	PRINT*,'STOP dans convadj'
223	PRINT*,'ngrid =',ngrid
224	PRINT*,'klon =',klon
225	ENDIF
226	!
227	write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
228	do ig=1,klon
229	if (prt_level.ge.20) then
230	print*,'th_main ig f0',ig,f0(ig)
231	endif
232	f0(ig)=max(f0(ig),1.e-2)
233	zmax0(ig)=max(zmax0(ig),40.)
234	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
235	enddo
236
237	!-----------------------------------------------------------------------
238	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
239	! --------------------------------------------------------------------
240	!
241	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
242	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
243
244	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
245
246	!------------------------------------------------------------------------
247	! --------------------
248	!
249	!
250	! + + + + + + + + + + +
251	!
252	!
253	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
254	! wh,wt,wo ...
255	!
256	! + + + + + + + + + + + zh,zu,zv,zo,rho
257	!
258	!
259	! -------------------- zlev(1)
260	! \\\\\\\\\\\\\\\\\\\\
261	!
262	!
263
264	!-----------------------------------------------------------------------
265	! Calcul des altitudes des couches
266	!-----------------------------------------------------------------------
267
268	do l=2,nlay
269	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
270	enddo
271	zlev(:,1)=0.
272	zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
273	do l=1,nlay
274	zlay(:,l)=pphi(:,l)/RG
275	enddo
276	!calcul de l epaisseur des couches
277	do l=1,nlay
278	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
279	enddo
280
281	! print*,'2 OK convect8'
282	!-----------------------------------------------------------------------
283	! Calcul des densites
284	!-----------------------------------------------------------------------
285
286	do l=1,nlay
287	rho(:,l)=pplay(:,l)/(zpspsk(:,l)RDztv(:,l))
288	enddo
289
290	!IM
291	if (prt_level.ge.10)write(lunout,*) &
292	& 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
293	rhobarz(:,1)=rho(:,1)
294
295	do l=2,nlay
296	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
297	enddo
298
299	!calcul de la masse
300	do l=1,nlay
301	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
302	enddo
303
304	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
305
306	!------------------------------------------------------------------
307	!
308	! /\|\
309	! -------- \| F_k+1 -------
310	! ----> D_k
311	! /\|\ <---- E_k , A_k
312	! -------- \| F_k ---------
313	! ----> D_k-1
314	! <---- E_k-1 , A_k-1
315	!
316	!
317	!
318	!
319	!
320	! ---------------------------
321	!
322	! ----- F_lmax+1=0 ---------- \
323	! lmax (zmax) \|
324	! --------------------------- \|
325	! \|
326	! --------------------------- \|
327	! \|
328	! --------------------------- \|
329	! \|
330	! --------------------------- \|
331	! \|
332	! --------------------------- \|
333	! \| E
334	! --------------------------- \| D
335	! \|
336	! --------------------------- \|
337	! \|
338	! --------------------------- \ \|
339	! lalim \| \|
340	! --------------------------- \| \|
341	! \| \|
342	! --------------------------- \| \|
343	! \| A \|
344	! --------------------------- \| \|
345	! \| \|
346	! --------------------------- \| \|
347	! lmin (=1 pour le moment) \| \|
348	! ----- F_lmin=0 ------------ / /
349	!
350	! ---------------------------
351	! //////////////////////////
352	!
353	!
354	!=============================================================================
355	! Calculs initiaux ne faisant pas intervenir les changements de phase
356	!=============================================================================
357
358	!------------------------------------------------------------------
359	! 1. alim_star est le profil vertical de l'alimentation a la base du
360	! panache thermique, calcule a partir de la flotabilite de l'air sec
361	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
362	!------------------------------------------------------------------
363	!
364	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
365	lmin=1
366
367	!-----------------------------------------------------------------------------
368	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
369	! panache sec conservatif (e=d=0) alimente selon alim_star
370	! Il s'agit d'un calcul de type CAPE
371	! zmax_sec est utilise pour determiner la geometrie du thermique.
372	!------------------------------------------------------------------------------
373	!---------------------------------------------------------------------------------
374	!calcul du melange et des variables dans le thermique
375	!--------------------------------------------------------------------------------
376	!
377	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
378	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
379	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
380	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
381	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
382	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
383	& ,lev_out,lunout1,igout)
384	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
385
386	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
387	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
388
389	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
390	if (prt_level.ge.10) then
391	write(lunout1,*) 'Dans thermcell_main 2'
392	write(lunout1,*) 'lmin ',lmin(igout)
393	write(lunout1,*) 'lalim ',lalim(igout)
394	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
395	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
396	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
397	endif
398
399	!-------------------------------------------------------------------------------
400	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
401	!-------------------------------------------------------------------------------
402	!
403	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
404	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
405
406
407
408	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
409	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
410	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
411	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
412
413	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
414
415	!-------------------------------------------------------------------------------
416	! Fermeture,determination de f
417	!-------------------------------------------------------------------------------
418	!
419	!
420	write(lunout,*)'THERM NOUVEAU RIO2009 31B'
421	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
422	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
423
424	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
425	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
426
427	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
428	if (prt_level.ge.10) then
429	write(lunout1,*) 'Dans thermcell_main 1b'
430	write(lunout1,*) 'lmin ',lmin(igout)
431	write(lunout1,*) 'lalim ',lalim(igout)
432	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
433	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
434	& ,l=1,lalim(igout)+4)
435	endif
436
437
438
439	print*,'THERM 26JJJ'
440
441	! Choix de la fonction d'alimentation utilisee pour la fermeture.
442	! Apparemment sans importance
443	alim_star_clos(:,:)=alim_star(:,:)
444	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
445
446	! Appel avec la version seche
447	CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
448	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
449
450	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
451	! Appel avec les zmax et wmax tenant compte de la condensation
452	! Semble moins bien marcher
453	! CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
454	! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
455	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
456
457	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
458
459	if (tau_thermals>1.) then
460	lambda=exp(-ptimestep/tau_thermals)
461	f0=(1.-lambda)f+lambdaf0
462	else
463	f0=f
464	endif
465
466	! Test valable seulement en 1D mais pas genant
467	if (.not. (f0(1).ge.0.) ) then
468	stop'Dans thermcell_main'
469	endif
470
471	!-------------------------------------------------------------------------------
472	!deduction des flux
473	!-------------------------------------------------------------------------------
474
475	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
476	& lalim,lmax,alim_star, &
477	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
478	& detr,zqla,lev_out,lunout1,igout)
479	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
480
481	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
482	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
483	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
484
485	!------------------------------------------------------------------
486	! On ne prend pas directement les profils issus des calculs precedents
487	! mais on s'autorise genereusement une relaxation vers ceci avec
488	! une constante de temps tau_thermals (typiquement 1800s).
489	!------------------------------------------------------------------
490
491	if (tau_thermals>1.) then
492	lambda=exp(-ptimestep/tau_thermals)
493	fm0=(1.-lambda)fm+lambdafm0
494	entr0=(1.-lambda)entr+lambdaentr0
495	detr0=(1.-lambda)detr+lambdadetr0
496	else
497	fm0=fm
498	entr0=entr
499	detr0=detr
500	endif
501
502	!c------------------------------------------------------------------
503	! calcul du transport vertical
504	!------------------------------------------------------------------
505
506	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
507	& zthl,zdthladj,zta,lev_out)
508	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
509	& po,pdoadj,zoa,lev_out)
510
511	!------------------------------------------------------------------
512	! Calcul de la fraction de l'ascendance
513	!------------------------------------------------------------------
514	do ig=1,klon
515	fraca(ig,1)=0.
516	fraca(ig,nlay+1)=0.
517	enddo
518	do l=2,nlay
519	do ig=1,klon
520	if (zw2(ig,l).gt.1.e-10) then
521	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
522	else
523	fraca(ig,l)=0.
524	endif
525	enddo
526	enddo
527
528	!------------------------------------------------------------------
529	! calcul du transport vertical du moment horizontal
530	!------------------------------------------------------------------
531
532	!IM 090508
533	if (1.eq.1) then
534	!IM 070508 vers. _dq
535	! if (1.eq.0) then
536
537
538	! Calcul du transport de V tenant compte d'echange par gradient
539	! de pression horizontal avec l'environnement
540
541	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
542	& ,fraca,zmax &
543	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
544
545	else
546
547	! calcul purement conservatif pour le transport de V
548	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
549	& ,zu,pduadj,zua,lev_out)
550	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
551	& ,zv,pdvadj,zva,lev_out)
552	endif
553
554	! print*,'13 OK convect8'
555	do l=1,nlay
556	do ig=1,ngrid
557	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
558	enddo
559	enddo
560
561	if (prt_level.ge.1) print*,'14 OK convect8'
562	!------------------------------------------------------------------
563	! Calculs de diagnostiques pour les sorties
564	!------------------------------------------------------------------
565	!calcul de fraca pour les sorties
566
567	if (sorties) then
568	if (prt_level.ge.1) print*,'14a OK convect8'
569	! calcul du niveau de condensation
570	! initialisation
571	do ig=1,ngrid
572	nivcon(ig)=0
573	zcon(ig)=0.
574	enddo
575	!nouveau calcul
576	do ig=1,ngrid
577	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
578	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
579	enddo
580	do k=1,nlay
581	do ig=1,ngrid
582	if ((pcon(ig).le.pplay(ig,k)) &
583	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
584	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
585	endif
586	enddo
587	enddo
588	if (prt_level.ge.1) print*,'14b OK convect8'
589	do k=nlay,1,-1
590	do ig=1,ngrid
591	if (zqla(ig,k).gt.1e-10) then
592	nivcon(ig)=k
593	zcon(ig)=zlev(ig,k)
594	endif
595	enddo
596	enddo
597	if (prt_level.ge.1) print*,'14c OK convect8'
598	!calcul des moments
599	!initialisation
600	do l=1,nlay
601	do ig=1,ngrid
602	q2(ig,l)=0.
603	wth2(ig,l)=0.
604	wth3(ig,l)=0.
605	ratqscth(ig,l)=0.
606	ratqsdiff(ig,l)=0.
607	enddo
608	enddo
609	if (prt_level.ge.1) print*,'14d OK convect8'
610	if (prt_level.ge.10)write(lunout,*) &
611	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
612	do l=1,nlay
613	do ig=1,ngrid
614	zf=fraca(ig,l)
615	zf2=zf/(1.-zf)
616	!
617	if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
618	!
619	if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
620	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
621	if(zw2(ig,l).gt.1.e-10) then
622	wth2(ig,l)=zf2(zw2(ig,l))*2
623	else
624	wth2(ig,l)=0.
625	endif
626	! print*,'wth2=',wth2(ig,l)
627	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
628	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
629	if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
630	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
631	!test: on calcul q2/po=ratqsc
632	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
633	enddo
634	enddo
635	!calcul des flux: q, thetal et thetav
636	do l=1,nlay
637	do ig=1,ngrid
638	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
639	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
640	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
641	enddo
642	enddo
643	!calcul de ale_bl et alp_bl
644	!pour le calcul d'une valeur integree entre la surface et lmax
645	do ig=1,ngrid
646	alp_int(ig)=0.
647	ale_int(ig)=0.
648	n_int(ig)=0
649	enddo
650	!
651	do l=1,nlay
652	do ig=1,ngrid
653	if(l.LE.lmax(ig)) THEN
654	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)
655	ale_int(ig)=ale_int(ig)+0.5zw2(ig,l)*2
656	n_int(ig)=n_int(ig)+1
657	endif
658	enddo
659	enddo
660	! print*,'avant calcul ale et alp'
661	!calcul de ALE et ALP pour la convection
662	do ig=1,ngrid
663	! Alp_bl(ig)=0.5rhobarz(ig,lmix_bis(ig))wth3(ig,lmix(ig))
664	! Alp_bl(ig)=0.5rhobarz(ig,nivcon(ig))wth3(ig,nivcon(ig))
665	! Alp_bl(ig)=0.5rhobarz(ig,lmix(ig))wth3(ig,lmix(ig))
666	! & *0.1
667	!valeur integree de alp_bl * 0.5:
668	if (n_int(ig).gt.0) then
669	Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
670	! if (Alp_bl(ig).lt.0.) then
671	! Alp_bl(ig)=0.
672	endif
673	! endif
674	! write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
675	! s wth3(ig,nivcon(ig)),Alp_bl(ig)
676	! write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
677	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
678	! if (nivcon(ig).eq.1) then
679	! Ale_bl(ig)=0.
680	! else
681	!valeur max de ale_bl:
682	Ale_bl(ig)=0.5zw2(ig,lmix(ig))*2
683	! & /2.
684	! & *0.1
685	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
686	! if (n_int(ig).gt.0) then
687	! Ale_bl(ig)=ale_int(ig)/n_int(ig)
688	! Ale_bl(ig)=4.
689	! endif
690	! endif
691	! Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
692	! Ale_bl(ig)=wth2(ig,nivcon(ig))
693	! write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
694	enddo
695	!test:calcul de la ponderation des couches pour KE
696	!initialisations
697	! print*,'ponderation'
698	do ig=1,ngrid
699	fm_tot(ig)=0.
700	enddo
701	do ig=1,ngrid
702	do k=1,klev
703	wght_th(ig,k)=1.
704	enddo
705	enddo
706	do ig=1,ngrid
707	! lalim_conv(ig)=lmix_bis(ig)
708	!la hauteur de la couche alim_conv = hauteur couche alim_therm
709	lalim_conv(ig)=lalim(ig)
710	! zentr(ig)=zlev(ig,lalim(ig))
711	enddo
712	do ig=1,ngrid
713	do k=1,lalim_conv(ig)
714	fm_tot(ig)=fm_tot(ig)+fm(ig,k)
715	enddo
716	enddo
717	do ig=1,ngrid
718	do k=1,lalim_conv(ig)
719	if (fm_tot(ig).gt.1.e-10) then
720	! wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
721	endif
722	!on pondere chaque couche par a*
723	if (alim_star(ig,k).gt.1.e-10) then
724	wght_th(ig,k)=alim_star(ig,k)
725	else
726	wght_th(ig,k)=1.
727	endif
728	enddo
729	enddo
730	! print*,'apres wght_th'
731	!test pour prolonger la convection
732	do ig=1,ngrid
733	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
734	if ((alim_star(ig,1).lt.1.e-10)) then
735	lalim_conv(ig)=1
736	wght_th(ig,1)=1.
737	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
738	endif
739	enddo
740
741	!calcul du ratqscdiff
742	if (prt_level.ge.1) print*,'14e OK convect8'
743	var=0.
744	vardiff=0.
745	ratqsdiff(:,:)=0.
746	do ig=1,ngrid
747	do l=1,lalim(ig)
748	var=var+alim_star(ig,l)zqta(ig,l)1000.
749	enddo
750	enddo
751	if (prt_level.ge.1) print*,'14f OK convect8'
752	do ig=1,ngrid
753	do l=1,lalim(ig)
754	zf=fraca(ig,l)
755	zf2=zf/(1.-zf)
756	vardiff=vardiff+alim_star(ig,l) &
757	& (zqta(ig,l)1000.-var)**2
758	! ratqsdiff=ratqsdiff+alim_star(ig,l)*
759	! s (zqta(ig,l)1000.-po(ig,l)1000.)**2
760	enddo
761	enddo
762	if (prt_level.ge.1) print*,'14g OK convect8'
763	do l=1,nlay
764	do ig=1,ngrid
765	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
766	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
767	enddo
768	enddo
769	!--------------------------------------------------------------------
770	!
771	!ecriture des fichiers sortie
772	! print*,'15 OK convect8'
773
774	#ifdef wrgrads_thermcell
775	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
776	#include "thermcell_out3d.h"
777	#endif
778
779	endif
780
781	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
782
783	! if(icount.eq.501) stop'au pas 301 dans thermcell_main'
784	return
785	end
786
787	!-----------------------------------------------------------------------------
788
789	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
790	IMPLICIT NONE
791	#include "iniprint.h"
792
793	integer i, k, klon,klev
794	real pplev(klon,klev+1),pplay(klon,klev)
795	real ztv(klon,klev)
796	real po(klon,klev)
797	real ztva(klon,klev)
798	real zqla(klon,klev)
799	real f_star(klon,klev)
800	real zw2(klon,klev)
801	integer long(klon)
802	real seuil
803	character*21 comment
804
805	if (prt_level.ge.1) THEN
806	print*,'WARNING !!! TEST ',comment
807	endif
808	return
809
810	! test sur la hauteur des thermiques ...
811	do i=1,klon
812	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
813	if (prt_level.ge.10) then
814	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
815	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
816	do k=1,klev
817	write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000po(i,k),ztva(i,k),1000zqla(i,k),f_star(i,k),zw2(i,k)
818	enddo
819	! stop
820	endif
821	enddo
822
823
824	return
825	end
826

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

Download in other formats:

Original Format