Context Navigation

thermcell_main.F90 @ 2101

Last change on this file since 2101 was 2101, checked in by aboissinot, 6 years ago
Fix a bug in thermcell_alim.F90 where vertical and horizontal loops must be inverted. In thermal plume model, arrays size declaration is standardised (no longer done thanks to dimphy module but by way of arguments). Clean up some thermal plume model routines (remove uselesss variables...)
File size: 37.4 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),env_tke_max0(ngrid)
216	real n2(ngrid),s2(ngrid)
217	real ale_bl_stat(ngrid)
218	real therm_tke_max(ngrid,nlay),env_tke_max(ngrid,nlay)
219	real alp_bl_det(ngrid),alp_bl_fluct_m(ngrid),alp_bl_fluct_tke(ngrid),alp_bl_conv(ngrid),alp_bl_stat(ngrid)
220	!------Local
221	integer nsrf
222	real rhobarz0(ngrid) ! Densite au LCL
223	logical ok_lcl(ngrid) ! Existence du LCL des thermiques
224	integer klcl(ngrid) ! Niveau du LCL
225	real interp(ngrid) ! Coef d'interpolation pour le LCL
226	!------Triggering
227	real Su ! Surface unite: celle d'un updraft elementaire
228	parameter(Su=4e4)
229	real hcoef ! Coefficient directeur pour le calcul de s2
230	parameter(hcoef=1)
231	real hmincoef ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2
232	parameter(hmincoef=0.3)
233	real eps1 ! Fraction de surface occupee par la population 1 : eps1=n1s1/(fraca0Sd)
234	parameter(eps1=0.3)
235	real hmin(ngrid) ! Ordonnee a l'origine pour le calcul de s2
236	real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
237	real zmax_moy_coef
238	parameter(zmax_moy_coef=0.33)
239	real depth(ngrid) ! Epaisseur moyenne du cumulus
240	real w_max(ngrid) ! Vitesse max statistique
241	real s_max(ngrid)
242	!------Closure
243	real pbl_tke_max(ngrid,nlay) ! Profil de TKE moyenne
244	real pbl_tke_max0(ngrid) ! TKE moyenne au LCL
245	real w_ls(ngrid,nlay) ! Vitesse verticale grande echelle (m/s)
246	real coef_m ! On considere un rendement pour alp_bl_fluct_m
247	parameter(coef_m=1.)
248	real coef_tke ! On considere un rendement pour alp_bl_fluct_tke
249	parameter(coef_tke=1.)
250
251	!!! fin nrlmd le 10/04/2012
252
253	! Nouvelles variables pour la convection
254	real Ale_bl(ngrid)
255	real Alp_bl(ngrid)
256	real alp_int(ngrid),dp_int(ngrid),zdp
257	real ale_int(ngrid)
258	integer n_int(ngrid)
259	real fm_tot(ngrid)
260	real wght_th(ngrid,nlay)
261	integer lalim_conv(ngrid)
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	ELSE
530	print *, 'f0 =', f0(1)
531	ENDIF
532
533	!------------------------------------------------------------------------------
534	! Mass fluxes
535	!------------------------------------------------------------------------------
536
537	CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
538	& lalim,lmin,lmax,alim_star,entr_star,detr_star, &
539	& f,rhobarz,zlev,zw2,fm,entr,detr,zqla, &
540	& lev_out,lunout1,igout)
541
542	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
543	! CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
544
545	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
546	! On ne prend pas directement les profils issus des calculs precedents mais on
547	! s'autorise genereusement une relaxation vers ceci avec une constante de temps
548	! tau_thermals (typiquement 1800s).
549	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
550
551	IF (tau_thermals>1.) THEN
552	lambda = exp(-ptimestep/tau_thermals)
553	fm0 = (1.-lambda) * fm + lambda * fm0
554	entr0 = (1.-lambda) * entr + lambda * entr0
555	detr0 = (1.-lambda) * detr + lambda * detr0
556	ELSE
557	fm0(:,:) = fm(:,:)
558	entr0(:,:) = entr(:,:)
559	detr0(:,:) = detr(:,:)
560	ENDIF
561
562	!------------------------------------------------------------------------------
563	! Updraft fraction
564	!------------------------------------------------------------------------------
565
566	DO ig=1,ngrid
567	fraca(ig,1) = 0.
568	fraca(ig,nlay+1) = 0.
569	ENDDO
570
571	DO l=2,nlay
572	DO ig=1,ngrid
573	IF (zw2(ig,l).gt.0.) THEN
574	fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
575	ELSE
576	fraca(ig,l) = 0.
577	ENDIF
578	ENDDO
579	ENDDO
580
581	!==============================================================================
582	! Transport vertical
583	!==============================================================================
584
585	!------------------------------------------------------------------------------
586	! Calcul du transport vertical (de qt et tp)
587	!------------------------------------------------------------------------------
588
589	CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
590	& zthl,zdthladj,zta,lmin,lev_out)
591
592	CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
593	& po,pdoadj,zoa,lmin,lev_out)
594
595	DO l=1,nlay
596	DO ig=1,ngrid
597	pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l)
598	ENDDO
599	ENDDO
600
601	!------------------------------------------------------------------------------
602	! Calcul du transport vertical du moment horizontal
603	!------------------------------------------------------------------------------
604
605	IF (dvimpl.eq.0) THEN
606	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
607	! Calcul du transport de V tenant compte d'echange par gradient
608	! de pression horizontal avec l'environnement
609	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
610	CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca, &
611	& zmax,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
612	ELSE
613	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
614	! Calcul purement conservatif pour le transport de V=(u,v)
615	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
616	CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
617	& zu,pduadj,zua,lmin,lev_out)
618
619	CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
620	& zv,pdvadj,zva,lmin,lev_out)
621	ENDIF
622
623	!==============================================================================
624	! Calculs de diagnostiques pour les sorties
625	!==============================================================================
626
627	IF (sorties) THEN
628
629	!------------------------------------------------------------------------------
630	! Calcul du niveau de condensation
631	!------------------------------------------------------------------------------
632
633	if (prt_level.ge.1) print*,'14a OK convect8'
634
635	do ig=1,ngrid
636	nivcon(ig) = 0
637	zcon(ig) = 0.
638	enddo
639	!nouveau calcul
640	do ig=1,ngrid
641	CHI = zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
642	pcon(ig) = pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
643	enddo
644
645	!IM do k=1,nlay
646	do k=1,nlay-1
647	do ig=1,ngrid
648	if ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then
649	zcon2(ig) = zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
650	endif
651	enddo
652	enddo
653
654	ierr = 0
655
656	do ig=1,ngrid
657	if (pcon(ig).le.pplay(ig,nlay)) then
658	zcon2(ig) = zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
659	ierr = 1
660	endif
661	enddo
662
663	if (ierr==1) then
664	write(,) 'ERROR: thermal plumes rise the model top'
665	CALL abort
666	endif
667
668	if (prt_level.ge.1) print*,'14b OK convect8'
669
670	do k=nlay,1,-1
671	do ig=1,ngrid
672	if (zqla(ig,k).gt.1e-10) then
673	nivcon(ig) = k
674	zcon(ig) = zlev(ig,k)
675	endif
676	enddo
677	enddo
678
679	if (prt_level.ge.1) print*,'14c OK convect8'
680
681	!------------------------------------------------------------------------------
682	! Calcul des moments
683	!------------------------------------------------------------------------------
684
685	do l=1,nlay
686	do ig=1,ngrid
687	q2(ig,l) = 0.
688	wth2(ig,l) = 0.
689	wth3(ig,l) = 0.
690	ratqscth(ig,l) = 0.
691	ratqsdiff(ig,l) = 0.
692	enddo
693	enddo
694
695	if (prt_level.ge.1) print*,'14d OK convect8'
696
697	do l=1,nlay
698	do ig=1,ngrid
699	zf = fraca(ig,l)
700	zf2 = zf/(1.-zf)
701	thetath2(ig,l) = zf2(ztla(ig,l)-zthl(ig,l))*2
702
703	if (zw2(ig,l).gt.1.e-10) then
704	wth2(ig,l) = zf2(zw2(ig,l))*2
705	else
706	wth2(ig,l) = 0.
707	endif
708
709	wth3(ig,l) = zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
710	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
711	q2(ig,l) = zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
712	!test: on calcul q2/po=ratqsc
713	ratqscth(ig,l) = sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
714	enddo
715	enddo
716
717	!------------------------------------------------------------------------------
718	! Calcul des flux: q, thetal et thetav
719	!------------------------------------------------------------------------------
720
721	do l=1,nlay
722	do ig=1,ngrid
723	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
724	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
725	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
726	enddo
727	enddo
728
729	CALL thermcell_alp(ngrid,nlay,ptimestep, &
730	& pplay,pplev, &
731	& fm0,entr0,lmax, &
732	& Ale_bl,Alp_bl,lalim_conv,wght_th, &
733	& zw2,fraca, &
734	& pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax, &
735	& pbl_tke,pctsrf,omega,airephy, &
736	& zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0, &
737	& n2,s2,ale_bl_stat, &
738	& therm_tke_max,env_tke_max, &
739	& alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke, &
740	& alp_bl_conv,alp_bl_stat)
741
742	!------------------------------------------------------------------------------
743	! Calcul du ratqscdiff
744	!------------------------------------------------------------------------------
745
746	var = 0.
747	vardiff = 0.
748	ratqsdiff(:,:) = 0.
749
750	DO l=1,nlay
751	DO ig=1,ngrid
752	IF (l<=lalim(ig)) THEN
753	var = var + alim_star(ig,l) * zqta(ig,l) * 1000.
754	ENDIF
755	ENDDO
756	ENDDO
757
758	if (prt_level.ge.1) print*,'14f OK convect8'
759
760	DO l=1,nlay
761	DO ig=1,ngrid
762	IF (l<=lalim(ig)) THEN
763	zf = fraca(ig,l)
764	zf2 = zf / (1.-zf)
765	vardiff = vardiff + alim_star(ig,l) * (zqta(ig,l) * 1000. - var)**2
766	ENDIF
767	ENDDO
768	ENDDO
769
770	IF (prt_level.ge.1) print*,'14g OK convect8'
771
772	DO l=1,nlay
773	DO ig=1,ngrid
774	ratqsdiff(ig,l) = sqrt(vardiff) / (po(ig,l) * 1000.)
775	ENDDO
776	ENDDO
777
778	ENDIF
779
780
781	RETURN
782	END
783
784
785	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
786	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
787	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
788	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
789	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
790	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
791
792
793	SUBROUTINE test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po, &
794	ztva,zqla,f_star,zw2,comment)
795
796
797	USE print_control_mod, ONLY: prt_level
798
799	IMPLICIT NONE
800
801
802	!==============================================================================
803	! Declaration
804	!==============================================================================
805
806	! inputs:
807	! -------
808
809	INTEGER ngrid
810	INTEGER nlay
811	INTEGER long(ngrid)
812
813	REAL pplev(ngrid,nlay+1)
814	REAL pplay(ngrid,nlay)
815	REAL ztv(ngrid,nlay)
816	REAL po(ngrid,nlay)
817	REAL ztva(ngrid,nlay)
818	REAL zqla(ngrid,nlay)
819	REAL f_star(ngrid,nlay)
820	REAL zw2(ngrid,nlay)
821	REAL seuil
822
823	CHARACTER*21 comment
824
825	! local:
826	! ------
827
828	INTEGER i, k
829
830	!==============================================================================
831	! Test
832	!==============================================================================
833
834	IF (prt_level.ge.1) THEN
835	write(,) 'WARNING: in test, ', comment
836	ENDIF
837
838	return
839
840	! test sur la hauteur des thermiques ...
841	do i=1,ngrid
842	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
843	if (prt_level.ge.10) then
844	print *, 'WARNING ',comment,' au point ',i,' K= ',long(i)
845	print , ' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
846	do k=1,nlay
847	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)
848	enddo
849	endif
850	enddo
851
852
853	RETURN
854	END
855
856
857	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
858	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
859	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
860	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
861	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
862	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
863
864
865	SUBROUTINE thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
866	rg,pplev,therm_tke_max)
867
868
869	!==============================================================================
870	! Calcul du transport verticale dans la couche limite en presence
871	! de "thermiques" explicitement representes
872	! calcul du dq/dt une fois qu'on connait les ascendances
873	!
874	! Transport de la TKE par le thermique moyen pour la fermeture en ALP
875	! On transporte pbl_tke pour donner therm_tke
876	!==============================================================================
877
878	USE print_control_mod, ONLY: prt_level
879
880	IMPLICIT NONE
881
882
883	!==============================================================================
884	! Declarations
885	!==============================================================================
886
887	integer ngrid,nlay,nsrf
888
889	real ptimestep
890	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
891	real entr0(ngrid,nlay),rg
892	real therm_tke_max(ngrid,nlay)
893	real detr0(ngrid,nlay)
894
895	real masse(ngrid,nlay),fm(ngrid,nlay+1)
896	real entr(ngrid,nlay)
897	real q(ngrid,nlay)
898
899	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
900
901	real zzm
902
903	integer ig,k
904	integer isrf
905
906	!------------------------------------------------------------------------------
907	! Calcul du detrainement
908	!------------------------------------------------------------------------------
909
910	do k=1,nlay
911	detr0(:,k) = fm0(:,k) - fm0(:,k+1) + entr0(:,k)
912	masse0(:,k) = (pplev(:,k) - pplev(:,k+1)) / RG
913	enddo
914
915	!------------------------------------------------------------------------------
916	! Decalage vertical des entrainements et detrainements.
917	!------------------------------------------------------------------------------
918
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
924	do k=1,nlay-1
925	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
926	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
927	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
928	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
929	enddo
930
931	fm(:,nlay+1)=0.
932
933	!!! nrlmd le 16/09/2010
934	! calcul de la valeur dans les ascendances
935	! do ig=1,ngrid
936	! qa(ig,1) = q(ig,1)
937	! enddo
938
939	q(:,:)=therm_tke_max(:,:)
940
941	do ig=1,ngrid
942	qa(ig,1)=q(ig,1)
943	enddo
944
945	if (1==1) then
946	do k=2,nlay
947	do ig=1,ngrid
948	if ((fm(ig,k+1)+detr(ig,k))ptimestep.gt.1.e-5masse(ig,k)) then
949	qa(ig,k) = (fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
950	& / (fm(ig,k+1)+detr(ig,k))
951	else
952	qa(ig,k)=q(ig,k)
953	endif
954
955	! if (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'
956	! if (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'
957	enddo
958	enddo
959
960	!------------------------------------------------------------------------------
961	! Calcul du flux subsident
962	!------------------------------------------------------------------------------
963
964	do k=2,nlay
965	do ig=1,ngrid
966	wqd(ig,k) = fm(ig,k) * q(ig,k)
967	if (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'
968	enddo
969	enddo
970
971	do ig=1,ngrid
972	wqd(ig,1) = 0.
973	wqd(ig,nlay+1) = 0.
974	enddo
975
976	!------------------------------------------------------------------------------
977	! Calcul des tendances
978	!------------------------------------------------------------------------------
979
980	do k=1,nlay
981	do ig=1,ngrid
982	q(ig,k) = q(ig,k) + ptimestep / masse(ig,k) &
983	& * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k) &
984	& - wqd(ig,k) + wqd(ig,k+1))
985	enddo
986	enddo
987	endif
988
989	therm_tke_max(:,:) = q(:,:)
990
991
992	RETURN
993	END
994

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

Download in other formats:

Original Format