Context Navigation

thermcell_main.F90 @ 1015

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

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

Download in other formats:

Original Format