Context Navigation

thermcell_main.F90 @ 3722

Last change on this file since 3722 was 1399, checked in by idelkadi, 14 years ago

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

physiq.F : appel modifie a calltherm et fisrtilp
thermcell_main.F90 : changement dans les flag_thermals_ed
thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

cloudth.F90 : PDF d'eau bi-gaussiennes

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

File size: 27.0 KB

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

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

Download in other formats:

Original Format