Context Navigation

thermcell_main.F90 @ 5443

Last change on this file since 5443 was 2387, checked in by fhourdin, 9 years ago
Poursuite du travail sur les thermiques. Extraction d'une routine calculant l'alimentation laterale à la base. + une correction de Jean-Yves sur des initialisations maquantes dans thermcell_main (mais qui sont maintenant dans thermcell_alp).
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 32.0 KB

Line
1	!
2	! $Id: thermcell_main.F90 2387 2015-11-07 09:27:40Z evignon $
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	& ,ztva )
22
23	USE dimphy
24	USE ioipsl
25	USE indice_sol_mod
26	USE print_control_mod, ONLY: lunout,prt_level
27	IMPLICIT NONE
28
29	!=======================================================================
30	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
31	! Version du 09.02.07
32	! Calcul du transport vertical dans la couche limite en presence
33	! de "thermiques" explicitement representes avec processus nuageux
34	!
35	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
36	!
37	! le thermique est suppose homogene et dissipe par melange avec
38	! son environnement. la longueur l_mix controle l'efficacite du
39	! melange
40	!
41	! Le calcul du transport des differentes especes se fait en prenant
42	! en compte:
43	! 1. un flux de masse montant
44	! 2. un flux de masse descendant
45	! 3. un entrainement
46	! 4. un detrainement
47	!
48	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
49	! Introduction of an implicit computation of vertical advection in
50	! the environment of thermal plumes in thermcell_dq
51	! impl = 0 : explicit, 1 : implicit, -1 : old version
52	! controled by iflag_thermals =
53	! 15, 16 run with impl=-1 : numerical convergence with NPv3
54	! 17, 18 run with impl=1 : more stable
55	! 15 and 17 correspond to the activation of the stratocumulus "bidouille"
56	!
57	!=======================================================================
58
59
60	!-----------------------------------------------------------------------
61	! declarations:
62	! -------------
63
64	#include "YOMCST.h"
65	#include "YOETHF.h"
66	#include "FCTTRE.h"
67	#include "thermcell.h"
68
69	! arguments:
70	! ----------
71
72	!IM 140508
73	INTEGER itap
74
75	INTEGER ngrid,nlay
76	real ptimestep
77	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
78	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
79	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
80	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
81	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
82	real pphi(ngrid,nlay)
83	LOGICAL debut
84
85	! local:
86	! ------
87
88	integer icount
89
90	integer, save :: dvdq=1,dqimpl=-1
91	!$OMP THREADPRIVATE(dvdq,dqimpl)
92	data icount/0/
93	save icount
94	!$OMP THREADPRIVATE(icount)
95
96	integer,save :: igout=1
97	!$OMP THREADPRIVATE(igout)
98	integer,save :: lunout1=6
99	!$OMP THREADPRIVATE(lunout1)
100	integer,save :: lev_out=10
101	!$OMP THREADPRIVATE(lev_out)
102
103	REAL susqr2pi, Reuler
104
105	INTEGER ig,k,l,ll,ierr
106	real zsortie1d(klon)
107	INTEGER lmax(klon),lmin(klon),lalim(klon)
108	INTEGER lmix(klon)
109	INTEGER lmix_bis(klon)
110	real linter(klon)
111	real zmix(klon)
112	real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
113	! real fraca(klon,klev)
114
115	real zmax_sec(klon)
116	!on garde le zmax du pas de temps precedent
117	real zmax0(klon)
118	!FH/IM save zmax0
119
120	real lambda
121
122	real zlev(klon,klev+1),zlay(klon,klev)
123	real deltaz(klon,klev)
124	REAL zh(klon,klev)
125	real zthl(klon,klev),zdthladj(klon,klev)
126	REAL ztv(klon,klev)
127	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
128	real zl(klon,klev)
129	real zsortie(klon,klev)
130	real zva(klon,klev)
131	real zua(klon,klev)
132	real zoa(klon,klev)
133
134	real zta(klon,klev)
135	real zha(klon,klev)
136	real fraca(klon,klev+1)
137	real zf,zf2
138	real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
139	real q2(klon,klev)
140	! FH probleme de dimensionnement avec l'allocation dynamique
141	! common/comtherm/thetath2,wth2
142	real wq(klon,klev)
143	real wthl(klon,klev)
144	real wthv(klon,klev)
145
146	real ratqscth(klon,klev)
147	real var
148	real vardiff
149	real ratqsdiff(klon,klev)
150
151	logical sorties
152	real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
153	real zpspsk(klon,klev)
154
155	real wmax(klon)
156	real wmax_tmp(klon)
157	real wmax_sec(klon)
158	real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
159	real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
160
161	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
162	!niveau de condensation
163	integer nivcon(klon)
164	real zcon(klon)
165	REAL CHI
166	real zcon2(klon)
167	real pcon(klon)
168	real zqsat(klon,klev)
169	real zqsatth(klon,klev)
170
171	real f_star(klon,klev+1),entr_star(klon,klev)
172	real detr_star(klon,klev)
173	real alim_star_tot(klon)
174	real alim_star(klon,klev)
175	real alim_star_clos(klon,klev)
176	real f(klon), f0(klon)
177	!FH/IM save f0
178	real zlevinter(klon)
179	real seuil
180	real csc(klon,klev)
181
182	!!! nrlmd le 10/04/2012
183
184	!------Entrées
185	real pbl_tke(klon,klev+1,nbsrf)
186	real pctsrf(klon,nbsrf)
187	real omega(klon,klev)
188	real airephy(klon)
189	!------Sorties
190	real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
191	real therm_tke_max0(klon),env_tke_max0(klon)
192	real n2(klon),s2(klon)
193	real ale_bl_stat(klon)
194	real therm_tke_max(klon,klev),env_tke_max(klon,klev)
195	real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
196	!------Local
197	integer nsrf
198	real rhobarz0(klon) ! Densité au LCL
199	logical ok_lcl(klon) ! Existence du LCL des thermiques
200	integer klcl(klon) ! Niveau du LCL
201	real interp(klon) ! Coef d'interpolation pour le LCL
202	!--Triggering
203	real Su ! Surface unité: celle d'un updraft élémentaire
204	parameter(Su=4e4)
205	real hcoef ! Coefficient directeur pour le calcul de s2
206	parameter(hcoef=1)
207	real hmincoef ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
208	parameter(hmincoef=0.3)
209	real eps1 ! Fraction de surface occupée par la population 1 : eps1=n1s1/(fraca0Sd)
210	parameter(eps1=0.3)
211	real hmin(ngrid) ! Ordonnée à l'origine pour le calcul de s2
212	real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
213	real zmax_moy_coef
214	parameter(zmax_moy_coef=0.33)
215	real depth(klon) ! Epaisseur moyenne du cumulus
216	real w_max(klon) ! Vitesse max statistique
217	real s_max(klon)
218	!--Closure
219	real pbl_tke_max(klon,klev) ! Profil de TKE moyenne
220	real pbl_tke_max0(klon) ! TKE moyenne au LCL
221	real w_ls(klon,klev) ! Vitesse verticale grande échelle (m/s)
222	real coef_m ! On considère un rendement pour alp_bl_fluct_m
223	parameter(coef_m=1.)
224	real coef_tke ! On considère un rendement pour alp_bl_fluct_tke
225	parameter(coef_tke=1.)
226
227	!!! fin nrlmd le 10/04/2012
228
229	!
230	!nouvelles variables pour la convection
231	real Ale_bl(klon)
232	real Alp_bl(klon)
233	real alp_int(klon),dp_int(klon),zdp
234	real ale_int(klon)
235	integer n_int(klon)
236	real fm_tot(klon)
237	real wght_th(klon,klev)
238	integer lalim_conv(klon)
239	!v1d logical therm
240	!v1d save therm
241
242	character*2 str2
243	character*10 str10
244
245	character (len=20) :: modname='thermcell_main'
246	character (len=80) :: abort_message
247
248	EXTERNAL SCOPY
249	!
250
251	!-----------------------------------------------------------------------
252	! initialisation:
253	! ---------------
254	!
255
256	seuil=0.25
257
258	if (debut) then
259	if (iflag_thermals==15.or.iflag_thermals==16) then
260	dvdq=0
261	dqimpl=-1
262	else
263	dvdq=1
264	dqimpl=1
265	endif
266
267	fm0=0.
268	entr0=0.
269	detr0=0.
270	endif
271	fm=0. ; entr=0. ; detr=0.
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
514	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
515	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
516
517	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
518	if (prt_level.ge.10) then
519	write(lunout1,*) 'Dans thermcell_main 1b'
520	write(lunout1,*) 'lmin ',lmin(igout)
521	write(lunout1,*) 'lalim ',lalim(igout)
522	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
523	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
524	& ,l=1,lalim(igout)+4)
525	endif
526
527
528
529
530	! Choix de la fonction d'alimentation utilisee pour la fermeture.
531	! Apparemment sans importance
532	alim_star_clos(:,:)=alim_star(:,:)
533	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
534	!
535	!CR Appel de la fermeture seche
536	if (iflag_thermals_closure.eq.1) then
537
538	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
539	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
540
541	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
542	! Appel avec les zmax et wmax tenant compte de la condensation
543	! Semble moins bien marcher
544	else if (iflag_thermals_closure.eq.2) then
545
546	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
547	& zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
548
549	endif
550
551	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
552
553	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
554
555	if (tau_thermals>1.) then
556	lambda=exp(-ptimestep/tau_thermals)
557	f0=(1.-lambda)f+lambdaf0
558	else
559	f0=f
560	endif
561
562	! Test valable seulement en 1D mais pas genant
563	if (.not. (f0(1).ge.0.) ) then
564	abort_message = '.not. (f0(1).ge.0.)'
565	CALL abort_physic (modname,abort_message,1)
566	endif
567
568	!-------------------------------------------------------------------------------
569	!deduction des flux
570	!-------------------------------------------------------------------------------
571
572	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
573	& lalim,lmax,alim_star, &
574	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
575	& detr,zqla,lev_out,lunout1,igout)
576	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
577
578	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
579	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
580	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
581
582	!------------------------------------------------------------------
583	! On ne prend pas directement les profils issus des calculs precedents
584	! mais on s'autorise genereusement une relaxation vers ceci avec
585	! une constante de temps tau_thermals (typiquement 1800s).
586	!------------------------------------------------------------------
587
588	if (tau_thermals>1.) then
589	lambda=exp(-ptimestep/tau_thermals)
590	fm0=(1.-lambda)fm+lambdafm0
591	entr0=(1.-lambda)entr+lambdaentr0
592	detr0=(1.-lambda)detr+lambdadetr0
593	else
594	fm0=fm
595	entr0=entr
596	detr0=detr
597	endif
598
599	!c------------------------------------------------------------------
600	! calcul du transport vertical
601	!------------------------------------------------------------------
602
603	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
604	& zthl,zdthladj,zta,lev_out)
605	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
606	& po,pdoadj,zoa,lev_out)
607
608	!------------------------------------------------------------------
609	! Calcul de la fraction de l'ascendance
610	!------------------------------------------------------------------
611	do ig=1,klon
612	fraca(ig,1)=0.
613	fraca(ig,nlay+1)=0.
614	enddo
615	do l=2,nlay
616	do ig=1,klon
617	if (zw2(ig,l).gt.1.e-10) then
618	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
619	else
620	fraca(ig,l)=0.
621	endif
622	enddo
623	enddo
624
625	!------------------------------------------------------------------
626	! calcul du transport vertical du moment horizontal
627	!------------------------------------------------------------------
628
629	!IM 090508
630	if (dvdq == 0 ) then
631
632	! Calcul du transport de V tenant compte d'echange par gradient
633	! de pression horizontal avec l'environnement
634
635	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
636	! & ,fraca*dvdq,zmax &
637	& ,fraca,zmax &
638	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
639
640	else
641
642	! calcul purement conservatif pour le transport de V
643	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
644	& ,zu,pduadj,zua,lev_out)
645	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
646	& ,zv,pdvadj,zva,lev_out)
647
648	endif
649
650	! print*,'13 OK convect8'
651	do l=1,nlay
652	do ig=1,ngrid
653	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
654	enddo
655	enddo
656
657	if (prt_level.ge.1) print*,'14 OK convect8'
658	!------------------------------------------------------------------
659	! Calculs de diagnostiques pour les sorties
660	!------------------------------------------------------------------
661	!calcul de fraca pour les sorties
662
663	if (sorties) then
664	if (prt_level.ge.1) print*,'14a OK convect8'
665	! calcul du niveau de condensation
666	! initialisation
667	do ig=1,ngrid
668	nivcon(ig)=0
669	zcon(ig)=0.
670	enddo
671	!nouveau calcul
672	do ig=1,ngrid
673	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
674	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
675	enddo
676	!IM do k=1,nlay
677	do k=1,nlay-1
678	do ig=1,ngrid
679	if ((pcon(ig).le.pplay(ig,k)) &
680	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
681	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
682	endif
683	enddo
684	enddo
685	!IM
686	ierr=0
687	do ig=1,ngrid
688	if (pcon(ig).le.pplay(ig,nlay)) then
689	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
690	ierr=1
691	endif
692	enddo
693	if (ierr==1) then
694	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
695	CALL abort_physic (modname,abort_message,1)
696	endif
697
698	if (prt_level.ge.1) print*,'14b OK convect8'
699	do k=nlay,1,-1
700	do ig=1,ngrid
701	if (zqla(ig,k).gt.1e-10) then
702	nivcon(ig)=k
703	zcon(ig)=zlev(ig,k)
704	endif
705	enddo
706	enddo
707	if (prt_level.ge.1) print*,'14c OK convect8'
708	!calcul des moments
709	!initialisation
710	do l=1,nlay
711	do ig=1,ngrid
712	q2(ig,l)=0.
713	wth2(ig,l)=0.
714	wth3(ig,l)=0.
715	ratqscth(ig,l)=0.
716	ratqsdiff(ig,l)=0.
717	enddo
718	enddo
719	if (prt_level.ge.1) print*,'14d OK convect8'
720	if (prt_level.ge.10)write(lunout,*) &
721	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
722	do l=1,nlay
723	do ig=1,ngrid
724	zf=fraca(ig,l)
725	zf2=zf/(1.-zf)
726	!
727	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
728	if(zw2(ig,l).gt.1.e-10) then
729	wth2(ig,l)=zf2(zw2(ig,l))*2
730	else
731	wth2(ig,l)=0.
732	endif
733	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
734	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
735	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
736	!test: on calcul q2/po=ratqsc
737	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
738	enddo
739	enddo
740	!calcul des flux: q, thetal et thetav
741	do l=1,nlay
742	do ig=1,ngrid
743	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
744	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
745	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
746	enddo
747	enddo
748	!
749	! $Id: thermcell_main.F90 2387 2015-11-07 09:27:40Z evignon $
750	!
751	CALL thermcell_alp(ngrid,nlay,ptimestep &
752	& ,pplay,pplev &
753	& ,fm0,entr0,lmax &
754	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
755	& ,zw2,fraca &
756	!!! necessire en plus
757	& ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
758	!!! nrlmd le 10/04/2012
759	& ,pbl_tke,pctsrf,omega,airephy &
760	& ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
761	& ,n2,s2,ale_bl_stat &
762	& ,therm_tke_max,env_tke_max &
763	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
764	& ,alp_bl_conv,alp_bl_stat &
765	!!! fin nrlmd le 10/04/2012
766	& )
767
768
769
770	!calcul du ratqscdiff
771	if (prt_level.ge.1) print*,'14e OK convect8'
772	var=0.
773	vardiff=0.
774	ratqsdiff(:,:)=0.
775
776	do l=1,klev
777	do ig=1,ngrid
778	if (l<=lalim(ig)) then
779	var=var+alim_star(ig,l)zqta(ig,l)1000.
780	endif
781	enddo
782	enddo
783
784	if (prt_level.ge.1) print*,'14f OK convect8'
785
786	do l=1,klev
787	do ig=1,ngrid
788	if (l<=lalim(ig)) then
789	zf=fraca(ig,l)
790	zf2=zf/(1.-zf)
791	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
792	endif
793	enddo
794	enddo
795
796	if (prt_level.ge.1) print*,'14g OK convect8'
797	do l=1,nlay
798	do ig=1,ngrid
799	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
800	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
801	enddo
802	enddo
803	!--------------------------------------------------------------------
804	!
805	!ecriture des fichiers sortie
806	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
807
808	endif
809
810	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
811
812	return
813	end
814
815	!-----------------------------------------------------------------------------
816
817	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
818	USE print_control_mod, ONLY: prt_level
819	IMPLICIT NONE
820
821	integer i, k, klon,klev
822	real pplev(klon,klev+1),pplay(klon,klev)
823	real ztv(klon,klev)
824	real po(klon,klev)
825	real ztva(klon,klev)
826	real zqla(klon,klev)
827	real f_star(klon,klev)
828	real zw2(klon,klev)
829	integer long(klon)
830	real seuil
831	character*21 comment
832
833	if (prt_level.ge.1) THEN
834	print*,'WARNING !!! TEST ',comment
835	endif
836	return
837
838	! test sur la hauteur des thermiques ...
839	do i=1,klon
840	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
841	if (prt_level.ge.10) then
842	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
843	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
844	do k=1,klev
845	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)
846	enddo
847	endif
848	enddo
849
850
851	return
852	end
853
854	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
855	! On transporte pbl_tke pour donner therm_tke
856	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
857	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
858	& rg,pplev,therm_tke_max)
859	USE print_control_mod, ONLY: prt_level
860	implicit none
861
862	!=======================================================================
863	!
864	! Calcul du transport verticale dans la couche limite en presence
865	! de "thermiques" explicitement representes
866	! calcul du dq/dt une fois qu'on connait les ascendances
867	!
868	!=======================================================================
869
870	integer ngrid,nlay,nsrf
871
872	real ptimestep
873	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
874	real entr0(ngrid,nlay),rg
875	real therm_tke_max(ngrid,nlay)
876	real detr0(ngrid,nlay)
877
878
879	real masse(ngrid,nlay),fm(ngrid,nlay+1)
880	real entr(ngrid,nlay)
881	real q(ngrid,nlay)
882	integer lev_out ! niveau pour les print
883
884	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
885
886	real zzm
887
888	integer ig,k
889	integer isrf
890
891
892	lev_out=0
893
894
895	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
896
897	! calcul du detrainement
898	do k=1,nlay
899	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
900	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
901	enddo
902
903
904	! Decalage vertical des entrainements et detrainements.
905	masse(:,1)=0.5*masse0(:,1)
906	entr(:,1)=0.5*entr0(:,1)
907	detr(:,1)=0.5*detr0(:,1)
908	fm(:,1)=0.
909	do k=1,nlay-1
910	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
911	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
912	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
913	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
914	enddo
915	fm(:,nlay+1)=0.
916
917	!!! nrlmd le 16/09/2010
918	! calcul de la valeur dans les ascendances
919	! do ig=1,ngrid
920	! qa(ig,1)=q(ig,1)
921	! enddo
922	!!!
923
924	!do isrf=1,nsrf
925
926	! q(:,:)=therm_tke(:,:,isrf)
927	q(:,:)=therm_tke_max(:,:)
928	!!! nrlmd le 16/09/2010
929	do ig=1,ngrid
930	qa(ig,1)=q(ig,1)
931	enddo
932	!!!
933
934	if (1==1) then
935	do k=2,nlay
936	do ig=1,ngrid
937	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
938	& 1.e-5*masse(ig,k)) then
939	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
940	& /(fm(ig,k+1)+detr(ig,k))
941	else
942	qa(ig,k)=q(ig,k)
943	endif
944	if (qa(ig,k).lt.0.) then
945	! print*,'qa<0!!!'
946	endif
947	if (q(ig,k).lt.0.) then
948	! print*,'q<0!!!'
949	endif
950	enddo
951	enddo
952
953	! Calcul du flux subsident
954
955	do k=2,nlay
956	do ig=1,ngrid
957	wqd(ig,k)=fm(ig,k)*q(ig,k)
958	if (wqd(ig,k).lt.0.) then
959	! print*,'wqd<0!!!'
960	endif
961	enddo
962	enddo
963	do ig=1,ngrid
964	wqd(ig,1)=0.
965	wqd(ig,nlay+1)=0.
966	enddo
967
968	! Calcul des tendances
969	do k=1,nlay
970	do ig=1,ngrid
971	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
972	& -wqd(ig,k)+wqd(ig,k+1)) &
973	& *ptimestep/masse(ig,k)
974	enddo
975	enddo
976
977	endif
978
979	therm_tke_max(:,:)=q(:,:)
980
981	return
982	!!! fin nrlmd le 10/04/2012
983	end
984

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

Download in other formats:

Original Format