Context Navigation

thermcell_main_mars.F90 @ 1032

Last change on this file since 1032 was 1032, checked in by aslmd, 11 years ago
LMDZ.MARS cleaned and commented version of the thermal plume model with automatic arrays.
File size: 42.6 KB

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

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

Download in other formats:

Original Format