Context Navigation

thermcell_main.F90 @ 5441

Last change on this file since 5441 was 3451, checked in by fhourdin, 6 years ago

Saving old version of thermcell_plume in thermcell_plume_6A.F90
It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
to the 5B and 6A versions used for CMIP5 and CMIP6.
The latest was previously named thermcellV1_plume.
The new thermcell_plume is a clean version (removing obsolete
options) of thermcell_plume_6A.
The 3 versions are controled by
flag_thermals_ed <= 9 thermcell_plume_6A

<= 19 thermcell_plume_5B
else thermcell_plume (default 20 for convergence with 6A)

Fredho

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.7 KB

Line
1	!
2	! $Id: thermcell_main.F90 3451 2019-01-27 11:07:30Z fhourdin $
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
443	!=====================================================================
444	! Old version of thermcell_plume in thermcell_plume_6A.F90
445	! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
446	! to the 5B and 6A versions used for CMIP5 and CMIP6.
447	! The latest was previously named thermcellV1_plume.
448	! The new thermcell_plume is a clean version (removing obsolete
449	! options) of thermcell_plume_6A.
450	! The 3 versions are controled by
451	! flag_thermals_ed <= 9 thermcell_plume_6A
452	! <= 19 thermcell_plume_5B
453	! else thermcell_plume (default 20 for convergence with 6A)
454	! Fredho
455	!=====================================================================
456
457	if (iflag_thermals_ed<=9) then
458	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
459	CALL thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
460	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
461	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
462	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
463	& ,lev_out,lunout1,igout)
464
465	elseif (iflag_thermals_ed<=19) then
466	! print*,'THERM RIO et al 2010, version d Arnaud'
467	CALL thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
468	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
469	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
470	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
471	& ,lev_out,lunout1,igout)
472	else
473	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
474	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
475	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
476	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
477	& ,lev_out,lunout1,igout)
478	endif
479
480	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
481
482	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
483	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
484
485	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
486	if (prt_level.ge.10) then
487	write(lunout1,*) 'Dans thermcell_main 2'
488	write(lunout1,*) 'lmin ',lmin(igout)
489	write(lunout1,*) 'lalim ',lalim(igout)
490	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
491	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
492	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
493	endif
494
495	!-------------------------------------------------------------------------------
496	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
497	!-------------------------------------------------------------------------------
498	!
499	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
500	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
501	! Attention, w2 est transforme en sa racine carree dans cette routine
502	! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
503	wmax_tmp=0.
504	do l=1,nlay
505	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
506	enddo
507	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
508
509
510
511	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
512	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
513	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
514	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
515
516	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
517
518	!-------------------------------------------------------------------------------
519	! Fermeture,determination de f
520	!-------------------------------------------------------------------------------
521	!
522	!
523	!! write(lunout,*)'THERM NOUVEAU XXXXX'
524	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
525	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
526
527
528	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
529	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
530
531	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
532	if (prt_level.ge.10) then
533	write(lunout1,*) 'Dans thermcell_main 1b'
534	write(lunout1,*) 'lmin ',lmin(igout)
535	write(lunout1,*) 'lalim ',lalim(igout)
536	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
537	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
538	& ,l=1,lalim(igout)+4)
539	endif
540
541
542
543
544	! Choix de la fonction d'alimentation utilisee pour la fermeture.
545	! Apparemment sans importance
546	alim_star_clos(:,:)=alim_star(:,:)
547	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
548	!
549	!CR Appel de la fermeture seche
550	if (iflag_thermals_closure.eq.1) then
551
552	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
553	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
554
555	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
556	! Appel avec les zmax et wmax tenant compte de la condensation
557	! Semble moins bien marcher
558	else if (iflag_thermals_closure.eq.2) then
559
560	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
561	& zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
562
563	endif
564
565	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
566
567	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
568
569	if (tau_thermals>1.) then
570	lambda=exp(-ptimestep/tau_thermals)
571	f0=(1.-lambda)f+lambdaf0
572	else
573	f0=f
574	endif
575
576	! Test valable seulement en 1D mais pas genant
577	if (.not. (f0(1).ge.0.) ) then
578	abort_message = '.not. (f0(1).ge.0.)'
579	CALL abort_physic (modname,abort_message,1)
580	endif
581
582	!-------------------------------------------------------------------------------
583	!deduction des flux
584	!-------------------------------------------------------------------------------
585
586	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
587	& lalim,lmax,alim_star, &
588	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
589	& detr,zqla,lev_out,lunout1,igout)
590	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
591
592	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
593	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
594	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
595
596	!------------------------------------------------------------------
597	! On ne prend pas directement les profils issus des calculs precedents
598	! mais on s'autorise genereusement une relaxation vers ceci avec
599	! une constante de temps tau_thermals (typiquement 1800s).
600	!------------------------------------------------------------------
601
602	if (tau_thermals>1.) then
603	lambda=exp(-ptimestep/tau_thermals)
604	fm0=(1.-lambda)fm+lambdafm0
605	entr0=(1.-lambda)entr+lambdaentr0
606	detr0=(1.-lambda)detr+lambdadetr0
607	else
608	fm0=fm
609	entr0=entr
610	detr0=detr
611	endif
612
613	!c------------------------------------------------------------------
614	! calcul du transport vertical
615	!------------------------------------------------------------------
616
617	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
618	& zthl,zdthladj,zta,lev_out)
619	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
620	& po,pdoadj,zoa,lev_out)
621
622	!------------------------------------------------------------------
623	! Calcul de la fraction de l'ascendance
624	!------------------------------------------------------------------
625	do ig=1,klon
626	fraca(ig,1)=0.
627	fraca(ig,nlay+1)=0.
628	enddo
629	do l=2,nlay
630	do ig=1,klon
631	if (zw2(ig,l).gt.1.e-10) then
632	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
633	else
634	fraca(ig,l)=0.
635	endif
636	enddo
637	enddo
638
639	!------------------------------------------------------------------
640	! calcul du transport vertical du moment horizontal
641	!------------------------------------------------------------------
642
643	!IM 090508
644	if (dvdq == 0 ) then
645
646	! Calcul du transport de V tenant compte d'echange par gradient
647	! de pression horizontal avec l'environnement
648
649	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
650	! & ,fraca*dvdq,zmax &
651	& ,fraca,zmax &
652	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
653
654	else
655
656	! calcul purement conservatif pour le transport de V
657	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
658	& ,zu,pduadj,zua,lev_out)
659	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
660	& ,zv,pdvadj,zva,lev_out)
661
662	endif
663
664	! print*,'13 OK convect8'
665	do l=1,nlay
666	do ig=1,ngrid
667	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
668	enddo
669	enddo
670
671	if (prt_level.ge.1) print*,'14 OK convect8'
672	!------------------------------------------------------------------
673	! Calculs de diagnostiques pour les sorties
674	!------------------------------------------------------------------
675	!calcul de fraca pour les sorties
676
677	if (sorties) then
678	if (prt_level.ge.1) print*,'14a OK convect8'
679	! calcul du niveau de condensation
680	! initialisation
681	do ig=1,ngrid
682	nivcon(ig)=0
683	zcon(ig)=0.
684	enddo
685	!nouveau calcul
686	do ig=1,ngrid
687	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
688	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
689	enddo
690	!IM do k=1,nlay
691	do k=1,nlay-1
692	do ig=1,ngrid
693	if ((pcon(ig).le.pplay(ig,k)) &
694	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
695	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
696	endif
697	enddo
698	enddo
699	!IM
700	ierr=0
701	do ig=1,ngrid
702	if (pcon(ig).le.pplay(ig,nlay)) then
703	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
704	ierr=1
705	endif
706	enddo
707	if (ierr==1) then
708	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
709	CALL abort_physic (modname,abort_message,1)
710	endif
711
712	if (prt_level.ge.1) print*,'14b OK convect8'
713	do k=nlay,1,-1
714	do ig=1,ngrid
715	if (zqla(ig,k).gt.1e-10) then
716	nivcon(ig)=k
717	zcon(ig)=zlev(ig,k)
718	endif
719	enddo
720	enddo
721	if (prt_level.ge.1) print*,'14c OK convect8'
722	!calcul des moments
723	!initialisation
724	do l=1,nlay
725	do ig=1,ngrid
726	q2(ig,l)=0.
727	wth2(ig,l)=0.
728	wth3(ig,l)=0.
729	ratqscth(ig,l)=0.
730	ratqsdiff(ig,l)=0.
731	enddo
732	enddo
733	if (prt_level.ge.1) print*,'14d OK convect8'
734	if (prt_level.ge.10)write(lunout,*) &
735	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
736	do l=1,nlay
737	do ig=1,ngrid
738	zf=fraca(ig,l)
739	zf2=zf/(1.-zf)
740	!
741	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
742	if(zw2(ig,l).gt.1.e-10) then
743	wth2(ig,l)=zf2(zw2(ig,l))*2
744	else
745	wth2(ig,l)=0.
746	endif
747	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
748	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
749	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
750	!test: on calcul q2/po=ratqsc
751	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
752	enddo
753	enddo
754	!calcul des flux: q, thetal et thetav
755	do l=1,nlay
756	do ig=1,ngrid
757	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
758	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
759	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
760	enddo
761	enddo
762	!
763	! $Id: thermcell_main.F90 3451 2019-01-27 11:07:30Z fhourdin $
764	!
765	CALL thermcell_alp(ngrid,nlay,ptimestep &
766	& ,pplay,pplev &
767	& ,fm0,entr0,lmax &
768	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
769	& ,zw2,fraca &
770	!!! necessire en plus
771	& ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
772	!!! nrlmd le 10/04/2012
773	& ,pbl_tke,pctsrf,omega,airephy &
774	& ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
775	& ,n2,s2,ale_bl_stat &
776	& ,therm_tke_max,env_tke_max &
777	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
778	& ,alp_bl_conv,alp_bl_stat &
779	!!! fin nrlmd le 10/04/2012
780	& )
781
782
783
784	!calcul du ratqscdiff
785	if (prt_level.ge.1) print*,'14e OK convect8'
786	var=0.
787	vardiff=0.
788	ratqsdiff(:,:)=0.
789
790	do l=1,klev
791	do ig=1,ngrid
792	if (l<=lalim(ig)) then
793	var=var+alim_star(ig,l)zqta(ig,l)1000.
794	endif
795	enddo
796	enddo
797
798	if (prt_level.ge.1) print*,'14f OK convect8'
799
800	do l=1,klev
801	do ig=1,ngrid
802	if (l<=lalim(ig)) then
803	zf=fraca(ig,l)
804	zf2=zf/(1.-zf)
805	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
806	endif
807	enddo
808	enddo
809
810	if (prt_level.ge.1) print*,'14g OK convect8'
811	do l=1,nlay
812	do ig=1,ngrid
813	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
814	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
815	enddo
816	enddo
817	!--------------------------------------------------------------------
818	!
819	!ecriture des fichiers sortie
820	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
821
822	endif
823
824	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
825
826	return
827	end
828
829	!-----------------------------------------------------------------------------
830
831	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
832	USE print_control_mod, ONLY: prt_level
833	IMPLICIT NONE
834
835	integer i, k, klon,klev
836	real pplev(klon,klev+1),pplay(klon,klev)
837	real ztv(klon,klev)
838	real po(klon,klev)
839	real ztva(klon,klev)
840	real zqla(klon,klev)
841	real f_star(klon,klev)
842	real zw2(klon,klev)
843	integer long(klon)
844	real seuil
845	character*21 comment
846
847	if (prt_level.ge.1) THEN
848	print*,'WARNING !!! TEST ',comment
849	endif
850	return
851
852	! test sur la hauteur des thermiques ...
853	do i=1,klon
854	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
855	if (prt_level.ge.10) then
856	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
857	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
858	do k=1,klev
859	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)
860	enddo
861	endif
862	enddo
863
864
865	return
866	end
867
868	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
869	! On transporte pbl_tke pour donner therm_tke
870	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
871	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
872	& rg,pplev,therm_tke_max)
873	USE print_control_mod, ONLY: prt_level
874	implicit none
875
876	!=======================================================================
877	!
878	! Calcul du transport verticale dans la couche limite en presence
879	! de "thermiques" explicitement representes
880	! calcul du dq/dt une fois qu'on connait les ascendances
881	!
882	!=======================================================================
883
884	integer ngrid,nlay,nsrf
885
886	real ptimestep
887	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
888	real entr0(ngrid,nlay),rg
889	real therm_tke_max(ngrid,nlay)
890	real detr0(ngrid,nlay)
891
892
893	real masse(ngrid,nlay),fm(ngrid,nlay+1)
894	real entr(ngrid,nlay)
895	real q(ngrid,nlay)
896	integer lev_out ! niveau pour les print
897
898	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
899
900	real zzm
901
902	integer ig,k
903	integer isrf
904
905
906	lev_out=0
907
908
909	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
910
911	! calcul du detrainement
912	do k=1,nlay
913	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
914	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
915	enddo
916
917
918	! Decalage vertical des entrainements et detrainements.
919	masse(:,1)=0.5*masse0(:,1)
920	entr(:,1)=0.5*entr0(:,1)
921	detr(:,1)=0.5*detr0(:,1)
922	fm(:,1)=0.
923	do k=1,nlay-1
924	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
925	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
926	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
927	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
928	enddo
929	fm(:,nlay+1)=0.
930
931	!!! nrlmd le 16/09/2010
932	! calcul de la valeur dans les ascendances
933	! do ig=1,ngrid
934	! qa(ig,1)=q(ig,1)
935	! enddo
936	!!!
937
938	!do isrf=1,nsrf
939
940	! q(:,:)=therm_tke(:,:,isrf)
941	q(:,:)=therm_tke_max(:,:)
942	!!! nrlmd le 16/09/2010
943	do ig=1,ngrid
944	qa(ig,1)=q(ig,1)
945	enddo
946	!!!
947
948	if (1==1) then
949	do k=2,nlay
950	do ig=1,ngrid
951	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
952	& 1.e-5*masse(ig,k)) then
953	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
954	& /(fm(ig,k+1)+detr(ig,k))
955	else
956	qa(ig,k)=q(ig,k)
957	endif
958	if (qa(ig,k).lt.0.) then
959	! print*,'qa<0!!!'
960	endif
961	if (q(ig,k).lt.0.) then
962	! print*,'q<0!!!'
963	endif
964	enddo
965	enddo
966
967	! Calcul du flux subsident
968
969	do k=2,nlay
970	do ig=1,ngrid
971	wqd(ig,k)=fm(ig,k)*q(ig,k)
972	if (wqd(ig,k).lt.0.) then
973	! print*,'wqd<0!!!'
974	endif
975	enddo
976	enddo
977	do ig=1,ngrid
978	wqd(ig,1)=0.
979	wqd(ig,nlay+1)=0.
980	enddo
981
982	! Calcul des tendances
983	do k=1,nlay
984	do ig=1,ngrid
985	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
986	& -wqd(ig,k)+wqd(ig,k+1)) &
987	& *ptimestep/masse(ig,k)
988	enddo
989	enddo
990
991	endif
992
993	therm_tke_max(:,:)=q(:,:)
994
995	return
996	!!! fin nrlmd le 10/04/2012
997	end
998

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

Download in other formats:

Original Format