Context Navigation

thermcell_main.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 27.0 KB

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