Context Navigation

thermcell_main_mars.F90 @ 3026

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

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

Download in other formats:

Original Format