Context Navigation

thermcell_main.F90 @ 1387

Last change on this file since 1387 was 1146, checked in by Laurent Fairhead, 15 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.0 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	if (prt_level.ge.10)write(lunout,*) &
235	& 'WARNING thermcell_main f0=max(f0,1.e-2)'
236	do ig=1,klon
237	if (prt_level.ge.20) then
238	print*,'th_main ig f0',ig,f0(ig)
239	endif
240	f0(ig)=max(f0(ig),1.e-2)
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 à la base du
367	! panache thermique, calculé à partir de la flotabilité 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	CALL thermcell_init(ngrid,nlay,ztv,zlay,zlev, &
373	& lalim,lmin,alim_star,alim_star_tot,lev_out)
374
375	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin ')
376	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
377
378
379	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
380	if (prt_level.ge.10) then
381	write(lunout1,*) 'Dans thermcell_main 1'
382	write(lunout1,*) 'lmin ',lmin(igout)
383	write(lunout1,*) 'lalim ',lalim(igout)
384	write(lunout1,*) ' ig l alim_star thetav'
385	write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
386	& ,ztv(igout,l),l=1,lalim(igout)+4)
387	endif
388
389	!v1d do ig=1,klon
390	!v1d if (alim_star(ig,1).gt.1.e-10) then
391	!v1d therm=.true.
392	!v1d endif
393	!v1d enddo
394	!-----------------------------------------------------------------------------
395	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
396	! panache sec conservatif (e=d=0) alimente selon alim_star
397	! Il s'agit d'un calcul de type CAPE
398	! zmax_sec est utilisé pour déterminer la géométrie du thermique.
399	!------------------------------------------------------------------------------
400	!
401	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
402	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
403
404	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
405	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
406
407	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
408	if (prt_level.ge.10) then
409	write(lunout1,*) 'Dans thermcell_main 1b'
410	write(lunout1,*) 'lmin ',lmin(igout)
411	write(lunout1,*) 'lalim ',lalim(igout)
412	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
413	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
414	& ,l=1,lalim(igout)+4)
415	endif
416
417
418
419	!---------------------------------------------------------------------------------
420	!calcul du melange et des variables dans le thermique
421	!--------------------------------------------------------------------------------
422	!
423	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
424	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
425	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
426	& zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot, &
427	& lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &
428	& ztla,zqla,zqta,zha,zw2,zw_est,zqsatth,lmix,lmix_bis,linter &
429	& ,lev_out,lunout1,igout)
430	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
431
432	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
433	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
434
435	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
436	if (prt_level.ge.10) then
437	write(lunout1,*) 'Dans thermcell_main 2'
438	write(lunout1,*) 'lmin ',lmin(igout)
439	write(lunout1,*) 'lalim ',lalim(igout)
440	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
441	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
442	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
443	endif
444
445	!-------------------------------------------------------------------------------
446	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
447	!-------------------------------------------------------------------------------
448	!
449	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
450	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
451
452
453	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
454	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
455	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
456	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
457
458	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
459
460	!-------------------------------------------------------------------------------
461	! Fermeture,determination de f
462	!-------------------------------------------------------------------------------
463	!
464	!avant closure: on redéfinit lalim, alim_star_tot et alim_star
465	! do ig=1,klon
466	! do l=2,lalim(ig)
467	! alim_star(ig,l)=entr_star(ig,l)
468	! entr_star(ig,l)=0.
469	! enddo
470	! enddo
471
472	CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
473	& zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
474
475	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
476
477	if (tau_thermals>1.) then
478	lambda=exp(-ptimestep/tau_thermals)
479	f0=(1.-lambda)f+lambdaf0
480	else
481	f0=f
482	endif
483
484	! Test valable seulement en 1D mais pas genant
485	if (.not. (f0(1).ge.0.) ) then
486	stop 'Dans thermcell_main'
487	endif
488
489	!-------------------------------------------------------------------------------
490	!deduction des flux
491	!-------------------------------------------------------------------------------
492
493	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
494	& lalim,lmax,alim_star, &
495	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
496	& detr,zqla,lev_out,lunout1,igout)
497	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
498
499	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
500	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
501	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
502
503	!------------------------------------------------------------------
504	! On ne prend pas directement les profils issus des calculs precedents
505	! mais on s'autorise genereusement une relaxation vers ceci avec
506	! une constante de temps tau_thermals (typiquement 1800s).
507	!------------------------------------------------------------------
508
509	if (tau_thermals>1.) then
510	lambda=exp(-ptimestep/tau_thermals)
511	fm0=(1.-lambda)fm+lambdafm0
512	entr0=(1.-lambda)entr+lambdaentr0
513	! detr0=(1.-lambda)detr+lambdadetr0
514	else
515	fm0=fm
516	entr0=entr
517	detr0=detr
518	endif
519
520	!c------------------------------------------------------------------
521	! calcul du transport vertical
522	!------------------------------------------------------------------
523
524	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
525	& zthl,zdthladj,zta,lev_out)
526	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
527	& po,pdoadj,zoa,lev_out)
528
529	!------------------------------------------------------------------
530	! Calcul de la fraction de l'ascendance
531	!------------------------------------------------------------------
532	do ig=1,klon
533	fraca(ig,1)=0.
534	fraca(ig,nlay+1)=0.
535	enddo
536	do l=2,nlay
537	do ig=1,klon
538	if (zw2(ig,l).gt.1.e-10) then
539	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
540	else
541	fraca(ig,l)=0.
542	endif
543	enddo
544	enddo
545
546	!------------------------------------------------------------------
547	! calcul du transport vertical du moment horizontal
548	!------------------------------------------------------------------
549
550	!IM 090508
551	if (1.eq.1) then
552	!IM 070508 vers. _dq
553	! if (1.eq.0) then
554
555
556	! Calcul du transport de V tenant compte d'echange par gradient
557	! de pression horizontal avec l'environnement
558
559	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
560	& ,fraca,zmax &
561	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
562	!IM 050508 & ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
563	else
564
565	! calcul purement conservatif pour le transport de V
566	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
567	& ,zu,pduadj,zua,lev_out)
568	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
569	& ,zv,pdvadj,zva,lev_out)
570	endif
571
572	! print*,'13 OK convect8'
573	do l=1,nlay
574	do ig=1,ngrid
575	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
576	enddo
577	enddo
578
579	if (prt_level.ge.1) print*,'14 OK convect8'
580	!------------------------------------------------------------------
581	! Calculs de diagnostiques pour les sorties
582	!------------------------------------------------------------------
583	!calcul de fraca pour les sorties
584
585	if (sorties) then
586	if (prt_level.ge.1) print*,'14a OK convect8'
587	! calcul du niveau de condensation
588	! initialisation
589	do ig=1,ngrid
590	nivcon(ig)=0
591	zcon(ig)=0.
592	enddo
593	!nouveau calcul
594	do ig=1,ngrid
595	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
596	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
597	enddo
598	do k=1,nlay
599	do ig=1,ngrid
600	if ((pcon(ig).le.pplay(ig,k)) &
601	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
602	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
603	endif
604	enddo
605	enddo
606	if (prt_level.ge.1) print*,'14b OK convect8'
607	do k=nlay,1,-1
608	do ig=1,ngrid
609	if (zqla(ig,k).gt.1e-10) then
610	nivcon(ig)=k
611	zcon(ig)=zlev(ig,k)
612	endif
613	enddo
614	enddo
615	if (prt_level.ge.1) print*,'14c OK convect8'
616	!calcul des moments
617	!initialisation
618	do l=1,nlay
619	do ig=1,ngrid
620	q2(ig,l)=0.
621	wth2(ig,l)=0.
622	wth3(ig,l)=0.
623	ratqscth(ig,l)=0.
624	ratqsdiff(ig,l)=0.
625	enddo
626	enddo
627	if (prt_level.ge.1) print*,'14d OK convect8'
628	if (prt_level.ge.10)write(lunout,*) &
629	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
630	do l=1,nlay
631	do ig=1,ngrid
632	zf=fraca(ig,l)
633	zf2=zf/(1.-zf)
634	!
635	if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
636	!
637	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)
638	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))*2
639	if(zw2(ig,l).gt.1.e-10) then
640	wth2(ig,l)=zf2(zw2(ig,l))*2
641	else
642	wth2(ig,l)=0.
643	endif
644	! print*,'wth2=',wth2(ig,l)
645	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
646	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
647	if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
648	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
649	!test: on calcul q2/po=ratqsc
650	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
651	enddo
652	enddo
653	!calcul de ale_bl et alp_bl
654	!pour le calcul d'une valeur intégrée entre la surface et lmax
655	do ig=1,ngrid
656	alp_int(ig)=0.
657	ale_int(ig)=0.
658	n_int(ig)=0
659	enddo
660	!
661	do l=1,nlay
662	do ig=1,ngrid
663	if(l.LE.lmax(ig)) THEN
664	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)
665	ale_int(ig)=ale_int(ig)+0.5zw2(ig,l)*2
666	n_int(ig)=n_int(ig)+1
667	endif
668	enddo
669	enddo
670	! print*,'avant calcul ale et alp'
671	!calcul de ALE et ALP pour la convection
672	do ig=1,ngrid
673	! Alp_bl(ig)=0.5rhobarz(ig,lmix_bis(ig))wth3(ig,lmix(ig))
674	! Alp_bl(ig)=0.5rhobarz(ig,nivcon(ig))wth3(ig,nivcon(ig))
675	! Alp_bl(ig)=0.5rhobarz(ig,lmix(ig))wth3(ig,lmix(ig))
676	! & *0.1
677	!valeur integree de alp_bl * 0.5:
678	if (n_int(ig).gt.0) then
679	Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
680	! if (Alp_bl(ig).lt.0.) then
681	! Alp_bl(ig)=0.
682	endif
683	! endif
684	! write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
685	! s wth3(ig,nivcon(ig)),Alp_bl(ig)
686	! write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
687	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
688	! if (nivcon(ig).eq.1) then
689	! Ale_bl(ig)=0.
690	! else
691	!valeur max de ale_bl:
692	Ale_bl(ig)=0.5zw2(ig,lmix(ig))*2
693	! & /2.
694	! & *0.1
695	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
696	! if (n_int(ig).gt.0) then
697	! Ale_bl(ig)=ale_int(ig)/n_int(ig)
698	! Ale_bl(ig)=4.
699	! endif
700	! endif
701	! Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
702	! Ale_bl(ig)=wth2(ig,nivcon(ig))
703	! write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
704	enddo
705	!test:calcul de la ponderation des couches pour KE
706	!initialisations
707	! print*,'ponderation'
708	do ig=1,ngrid
709	fm_tot(ig)=0.
710	enddo
711	do ig=1,ngrid
712	do k=1,klev
713	wght_th(ig,k)=1.
714	enddo
715	enddo
716	do ig=1,ngrid
717	! lalim_conv(ig)=lmix_bis(ig)
718	!la hauteur de la couche alim_conv = hauteur couche alim_therm
719	lalim_conv(ig)=lalim(ig)
720	! zentr(ig)=zlev(ig,lalim(ig))
721	enddo
722	do ig=1,ngrid
723	do k=1,lalim_conv(ig)
724	fm_tot(ig)=fm_tot(ig)+fm(ig,k)
725	enddo
726	enddo
727	do ig=1,ngrid
728	do k=1,lalim_conv(ig)
729	if (fm_tot(ig).gt.1.e-10) then
730	! wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
731	endif
732	!on pondere chaque couche par a*
733	if (alim_star(ig,k).gt.1.e-10) then
734	wght_th(ig,k)=alim_star(ig,k)
735	else
736	wght_th(ig,k)=1.
737	endif
738	enddo
739	enddo
740	! print*,'apres wght_th'
741	!test pour prolonger la convection
742	do ig=1,ngrid
743	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
744	if ((alim_star(ig,1).lt.1.e-10)) then
745	lalim_conv(ig)=1
746	wght_th(ig,1)=1.
747	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
748	endif
749	enddo
750
751	!calcul du ratqscdiff
752	if (prt_level.ge.1) print*,'14e OK convect8'
753	var=0.
754	vardiff=0.
755	ratqsdiff(:,:)=0.
756	do ig=1,ngrid
757	do l=1,lalim(ig)
758	var=var+alim_star(ig,l)zqta(ig,l)1000.
759	enddo
760	enddo
761	if (prt_level.ge.1) print*,'14f OK convect8'
762	do ig=1,ngrid
763	do l=1,lalim(ig)
764	zf=fraca(ig,l)
765	zf2=zf/(1.-zf)
766	vardiff=vardiff+alim_star(ig,l) &
767	& (zqta(ig,l)1000.-var)**2
768	! ratqsdiff=ratqsdiff+alim_star(ig,l)*
769	! s (zqta(ig,l)1000.-po(ig,l)1000.)**2
770	enddo
771	enddo
772	if (prt_level.ge.1) print*,'14g OK convect8'
773	do l=1,nlay
774	do ig=1,ngrid
775	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
776	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
777	enddo
778	enddo
779	!--------------------------------------------------------------------
780	!
781	!ecriture des fichiers sortie
782	! print*,'15 OK convect8'
783
784	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
785	#ifdef wrgrads_thermcell
786	#include "thermcell_out3d.h"
787	#endif
788
789	endif
790
791	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
792
793	! if(icount.eq.501) stop'au pas 301 dans thermcell_main'
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	! stop
830	endif
831	enddo
832
833
834	return
835	end
836

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

Download in other formats:

Original Format