Context Navigation

thermcell_main.F90 @ 1645

Last change on this file since 1645 was 1638, checked in by idelkadi, 12 years ago
Introduction du declenchement stochastique de la convection
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 40.1 KB

Line
1	!
2	! $Id: thermcell_main.F90 1638 2012-07-23 11:11:05Z lguez $
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	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
11	& ,zmax0, f0,zw2,fraca,ztv &
12	& ,zpspsk,ztla,zthl &
13	!!! nrlmd le 10/04/2012
14	& ,pbl_tke,pctsrf,omega,airephy &
15	& ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
16	& ,n2,s2,ale_bl_stat &
17	& ,therm_tke_max,env_tke_max &
18	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
19	& ,alp_bl_conv,alp_bl_stat &
20	!!! fin nrlmd le 10/04/2012
21	& )
22
23	USE dimphy
24	USE comgeomphy , ONLY:rlond,rlatd
25	IMPLICIT NONE
26
27	!=======================================================================
28	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
29	! Version du 09.02.07
30	! Calcul du transport vertical dans la couche limite en presence
31	! de "thermiques" explicitement representes avec processus nuageux
32	!
33	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
34	!
35	! le thermique est suppose homogene et dissipe par melange avec
36	! son environnement. la longueur l_mix controle l'efficacite du
37	! melange
38	!
39	! Le calcul du transport des differentes especes se fait en prenant
40	! en compte:
41	! 1. un flux de masse montant
42	! 2. un flux de masse descendant
43	! 3. un entrainement
44	! 4. un detrainement
45	!
46	!=======================================================================
47
48	!-----------------------------------------------------------------------
49	! declarations:
50	! -------------
51
52	#include "dimensions.h"
53	#include "YOMCST.h"
54	#include "YOETHF.h"
55	#include "FCTTRE.h"
56	#include "iniprint.h"
57	#include "thermcell.h"
58	!!! nrlmd le 10/04/2012
59	#include "indicesol.h"
60	!!! fin nrlmd le 10/04/2012
61
62	! arguments:
63	! ----------
64
65	!IM 140508
66	INTEGER itap
67
68	INTEGER ngrid,nlay
69	real ptimestep
70	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
71	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
72	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
73	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
74	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
75	real pphi(ngrid,nlay)
76
77	! local:
78	! ------
79
80	integer icount
81	data icount/0/
82	save icount
83	!$OMP THREADPRIVATE(icount)
84
85	integer,save :: igout=1
86	!$OMP THREADPRIVATE(igout)
87	integer,save :: lunout1=6
88	!$OMP THREADPRIVATE(lunout1)
89	integer,save :: lev_out=10
90	!$OMP THREADPRIVATE(lev_out)
91
92	REAL susqr2pi, Reuler
93
94	INTEGER ig,k,l,ll,ierr
95	real zsortie1d(klon)
96	INTEGER lmax(klon),lmin(klon),lalim(klon)
97	INTEGER lmix(klon)
98	INTEGER lmix_bis(klon)
99	real linter(klon)
100	real zmix(klon)
101	real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
102	! real fraca(klon,klev)
103
104	real zmax_sec(klon)
105	!on garde le zmax du pas de temps precedent
106	real zmax0(klon)
107	!FH/IM save zmax0
108
109	real lambda
110
111	real zlev(klon,klev+1),zlay(klon,klev)
112	real deltaz(klon,klev)
113	REAL zh(klon,klev)
114	real zthl(klon,klev),zdthladj(klon,klev)
115	REAL ztv(klon,klev)
116	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
117	real zl(klon,klev)
118	real zsortie(klon,klev)
119	real zva(klon,klev)
120	real zua(klon,klev)
121	real zoa(klon,klev)
122
123	real zta(klon,klev)
124	real zha(klon,klev)
125	real fraca(klon,klev+1)
126	real zf,zf2
127	real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
128	real q2(klon,klev)
129	! FH probleme de dimensionnement avec l'allocation dynamique
130	! common/comtherm/thetath2,wth2
131	real wq(klon,klev)
132	real wthl(klon,klev)
133	real wthv(klon,klev)
134
135	real ratqscth(klon,klev)
136	real var
137	real vardiff
138	real ratqsdiff(klon,klev)
139
140	logical sorties
141	real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
142	real zpspsk(klon,klev)
143
144	real wmax(klon)
145	real wmax_tmp(klon)
146	real wmax_sec(klon)
147	real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
148	real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
149
150	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
151	!niveau de condensation
152	integer nivcon(klon)
153	real zcon(klon)
154	REAL CHI
155	real zcon2(klon)
156	real pcon(klon)
157	real zqsat(klon,klev)
158	real zqsatth(klon,klev)
159
160	real f_star(klon,klev+1),entr_star(klon,klev)
161	real detr_star(klon,klev)
162	real alim_star_tot(klon)
163	real alim_star(klon,klev)
164	real alim_star_clos(klon,klev)
165	real f(klon), f0(klon)
166	!FH/IM save f0
167	real zlevinter(klon)
168	logical debut
169	real seuil
170	real csc(klon,klev)
171
172	!!! nrlmd le 10/04/2012
173
174	!------Entrées
175	real pbl_tke(klon,klev+1,nbsrf)
176	real pctsrf(klon,nbsrf)
177	real omega(klon,klev)
178	real airephy(klon)
179	!------Sorties
180	real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
181	real therm_tke_max0(klon),env_tke_max0(klon)
182	real n2(klon),s2(klon)
183	real ale_bl_stat(klon)
184	real therm_tke_max(klon,klev),env_tke_max(klon,klev)
185	real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
186	!------Local
187	integer nsrf
188	real rhobarz0(klon) ! Densité au LCL
189	logical ok_lcl(klon) ! Existence du LCL des thermiques
190	integer klcl(klon) ! Niveau du LCL
191	real interp(klon) ! Coef d'interpolation pour le LCL
192	!--Triggering
193	real Su ! Surface unité: celle d'un updraft élémentaire
194	parameter(Su=4e4)
195	real hcoef ! Coefficient directeur pour le calcul de s2
196	parameter(hcoef=1)
197	real hmincoef ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
198	parameter(hmincoef=0.3)
199	real eps1 ! Fraction de surface occupée par la population 1 : eps1=n1s1/(fraca0Sd)
200	parameter(eps1=0.3)
201	real hmin(ngrid) ! Ordonnée à l'origine pour le calcul de s2
202	real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
203	real zmax_moy_coef
204	parameter(zmax_moy_coef=0.33)
205	real depth(klon) ! Epaisseur moyenne du cumulus
206	real w_max(klon) ! Vitesse max statistique
207	real s_max(klon)
208	!--Closure
209	real pbl_tke_max(klon,klev) ! Profil de TKE moyenne
210	real pbl_tke_max0(klon) ! TKE moyenne au LCL
211	real w_ls(klon,klev) ! Vitesse verticale grande échelle (m/s)
212	real coef_m ! On considère un rendement pour alp_bl_fluct_m
213	parameter(coef_m=1.)
214	real coef_tke ! On considère un rendement pour alp_bl_fluct_tke
215	parameter(coef_tke=1.)
216
217	!!! fin nrlmd le 10/04/2012
218
219	!
220	!nouvelles variables pour la convection
221	real Ale_bl(klon)
222	real Alp_bl(klon)
223	real alp_int(klon),dp_int(klon),zdp
224	real ale_int(klon)
225	integer n_int(klon)
226	real fm_tot(klon)
227	real wght_th(klon,klev)
228	integer lalim_conv(klon)
229	!v1d logical therm
230	!v1d save therm
231
232	character*2 str2
233	character*10 str10
234
235	character (len=20) :: modname='thermcell_main'
236	character (len=80) :: abort_message
237
238	EXTERNAL SCOPY
239	!
240
241	!-----------------------------------------------------------------------
242	! initialisation:
243	! ---------------
244	!
245
246	seuil=0.25
247
248	if (debut) then
249	fm0=0.
250	entr0=0.
251	detr0=0.
252
253
254	#undef wrgrads_thermcell
255	#ifdef wrgrads_thermcell
256	! Initialisation des sorties grads pour les thermiques.
257	! Pour l'instant en 1D sur le point igout.
258	! Utilise par thermcell_out3d.h
259	str10='therm'
260	call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
261	& rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., &
262	& ptimestep,str10,'therm ')
263	#endif
264
265
266
267	endif
268
269	fm=0. ; entr=0. ; detr=0.
270
271
272	icount=icount+1
273
274	!IM 090508 beg
275	!print*,'====================================================================='
276	!print*,'====================================================================='
277	!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
278	!print*,'====================================================================='
279	!print*,'====================================================================='
280	!IM 090508 end
281
282	if (prt_level.ge.1) print*,'thermcell_main V4'
283
284	sorties=.true.
285	IF(ngrid.NE.klon) THEN
286	PRINT*
287	PRINT*,'STOP dans convadj'
288	PRINT*,'ngrid =',ngrid
289	PRINT*,'klon =',klon
290	ENDIF
291	!
292	! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
293	do ig=1,klon
294	f0(ig)=max(f0(ig),1.e-2)
295	zmax0(ig)=max(zmax0(ig),40.)
296	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
297	enddo
298
299	if (prt_level.ge.20) then
300	do ig=1,ngrid
301	print*,'th_main ig f0',ig,f0(ig)
302	enddo
303	endif
304	!-----------------------------------------------------------------------
305	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
306	! --------------------------------------------------------------------
307	!
308	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
309	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
310
311	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
312
313	!------------------------------------------------------------------------
314	! --------------------
315	!
316	!
317	! + + + + + + + + + + +
318	!
319	!
320	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
321	! wh,wt,wo ...
322	!
323	! + + + + + + + + + + + zh,zu,zv,zo,rho
324	!
325	!
326	! -------------------- zlev(1)
327	! \\\\\\\\\\\\\\\\\\\\
328	!
329	!
330
331	!-----------------------------------------------------------------------
332	! Calcul des altitudes des couches
333	!-----------------------------------------------------------------------
334
335	do l=2,nlay
336	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
337	enddo
338	zlev(:,1)=0.
339	zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
340	do l=1,nlay
341	zlay(:,l)=pphi(:,l)/RG
342	enddo
343	!calcul de l epaisseur des couches
344	do l=1,nlay
345	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
346	enddo
347
348	! print*,'2 OK convect8'
349	!-----------------------------------------------------------------------
350	! Calcul des densites
351	!-----------------------------------------------------------------------
352
353	rho(:,:)=pplay(:,:)/(zpspsk(:,:)RDztv(:,:))
354
355	if (prt_level.ge.10)write(lunout,*) &
356	& 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
357	rhobarz(:,1)=rho(:,1)
358
359	do l=2,nlay
360	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
361	enddo
362
363	!calcul de la masse
364	do l=1,nlay
365	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
366	enddo
367
368	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
369
370	!------------------------------------------------------------------
371	!
372	! /\|\
373	! -------- \| F_k+1 -------
374	! ----> D_k
375	! /\|\ <---- E_k , A_k
376	! -------- \| F_k ---------
377	! ----> D_k-1
378	! <---- E_k-1 , A_k-1
379	!
380	!
381	!
382	!
383	!
384	! ---------------------------
385	!
386	! ----- F_lmax+1=0 ---------- \
387	! lmax (zmax) \|
388	! --------------------------- \|
389	! \|
390	! --------------------------- \|
391	! \|
392	! --------------------------- \|
393	! \|
394	! --------------------------- \|
395	! \|
396	! --------------------------- \|
397	! \| E
398	! --------------------------- \| D
399	! \|
400	! --------------------------- \|
401	! \|
402	! --------------------------- \ \|
403	! lalim \| \|
404	! --------------------------- \| \|
405	! \| \|
406	! --------------------------- \| \|
407	! \| A \|
408	! --------------------------- \| \|
409	! \| \|
410	! --------------------------- \| \|
411	! lmin (=1 pour le moment) \| \|
412	! ----- F_lmin=0 ------------ / /
413	!
414	! ---------------------------
415	! //////////////////////////
416	!
417	!
418	!=============================================================================
419	! Calculs initiaux ne faisant pas intervenir les changements de phase
420	!=============================================================================
421
422	!------------------------------------------------------------------
423	! 1. alim_star est le profil vertical de l'alimentation a la base du
424	! panache thermique, calcule a partir de la flotabilite de l'air sec
425	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
426	!------------------------------------------------------------------
427	!
428	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
429	lmin=1
430
431	!-----------------------------------------------------------------------------
432	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
433	! panache sec conservatif (e=d=0) alimente selon alim_star
434	! Il s'agit d'un calcul de type CAPE
435	! zmax_sec est utilise pour determiner la geometrie du thermique.
436	!------------------------------------------------------------------------------
437	!---------------------------------------------------------------------------------
438	!calcul du melange et des variables dans le thermique
439	!--------------------------------------------------------------------------------
440	!
441	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
442	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
443
444	! Gestion temporaire de plusieurs appels à thermcell_plume au travers
445	! de la variable iflag_thermals
446
447	! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
448	if (iflag_thermals_ed<=9) then
449	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
450	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
451	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
452	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
453	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
454	& ,lev_out,lunout1,igout)
455
456	elseif (iflag_thermals_ed>9) then
457	! print*,'THERM RIO et al 2010, version d Arnaud'
458	CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
459	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
460	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
461	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
462	& ,lev_out,lunout1,igout)
463
464	endif
465
466	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
467
468	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
469	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
470
471	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
472	if (prt_level.ge.10) then
473	write(lunout1,*) 'Dans thermcell_main 2'
474	write(lunout1,*) 'lmin ',lmin(igout)
475	write(lunout1,*) 'lalim ',lalim(igout)
476	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
477	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
478	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
479	endif
480
481	!-------------------------------------------------------------------------------
482	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
483	!-------------------------------------------------------------------------------
484	!
485	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
486	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
487	! Attention, w2 est transforme en sa racine carree dans cette routine
488	! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
489	wmax_tmp=0.
490	do l=1,nlay
491	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
492	enddo
493	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
494
495
496
497	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
498	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
499	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
500	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
501
502	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
503
504	!-------------------------------------------------------------------------------
505	! Fermeture,determination de f
506	!-------------------------------------------------------------------------------
507	!
508	!
509	!! write(lunout,*)'THERM NOUVEAU XXXXX'
510	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
511	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
512
513	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
514	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
515
516	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
517	if (prt_level.ge.10) then
518	write(lunout1,*) 'Dans thermcell_main 1b'
519	write(lunout1,*) 'lmin ',lmin(igout)
520	write(lunout1,*) 'lalim ',lalim(igout)
521	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
522	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
523	& ,l=1,lalim(igout)+4)
524	endif
525
526
527
528
529	! Choix de la fonction d'alimentation utilisee pour la fermeture.
530	! Apparemment sans importance
531	alim_star_clos(:,:)=alim_star(:,:)
532	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
533
534	! Appel avec la version seche
535	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
536	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
537
538	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
539	! Appel avec les zmax et wmax tenant compte de la condensation
540	! Semble moins bien marcher
541	! CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
542	! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
543	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
544
545	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
546
547	if (tau_thermals>1.) then
548	lambda=exp(-ptimestep/tau_thermals)
549	f0=(1.-lambda)f+lambdaf0
550	else
551	f0=f
552	endif
553
554	! Test valable seulement en 1D mais pas genant
555	if (.not. (f0(1).ge.0.) ) then
556	abort_message = '.not. (f0(1).ge.0.)'
557	CALL abort_gcm (modname,abort_message,1)
558	endif
559
560	!-------------------------------------------------------------------------------
561	!deduction des flux
562	!-------------------------------------------------------------------------------
563
564	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
565	& lalim,lmax,alim_star, &
566	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
567	& detr,zqla,lev_out,lunout1,igout)
568	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
569
570	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
571	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
572	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
573
574	!------------------------------------------------------------------
575	! On ne prend pas directement les profils issus des calculs precedents
576	! mais on s'autorise genereusement une relaxation vers ceci avec
577	! une constante de temps tau_thermals (typiquement 1800s).
578	!------------------------------------------------------------------
579
580	if (tau_thermals>1.) then
581	lambda=exp(-ptimestep/tau_thermals)
582	fm0=(1.-lambda)fm+lambdafm0
583	entr0=(1.-lambda)entr+lambdaentr0
584	detr0=(1.-lambda)detr+lambdadetr0
585	else
586	fm0=fm
587	entr0=entr
588	detr0=detr
589	endif
590
591	!c------------------------------------------------------------------
592	! calcul du transport vertical
593	!------------------------------------------------------------------
594
595	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
596	& zthl,zdthladj,zta,lev_out)
597	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
598	& po,pdoadj,zoa,lev_out)
599
600	!------------------------------------------------------------------
601	! Calcul de la fraction de l'ascendance
602	!------------------------------------------------------------------
603	do ig=1,klon
604	fraca(ig,1)=0.
605	fraca(ig,nlay+1)=0.
606	enddo
607	do l=2,nlay
608	do ig=1,klon
609	if (zw2(ig,l).gt.1.e-10) then
610	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
611	else
612	fraca(ig,l)=0.
613	endif
614	enddo
615	enddo
616
617	!------------------------------------------------------------------
618	! calcul du transport vertical du moment horizontal
619	!------------------------------------------------------------------
620
621	!IM 090508
622	if (1.eq.1) then
623	!IM 070508 vers. _dq
624	! if (1.eq.0) then
625
626
627	! Calcul du transport de V tenant compte d'echange par gradient
628	! de pression horizontal avec l'environnement
629
630	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
631	& ,fraca,zmax &
632	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
633
634	else
635
636	! calcul purement conservatif pour le transport de V
637	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
638	& ,zu,pduadj,zua,lev_out)
639	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
640	& ,zv,pdvadj,zva,lev_out)
641	endif
642
643	! print*,'13 OK convect8'
644	do l=1,nlay
645	do ig=1,ngrid
646	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
647	enddo
648	enddo
649
650	if (prt_level.ge.1) print*,'14 OK convect8'
651	!------------------------------------------------------------------
652	! Calculs de diagnostiques pour les sorties
653	!------------------------------------------------------------------
654	!calcul de fraca pour les sorties
655
656	if (sorties) then
657	if (prt_level.ge.1) print*,'14a OK convect8'
658	! calcul du niveau de condensation
659	! initialisation
660	do ig=1,ngrid
661	nivcon(ig)=0
662	zcon(ig)=0.
663	enddo
664	!nouveau calcul
665	do ig=1,ngrid
666	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
667	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
668	enddo
669	!IM do k=1,nlay
670	do k=1,nlay-1
671	do ig=1,ngrid
672	if ((pcon(ig).le.pplay(ig,k)) &
673	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
674	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
675	endif
676	enddo
677	enddo
678	!IM
679	ierr=0
680	do ig=1,ngrid
681	if (pcon(ig).le.pplay(ig,nlay)) then
682	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
683	ierr=1
684	endif
685	enddo
686	if (ierr==1) then
687	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
688	CALL abort_gcm (modname,abort_message,1)
689	endif
690
691	if (prt_level.ge.1) print*,'14b OK convect8'
692	do k=nlay,1,-1
693	do ig=1,ngrid
694	if (zqla(ig,k).gt.1e-10) then
695	nivcon(ig)=k
696	zcon(ig)=zlev(ig,k)
697	endif
698	enddo
699	enddo
700	if (prt_level.ge.1) print*,'14c OK convect8'
701	!calcul des moments
702	!initialisation
703	do l=1,nlay
704	do ig=1,ngrid
705	q2(ig,l)=0.
706	wth2(ig,l)=0.
707	wth3(ig,l)=0.
708	ratqscth(ig,l)=0.
709	ratqsdiff(ig,l)=0.
710	enddo
711	enddo
712	if (prt_level.ge.1) print*,'14d OK convect8'
713	if (prt_level.ge.10)write(lunout,*) &
714	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
715	do l=1,nlay
716	do ig=1,ngrid
717	zf=fraca(ig,l)
718	zf2=zf/(1.-zf)
719	!
720	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
721	if(zw2(ig,l).gt.1.e-10) then
722	wth2(ig,l)=zf2(zw2(ig,l))*2
723	else
724	wth2(ig,l)=0.
725	endif
726	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
727	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
728	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
729	!test: on calcul q2/po=ratqsc
730	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
731	enddo
732	enddo
733	!calcul des flux: q, thetal et thetav
734	do l=1,nlay
735	do ig=1,ngrid
736	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
737	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
738	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
739	enddo
740	enddo
741	!
742
743	!!! nrlmd le 10/04/2012
744
745	!------------Test sur le LCL des thermiques
746	do ig=1,ngrid
747	ok_lcl(ig)=.false.
748	if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
749	enddo
750
751	!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
752	do l=1,nlay-1
753	do ig=1,ngrid
754	if (ok_lcl(ig)) then
755	if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
756	klcl(ig)=l
757	interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
758	endif
759	endif
760	enddo
761	enddo
762
763	!------------Hauteur des thermiques
764	!!jyg le 27/04/2012
765	!! do ig =1,ngrid
766	!! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
767	!! & -rhobarz(ig,klcl(ig)))*interp(ig)
768	!! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
769	!! zmax(ig)=pphi(ig,lmax(ig))/rg
770	!! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
771	!! enddo
772	do ig =1,ngrid
773	zmax(ig)=pphi(ig,lmax(ig))/rg
774	if (ok_lcl(ig)) then
775	rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
776	& -rhobarz(ig,klcl(ig)))*interp(ig)
777	zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
778	zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
779	else
780	rhobarz0(ig)=0.
781	zlcl(ig)=zmax(ig)
782	endif
783	enddo
784	!!jyg fin
785
786	!------------Calcul des propriétés du thermique au LCL
787	IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
788
789	!-----Initialisation de la TKE moyenne
790	do l=1,nlay
791	do ig=1,ngrid
792	pbl_tke_max(ig,l)=0.
793	enddo
794	enddo
795
796	!-----Calcul de la TKE moyenne
797	do nsrf=1,nbsrf
798	do l=1,nlay
799	do ig=1,ngrid
800	pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
801	enddo
802	enddo
803	enddo
804
805	!-----Initialisations des TKE dans et hors des thermiques
806	do l=1,nlay
807	do ig=1,ngrid
808	therm_tke_max(ig,l)=pbl_tke_max(ig,l)
809	env_tke_max(ig,l)=pbl_tke_max(ig,l)
810	enddo
811	enddo
812
813	!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
814	call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
815	& rg,pplev,therm_tke_max)
816	! print *,' thermcell_tke_transport -> ' !!jyg
817
818	!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
819	do l=1,nlay
820	do ig=1,ngrid
821	pbl_tke_max(ig,l)=fraca(ig,l)therm_tke_max(ig,l)+(1.-fraca(ig,l))env_tke_max(ig,l) ! Recalcul de TKE moyenne aprés transport de TKE_TH
822	env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l)) ! Recalcul de TKE dans l'environnement aprés transport de TKE_TH
823	w_ls(ig,l)=-1.omega(ig,l)/(RGrhobarz(ig,l)) ! Vitesse verticale de grande échelle
824	enddo
825	enddo
826	! print *,' apres w_ls = ' !!jyg
827
828	do ig=1,ngrid
829	if (ok_lcl(ig)) then
830	fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
831	& -fraca(ig,klcl(ig)))*interp(ig)
832	w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
833	& -zw2(ig,klcl(ig)))*interp(ig)
834	w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
835	& -w_ls(ig,klcl(ig)))*interp(ig)
836	therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
837	& +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
838	env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
839	& -env_tke_max(ig,klcl(ig)))*interp(ig)
840	pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
841	& -pbl_tke_max(ig,klcl(ig)))*interp(ig)
842	if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
843	if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
844	if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
845	else
846	fraca0(ig)=0.
847	w0(ig)=0.
848	!!jyg le 27/04/2012
849	!! zlcl(ig)=0.
850	!!
851	endif
852	enddo
853
854	ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
855	! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
856
857	!------------Triggering------------------
858	IF (iflag_trig_bl.ge.1) THEN
859
860	!-----Initialisations
861	depth(:)=0.
862	n2(:)=0.
863	s2(:)=0.
864	s_max(:)=0.
865
866	!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
867	do ig=1,ngrid
868	zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
869	depth(ig)=zmax_moy(ig)-zlcl(ig)
870	hmin(ig)=hmincoef*zlcl(ig)
871	if (depth(ig).ge.10.) then
872	s2(ig)=(hcoefdepth(ig)+hmin(ig))*2
873	n2(ig)=(1.-eps1)fraca0(ig)airephy(ig)/s2(ig)
874	!!
875	!!jyg le 27/04/2012
876	!! s_max(ig)=s2(ig)*log(n2(ig))
877	!! if (n2(ig) .lt. 1) s_max(ig)=0.
878	s_max(ig)=s2(ig)*log(max(n2(ig),1.))
879	!!fin jyg
880	else
881	s2(ig)=0.
882	n2(ig)=0.
883	s_max(ig)=0.
884	endif
885	enddo
886	! print *,'avant Calcul de Wmax ' !!jyg
887
888	!-----Calcul de Wmax et ALE_BL_STAT associée
889	!!jyg le 30/04/2012
890	!! do ig=1,ngrid
891	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
892	!! w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/su)-log(2.3.14)-log(2.log(s_max(ig)/su)-log(2.*3.14))))
893	!! ale_bl_stat(ig)=0.5w_max(ig)*2
894	!! else
895	!! w_max(ig)=0.
896	!! ale_bl_stat(ig)=0.
897	!! endif
898	!! enddo
899	susqr2pi=susqrt(2.Rpi)
900	Reuler=exp(1.)
901	do ig=1,ngrid
902	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
903	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
904	ale_bl_stat(ig)=0.5w_max(ig)*2
905	else
906	w_max(ig)=0.
907	ale_bl_stat(ig)=0.
908	endif
909	enddo
910
911	ENDIF ! iflag_trig_bl
912	! print *,'ENDIF iflag_trig_bl' !!jyg
913
914	!------------Closure------------------
915
916	IF (iflag_clos_bl.ge.1) THEN
917
918	!-----Calcul de ALP_BL_STAT
919	do ig=1,ngrid
920	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
921	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
922	& (w0(ig)**2)
923	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
924	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
925	if (iflag_clos_bl.ge.2) then
926	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
927	& (w0(ig)**2)
928	else
929	alp_bl_conv(ig)=0.
930	endif
931	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
932	enddo
933
934	!-----Sécurité ALP infinie
935	do ig=1,ngrid
936	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
937	enddo
938
939	ENDIF ! (iflag_clos_bl.ge.1)
940
941	!!! fin nrlmd le 10/04/2012
942
943	if (prt_level.ge.10) then
944	ig=igout
945	do l=1,nlay
946	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
947	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
948	enddo
949	endif
950
951	! print*,'avant calcul ale et alp'
952	!calcul de ALE et ALP pour la convection
953	Alp_bl(:)=0.
954	Ale_bl(:)=0.
955	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
956	do l=1,nlay
957	do ig=1,ngrid
958	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
959	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
960	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
961	enddo
962	enddo
963
964	!test:calcul de la ponderation des couches pour KE
965	!initialisations
966
967	fm_tot(:)=0.
968	wght_th(:,:)=1.
969	lalim_conv(:)=lalim(:)
970
971	do k=1,klev
972	do ig=1,ngrid
973	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
974	enddo
975	enddo
976
977	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
978	! plus petit que 1.e-10, on prend wght_th=1.
979	do k=1,klev
980	do ig=1,ngrid
981	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
982	wght_th(ig,k)=alim_star(ig,k)
983	endif
984	enddo
985	enddo
986
987	! print*,'apres wght_th'
988	!test pour prolonger la convection
989	do ig=1,ngrid
990	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
991	if ((alim_star(ig,1).lt.1.e-10)) then
992	lalim_conv(ig)=1
993	wght_th(ig,1)=1.
994	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
995	endif
996	enddo
997
998	!------------------------------------------------------------------------
999	! Modif CR/FH 20110310 : Alp integree sur la verticale.
1000	! Integrale verticale de ALP.
1001	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
1002	! couches
1003	!------------------------------------------------------------------------
1004
1005	alp_int(:)=0.
1006	dp_int(:)=0.
1007	do l=2,nlay
1008	do ig=1,ngrid
1009	if(l.LE.lmax(ig)) THEN
1010	zdp=pplay(ig,l-1)-pplay(ig,l)
1011	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
1012	dp_int(ig)=dp_int(ig)+zdp
1013	endif
1014	enddo
1015	enddo
1016
1017	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
1018	do ig=1,ngrid
1019	!valeur integree de alp_bl * 0.5:
1020	if (dp_int(ig)>0.) then
1021	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
1022	endif
1023	enddo!
1024	endif
1025
1026
1027	! Facteur multiplicatif sur Alp_bl
1028	Alp_bl(:)=alp_bl_k*Alp_bl(:)
1029
1030	!------------------------------------------------------------------------
1031
1032
1033	!calcul du ratqscdiff
1034	if (prt_level.ge.1) print*,'14e OK convect8'
1035	var=0.
1036	vardiff=0.
1037	ratqsdiff(:,:)=0.
1038
1039	do l=1,klev
1040	do ig=1,ngrid
1041	if (l<=lalim(ig)) then
1042	var=var+alim_star(ig,l)zqta(ig,l)1000.
1043	endif
1044	enddo
1045	enddo
1046
1047	if (prt_level.ge.1) print*,'14f OK convect8'
1048
1049	do l=1,klev
1050	do ig=1,ngrid
1051	if (l<=lalim(ig)) then
1052	zf=fraca(ig,l)
1053	zf2=zf/(1.-zf)
1054	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
1055	endif
1056	enddo
1057	enddo
1058
1059	if (prt_level.ge.1) print*,'14g OK convect8'
1060	do l=1,nlay
1061	do ig=1,ngrid
1062	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
1063	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
1064	enddo
1065	enddo
1066	!--------------------------------------------------------------------
1067	!
1068	!ecriture des fichiers sortie
1069	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
1070
1071	#ifdef wrgrads_thermcell
1072	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
1073	#include "thermcell_out3d.h"
1074	#endif
1075
1076	endif
1077
1078	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
1079
1080	return
1081	end
1082
1083	!-----------------------------------------------------------------------------
1084
1085	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
1086	IMPLICIT NONE
1087	#include "iniprint.h"
1088
1089	integer i, k, klon,klev
1090	real pplev(klon,klev+1),pplay(klon,klev)
1091	real ztv(klon,klev)
1092	real po(klon,klev)
1093	real ztva(klon,klev)
1094	real zqla(klon,klev)
1095	real f_star(klon,klev)
1096	real zw2(klon,klev)
1097	integer long(klon)
1098	real seuil
1099	character*21 comment
1100
1101	if (prt_level.ge.1) THEN
1102	print*,'WARNING !!! TEST ',comment
1103	endif
1104	return
1105
1106	! test sur la hauteur des thermiques ...
1107	do i=1,klon
1108	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
1109	if (prt_level.ge.10) then
1110	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
1111	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
1112	do k=1,klev
1113	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)
1114	enddo
1115	endif
1116	enddo
1117
1118
1119	return
1120	end
1121
1122	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
1123	! On transporte pbl_tke pour donner therm_tke
1124	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
1125	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
1126	& rg,pplev,therm_tke_max)
1127	implicit none
1128
1129	#include "iniprint.h"
1130	!=======================================================================
1131	!
1132	! Calcul du transport verticale dans la couche limite en presence
1133	! de "thermiques" explicitement representes
1134	! calcul du dq/dt une fois qu'on connait les ascendances
1135	!
1136	!=======================================================================
1137
1138	integer ngrid,nlay,nsrf
1139
1140	real ptimestep
1141	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
1142	real entr0(ngrid,nlay),rg
1143	real therm_tke_max(ngrid,nlay)
1144	real detr0(ngrid,nlay)
1145
1146
1147	real masse(ngrid,nlay),fm(ngrid,nlay+1)
1148	real entr(ngrid,nlay)
1149	real q(ngrid,nlay)
1150	integer lev_out ! niveau pour les print
1151
1152	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
1153
1154	real zzm
1155
1156	integer ig,k
1157	integer isrf
1158
1159
1160	lev_out=0
1161
1162
1163	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
1164
1165	! calcul du detrainement
1166	do k=1,nlay
1167	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
1168	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
1169	enddo
1170
1171
1172	! Decalage vertical des entrainements et detrainements.
1173	masse(:,1)=0.5*masse0(:,1)
1174	entr(:,1)=0.5*entr0(:,1)
1175	detr(:,1)=0.5*detr0(:,1)
1176	fm(:,1)=0.
1177	do k=1,nlay-1
1178	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
1179	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
1180	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
1181	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
1182	enddo
1183	fm(:,nlay+1)=0.
1184
1185	!!! nrlmd le 16/09/2010
1186	! calcul de la valeur dans les ascendances
1187	! do ig=1,ngrid
1188	! qa(ig,1)=q(ig,1)
1189	! enddo
1190	!!!
1191
1192	!do isrf=1,nsrf
1193
1194	! q(:,:)=therm_tke(:,:,isrf)
1195	q(:,:)=therm_tke_max(:,:)
1196	!!! nrlmd le 16/09/2010
1197	do ig=1,ngrid
1198	qa(ig,1)=q(ig,1)
1199	enddo
1200	!!!
1201
1202	if (1==1) then
1203	do k=2,nlay
1204	do ig=1,ngrid
1205	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
1206	& 1.e-5*masse(ig,k)) then
1207	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
1208	& /(fm(ig,k+1)+detr(ig,k))
1209	else
1210	qa(ig,k)=q(ig,k)
1211	endif
1212	if (qa(ig,k).lt.0.) then
1213	! print*,'qa<0!!!'
1214	endif
1215	if (q(ig,k).lt.0.) then
1216	! print*,'q<0!!!'
1217	endif
1218	enddo
1219	enddo
1220
1221	! Calcul du flux subsident
1222
1223	do k=2,nlay
1224	do ig=1,ngrid
1225	wqd(ig,k)=fm(ig,k)*q(ig,k)
1226	if (wqd(ig,k).lt.0.) then
1227	! print*,'wqd<0!!!'
1228	endif
1229	enddo
1230	enddo
1231	do ig=1,ngrid
1232	wqd(ig,1)=0.
1233	wqd(ig,nlay+1)=0.
1234	enddo
1235
1236	! Calcul des tendances
1237	do k=1,nlay
1238	do ig=1,ngrid
1239	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
1240	& -wqd(ig,k)+wqd(ig,k+1)) &
1241	& *ptimestep/masse(ig,k)
1242	enddo
1243	enddo
1244
1245	endif
1246
1247	therm_tke_max(:,:)=q(:,:)
1248
1249	return
1250	!!! fin nrlmd le 10/04/2012
1251	end
1252

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

Download in other formats:

Original Format