Context Navigation

thermcell_main.F90 @ 2078

Last change on this file since 2078 was 2072, checked in by aslmd, 7 years ago
ioipsl is not used in thermcell_main et thermcell_plume
File size: 36.8 KB

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

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

Download in other formats:

Original Format