Context Navigation

thermcell_main_mars.F90 @ 1130

Last change on this file since 1130 was 1130, checked in by emillour, 11 years ago

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
dyn3d/control.h is now module control_mod.F90
rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

File size: 43.4 KB

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

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

Download in other formats:

Original Format