Context Navigation

thermcell_main.F90 @ 2104

Last change on this file since 2104 was 2104, checked in by aboissinot, 6 years ago
Fix an oversight in thermcell_main
File size: 37.3 KB

Line
1	!
2	!
3	!
4	SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep, &
5	pplay,pplev,pphi,firstcall, &
6	pu,pv,pt,po, &
7	pduadj,pdvadj,pdtadj,pdoadj, &
8	f0,fm0,entr0,detr0, &
9	zqta,zqla,ztv,ztva,ztla,zthl,zqsatth, &
10	zw2,fraca, &
11	lmin,lmix,lalim,lmax, &
12	zpspsk,ratqscth,ratqsdiff, &
13	Ale_bl,Alp_bl,lalim_conv,wght_th, &
14	!!! nrlmd le 10/04/2012
15	pbl_tke,pctsrf,omega,airephy, &
16	zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0, &
17	n2,s2,ale_bl_stat, &
18	therm_tke_max,env_tke_max, &
19	alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke, &
20	alp_bl_conv,alp_bl_stat)
21	!!! fin nrlmd le 10/04/2012
22
23
24	!==============================================================================
25	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
26	! Version du 09.02.07
27	! Calcul du transport vertical dans la couche limite en presence
28	! de "thermiques" explicitement representes avec processus nuageux
29	!
30	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
31	!
32	! le thermique est suppose homogene et dissipe par melange avec
33	! son environnement. la longueur l_mix controle l'efficacite du
34	! melange
35	!
36	! Le calcul du transport des differentes especes se fait en prenant
37	! en compte:
38	! 1. un flux de masse montant
39	! 2. un flux de masse descendant
40	! 3. un entrainement
41	! 4. un detrainement
42	!
43	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
44	! Introduction of an implicit computation of vertical advection in
45	! the environment of thermal plumes in thermcell_dq
46	! impl = 0 : explicit, 1 : implicit, -1 : old version
47	! controled by iflag_thermals =
48	! 15, 16 run with impl=-1 : numerical convergence with NPv3
49	! 17, 18 run with impl=1 : more stable
50	! 15 and 17 correspond to the activation of the stratocumulus "bidouille"
51	!
52	!==============================================================================
53
54	USE thermcell_mod
55	USE print_control_mod, ONLY: lunout, prt_level
56
57	IMPLICIT NONE
58
59
60	!==============================================================================
61	! Declaration
62	!==============================================================================
63
64	! inputs:
65	! -------
66
67	INTEGER itap
68	INTEGER ngrid
69	INTEGER nlay
70
71	REAL ptimestep
72	REAL pplay(ngrid,nlay)
73	REAL pplev(ngrid,nlay+1)
74	REAL pphi(ngrid,nlay)
75
76	REAL pt(ngrid,nlay) !
77	REAL pu(ngrid,nlay) !
78	REAL pv(ngrid,nlay) !
79	REAL po(ngrid,nlay) !
80
81	LOGICAL firstcall
82
83	! outputs:
84	! --------
85
86	REAL pdtadj(ngrid,nlay) ! t convective variations
87	REAL pduadj(ngrid,nlay) ! u convective variations
88	REAL pdvadj(ngrid,nlay) ! v convective variations
89	REAL pdoadj(ngrid,nlay) ! water convective variations
90
91	REAL fm0(ngrid,nlay+1) ! mass flux (after possible time relaxation)
92	REAL entr0(ngrid,nlay) ! entrainment (after possible time relaxation)
93	REAL detr0(ngrid,nlay) ! detrainment (after possible time relaxation)
94	REAL f0(ngrid) ! mass flux norm (after possible time relaxation)
95
96	! local:
97	! ------
98
99	INTEGER, save :: dvimpl=1
100	!$OMP THREADPRIVATE(dvimpl)
101
102	INTEGER, save :: dqimpl=-1
103	!$OMP THREADPRIVATE(dqimpl)
104
105	INTEGER, SAVE :: igout=1
106	!$OMP THREADPRIVATE(igout)
107
108	INTEGER, SAVE :: lunout1=6
109	!$OMP THREADPRIVATE(lunout1)
110
111	INTEGER, SAVE :: lev_out=10
112	!$OMP THREADPRIVATE(lev_out)
113
114	INTEGER ig,k,l,ll,ierr
115	INTEGER lmix_bis(ngrid)
116	INTEGER lmax(ngrid) !
117	INTEGER lmin(ngrid) !
118	INTEGER lalim(ngrid) !
119	INTEGER lmix(ngrid) !
120
121	REAL linter(ngrid)
122	REAL zmix(ngrid)
123	REAL zmax(ngrid)
124	REAL zw2(ngrid,nlay+1)
125	REAL zw_est(ngrid,nlay+1)
126	REAL zmax_sec(ngrid)
127
128	REAL zlay(ngrid,nlay) ! layers altitude
129	REAL zlev(ngrid,nlay+1) ! levels altitude
130	REAL rho(ngrid,nlay) ! layers density
131	REAL rhobarz(ngrid,nlay) ! levels density
132	REAL deltaz(ngrid,nlay) ! layers heigth
133	REAL masse(ngrid,nlay) ! layers mass
134	REAL zpspsk(ngrid,nlay) ! Exner function
135
136	REAL zu(ngrid,nlay) ! environment zonal wind
137	REAL zv(ngrid,nlay) ! environment meridional wind
138	REAL zo(ngrid,nlay) ! environment water tracer
139	REAL zva(ngrid,nlay) ! plume zonal wind
140	REAL zua(ngrid,nlay) ! plume meridional wind
141	REAL zoa(ngrid,nlay) ! plume water tracer
142
143	REAL zt(ngrid,nlay) ! T environment
144	REAL zh(ngrid,nlay) ! T,TP environment
145	REAL zthl(ngrid,nlay) ! TP environment
146	REAL ztv(ngrid,nlay) ! TPV environment ?
147	REAL zl(ngrid,nlay) ! ql environment
148
149	REAL zta(ngrid,nlay) !
150	REAL zha(ngrid,nlay) ! TRPV plume
151	REAL ztla(ngrid,nlay) ! TP plume
152	REAL ztva(ngrid,nlay) ! TRPV plume
153	REAL ztva_est(ngrid,nlay) ! TRPV plume (temporary)
154	REAL zqla(ngrid,nlay) ! qv plume
155	REAL zqta(ngrid,nlay) ! qt plume
156
157	REAL wmax(ngrid) ! maximal vertical speed
158	REAL wmax_tmp(ngrid) ! temporary maximal vertical speed
159	REAL wmax_sec(ngrid) ! maximal vertical speed if dry case
160
161	REAL fraca(ngrid,nlay+1) ! updraft fraction
162	REAL f_star(ngrid,nlay+1) ! normalized mass flux
163	REAL entr_star(ngrid,nlay) ! normalized entrainment
164	REAL detr_star(ngrid,nlay) ! normalized detrainment
165	REAL alim_star_tot(ngrid) ! integrated alimentation
166	REAL alim_star(ngrid,nlay) ! normalized alimentation
167	REAL alim_star_clos(ngrid,nlay) ! closure alimentation
168
169	REAL fm(ngrid,nlay+1) ! mass flux
170	REAL entr(ngrid,nlay) ! entrainment
171	REAL detr(ngrid,nlay) ! detrainment
172	REAL f(ngrid) ! mass flux norm
173
174	REAL zdthladj(ngrid,nlay) !
175	REAL lambda ! time relaxation coefficent
176
177	REAL zsortie(ngrid,nlay)
178	REAL zsortie1d(ngrid)
179	REAL susqr2pi, Reuler
180	REAL zf
181	REAL zf2
182	REAL thetath2(ngrid,nlay)
183	REAL wth2(ngrid,nlay)
184	REAL wth3(ngrid,nlay)
185	REAL q2(ngrid,nlay)
186	! FH : probleme de dimensionnement avec l'allocation dynamique
187	! common/comtherm/thetath2,wth2
188	real wq(ngrid,nlay)
189	real wthl(ngrid,nlay)
190	real wthv(ngrid,nlay)
191	real ratqscth(ngrid,nlay)
192	real var
193	real vardiff
194	real ratqsdiff(ngrid,nlay)
195	! niveau de condensation
196	integer nivcon(ngrid)
197	real zcon(ngrid)
198	real CHI
199	real zcon2(ngrid)
200	real pcon(ngrid)
201	real zqsat(ngrid,nlay)
202	real zqsatth(ngrid,nlay)
203	real zlevinter(ngrid)
204	real seuil
205
206	!!! nrlmd le 10/04/2012
207
208	!------Entrees
209	real pbl_tke(ngrid,nlay+1,nbsrf)
210	real pctsrf(ngrid,nbsrf)
211	real omega(ngrid,nlay)
212	real airephy(ngrid)
213	!------Sorties
214	real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid)
215	real therm_tke_max0(ngrid)
216	real env_tke_max0(ngrid)
217	real n2(ngrid),s2(ngrid)
218	real ale_bl_stat(ngrid)
219	real therm_tke_max(ngrid,nlay)
220	real env_tke_max(ngrid,nlay)
221	real alp_bl_det(ngrid)
222	real alp_bl_fluct_m(ngrid)
223	real alp_bl_fluct_tke(ngrid)
224	real alp_bl_conv(ngrid)
225	real alp_bl_stat(ngrid)
226	!------Local
227	integer nsrf
228	real rhobarz0(ngrid) ! Densite au LCL
229	logical ok_lcl(ngrid) ! Existence du LCL des thermiques
230	integer klcl(ngrid) ! Niveau du LCL
231	real interp(ngrid) ! Coef d'interpolation pour le LCL
232	!------Triggering
233	real,parameter :: Su=4e4 ! Surface unite: celle d'un updraft elementaire
234	real,parameter :: hcoef=1. ! Coefficient directeur pour le calcul de s2
235	real,parameter :: hmincoef=0.3 ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 (hmincoef=0.3)
236	real,parameter :: eps1=0.3 ! Fraction de surface occupee par la population 1 : eps1=n1s1/(fraca0Sd)
237	real hmin(ngrid) ! Ordonnee a l'origine pour le calcul de s2
238	real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
239	real,parameter :: zmax_moy_coef=0.33
240	real depth(ngrid) ! Epaisseur moyenne du cumulus
241	real w_max(ngrid) ! Vitesse max statistique
242	real s_max(ngrid)
243	!------Closure
244	real pbl_tke_max(ngrid,nlay) ! Profil de TKE moyenne
245	real pbl_tke_max0(ngrid) ! TKE moyenne au LCL
246	real w_ls(ngrid,nlay) ! Vitesse verticale grande echelle (m/s)
247	real,parameter :: coef_m=1. ! On considere un rendement pour alp_bl_fluct_m
248	real,parameter :: coef_tke=1. ! On considere un rendement pour alp_bl_fluct_tke
249
250	!!! fin nrlmd le 10/04/2012
251
252	! Nouvelles variables pour la convection
253	integer lalim_conv(ngrid)
254	integer n_int(ngrid)
255	real Ale_bl(ngrid)
256	real Alp_bl(ngrid)
257	real alp_int(ngrid)
258	real dp_int(ngrid),zdp
259	real ale_int(ngrid)
260	real fm_tot(ngrid)
261	real wght_th(ngrid,nlay)
262
263	CHARACTER*2 str2
264	CHARACTER*10 str10
265
266	CHARACTER (len=20) :: modname='thermcell_main'
267	CHARACTER (len=80) :: abort_message
268
269	EXTERNAL SCOPY
270
271	!==============================================================================
272	! Initialization
273	!==============================================================================
274
275	seuil = 0.25
276
277	IF (firstcall) THEN
278	IF (iflag_thermals==15.or.iflag_thermals==16) THEN
279	dvimpl = 0
280	dqimpl = -1
281	ELSE
282	dvimpl = 1
283	dqimpl = 1
284	ENDIF
285
286	fm0(:,:) = 0.
287	entr0(:,:) = 0.
288	detr0(:,:) = 0.
289	ENDIF
290
291	fm(:,:) = 0.
292	entr(:,:) = 0.
293	detr(:,:) = 0.
294	f(:) = 0.
295
296	DO ig=1,ngrid
297	f0(ig) = max(f0(ig), 1.e-2)
298	ENDDO
299
300	IF (prt_level.ge.20) then
301	DO ig=1,ngrid
302	print *, 'ig,f0', ig, f0(ig)
303	ENDDO
304	ENDIF
305
306	wmax_tmp(:) = 0.
307
308	!------------------------------------------------------------------------------
309	! Calcul de T,q,ql a partir de Tl et qt dans l environnement
310	!------------------------------------------------------------------------------
311
312	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
313	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
314
315	!------------------------------------------------------------------------------
316	!
317	! --------------------
318	!
319	!
320	! + + + + + + + + + + +
321	!
322	!
323	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
324	! wh,wt,wo ...
325	!
326	! + + + + + + + + + + + zh,zu,zv,zo,rho
327	!
328	!
329	! -------------------- zlev(1)
330	! \\\\\\\\\\\\\\\\\\\\
331	!
332	!------------------------------------------------------------------------------
333	! Calcul des altitudes des couches
334	!------------------------------------------------------------------------------
335
336	DO l=2,nlay
337	zlev(:,l) = 0.5 * (pphi(:,l) + pphi(:,l-1)) / RG
338	ENDDO
339
340	zlev(:,1) = 0.
341	zlev(:,nlay+1) = (2. * pphi(:,nlay) - pphi(:,nlay-1)) / RG
342
343	DO l=1,nlay
344	zlay(:,l) = pphi(:,l)/RG
345	ENDDO
346
347	!------------------------------------------------------------------------------
348	! Calcul de l'epaisseur des couches
349	!------------------------------------------------------------------------------
350
351	DO l=1,nlay
352	deltaz(:,l) = zlev(:,l+1)-zlev(:,l)
353	ENDDO
354
355	!------------------------------------------------------------------------------
356	! Calcul des densites
357	!------------------------------------------------------------------------------
358
359	rho(:,:) = pplay(:,:) / (zpspsk(:,:) * RD * ztv(:,:))
360
361	IF (prt_level.ge.10) THEN
362	write(lunout,*) 'WARNING: thermcell_main rhobarz(:,1)=rho(:,1)'
363	ENDIF
364
365	rhobarz(:,1) = rho(:,1)
366
367	DO l=2,nlay
368	rhobarz(:,l) = 0.5 * (rho(:,l) + rho(:,l-1))
369	ENDDO
370
371	!------------------------------------------------------------------------------
372	! Calcul de la masse
373	!------------------------------------------------------------------------------
374
375	DO l=1,nlay
376	masse(:,l) = (pplev(:,l) - pplev(:,l+1)) / RG
377	ENDDO
378
379	IF (prt_level.ge.1) print *, 'thermcell_main apres initialisation'
380
381	!------------------------------------------------------------------------------
382	!
383	! /\|\
384	! -------- \| F_k+1 -------
385	! ----> D_k
386	! /\|\ <---- E_k , A_k
387	! -------- \| F_k ---------
388	! ----> D_k-1
389	! <---- E_k-1 , A_k-1
390	!
391	!
392	!
393	!
394	!
395	! ---------------------------
396	!
397	! ----- F_lmax+1=0 ---------- \
398	! lmax (zmax) \|
399	! --------------------------- \|
400	! \|
401	! --------------------------- \|
402	! \|
403	! --------------------------- \|
404	! \|
405	! --------------------------- \|
406	! \|
407	! --------------------------- \|
408	! \| E
409	! --------------------------- \| D
410	! \|
411	! --------------------------- \|
412	! \|
413	! --------------------------- \ \|
414	! lalim \| \|
415	! --------------------------- \| \|
416	! \| \|
417	! --------------------------- \| \|
418	! \| A \|
419	! --------------------------- \| \|
420	! \| \|
421	! --------------------------- \| \|
422	! lmin \| \|
423	! ----- F_lmin=0 ------------ / /
424	!
425	! ---------------------------
426	! ////////////////////////////
427	!
428	!------------------------------------------------------------------------------
429
430	!==============================================================================
431	! Calculs initiaux ne faisant pas intervenir les changements de phase
432	!==============================================================================
433
434	!------------------------------------------------------------------------------
435	! 1. alim_star est le profil vertical de l'alimentation a la base du
436	! panache thermique, calcule a partir de la flotabilite de l'air sec
437	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
438	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
439	! panache sec conservatif (e=d=0) alimente selon alim_star
440	! Il s'agit d'un calcul de type CAPE
441	! zmax_sec est utilise pour determiner la geometrie du thermique.
442	!------------------------------------------------------------------------------
443
444	entr_star(:,:) = 0.
445	detr_star(:,:) = 0.
446	alim_star(:,:) = 0.
447
448	alim_star_tot(:) = 0.
449
450	lmin(:) = 1
451
452	!==============================================================================
453	! Calcul du melange et des variables dans le thermique
454	!==============================================================================
455
456	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl, &
457	& po,zl,rhobarz,zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
458	& lalim,f0,detr_star,entr_star,f_star,ztva, &
459	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis, &
460	& lmin,lev_out,lunout1,igout)
461
462	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
463	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
464
465	!==============================================================================
466	! Thermal plumes characteristics: zmax, zmix, wmax
467	!==============================================================================
468
469	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
470	& zlev,lmax,zmax,zmix,wmax,f_star,lev_out)
471
472	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
473	! AB : WARNING: zw2 became its square root in thermcell_height!
474	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
475
476	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
477	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
478	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
479	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
480
481	!==============================================================================
482	! Closure and mass fluxes computation
483	!==============================================================================
484
485	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
486	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
487
488	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
489	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
490
491	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
492	! Choix de la fonction d'alimentation utilisee pour la fermeture.
493	! Apparemment sans importance
494	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
495	alim_star_clos(:,:) = alim_star(:,:)
496	alim_star_clos(:,:) = entr_star(:,:) + alim_star(:,:)
497
498	!------------------------------------------------------------------------------
499	! Closure (dry if iflag_thermals_closure=1, moist if iflag_thermals_closure=2)
500	!------------------------------------------------------------------------------
501
502	IF (iflag_thermals_closure.eq.1) THEN
503	CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, &
504	& lalim,alim_star_clos,f_star, &
505	& zmax_sec,wmax_sec,f,lev_out)
506	ELSEIF (iflag_thermals_closure.eq.2) THEN
507	CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, &
508	& lalim,alim_star,f_star, &
509	& zmax,wmax,f,lev_out)
510	ELSE
511	print *, 'ERROR: no closure had been selected!'
512	call abort
513	ENDIF
514
515	IF (tau_thermals>1.) THEN
516	lambda = exp(-ptimestep/tau_thermals)
517	f0(:) = (1.-lambda) * f(:) + lambda * f0(:)
518	ELSE
519	f0(:) = f(:)
520	ENDIF
521
522	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
523	! Test valable seulement en 1D mais pas genant
524	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
525	IF (.not. (f0(1).ge.0.) ) THEN
526	abort_message = '.not. (f0(1).ge.0.)'
527	print *, 'f0 =', f0(1)
528	CALL abort_physic(modname,abort_message,1)
529	ENDIF
530
531	!------------------------------------------------------------------------------
532	! Mass fluxes
533	!------------------------------------------------------------------------------
534
535	CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
536	& lalim,lmin,lmax,alim_star,entr_star,detr_star, &
537	& f,rhobarz,zlev,zw2,fm,entr,detr,zqla, &
538	& lev_out,lunout1,igout)
539
540	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
541	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
542
543	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
544	! On ne prend pas directement les profils issus des calculs precedents mais on
545	! s'autorise genereusement une relaxation vers ceci avec une constante de temps
546	! tau_thermals (typiquement 1800s).
547	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
548
549	IF (tau_thermals>1.) THEN
550	lambda = exp(-ptimestep/tau_thermals)
551	fm0 = (1.-lambda) * fm + lambda * fm0
552	entr0 = (1.-lambda) * entr + lambda * entr0
553	detr0 = (1.-lambda) * detr + lambda * detr0
554	ELSE
555	fm0(:,:) = fm(:,:)
556	entr0(:,:) = entr(:,:)
557	detr0(:,:) = detr(:,:)
558	ENDIF
559
560	!------------------------------------------------------------------------------
561	! Updraft fraction
562	!------------------------------------------------------------------------------
563
564	DO ig=1,ngrid
565	fraca(ig,1) = 0.
566	fraca(ig,nlay+1) = 0.
567	ENDDO
568
569	DO l=2,nlay
570	DO ig=1,ngrid
571	IF (zw2(ig,l).gt.0.) THEN
572	fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
573	ELSE
574	fraca(ig,l) = 0.
575	ENDIF
576	ENDDO
577	ENDDO
578
579	!==============================================================================
580	! Transport vertical
581	!==============================================================================
582
583	!------------------------------------------------------------------------------
584	! Calcul du transport vertical (de qt et tp)
585	!------------------------------------------------------------------------------
586
587	CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
588	& zthl,zdthladj,zta,lmin,lev_out)
589
590	CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
591	& po,pdoadj,zoa,lmin,lev_out)
592
593	DO l=1,nlay
594	DO ig=1,ngrid
595	pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l)
596	ENDDO
597	ENDDO
598
599	!------------------------------------------------------------------------------
600	! Calcul du transport vertical du moment horizontal
601	!------------------------------------------------------------------------------
602
603	IF (dvimpl.eq.0) THEN
604	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
605	! Calcul du transport de V tenant compte d'echange par gradient
606	! de pression horizontal avec l'environnement
607	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
608	CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca, &
609	& zmax,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
610	ELSE
611	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
612	! Calcul purement conservatif pour le transport de V=(u,v)
613	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
614	CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
615	& zu,pduadj,zua,lmin,lev_out)
616
617	CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
618	& zv,pdvadj,zva,lmin,lev_out)
619	ENDIF
620
621	!==============================================================================
622	! Calculs de diagnostiques pour les sorties
623	!==============================================================================
624
625	IF (sorties) THEN
626
627	!------------------------------------------------------------------------------
628	! Calcul du niveau de condensation
629	!------------------------------------------------------------------------------
630
631	DO ig=1,ngrid
632	nivcon(ig) = 0
633	zcon(ig) = 0.
634	ENDDO
635	!nouveau calcul
636	do ig=1,ngrid
637	CHI = zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
638	pcon(ig) = pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
639	ENDDO
640
641	!IM do k=1,nlay
642	do k=1,nlay-1
643	do ig=1,ngrid
644	IF ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then
645	zcon2(ig) = zlay(ig,k) - (pcon(ig) - pplay(ig,k)) &
646	/ (RG * rho(ig,k)) / 100.
647	ENDIF
648	ENDDO
649	ENDDO
650
651	ierr = 0
652
653	do ig=1,ngrid
654	IF (pcon(ig).le.pplay(ig,nlay)) then
655	zcon2(ig) = zlay(ig,nlay) - (pcon(ig) - pplay(ig,nlay)) &
656	/ (RG * rho(ig,nlay)) / 100.
657	ierr = 1
658	ENDIF
659	ENDDO
660
661	IF (ierr==1) then
662	write(,) 'ERROR: thermal plumes rise the model top'
663	CALL abort
664	ENDIF
665
666	IF (prt_level.ge.1) print*,'14b OK convect8'
667
668	do k=nlay,1,-1
669	do ig=1,ngrid
670	IF (zqla(ig,k).gt.1e-10) then
671	nivcon(ig) = k
672	zcon(ig) = zlev(ig,k)
673	ENDIF
674	ENDDO
675	ENDDO
676
677	IF (prt_level.ge.1) print*,'14c OK convect8'
678
679	!------------------------------------------------------------------------------
680	! Calcul des moments
681	!------------------------------------------------------------------------------
682
683	do l=1,nlay
684	do ig=1,ngrid
685	q2(ig,l) = 0.
686	wth2(ig,l) = 0.
687	wth3(ig,l) = 0.
688	ratqscth(ig,l) = 0.
689	ratqsdiff(ig,l) = 0.
690	ENDDO
691	ENDDO
692
693	IF (prt_level.ge.1) print*,'14d OK convect8'
694
695	do l=1,nlay
696	do ig=1,ngrid
697	zf = fraca(ig,l)
698	zf2 = zf/(1.-zf)
699	thetath2(ig,l) = zf2(ztla(ig,l)-zthl(ig,l))*2
700
701	IF (zw2(ig,l).gt.1.e-10) then
702	wth2(ig,l) = zf2(zw2(ig,l))*2
703	else
704	wth2(ig,l) = 0.
705	ENDIF
706
707	wth3(ig,l) = zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
708	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
709	q2(ig,l) = zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
710	!test: on calcul q2/po=ratqsc
711	ratqscth(ig,l) = sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
712	ENDDO
713	ENDDO
714
715	!------------------------------------------------------------------------------
716	! Calcul des flux: q, thetal et thetav
717	!------------------------------------------------------------------------------
718
719	do l=1,nlay
720	do ig=1,ngrid
721	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
722	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
723	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
724	ENDDO
725	ENDDO
726
727	CALL thermcell_alp(ngrid,nlay,ptimestep, &
728	& pplay,pplev, &
729	& fm0,entr0,lmax, &
730	& Ale_bl,Alp_bl,lalim_conv,wght_th, &
731	& zw2,fraca,pcon, &
732	& rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax, &
733	& pbl_tke,pctsrf,omega,airephy, &
734	& zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,&
735	& n2,s2,ale_bl_stat, &
736	& therm_tke_max,env_tke_max, &
737	& alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke, &
738	& alp_bl_conv,alp_bl_stat)
739
740	!------------------------------------------------------------------------------
741	! Calcul du ratqscdiff
742	!------------------------------------------------------------------------------
743
744	var = 0.
745	vardiff = 0.
746	ratqsdiff(:,:) = 0.
747
748	DO l=1,nlay
749	DO ig=1,ngrid
750	IF (l<=lalim(ig)) THEN
751	var = var + alim_star(ig,l) * zqta(ig,l) * 1000.
752	ENDIF
753	ENDDO
754	ENDDO
755
756	IF (prt_level.ge.1) print*,'14f OK convect8'
757
758	DO l=1,nlay
759	DO ig=1,ngrid
760	IF (l<=lalim(ig)) THEN
761	zf = fraca(ig,l)
762	zf2 = zf / (1.-zf)
763	vardiff = vardiff + alim_star(ig,l) &
764	* (zqta(ig,l) * 1000. - var)**2
765	ENDIF
766	ENDDO
767	ENDDO
768
769	IF (prt_level.ge.1) print*,'14g OK convect8'
770
771	DO l=1,nlay
772	DO ig=1,ngrid
773	ratqsdiff(ig,l) = sqrt(vardiff) / (po(ig,l) * 1000.)
774	ENDDO
775	ENDDO
776
777	ENDIF
778
779
780	RETURN
781	END
782
783
784	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
785	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
786	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
787	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
788	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
789	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
790
791
792	SUBROUTINE test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po, &
793	ztva,zqla,f_star,zw2,comment)
794
795
796	USE print_control_mod, ONLY: prt_level
797
798	IMPLICIT NONE
799
800
801	!==============================================================================
802	! Declaration
803	!==============================================================================
804
805	! inputs:
806	! -------
807
808	INTEGER ngrid
809	INTEGER nlay
810	INTEGER long(ngrid)
811
812	REAL pplev(ngrid,nlay+1)
813	REAL pplay(ngrid,nlay)
814	REAL ztv(ngrid,nlay)
815	REAL po(ngrid,nlay)
816	REAL ztva(ngrid,nlay)
817	REAL zqla(ngrid,nlay)
818	REAL f_star(ngrid,nlay)
819	REAL zw2(ngrid,nlay)
820	REAL seuil
821
822	CHARACTER*21 comment
823
824	! local:
825	! ------
826
827	INTEGER i, k
828
829	!==============================================================================
830	! Test
831	!==============================================================================
832
833	IF (prt_level.ge.1) THEN
834	write(,) 'WARNING: in test, ', comment
835	ENDIF
836
837	return
838
839	! test sur la hauteur des thermiques ...
840	do i=1,ngrid
841	!IMtemp IF (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
842	IF (prt_level.ge.10) then
843	print *, 'WARNING ',comment,' au point ',i,' K= ',long(i)
844	print , ' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
845	do k=1,nlay
846	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)
847	ENDDO
848	ENDIF
849	ENDDO
850
851
852	RETURN
853	END
854
855
856	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
857	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
858	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
859	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
860	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
861	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
862
863
864	SUBROUTINE thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
865	rg,pplev,therm_tke_max)
866
867
868	!==============================================================================
869	! Calcul du transport verticale dans la couche limite en presence
870	! de "thermiques" explicitement representes
871	! calcul du dq/dt une fois qu'on connait les ascendances
872	!
873	! Transport de la TKE par le thermique moyen pour la fermeture en ALP
874	! On transporte pbl_tke pour donner therm_tke
875	!==============================================================================
876
877	USE print_control_mod, ONLY: prt_level
878
879	IMPLICIT NONE
880
881
882	!==============================================================================
883	! Declarations
884	!==============================================================================
885
886	integer ngrid,nlay,nsrf
887
888	real ptimestep
889	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
890	real entr0(ngrid,nlay),rg
891	real therm_tke_max(ngrid,nlay)
892	real detr0(ngrid,nlay)
893
894	real masse(ngrid,nlay),fm(ngrid,nlay+1)
895	real entr(ngrid,nlay)
896	real q(ngrid,nlay)
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	! Calcul du detrainement
907	!------------------------------------------------------------------------------
908
909	do k=1,nlay
910	detr0(:,k) = fm0(:,k) - fm0(:,k+1) + entr0(:,k)
911	masse0(:,k) = (pplev(:,k) - pplev(:,k+1)) / RG
912	ENDDO
913
914	!------------------------------------------------------------------------------
915	! Decalage vertical des entrainements et detrainements.
916	!------------------------------------------------------------------------------
917
918	masse(:,1)=0.5*masse0(:,1)
919	entr(:,1)=0.5*entr0(:,1)
920	detr(:,1)=0.5*detr0(:,1)
921	fm(:,1)=0.
922
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
930	fm(:,nlay+1)=0.
931
932	!!! nrlmd le 16/09/2010
933	! calcul de la valeur dans les ascendances
934	! do ig=1,ngrid
935	! qa(ig,1) = q(ig,1)
936	! ENDDO
937
938	q(:,:)=therm_tke_max(:,:)
939
940	do ig=1,ngrid
941	qa(ig,1)=q(ig,1)
942	ENDDO
943
944	IF (1==1) then
945	do k=2,nlay
946	do ig=1,ngrid
947	IF ((fm(ig,k+1)+detr(ig,k))ptimestep.gt.1.e-5masse(ig,k)) then
948	qa(ig,k) = (fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
949	& / (fm(ig,k+1)+detr(ig,k))
950	else
951	qa(ig,k)=q(ig,k)
952	ENDIF
953
954	! IF (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'
955	! IF (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'
956	ENDDO
957	ENDDO
958
959	!------------------------------------------------------------------------------
960	! Calcul du flux subsident
961	!------------------------------------------------------------------------------
962
963	do k=2,nlay
964	do ig=1,ngrid
965	wqd(ig,k) = fm(ig,k) * q(ig,k)
966	IF (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'
967	ENDDO
968	ENDDO
969
970	do ig=1,ngrid
971	wqd(ig,1) = 0.
972	wqd(ig,nlay+1) = 0.
973	ENDDO
974
975	!------------------------------------------------------------------------------
976	! Calcul des tendances
977	!------------------------------------------------------------------------------
978
979	do k=1,nlay
980	do ig=1,ngrid
981	q(ig,k) = q(ig,k) + ptimestep / masse(ig,k) &
982	& * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k) &
983	& - wqd(ig,k) + wqd(ig,k+1))
984	ENDDO
985	ENDDO
986	ENDIF
987
988	therm_tke_max(:,:) = q(:,:)
989
990
991	RETURN
992	END
993

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

Download in other formats:

Original Format