Context Navigation

thermcell_main.F90 @ 1123

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

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

Download in other formats:

Original Format