Context Navigation

thermcell_main_mars.F90 @ 2616

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

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

Download in other formats:

Original Format