Context Navigation

thermcell_main.F90 @ 5425

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

Modifications to thermals for TKE transport

Modifications aux thermiques pour le transport de la TKE

pbl_surface_mode.F90 : ok_flux_surf=.false. seulement pour klon>1
physiq.F : option iflag_pbl=10 pour transporter la TKE avec les thermiques.
calltherm.F90 : passage de iflag_thermals_ed en argument pour thermcell_plume
thermcell_main.F90 : Appel a plusieurs version de thermcell_plume en option
thermcell_plume.F90 : plusieurs versions dans le meme fichier (temporaire)
thermcell_height.F90 : verrue pour les cas ou les thermiques montent tout

en haut

yamada4 : inclusion de la diffusion verticale en option iflag_pbl=9

+ variables anciennement common, puis save/allocatable, remises en local

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 27.3 KB

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

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

Download in other formats:

Original Format