Context Navigation

thermcell_main_mars.F90 @ 1088

Last change on this file since 1088 was 1033, checked in by aslmd, 12 years ago
LMDZ.MARS additional comments of thermal plume model. moved tuning variables in comtherm_h
File size: 42.7 KB

Line
1	!=======================================================================
2	! THERMCELL_MAIN_MARS
3	!=======================================================================
4	! This routine is called by calltherm_interface and is inside a sub-timestep
5	! loop. It computes thermals properties from parametrized entrainment and
6	! detrainment rate as well as the source profile.
7	! Mass flux are then computed and temperature and CO2 MMR are transported.
8	!=======================================================================
9	! Author : A. Colaitis 2011-01-05 (with updates 2011-2013)
10	! after C. Rio and F. Hourdin
11	! Institution : Laboratoire de Meteorologie Dynamique (LMD) Paris, France
12	! -----------------------------------------------------------------------
13	! Corresponding author : A. Spiga aymeric.spiga_AT_upmc.fr
14	! -----------------------------------------------------------------------
15	! ASSOCIATED FILES
16	! --> calltherm_interface.F90
17	! --> thermcell_dqup.F90
18	! --> comtherm_h.F90
19	!=======================================================================
20	! Reference paper:
21	! A. Colaïtis, A. Spiga, F. Hourdin, C. Rio, F. Forget, and E. Millour.
22	! A thermal plume model for the Martian convective boundary layer.
23	! Journal of Geophysical Research (Planets), 118:1468-1487, July 2013.
24	! http://dx.doi.org/10.1002/jgre.20104
25	! http://arxiv.org/abs/1306.6215
26	! -----------------------------------------------------------------------
27	! Reference paper for terrestrial plume model:
28	! C. Rio and F. Hourdin.
29	! A thermal plume model for the convective boundary layer : Representation of cumulus clouds.
30	! Journal of the Atmospheric Sciences, 65:407-425, 2008.
31	! -----------------------------------------------------------------------
32
33	SUBROUTINE thermcell_main_mars(ngrid,nlayer,nq &
34	& ,tracer,igcm_co2 &
35	& ,ptimestep &
36	& ,pplay,pplev,pphi,zlev,zlay &
37	& ,pu,pv,pt,pq,pq2 &
38	& ,pdtadj,pdqadj &
39	& ,fm,entr,detr,lmax,zmax &
40	& ,zw2,fraca &
41	& ,zpopsk,heatFlux,heatFlux_down &
42	& ,buoyancyOut, buoyancyEst)
43
44	USE comtherm_h
45
46	IMPLICIT NONE
47
48	!=======================================================================
49
50	! SHARED VARIABLES. This needs adaptations in another climate model.
51	#include "comcstfi.h" !contains physical constant values such as
52	! "g" : gravitational acceleration (m.s-2)
53	! "r" : recuced gas constant (J.K-1.mol-1)
54
55	! ============== INPUTS ==============
56
57	INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
58	INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
59	INTEGER, INTENT(IN) :: nq ! number of tracer species
60	LOGICAL, INTENT(IN) :: tracer ! =.true. if tracers are present and to be transported
61	INTEGER, INTENT(IN) :: igcm_co2 ! index of the CO2 tracer in mixing ratio array
62	! --> 0 if no tracer is CO2 (or no tracer at all)
63	! --> this prepares special treatment for polar night mixing
64	REAL, INTENT(IN) :: ptimestep !subtimestep (s)
65	REAL, INTENT(IN) :: pt(ngrid,nlayer) !temperature (K)
66	REAL, INTENT(IN) :: pu(ngrid,nlayer) !u component of the wind (ms-1)
67	REAL, INTENT(IN) :: pv(ngrid,nlayer) !v component of the wind (ms-1)
68	REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) !tracer concentration (kg/kg)
69	REAL, INTENT(IN) :: pq2(ngrid,nlayer) ! Turbulent Kinetic Energy
70	REAL, INTENT(IN) :: pplay(ngrid,nlayer) !Pressure at the middle of the layers (Pa)
71	REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) !intermediate pressure levels (Pa)
72	REAL, INTENT(IN) :: pphi(ngrid,nlayer) !Geopotential at the middle of the layers (m2s-2)
73	REAL, INTENT(IN) :: zlay(ngrid,nlayer) ! altitude at the middle of the layers
74	REAL, INTENT(IN) :: zlev(ngrid,nlayer+1) ! altitude at layer boundaries
75
76	! ============== OUTPUTS ==============
77
78	! TEMPERATURE
79	REAL, INTENT(OUT) :: pdtadj(ngrid,nlayer) !temperature change from thermals dT/dt (K/s)
80
81	! DIAGNOSTICS
82	REAL, INTENT(OUT) :: zw2(ngrid,nlayer+1) ! vertical velocity (m/s)
83	REAL, INTENT(OUT) :: heatFlux(ngrid,nlayer) ! interface heatflux
84	REAL, INTENT(OUT) :: heatFlux_down(ngrid,nlayer) ! interface heat flux from downdraft
85
86	! ============== LOCAL ================
87	REAL :: pdqadj(ngrid,nlayer,nq) !tracer change from thermals dq/dt, only for CO2 (the rest can be advected outside of the loop)
88
89	! dummy variables when output not needed :
90
91	REAL :: buoyancyOut(ngrid,nlayer) ! interlayer buoyancy term
92	REAL :: buoyancyEst(ngrid,nlayer) ! interlayer estimated buoyancy term
93
94	! ============== LOCAL ================
95
96	INTEGER ig,k,l,ll,iq
97	INTEGER lmax(ngrid),lmin(ngrid),lalim(ngrid)
98	REAL zmax(ngrid)
99	REAL ztva(ngrid,nlayer),zw_est(ngrid,nlayer+1),ztva_est(ngrid,nlayer)
100	REAL zh(ngrid,nlayer)
101	REAL zdthladj(ngrid,nlayer)
102	REAL zdthladj_down(ngrid,nlayer)
103	REAL ztvd(ngrid,nlayer)
104	REAL ztv(ngrid,nlayer)
105	REAL zu(ngrid,nlayer),zv(ngrid,nlayer),zo(ngrid,nlayer)
106	REAL zva(ngrid,nlayer)
107	REAL zua(ngrid,nlayer)
108
109	REAL zta(ngrid,nlayer)
110	REAL fraca(ngrid,nlayer+1)
111	REAL q2(ngrid,nlayer)
112	REAL rho(ngrid,nlayer),rhobarz(ngrid,nlayer),masse(ngrid,nlayer)
113	REAL zpopsk(ngrid,nlayer)
114
115	REAL wmax(ngrid)
116	REAL fm(ngrid,nlayer+1),entr(ngrid,nlayer),detr(ngrid,nlayer)
117
118	REAL fm_down(ngrid,nlayer+1)
119
120	REAL ztla(ngrid,nlayer)
121
122	REAL f_star(ngrid,nlayer+1),entr_star(ngrid,nlayer)
123	REAL detr_star(ngrid,nlayer)
124	REAL alim_star_tot(ngrid)
125	REAL alim_star(ngrid,nlayer)
126	REAL alim_star_clos(ngrid,nlayer)
127	REAL f(ngrid)
128
129	REAL detrmod(ngrid,nlayer)
130
131	REAL teta_th_int(ngrid,nlayer)
132	REAL teta_env_int(ngrid,nlayer)
133	REAL teta_down_int(ngrid,nlayer)
134
135	CHARACTER (LEN=80) :: abort_message
136	INTEGER ndt
137
138	! ============= PLUME VARIABLES ============
139
140	REAL w_est(ngrid,nlayer+1)
141	REAL wa_moy(ngrid,nlayer+1)
142	REAL wmaxa(ngrid)
143	REAL zdz,zbuoy(ngrid,nlayer),zw2m
144	LOGICAL activecell(ngrid),activetmp(ngrid)
145	INTEGER tic
146
147	! ==========================================
148
149	! ============= HEIGHT VARIABLES ===========
150
151	REAL num(ngrid)
152	REAL denom(ngrid)
153	REAL zlevinter(ngrid)
154	INTEGER zlmax
155
156	! =========================================
157
158	! ============= CLOSURE VARIABLES =========
159
160	REAL zdenom(ngrid)
161	REAL alim_star2(ngrid)
162	REAL alim_star_tot_clos(ngrid)
163	INTEGER llmax
164
165	! =========================================
166
167	! ============= FLUX2 VARIABLES ===========
168
169	INTEGER ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
170	INTEGER ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
171	REAL zfm
172	REAL f_old,ddd0,eee0,ddd,eee,zzz
173	REAL fomass_max,alphamax
174
175	! =========================================
176
177	! ============== Theta_M Variables ========
178
179	REAL m_co2, m_noco2, A , B
180	SAVE A, B
181	REAL zhc(ngrid,nlayer)
182	REAL ratiom(ngrid,nlayer)
183
184	! =========================================
185
186	!-----------------------------------------------------------------------
187	! initialization:
188	! ---------------
189
190	entr(:,:)=0. ! entrainment mass flux
191	detr(:,:)=0. ! detrainment mass flux
192	fm(:,:)=0. ! upward mass flux
193	zhc(:,:)=pt(:,:)/zpopsk(:,:) ! potential temperature
194	ndt=1
195
196	!.......................................................................
197	! Special treatment for co2:
198	!.......................................................................
199	! **********************************************************************
200	! In order to take into account the effect of vertical molar mass
201	! gradient on convection, we define modified theta that depends
202	! on the mass mixing ratio of Co2 in the cell.
203	! See for details:
204	!
205	! Forget, F. and Millour, E. et al. "Non condensable gas enrichment and depletion
206	! in the martian polar regions", third international workshop on the Mars Atmosphere:
207	! Modeling and Observations, 1447, 9106. year: 2008
208	!
209	! This is especially important for modelling polar convection.
210	! **********************************************************************
211	if (igcm_co2.ne.0) then
212
213	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
214	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
215	! Compute A and B coefficient use to compute
216	! mean molecular mass Mair defined by
217	! 1/Mair = q(igcm_co2)/m_co2 + (1-q(igcm_co2))/m_noco2
218	! 1/Mair = A*q(igcm_co2) + B
219	A =(1/m_co2 - 1/m_noco2)
220	B=1/m_noco2
221
222	! Special case if one of the tracers is CO2 gas
223	DO l=1,nlayer
224	DO ig=1,ngrid
225	ztv(ig,l) = zhc(ig,l)(Apq(ig,l,igcm_co2)+B)
226	ENDDO
227	ENDDO
228	else
229	ztv(:,:)=zhc(:,:)
230	end if
231
232	!------------------------------------------------------------------------
233	! where are the different quantities defined ?
234	!------------------------------------------------------------------------
235	! --------------------
236	!
237	!
238	! + + + + + + + + + + +
239	!
240	!
241	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
242	! wh,wt,wo ...
243	!
244	! + + + + + + + + + + + zh,zu,zv,zo,rho
245	!
246	!
247	! -------------------- zlev(1)
248	! \\\\\\\\\\\\\\\\\\\\
249	!
250	!
251
252	!-----------------------------------------------------------------------
253	! Densities at layer and layer interface (see above), mass:
254	!-----------------------------------------------------------------------
255
256	rho(:,:)=pplay(:,:)/(r*pt(:,:))
257
258	rhobarz(:,1)=rho(:,1)
259
260	do l=2,nlayer
261	rhobarz(:,l)=pplev(:,l)/(r0.5(pt(:,l)+pt(:,l-1)))
262	enddo
263
264	! mass computation
265	do l=1,nlayer
266	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
267	enddo
268
269
270	!-----------------------------------------------------------------
271	! Schematic representation of an updraft:
272	!------------------------------------------------------------------
273	!
274	! /\|\
275	! -------- \| F_k+1 -------
276	! ----> D_k
277	! /\|\ <---- E_k , A_k
278	! -------- \| F_k ---------
279	! ----> D_k-1
280	! <---- E_k-1 , A_k-1
281	!
282	!
283	! ---------------------------
284	!
285	! ----- F_lmax+1=0 ---------- \
286	! lmax (zmax) \|
287	! --------------------------- \|
288	! \|
289	! --------------------------- \|
290	! \|
291	! --------------------------- \|
292	! \|
293	! --------------------------- \|
294	! \|
295	! --------------------------- \|
296	! \| E
297	! --------------------------- \| D
298	! \|
299	! --------------------------- \|
300	! \|
301	! --------------------------- \ \|
302	! lalim \| \|
303	! --------------------------- \| \|
304	! \| \|
305	! --------------------------- \| \|
306	! \| A \|
307	! --------------------------- \| \|
308	! \| \|
309	! --------------------------- \| \|
310	! lmin (=1 pour le moment) \| \|
311	! ----- F_lmin=0 ------------ / /
312	!
313	! ---------------------------
314	! //////////////////////////
315	!
316
317	!=============================================================================
318	! Mars version: no phase change is considered, we use a "dry" definition
319	! for the potential temperature.
320	!=============================================================================
321
322	!------------------------------------------------------------------
323	! 1. alim_star is the source layer vertical profile in the lowest layers
324	! of the thermal plume. Computed from the air buoyancy
325	! 2. lmin and lalim are the indices of begining and end of source profile
326	!------------------------------------------------------------------
327	!
328	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
329	lmin=1
330
331	!-----------------------------------------------------------------------------
332	! 3. wmax and zmax are maximum vertical velocity and altitude of a
333	! conservative plume (entrainment = detrainment = 0) using only
334	! the source layer. This is a CAPE computation used for determining
335	! the closure mass flux.
336	!-----------------------------------------------------------------------------
337
338	! ===========================================================================
339	! ===================== PLUME ===============================================
340	! ===========================================================================
341
342	! Initialization
343	ztva(:,:)=ztv(:,:) ! temperature in the updraft = temperature of the env.
344	ztva_est(:,:)=ztva(:,:) ! estimated temp. in the updraft
345	ztla(:,:)=0. !intermediary variable
346	zdz=0. !layer thickness
347	zbuoy(:,:)=0. !buoyancy
348	w_est(:,:)=0. !estimated vertical velocity
349	f_star(:,:)=0. !non-dimensional upward mass flux f*
350	wa_moy(:,:)=0. !vertical velocity
351
352	! Some more initializations
353	wmaxa(:)=0.
354	lalim(:)=1
355
356	!-------------------------------------------------------------------------
357	! We consider as an activecell columns where the two first layers are
358	! convectively unstable
359	! When it is the case, we compute the source layer profile (alim_star)
360	! see paper appendix 4.1 for details on the source layer
361	!-------------------------------------------------------------------------
362
363	activecell(:)=ztv(:,1)>ztv(:,2)
364	do ig=1,ngrid
365	if (ztv(ig,1)>=(ztv(ig,2))) then
366	alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.) &
367	& *sqrt(zlev(ig,2))
368	lalim(ig)=2
369	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1)
370	endif
371	enddo
372
373	do l=2,nlayer-1
374	do ig=1,ngrid
375	if (ztv(ig,l)>(ztv(ig,l+1)) .and. ztv(ig,1)>=ztv(ig,l) &
376	& .and. (alim_star(ig,l-1).ne. 0.)) then
377	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
378	& *sqrt(zlev(ig,l+1))
379	lalim(ig)=l+1
380	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
381	endif
382	enddo
383	enddo
384	do l=1,nlayer
385	do ig=1,ngrid
386	if (alim_star_tot(ig) > 1.e-10 ) then
387	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
388	endif
389	enddo
390	enddo
391
392	alim_star_tot(:)=1.
393
394	! We compute the initial squared velocity (zw2) and non-dimensional upward mass flux
395	! (f_star) in the first and second layer from the source profile.
396
397	do ig=1,ngrid
398	if (activecell(ig)) then
399	ztla(ig,1)=ztv(ig,1)
400	f_star(ig,1)=0.
401	f_star(ig,2)=alim_star(ig,1)
402	zw2(ig,2)=2.g(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
403	& *(zlev(ig,2)-zlev(ig,1)) &
404	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) !0.4=von Karman constant
405	w_est(ig,2)=zw2(ig,2)
406	endif
407	enddo
408
409	!==============================================================================
410	!==============================================================================
411	!==============================================================================
412	! LOOP ON VERTICAL LEVELS
413	!==============================================================================
414	do l=2,nlayer-1
415	!==============================================================================
416	!==============================================================================
417	!==============================================================================
418
419
420	! is the thermal plume still active ?
421	do ig=1,ngrid
422	activecell(ig)=activecell(ig) &
423	& .and. zw2(ig,l)>1.e-10 &
424	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
425	enddo
426
427	!---------------------------------------------------------------------------
428	!
429	! .I. INITIALIZATION
430	!
431	! Computations of the temperature and buoyancy properties in layer l,
432	! without accounting for entrainment and detrainment. We are therefore
433	! assuming constant temperature in the updraft
434	!
435	! This computation yields an estimation of the buoyancy (zbuoy) and thereforce
436	! an estimation of the velocity squared (w_est)
437	!---------------------------------------------------------------------------
438
439	do ig=1,ngrid
440	if(activecell(ig)) then
441	ztva_est(ig,l)=ztla(ig,l-1)
442
443	zdz=zlev(ig,l+1)-zlev(ig,l)
444	zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
445
446	! Estimated vertical velocity squared
447	! (discretized version of equation 12 in paragraph 40 of paper)
448
449	if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
450	w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.zdza1zbuoy(ig,l)-2.zdzw_est(ig,l)b1 &
451	& -2.(1.-omega)zdzw_est(ig,l)ae(a1zbuoy(ig,l)/w_est(ig,l)-b1)**be)
452	else
453	w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.zdza1invzbuoy(ig,l)-2.zdzw_est(ig,l)b1inv)
454	endif
455	if (w_est(ig,l+1).lt.0.) then
456	w_est(ig,l+1)=zw2(ig,l)
457	endif
458	endif
459	enddo
460
461	!-------------------------------------------------
462	! Compute corresponding non-dimensional (ND) entrainment and detrainment rates
463	!-------------------------------------------------
464
465	do ig=1,ngrid
466	if (activecell(ig)) then
467
468	zw2m=w_est(ig,l+1)
469
470	if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
471
472	! ND entrainment rate, see equation 16 of paper (paragraph 43)
473
474	entr_star(ig,l)=f_star(ig,l)zdz &
475	& MAX(0.,ae(a1(zbuoy(ig,l)/zw2m)-b1)**be)
476	else
477	entr_star(ig,l)=0.
478	endif
479
480	if(zbuoy(ig,l) .gt. 0.) then
481	if(l .lt. lalim(ig)) then
482
483	detr_star(ig,l)=0.
484	else
485
486	! ND detrainment rate, see paragraph 44 of paper
487
488	detr_star(ig,l) = f_star(ig,l)zdz &
489	& ad
490
491	endif
492	else
493	detr_star(ig,l)=f_star(ig,l)zdz &
494	& MAX(ad,bd*zbuoy(ig,l)/zw2m)
495
496	endif
497
498	! If we are still in the source layer, we define the source layer entr. rate (alim_star) as the
499	! maximum between the source entrainment rate and the estimated entrainment rate.
500
501	if (l.lt.lalim(ig)) then
502	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
503	entr_star(ig,l)=0.
504	endif
505
506	! Compute the non-dimensional upward mass flux at layer l+1
507	! using equation 11 of appendix 4.2 in paper
508
509	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
510	& -detr_star(ig,l)
511
512	endif
513	enddo
514
515	! -----------------------------------------------------------------------------------
516	!
517	! .II. CONVERGENCE LOOP
518	!
519	! We have estimated a vertical velocity profile and refined the source layer profile
520	! We now conduct iterations to compute:
521	!
522	! - the temperature inside the updraft from the estimated entrainment/source, detrainment,
523	! and upward mass flux.
524	! - the buoyancy from the new temperature inside the updraft
525	! - the vertical velocity from the new buoyancy
526	! - the entr., detr. and upward mass flux from the new buoyancy and vertical velocity
527	!
528	! This loop (tic) converges quickly. We have hardcoded 6 iterations from empirical observations.
529	! Convergence occurs in 1 or 2 iterations in most cases.
530	! -----------------------------------------------------------------------------------
531
532	! -----------------------------------------------------------------------------------
533	! -----------------------------------------------------------------------------------
534	DO tic=0,5 ! internal convergence loop
535	! -----------------------------------------------------------------------------------
536	! -----------------------------------------------------------------------------------
537
538	! Is the cell still active ?
539	activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
540
541	! If the cell is active, compute temperature inside updraft
542	do ig=1,ngrid
543	if (activetmp(ig)) then
544
545	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
546	& (alim_star(ig,l)+entr_star(ig,l))*ztv(ig,l)) &
547	& /(f_star(ig,l+1)+detr_star(ig,l))
548
549	endif
550	enddo
551
552	! Is the cell still active with respect to temperature variations ?
553	activetmp(:)=activetmp(:).and.(abs(ztla(:,l)-ztva(:,l)).gt.0.01)
554
555	! Compute new buoyancu and vertical velocity
556	do ig=1,ngrid
557	if (activetmp(ig)) then
558	ztva(ig,l) = ztla(ig,l)
559	zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
560
561	! (discretized version of equation 12 in paragraph 40 of paper)
562	if (((a1*zbuoy(ig,l)/zw2(ig,l)-b1) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then
563	zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.zdza1*zbuoy(ig,l)- &
564	& 2.zdzzw2(ig,l)b1-2.(1.-omega)zdzzw2(ig,l)ae(a1zbuoy(ig,l)/zw2(ig,l)-b1)*be)
565	else
566	zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.zdza1invzbuoy(ig,l)-2.zdzzw2(ig,l)b1inv)
567	endif
568	endif
569	enddo
570
571	! ================ RECOMPUTE ENTR, DETR, and F FROM NEW W2 ===================
572	! ND entrainment rate, see equation 16 of paper (paragraph 43)
573	! ND detrainment rate, see paragraph 44 of paper
574
575	do ig=1,ngrid
576	if (activetmp(ig)) then
577
578	zw2m=zw2(ig,l+1)
579	if(zw2m .gt. 0) then
580	if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
581	entr_star(ig,l)=f_star(ig,l)zdz &
582	& MAX(0.,ae(a1(zbuoy(ig,l)/zw2m)-b1)**be)
583	else
584	entr_star(ig,l)=0.
585	endif
586
587	if(zbuoy(ig,l) .gt. 0.) then
588	if(l .lt. lalim(ig)) then
589
590	detr_star(ig,l)=0.
591
592	else
593	detr_star(ig,l) = f_star(ig,l)zdz &
594	& ad
595
596	endif
597	else
598	detr_star(ig,l)=f_star(ig,l)zdz &
599	& MAX(ad,bd*zbuoy(ig,l)/zw2m)
600
601	endif
602	else
603	entr_star(ig,l)=0.
604	detr_star(ig,l)=0.
605	endif
606
607	! If we are still in the source layer, we define the source layer entr. rate (alim_star) as the
608	! maximum between the source entrainment rate and the estimated entrainment rate.
609
610	if (l.lt.lalim(ig)) then
611	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
612	entr_star(ig,l)=0.
613	endif
614
615	! Compute the non-dimensional upward mass flux at layer l+1
616	! using equation 11 of appendix 4.2 in paper
617
618	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
619	& -detr_star(ig,l)
620
621	endif
622	enddo
623	! -----------------------------------------------------------------------------------
624	! -----------------------------------------------------------------------------------
625	ENDDO ! of internal convergence loop
626	! -----------------------------------------------------------------------------------
627	! -----------------------------------------------------------------------------------
628
629	!---------------------------------------------------------------------------
630	! Miscellaneous computations for height
631	!---------------------------------------------------------------------------
632
633	do ig=1,ngrid
634	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
635	IF (thermverbose) THEN
636	print*,'thermcell_plume, particular case in velocity profile'
637	ENDIF
638	zw2(ig,l+1)=0.
639	endif
640
641	if (zw2(ig,l+1).lt.0.) then
642	zw2(ig,l+1)=0.
643	endif
644	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
645
646	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
647	wmaxa(ig)=wa_moy(ig,l+1)
648	endif
649	enddo
650
651	!=========================================================================
652	!=========================================================================
653	!=========================================================================
654	! END OF THE LOOP ON VERTICAL LEVELS
655	enddo
656	!=========================================================================
657	!=========================================================================
658	!=========================================================================
659
660	! Recompute the source layer total entrainment alim_star_tot
661	! as alim_star may have been modified in the above loop. Renormalization of
662	! alim_star.
663
664	do ig=1,ngrid
665	alim_star_tot(ig)=0.
666	enddo
667	do ig=1,ngrid
668	do l=1,lalim(ig)-1
669	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
670	enddo
671	enddo
672
673	do l=1,nlayer
674	do ig=1,ngrid
675	if (alim_star_tot(ig) > 1.e-10 ) then
676	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
677	endif
678	enddo
679	enddo
680
681	! ===========================================================================
682	! ================= FIN PLUME ===============================================
683	! ===========================================================================
684
685	! ===========================================================================
686	! ================= HEIGHT ==================================================
687	! ===========================================================================
688
689	! WARNING, W2 (squared velocity) IS TRANSFORMED IN ITS SQUARE ROOT HERE
690
691	!-------------------------------------------------------------------------------
692	! Computations of the thermal height zmax and maximum vertical velocity wmax
693	!-------------------------------------------------------------------------------
694
695	! Index of the thermal plume height
696	do ig=1,ngrid
697	lmax(ig)=lalim(ig)
698	enddo
699	do ig=1,ngrid
700	do l=nlayer,lalim(ig)+1,-1
701	if (zw2(ig,l).le.1.e-10) then
702	lmax(ig)=l-1
703	endif
704	enddo
705	enddo
706
707	! Particular case when the thermal reached the model top, which is not a good sign
708	do ig=1,ngrid
709	if ( zw2(ig,nlayer) > 1.e-10 ) then
710	print*,'WARNING !!!!! W2 non-zero in last layer'
711	lmax(ig)=nlayer
712	endif
713	enddo
714
715	! Maximum vertical velocity zw2
716	do ig=1,ngrid
717	wmax(ig)=0.
718	enddo
719
720	do l=1,nlayer
721	do ig=1,ngrid
722	if (l.le.lmax(ig)) then
723	if (zw2(ig,l).lt.0.)then
724	! print*,'pb2 zw2<0',zw2(ig,l)
725	zw2(ig,l)=0.
726	endif
727	zw2(ig,l)=sqrt(zw2(ig,l))
728	wmax(ig)=max(wmax(ig),zw2(ig,l))
729	else
730	zw2(ig,l)=0.
731	endif
732	enddo
733	enddo
734
735	! Height of the thermal plume, defined as the following:
736	! zmax=Integral[zw(z)dz]/Integral[w(z)*dz]
737	!
738	do ig=1,ngrid
739	zmax(ig)=0.
740	zlevinter(ig)=zlev(ig,1)
741	enddo
742
743	num(:)=0.
744	denom(:)=0.
745	do ig=1,ngrid
746	do l=1,nlayer
747	num(ig)=num(ig)+zw2(ig,l)zlev(ig,l)(zlev(ig,l+1)-zlev(ig,l))
748	denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
749	enddo
750	enddo
751	do ig=1,ngrid
752	if (denom(ig).gt.1.e-10) then
753	zmax(ig)=2.*num(ig)/denom(ig)
754	endif
755	enddo
756
757	! ===========================================================================
758	! ================= FIN HEIGHT ==============================================
759	! ===========================================================================
760
761	zlmax=MAXVAL(lmax(:))+2
762	if (zlmax .ge. nlayer) then
763	print*,'thermals have reached last layer of the model'
764	print*,'this is not good !'
765	endif
766
767	! alim_star_clos is the source profile used for closure. It consists of the
768	! modified source profile in the source layers, and the entrainment profile
769	! above it.
770
771	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
772
773	! ===========================================================================
774	! ============= CLOSURE =====================================================
775	! ===========================================================================
776
777	!-------------------------------------------------------------------------------
778	! Closure, determination of the upward mass flux
779	!-------------------------------------------------------------------------------
780	! Init.
781
782	alim_star2(:)=0.
783	alim_star_tot_clos(:)=0.
784	f(:)=0.
785
786	! llmax is the index of the heighest thermal in the simulation domain
787	llmax=1
788	do ig=1,ngrid
789	if (lalim(ig)>llmax) llmax=lalim(ig)
790	enddo
791
792	! Integral of a*2/(rho Delta z), see equation 13 of appendix 4.2 in paper
793
794	do k=1,llmax-1
795	do ig=1,ngrid
796	if (k<lalim(ig)) then
797	alim_star2(ig)=alim_star2(ig)+alim_star_clos(ig,k)*alim_star_clos(ig,k) &
798	& /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
799	alim_star_tot_clos(ig)=alim_star_tot_clos(ig)+alim_star_clos(ig,k)
800	endif
801	enddo
802	enddo
803
804	! Closure mass flux, equation 13 of appendix 4.2 in paper
805
806	do ig=1,ngrid
807	if (alim_star2(ig)>1.e-10) then
808	f(ig)=wmax(ig)*alim_star_tot_clos(ig)/ &
809	& (max(500.,zmax(ig))r_aspect_thermalsalim_star2(ig))
810
811	endif
812	enddo
813
814	! ===========================================================================
815	! ============= FIN CLOSURE =================================================
816	! ===========================================================================
817
818
819	! ===========================================================================
820	! ============= FLUX2 =======================================================
821	! ===========================================================================
822
823	!-------------------------------------------------------------------------------
824	! With the closure mass flux, we can compute the entrainment, detrainment and
825	! upward mass flux from the non-dimensional ones.
826	!-------------------------------------------------------------------------------
827
828	fomass_max=0.8 !maximum mass fraction of a cell that can go upward in an
829	! updraft
830	alphamax=0.5 !maximum updraft coverage in a cell
831
832
833	! these variables allow to follow corrections made to the mass flux when thermverbose=.true.
834	ncorecfm1=0
835	ncorecfm2=0
836	ncorecfm3=0
837	ncorecfm4=0
838	ncorecfm5=0
839	ncorecfm6=0
840	ncorecfm7=0
841	ncorecfm8=0
842	ncorecalpha=0
843
844	!-------------------------------------------------------------------------
845	! Multiply by the closure mass flux
846	!-------------------------------------------------------------------------
847
848	do l=1,zlmax
849	entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
850	detr(:,l)=f(:)*detr_star(:,l)
851	enddo
852
853	! Reconstruct the updraft mass flux everywhere
854
855	do l=1,zlmax
856	do ig=1,ngrid
857	if (l.lt.lmax(ig)) then
858	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
859	elseif(l.eq.lmax(ig)) then
860	fm(ig,l+1)=0.
861	detr(ig,l)=fm(ig,l)+entr(ig,l)
862	else
863	fm(ig,l+1)=0.
864	endif
865	enddo
866	enddo
867
868
869	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
870	!
871	! Now we will reconstruct once again the upward
872	! mass flux, but we will apply corrections
873	! in some cases. We can compare to the
874	! previously computed mass flux (above)
875	!
876	! This verification is done level by level
877	!
878	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
879
880	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
881
882	do l=1,zlmax !loop on the levels
883	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
884
885	! Upward mass flux at level l+1
886
887	do ig=1,ngrid
888	if (l.lt.lmax(ig)) then
889	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
890	elseif(l.eq.lmax(ig)) then
891	fm(ig,l+1)=0.
892	detr(ig,l)=fm(ig,l)+entr(ig,l)
893	else
894	fm(ig,l+1)=0.
895	endif
896	enddo
897
898
899	!-------------------------------------------------------------------------
900	! Upward mass flux should be positive
901	!-------------------------------------------------------------------------
902
903	do ig=1,ngrid
904
905	if (fm(ig,l+1).lt.0.) then
906	if((l+1) .eq. lmax(ig)) then
907	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
908	fm(ig,l+1)=0.
909	entr(ig,l+1)=0.
910	ncorecfm2=ncorecfm2+1
911	else
912	IF (thermverbose) THEN
913	print*,'fm(l+1)<0 : ig, l+1,lmax :',ig,l+1,lmax(ig),fm(ig,l+1)
914	ENDIF
915	ncorecfm1=ncorecfm1+1
916	fm(ig,l+1)=fm(ig,l)
917	detr(ig,l)=entr(ig,l)
918	endif
919	endif
920
921	enddo
922
923	!-------------------------------------------------------------------------
924	! Detrainment should be lower than upward mass flux
925	!-------------------------------------------------------------------------
926
927	do ig=1,ngrid
928	if (detr(ig,l).gt.fm(ig,l)) then
929	ncorecfm6=ncorecfm6+1
930	detr(ig,l)=fm(ig,l)
931	entr(ig,l)=fm(ig,l+1)
932
933	! When detrainment is stronger than upward mass flux, and we are above the
934	! thermal last level, the plume is stopped
935
936	if(l.gt.lmax(ig)) then
937	detr(ig,l)=0.
938	fm(ig,l+1)=0.
939	entr(ig,l)=0.
940	endif
941
942	endif
943
944	enddo
945
946	!-------------------------------------------------------------------------
947	! Check again for mass flux positivity
948	!-------------------------------------------------------------------------
949
950	do ig=1,ngrid
951	if (fm(ig,l+1).lt.0.) then
952	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
953	fm(ig,l+1)=0.
954	ncorecfm2=ncorecfm2+1
955	endif
956	enddo
957
958	!-----------------------------------------------------------------------
959	! Fractional coverage should be less than 1
960	!-----------------------------------------------------------------------
961
962	do ig=1,ngrid
963	if (zw2(ig,l+1).gt.1.e-10) then
964	zfm=rhobarz(ig,l+1)zw2(ig,l+1)alphamax
965	if ( fm(ig,l+1) .gt. zfm) then
966	f_old=fm(ig,l+1)
967	fm(ig,l+1)=zfm
968	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
969	ncorecalpha=ncorecalpha+1
970	endif
971	endif
972
973	enddo
974
975	enddo ! on vertical levels
976
977	!-----------------------------------------------------------------------
978	!
979	! We limit the total mass going from one level to the next, compared to the
980	! initial total mass fo the cell
981	!
982	!-----------------------------------------------------------------------
983
984	do l=1,zlmax
985	do ig=1,ngrid
986	eee0=entr(ig,l)
987	ddd0=detr(ig,l)
988	eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
989	ddd=detr(ig,l)-eee
990	if (eee.gt.0.) then
991	ncorecfm3=ncorecfm3+1
992	entr(ig,l)=entr(ig,l)-eee
993	if ( ddd.gt.0.) then
994	! The entrainment is too strong but we can compensate the excess by a detrainment decrease
995	detr(ig,l)=ddd
996	else
997	! The entrainment is too strong and we compensate the excess by a stronger entrainment
998	! in the layer above
999	if(l.eq.lmax(ig)) then
1000	detr(ig,l)=fm(ig,l)+entr(ig,l)
1001	else
1002	entr(ig,l+1)=entr(ig,l+1)-ddd
1003	detr(ig,l)=0.
1004	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
1005	detr(ig,l)=0.
1006	endif
1007	endif
1008	endif
1009	enddo
1010	enddo
1011
1012	! Check again that everything cancels at zmax
1013	do ig=1,ngrid
1014	fm(ig,lmax(ig)+1)=0.
1015	entr(ig,lmax(ig))=0.
1016	detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
1017	enddo
1018
1019	!-----------------------------------------------------------------------
1020	! Summary of the number of modifications that were necessary (if thermverbose=.true.
1021	! and only if there were a lot of them)
1022	!-----------------------------------------------------------------------
1023
1024	!IM 090508 beg
1025	IF (thermverbose) THEN
1026	if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngrid/4. ) then
1027	print*,'thermcell warning : large number of corrections'
1028	print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
1029	& ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
1030	& ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
1031	& ncorecfm6,'x fm6', &
1032	& ncorecfm7,'x fm7', &
1033	& ncorecfm8,'x fm8', &
1034	& ncorecalpha,'x alpha'
1035	endif
1036	ENDIF
1037
1038	! ===========================================================================
1039	! ============= FIN FLUX2 ===================================================
1040	! ===========================================================================
1041
1042
1043	! ===========================================================================
1044	! ============= TRANSPORT ===================================================
1045	! ===========================================================================
1046
1047	!------------------------------------------------------------------
1048	! vertical transport computation
1049	!------------------------------------------------------------------
1050
1051	! ------------------------------------------------------------------
1052	! IN THE UPDRAFT
1053	! ------------------------------------------------------------------
1054
1055	zdthladj(:,:)=0.
1056	! Based on equation 14 in appendix 4.2
1057
1058	do ig=1,ngrid
1059	if(lmax(ig) .gt. 1) then
1060	do k=1,lmax(ig)
1061	zdthladj(ig,k)=(1./masse(ig,k))(fm(ig,k+1)ztv(ig,k+1)- &
1062	& fm(ig,k)ztv(ig,k)+fm(ig,k)ztva(ig,k)-fm(ig,k+1)*ztva(ig,k+1))
1063	if (ztv(ig,k) + ptimestep*zdthladj(ig,k) .le. 0.) then
1064	IF (thermverbose) THEN
1065	print,'Teta<0 in thermcell_dTeta up: qenv .. dq : ', ztv(ig,k),ptimestepzdthladj(ig,k)
1066	ENDIF
1067	if(ztv(ig,k) .gt. 0.) then
1068	zdthladj(ig,k)=0.
1069	endif
1070	endif
1071	enddo
1072	endif
1073	enddo
1074
1075	! ------------------------------------------------------------------
1076	! DOWNDRAFT PARAMETERIZATION
1077	! ------------------------------------------------------------------
1078
1079	ztvd(:,:)=ztv(:,:)
1080	fm_down(:,:)=0.
1081	do ig=1,ngrid
1082	if (lmax(ig) .gt. 1) then
1083	do l=1,lmax(ig)
1084	if(zlay(ig,l) .le. zmax(ig)) then
1085
1086	! see equation 18 of paragraph 48 in paper
1087	fm_down(ig,l) =fm(ig,l)* &
1088	& max(fdfu,-4*max(0.,(zlay(ig,l)/zmax(ig)))-0.6)
1089	endif
1090
1091	if(zlay(ig,l) .le. zmax(ig)) then
1092	! see equation 19 of paragraph 49 in paper
1093	ztvd(ig,l)=min(ztv(ig,l),ztv(ig,l)*((zlay(ig,l)/zmax(ig))/400. + 0.997832))
1094	else
1095	ztvd(ig,l)=ztv(ig,l)
1096	endif
1097
1098	enddo
1099	endif
1100	enddo
1101
1102	! ------------------------------------------------------------------
1103	! TRANSPORT IN DOWNDRAFT
1104	! ------------------------------------------------------------------
1105
1106	zdthladj_down(:,:)=0.
1107
1108	do ig=1,ngrid
1109	if(lmax(ig) .gt. 1) then
1110	! No downdraft in the very-near surface layer, we begin at k=3
1111	! Based on equation 14 in appendix 4.2
1112
1113	do k=3,lmax(ig)
1114	zdthladj_down(ig,k)=(1./masse(ig,k))(fm_down(ig,k+1)ztv(ig,k+1)- &
1115	& fm_down(ig,k)ztv(ig,k)+fm_down(ig,k)ztvd(ig,k)-fm_down(ig,k+1)*ztvd(ig,k+1))
1116	if (ztv(ig,k) + ptimestep*zdthladj_down(ig,k) .le. 0.) then
1117	IF (thermverbose) THEN
1118	print,'q<0 in thermcell_dTeta down: qenv .. dq : ', ztv(ig,k),ptimestepzdthladj_down(ig,k)
1119	ENDIF
1120	if(ztv(ig,k) .gt. 0.) then
1121	zdthladj(ig,k)=0.
1122	endif
1123	endif
1124	enddo
1125	endif
1126	enddo
1127
1128	!------------------------------------------------------------------
1129	! Final fraction coverage with the clean upward mass flux, computed at interfaces
1130	!------------------------------------------------------------------
1131	fraca(:,:)=0.
1132	do l=2,zlmax
1133	do ig=1,ngrid
1134	if (zw2(ig,l).gt.1.e-10) then
1135	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
1136	else
1137	fraca(ig,l)=0.
1138	endif
1139	enddo
1140	enddo
1141
1142	!------------------------------------------------------------------
1143	! Transport of C02 Tracer
1144	!------------------------------------------------------------------
1145
1146	! We only transport co2 tracer because it is coupled to the scheme through theta_m
1147	! The rest is transported outside the sub-timestep loop
1148
1149	ratiom(:,:)=1.
1150
1151	if (igcm_co2.ne.0) then
1152	detrmod(:,:)=0.
1153	do k=1,zlmax
1154	do ig=1,ngrid
1155	detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
1156	& +entr(ig,k)
1157	if (detrmod(ig,k).lt.0.) then
1158	entr(ig,k)=entr(ig,k)-detrmod(ig,k)
1159	detrmod(ig,k)=0.
1160	endif
1161	enddo
1162	enddo
1163
1164	call thermcell_dqup(ngrid,nlayer,ptimestep &
1165	& ,fm,entr,detrmod, &
1166	& masse,pq(:,:,igcm_co2),pdqadj(:,:,igcm_co2),ndt,zlmax)
1167
1168	! Compute the ratio between theta and theta_m
1169
1170	do l=1,zlmax
1171	do ig=1,ngrid
1172	ratiom(ig,l)=1./(A(pq(ig,l,igcm_co2)+pdqadj(ig,l,igcm_co2)ptimestep)+B)
1173	enddo
1174	enddo
1175
1176	endif
1177
1178	!------------------------------------------------------------------
1179	! incrementation dt
1180	!------------------------------------------------------------------
1181
1182	pdtadj(:,:)=0.
1183	do l=1,zlmax
1184	do ig=1,ngrid
1185	pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))zpopsk(ig,l)ratiom(ig,l)
1186	enddo
1187	enddo
1188
1189	! ===========================================================================
1190	! ============= FIN TRANSPORT ===============================================
1191	! ===========================================================================
1192
1193
1194	!------------------------------------------------------------------
1195	! Diagnostics for outputs
1196	!------------------------------------------------------------------
1197	! We compute interface values for teta env and th. The last interface
1198	! value does not matter, as the mass flux is 0 there.
1199
1200
1201	do l=1,nlayer-1
1202	do ig=1,ngrid
1203	teta_th_int(ig,l)=0.5(ztva(ig,l+1)+ztva(ig,l))ratiom(ig,l)
1204	teta_down_int(ig,l) = 0.5(ztvd(ig,l+1)+ztvd(ig,l))ratiom(ig,l)
1205	teta_env_int(ig,l)=0.5(ztv(ig,l+1)+ztv(ig,l))ratiom(ig,l)
1206	enddo
1207	enddo
1208	do ig=1,ngrid
1209	teta_th_int(ig,nlayer)=teta_th_int(ig,nlayer-1)
1210	teta_env_int(ig,nlayer)=teta_env_int(ig,nlayer-1)
1211	teta_down_int(ig,nlayer)=teta_down_int(ig,nlayer-1)
1212	enddo
1213	heatFlux(:,:)=0.
1214	buoyancyOut(:,:)=0.
1215	buoyancyEst(:,:)=0.
1216	heatFlux_down(:,:)=0.
1217	do l=1,zlmax
1218	do ig=1,ngrid
1219	heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
1220	buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1221	buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1222	heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l)
1223	enddo
1224	enddo
1225
1226	return
1227	end

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

Download in other formats:

Original Format