Context Navigation

thermcell_main.F90 @ 982

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

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

Download in other formats:

Original Format