Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_camzm.F @ 6

Last change on this file since 6 was 1, checked in by lfita, 11 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 128.5 KB

Line
1	#define WRF_PORT
2
3	module module_cu_camzm
4
5	!Based on zm_conv.F90 from CAM
6
7	!---------------------------------------------------------------------------------
8	! Purpose:
9	!
10	! Interface from Zhang-McFarlane convection scheme, includes evaporation of convective
11	! precip from the ZM scheme
12	!
13	! Apr 2006: RBN: Code added to perform a dilute ascent for closure of the CM mass flux
14	! based on an entraining plume a la Raymond and Blythe (1992)
15	!
16	! Author: Byron Boville, from code in tphysbc
17	! Ported from CAM to WRF by William.Gustafson@pnl.gov, Nov. 2009
18	! updated to CESM1_0_1, Nov. 2010
19	!
20	!---------------------------------------------------------------------------------
21
22	use shr_kind_mod, only: r8 => shr_kind_r8
23	use cldwat, only: cldwat_fice
24	use physconst, only: cpair, epsilo, gravit, latice, latvap, tmelt, rair, &
25	cpwv, cpliq, rh2o
26	#ifndef WRF_PORT
27	use spmd_utils, only: masterproc
28	use ppgrid, only: pcols, pver, pverp
29
30	use abortutils, only: endrun
31	use cam_logfile, only: iulog
32
33	#else
34	use module_cam_support, only: masterproc, &
35	pcols, pver, pverp, &
36	endrun, &
37	iulog
38	use module_wrf_error
39	#endif
40
41	implicit none
42
43	save
44	private ! Make default type private to the module
45	!
46	! PUBLIC: interfaces
47	!
48	public zmconv_readnl ! read zmconv_nl namelist
49	public zm_convi ! ZM schemea
50	public zm_convr ! ZM schemea
51	public zm_conv_evap ! evaporation of precip from ZM schemea
52	public convtran ! convective transport
53	public momtran ! convective momentum transport
54
55	!
56	! Private data
57	!
58	real(r8), parameter :: unset_r8 = huge(1.0_r8)
59	real(r8) :: zmconv_c0_lnd = unset_r8
60	real(r8) :: zmconv_c0_ocn = unset_r8
61	real(r8) :: zmconv_ke = unset_r8
62	real(r8) rl ! wg latent heat of vaporization.
63	real(r8) cpres ! specific heat at constant pressure in j/kg-degk.
64	real(r8), parameter :: capelmt = 70._r8 ! threshold value for cape for deep convection.
65	real(r8) :: ke ! Tunable evaporation efficiency set from namelist input zmconv_ke
66	real(r8) :: c0_lnd ! set from namelist input zmconv_c0_lnd
67	real(r8) :: c0_ocn ! set from namelist input zmconv_c0_ocn
68	real(r8) tau ! convective time scale
69	real(r8),parameter :: a = 21.656_r8
70	real(r8),parameter :: b = 5418._r8
71	real(r8),parameter :: c1 = 6.112_r8
72	real(r8),parameter :: c2 = 17.67_r8
73	real(r8),parameter :: c3 = 243.5_r8
74	real(r8) :: tfreez
75	real(r8) :: eps1
76
77
78	logical :: no_deep_pbl ! default = .false.
79	! no_deep_pbl = .true. eliminates deep convection entirely within PBL
80
81
82	!moved from moistconvection.F90
83	real(r8) :: rgrav ! reciprocal of grav
84	real(r8) :: rgas ! gas constant for dry air
85	real(r8) :: grav ! = gravit
86	real(r8) :: cp ! = cpres = cpair
87
88	integer limcnv ! top interface level limit for convection
89	real(r8),parameter :: tiedke_add = 0.5_r8
90	contains
91
92	subroutine zmconv_readnl(nlfile)
93	#ifndef WRF_PORT
94	use namelist_utils, only: find_group_name
95	use units, only: getunit, freeunit
96	use mpishorthand
97	#endif
98
99	character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
100
101	#ifndef WRF_PORT
102	! Local variables
103	integer :: unitn, ierr
104	character(len=*), parameter :: subname = 'zmconv_readnl'
105
106	namelist /zmconv_nl/ zmconv_c0_lnd, zmconv_c0_ocn, zmconv_ke
107	!-----------------------------------------------------------------------------
108
109	if (masterproc) then
110	unitn = getunit()
111	open( unitn, file=trim(nlfile), status='old' )
112	call find_group_name(unitn, 'zmconv_nl', status=ierr)
113	if (ierr == 0) then
114	read(unitn, zmconv_nl, iostat=ierr)
115	if (ierr /= 0) then
116	call endrun(subname // ':: ERROR reading namelist')
117	end if
118	end if
119	close(unitn)
120	call freeunit(unitn)
121
122	! set local variables
123	c0_lnd = zmconv_c0_lnd
124	c0_ocn = zmconv_c0_ocn
125	ke = zmconv_ke
126
127	end if
128
129	#ifdef SPMD
130	! Broadcast namelist variables
131	call mpibcast(c0_lnd, 1, mpir8, 0, mpicom)
132	call mpibcast(c0_ocn, 1, mpir8, 0, mpicom)
133	call mpibcast(ke, 1, mpir8, 0, mpicom)
134	#endif
135
136	#else
137	! WRF_PORT currently uses hard-wired values for the namelist input. The
138	! values could easily be setup to come from the Registry in the future.
139	! The hard-wired values are the defaults for the fv core. They should be
140	! verified by somebody knowledgable on the matter.
141	c0_lnd = 0.0035d0
142	c0_ocn = 0.0035d0
143	ke = 1.0e-6
144	#endif
145
146	end subroutine zmconv_readnl
147
148
149	subroutine zm_convi(limcnv_in, no_deep_pbl_in)
150
151	#ifndef WRF_PORT
152	use dycore, only: dycore_is, get_resolution
153	#endif
154
155	integer, intent(in) :: limcnv_in ! top interface level limit for convection
156	logical, intent(in), optional :: no_deep_pbl_in ! no_deep_pbl = .true. eliminates ZM convection entirely within PBL
157
158	! local variables
159	character(len=32) :: hgrid ! horizontal grid specifier
160
161	! Initialization of ZM constants
162	limcnv = limcnv_in
163	tfreez = tmelt
164	eps1 = epsilo
165	rl = latvap
166	cpres = cpair
167	rgrav = 1.0_r8/gravit
168	rgas = rair
169	grav = gravit
170	cp = cpres
171
172	if ( present(no_deep_pbl_in) ) then
173	no_deep_pbl = no_deep_pbl_in
174	else
175	no_deep_pbl = .false.
176	endif
177
178	! tau=4800. were used in canadian climate center. however, in echam3 t42,
179	! convection is too weak, thus adjusted to 2400.
180
181	#ifndef WRF_PORT
182	hgrid = get_resolution()
183	#endif
184	tau = 3600._r8
185
186	if ( masterproc ) then
187	write(iulog,*) 'tuning parameters zm_convi: tau',tau
188	#ifdef WRF_PORT
189	call wrf_debug(100,iulog)
190	#endif
191	write(iulog,*) 'tuning parameters zm_convi: c0_lnd',c0_lnd, ', c0_ocn', c0_ocn
192	#ifdef WRF_PORT
193	call wrf_debug(100,iulog)
194	#endif
195	write(iulog,*) 'tuning parameters zm_convi: ke',ke
196	#ifdef WRF_PORT
197	call wrf_debug(100,iulog)
198	#endif
199	write(iulog,*) 'tuning parameters zm_convi: no_deep_pbl',no_deep_pbl
200	#ifdef WRF_PORT
201	call wrf_debug(100,iulog)
202	endif
203	#endif
204
205	if (masterproc) then
206	write(iulog,)'* ZM: DILUTE Buoyancy Calculation **'
207	#ifdef WRF_PORT
208	call wrf_debug(100,iulog)
209	#endif
210	end if
211
212	end subroutine zm_convi
213
214
215
216	subroutine zm_convr(lchnk ,ncol , &
217	t ,qh ,prec ,jctop ,jcbot , &
218	pblh ,zm ,geos ,zi ,qtnd , &
219	heat ,pap ,paph ,dpp , &
220	delt ,mcon ,cme ,cape , &
221	tpert ,dlf ,pflx ,zdu ,rprd , &
222	mu ,md ,du ,eu ,ed , &
223	dp ,dsubcld ,jt ,maxg ,ideep , &
224	lengath ,ql ,rliq ,landfrac)
225	!-----------------------------------------------------------------------
226	!
227	! Purpose:
228	! Main driver for zhang-mcfarlane convection scheme
229	!
230	! Method:
231	! performs deep convective adjustment based on mass-flux closure
232	! algorithm.
233	!
234	! Author:guang jun zhang, m.lazare, n.mcfarlane. CAM Contact: P. Rasch
235	!
236	! This is contributed code not fully standardized by the CAM core group.
237	! All variables have been typed, where most are identified in comments
238	! The current procedure will be reimplemented in a subsequent version
239	! of the CAM where it will include a more straightforward formulation
240	! and will make use of the standard CAM nomenclature
241	!
242	!-----------------------------------------------------------------------
243	#ifdef WRF_PORT
244	use module_cam_support, only: pcnst
245	#else
246	use constituents, only: pcnst
247	#endif
248
249	!
250	! ********************** index of variables ********************
251	!
252	! wg * alpha array of vertical differencing used (=1. for upstream).
253	! w * cape convective available potential energy.
254	! wg * capeg gathered convective available potential energy.
255	! c * capelmt threshold value for cape for deep convection.
256	! ic * cpres specific heat at constant pressure in j/kg-degk.
257	! i * dpp
258	! ic * delt length of model time-step in seconds.
259	! wg * dp layer thickness in mbs (between upper/lower interface).
260	! wg * dqdt mixing ratio tendency at gathered points.
261	! wg * dsdt dry static energy ("temp") tendency at gathered points.
262	! wg * dudt u-wind tendency at gathered points.
263	! wg * dvdt v-wind tendency at gathered points.
264	! wg * dsubcld layer thickness in mbs between lcl and maxi.
265	! ic * grav acceleration due to gravity in m/sec2.
266	! wg * du detrainment in updraft. specified in mid-layer
267	! wg * ed entrainment in downdraft.
268	! wg * eu entrainment in updraft.
269	! wg * hmn moist static energy.
270	! wg * hsat saturated moist static energy.
271	! w * ideep holds position of gathered points vs longitude index.
272	! ic * pver number of model levels.
273	! wg * j0 detrainment initiation level index.
274	! wg * jd downdraft initiation level index.
275	! ic * jlatpr gaussian latitude index for printing grids (if needed).
276	! wg * jt top level index of deep cumulus convection.
277	! w * lcl base level index of deep cumulus convection.
278	! wg * lclg gathered values of lcl.
279	! w * lel index of highest theoretical convective plume.
280	! wg * lelg gathered values of lel.
281	! w * lon index of onset level for deep convection.
282	! w * maxi index of level with largest moist static energy.
283	! wg * maxg gathered values of maxi.
284	! wg * mb cloud base mass flux.
285	! wg * mc net upward (scaled by mb) cloud mass flux.
286	! wg * md downward cloud mass flux (positive up).
287	! wg * mu upward cloud mass flux (positive up). specified
288	! at interface
289	! ic * msg number of missing moisture levels at the top of model.
290	! w * p grid slice of ambient mid-layer pressure in mbs.
291	! i * pblt row of pbl top indices.
292	! w * pcpdh scaled surface pressure.
293	! w * pf grid slice of ambient interface pressure in mbs.
294	! wg * pg grid slice of gathered values of p.
295	! w * q grid slice of mixing ratio.
296	! wg * qd grid slice of mixing ratio in downdraft.
297	! wg * qg grid slice of gathered values of q.
298	! i/o * qh grid slice of specific humidity.
299	! w * qh0 grid slice of initial specific humidity.
300	! wg * qhat grid slice of upper interface mixing ratio.
301	! wg * ql grid slice of cloud liquid water.
302	! wg * qs grid slice of saturation mixing ratio.
303	! w * qstp grid slice of parcel temp. saturation mixing ratio.
304	! wg * qstpg grid slice of gathered values of qstp.
305	! wg * qu grid slice of mixing ratio in updraft.
306	! ic * rgas dry air gas constant.
307	! wg * rl latent heat of vaporization.
308	! w * s grid slice of scaled dry static energy (t+gz/cp).
309	! wg * sd grid slice of dry static energy in downdraft.
310	! wg * sg grid slice of gathered values of s.
311	! wg * shat grid slice of upper interface dry static energy.
312	! wg * su grid slice of dry static energy in updraft.
313	! i/o * t
314	! o * jctop row of top-of-deep-convection indices passed out.
315	! O * jcbot row of base of cloud indices passed out.
316	! wg * tg grid slice of gathered values of t.
317	! w * tl row of parcel temperature at lcl.
318	! wg * tlg grid slice of gathered values of tl.
319	! w * tp grid slice of parcel temperatures.
320	! wg * tpg grid slice of gathered values of tp.
321	! i/o * u grid slice of u-wind (real).
322	! wg * ug grid slice of gathered values of u.
323	! i/o * utg grid slice of u-wind tendency (real).
324	! i/o * v grid slice of v-wind (real).
325	! w * va work array re-used by called subroutines.
326	! wg * vg grid slice of gathered values of v.
327	! i/o * vtg grid slice of v-wind tendency (real).
328	! i * w grid slice of diagnosed large-scale vertical velocity.
329	! w * z grid slice of ambient mid-layer height in metres.
330	! w * zf grid slice of ambient interface height in metres.
331	! wg * zfg grid slice of gathered values of zf.
332	! wg * zg grid slice of gathered values of z.
333	!
334	!-----------------------------------------------------------------------
335	!
336	! multi-level i/o fields:
337	! i => input arrays.
338	! i/o => input/output arrays.
339	! w => work arrays.
340	! wg => work arrays operating only on gathered points.
341	! ic => input data constants.
342	! c => data constants pertaining to subroutine itself.
343	!
344	! input arguments
345	!
346	integer, intent(in) :: lchnk ! chunk identifier
347	integer, intent(in) :: ncol ! number of atmospheric columns
348
349	real(r8), intent(in) :: t(pcols,pver) ! grid slice of temperature at mid-layer.
350	real(r8), intent(in) :: qh(pcols,pver,pcnst) ! grid slice of specific humidity.
351	real(r8), intent(in) :: pap(pcols,pver)
352	real(r8), intent(in) :: paph(pcols,pver+1)
353	real(r8), intent(in) :: dpp(pcols,pver) ! local sigma half-level thickness (i.e. dshj).
354	real(r8), intent(in) :: zm(pcols,pver)
355	real(r8), intent(in) :: geos(pcols)
356	real(r8), intent(in) :: zi(pcols,pver+1)
357	real(r8), intent(in) :: pblh(pcols)
358	real(r8), intent(in) :: tpert(pcols)
359	real(r8), intent(in) :: landfrac(pcols) ! RBN Landfrac
360	!
361	! output arguments
362	!
363	real(r8), intent(out) :: qtnd(pcols,pver) ! specific humidity tendency (kg/kg/s)
364	real(r8), intent(out) :: heat(pcols,pver) ! heating rate (dry static energy tendency, W/kg)
365	real(r8), intent(out) :: mcon(pcols,pverp)
366	real(r8), intent(out) :: dlf(pcols,pver) ! scattrd version of the detraining cld h2o tend
367	real(r8), intent(out) :: pflx(pcols,pverp) ! scattered precip flux at each level
368	real(r8), intent(out) :: cme(pcols,pver)
369	real(r8), intent(out) :: cape(pcols) ! w convective available potential energy.
370	real(r8), intent(out) :: zdu(pcols,pver)
371	real(r8), intent(out) :: rprd(pcols,pver) ! rain production rate
372	! move these vars from local storage to output so that convective
373	! transports can be done in outside of conv_cam.
374	real(r8), intent(out) :: mu(pcols,pver)
375	real(r8), intent(out) :: eu(pcols,pver)
376	real(r8), intent(out) :: du(pcols,pver)
377	real(r8), intent(out) :: md(pcols,pver)
378	real(r8), intent(out) :: ed(pcols,pver)
379	real(r8), intent(out) :: dp(pcols,pver) ! wg layer thickness in mbs (between upper/lower interface).
380	real(r8), intent(out) :: dsubcld(pcols) ! wg layer thickness in mbs between lcl and maxi.
381	real(r8), intent(out) :: jctop(pcols) ! o row of top-of-deep-convection indices passed out.
382	real(r8), intent(out) :: jcbot(pcols) ! o row of base of cloud indices passed out.
383	real(r8), intent(out) :: prec(pcols)
384	real(r8), intent(out) :: rliq(pcols) ! reserved liquid (not yet in cldliq) for energy integrals
385
386	real(r8) zs(pcols)
387	real(r8) dlg(pcols,pver) ! gathrd version of the detraining cld h2o tend
388	real(r8) pflxg(pcols,pverp) ! gather precip flux at each level
389	real(r8) cug(pcols,pver) ! gathered condensation rate
390	real(r8) evpg(pcols,pver) ! gathered evap rate of rain in downdraft
391	real(r8) mumax(pcols)
392	integer jt(pcols) ! wg top level index of deep cumulus convection.
393	integer maxg(pcols) ! wg gathered values of maxi.
394	integer ideep(pcols) ! w holds position of gathered points vs longitude index.
395	integer lengath
396	! diagnostic field used by chem/wetdep codes
397	real(r8) ql(pcols,pver) ! wg grid slice of cloud liquid water.
398	!
399	real(r8) pblt(pcols) ! i row of pbl top indices.
400
401	!
402	!-----------------------------------------------------------------------
403	!
404	! general work fields (local variables):
405	!
406	real(r8) q(pcols,pver) ! w grid slice of mixing ratio.
407	real(r8) p(pcols,pver) ! w grid slice of ambient mid-layer pressure in mbs.
408	real(r8) z(pcols,pver) ! w grid slice of ambient mid-layer height in metres.
409	real(r8) s(pcols,pver) ! w grid slice of scaled dry static energy (t+gz/cp).
410	real(r8) tp(pcols,pver) ! w grid slice of parcel temperatures.
411	real(r8) zf(pcols,pver+1) ! w grid slice of ambient interface height in metres.
412	real(r8) pf(pcols,pver+1) ! w grid slice of ambient interface pressure in mbs.
413	real(r8) qstp(pcols,pver) ! w grid slice of parcel temp. saturation mixing ratio.
414
415	real(r8) tl(pcols) ! w row of parcel temperature at lcl.
416
417	integer lcl(pcols) ! w base level index of deep cumulus convection.
418	integer lel(pcols) ! w index of highest theoretical convective plume.
419	integer lon(pcols) ! w index of onset level for deep convection.
420	integer maxi(pcols) ! w index of level with largest moist static energy.
421	integer index(pcols)
422	real(r8) precip
423	!
424	! gathered work fields:
425	!
426	real(r8) qg(pcols,pver) ! wg grid slice of gathered values of q.
427	real(r8) tg(pcols,pver) ! w grid slice of temperature at interface.
428	real(r8) pg(pcols,pver) ! wg grid slice of gathered values of p.
429	real(r8) zg(pcols,pver) ! wg grid slice of gathered values of z.
430	real(r8) sg(pcols,pver) ! wg grid slice of gathered values of s.
431	real(r8) tpg(pcols,pver) ! wg grid slice of gathered values of tp.
432	real(r8) zfg(pcols,pver+1) ! wg grid slice of gathered values of zf.
433	real(r8) qstpg(pcols,pver) ! wg grid slice of gathered values of qstp.
434	real(r8) ug(pcols,pver) ! wg grid slice of gathered values of u.
435	real(r8) vg(pcols,pver) ! wg grid slice of gathered values of v.
436	real(r8) cmeg(pcols,pver)
437
438	real(r8) rprdg(pcols,pver) ! wg gathered rain production rate
439	real(r8) capeg(pcols) ! wg gathered convective available potential energy.
440	real(r8) tlg(pcols) ! wg grid slice of gathered values of tl.
441	real(r8) landfracg(pcols) ! wg grid slice of landfrac
442
443	integer lclg(pcols) ! wg gathered values of lcl.
444	integer lelg(pcols)
445	!
446	! work fields arising from gathered calculations.
447	!
448	real(r8) dqdt(pcols,pver) ! wg mixing ratio tendency at gathered points.
449	real(r8) dsdt(pcols,pver) ! wg dry static energy ("temp") tendency at gathered points.
450	! real(r8) alpha(pcols,pver) ! array of vertical differencing used (=1. for upstream).
451	real(r8) sd(pcols,pver) ! wg grid slice of dry static energy in downdraft.
452	real(r8) qd(pcols,pver) ! wg grid slice of mixing ratio in downdraft.
453	real(r8) mc(pcols,pver) ! wg net upward (scaled by mb) cloud mass flux.
454	real(r8) qhat(pcols,pver) ! wg grid slice of upper interface mixing ratio.
455	real(r8) qu(pcols,pver) ! wg grid slice of mixing ratio in updraft.
456	real(r8) su(pcols,pver) ! wg grid slice of dry static energy in updraft.
457	real(r8) qs(pcols,pver) ! wg grid slice of saturation mixing ratio.
458	real(r8) shat(pcols,pver) ! wg grid slice of upper interface dry static energy.
459	real(r8) hmn(pcols,pver) ! wg moist static energy.
460	real(r8) hsat(pcols,pver) ! wg saturated moist static energy.
461	real(r8) qlg(pcols,pver)
462	real(r8) dudt(pcols,pver) ! wg u-wind tendency at gathered points.
463	real(r8) dvdt(pcols,pver) ! wg v-wind tendency at gathered points.
464	! real(r8) ud(pcols,pver)
465	! real(r8) vd(pcols,pver)
466
467	real(r8) mb(pcols) ! wg cloud base mass flux.
468
469	integer jlcl(pcols)
470	integer j0(pcols) ! wg detrainment initiation level index.
471	integer jd(pcols) ! wg downdraft initiation level index.
472
473	real(r8) delt ! length of model time-step in seconds.
474
475	integer i
476	integer ii
477	integer k
478	integer msg ! ic number of missing moisture levels at the top of model.
479	real(r8) qdifr
480	real(r8) sdifr
481	!
482	!--------------------------Data statements------------------------------
483	!
484	! Set internal variable "msg" (convection limit) to "limcnv-1"
485	!
486	msg = limcnv - 1
487	!
488	! initialize necessary arrays.
489	! zero out variables not used in cam
490	!
491	qtnd(:,:) = 0._r8
492	heat(:,:) = 0._r8
493	mcon(:,:) = 0._r8
494	rliq(:ncol) = 0._r8
495	!
496	! initialize convective tendencies
497	!
498	prec(:ncol) = 0._r8
499	do k = 1,pver
500	do i = 1,ncol
501	dqdt(i,k) = 0._r8
502	dsdt(i,k) = 0._r8
503	dudt(i,k) = 0._r8
504	dvdt(i,k) = 0._r8
505	pflx(i,k) = 0._r8
506	pflxg(i,k) = 0._r8
507	cme(i,k) = 0._r8
508	rprd(i,k) = 0._r8
509	zdu(i,k) = 0._r8
510	ql(i,k) = 0._r8
511	qlg(i,k) = 0._r8
512	dlf(i,k) = 0._r8
513	dlg(i,k) = 0._r8
514	end do
515	end do
516	do i = 1,ncol
517	pflx(i,pverp) = 0
518	pflxg(i,pverp) = 0
519	end do
520	!
521	do i = 1,ncol
522	pblt(i) = pver
523	dsubcld(i) = 0._r8
524
525	jctop(i) = pver
526	jcbot(i) = 1
527	end do
528	!
529	! calculate local pressure (mbs) and height (m) for both interface
530	! and mid-layer locations.
531	!
532	do i = 1,ncol
533	zs(i) = geos(i)*rgrav
534	pf(i,pver+1) = paph(i,pver+1)*0.01_r8
535	zf(i,pver+1) = zi(i,pver+1) + zs(i)
536	end do
537	do k = 1,pver
538	do i = 1,ncol
539	p(i,k) = pap(i,k)*0.01_r8
540	pf(i,k) = paph(i,k)*0.01_r8
541	z(i,k) = zm(i,k) + zs(i)
542	zf(i,k) = zi(i,k) + zs(i)
543	end do
544	end do
545	!
546	do k = pver - 1,msg + 1,-1
547	do i = 1,ncol
548	if (abs(z(i,k)-zs(i)-pblh(i)) < (zf(i,k)-zf(i,k+1))*0.5_r8) pblt(i) = k
549	end do
550	end do
551	!
552	! store incoming specific humidity field for subsequent calculation
553	! of precipitation (through change in storage).
554	! define dry static energy (normalized by cp).
555	!
556	do k = 1,pver
557	do i = 1,ncol
558	q(i,k) = qh(i,k,1)
559	s(i,k) = t(i,k) + (grav/cpres)*z(i,k)
560	tp(i,k)=0.0_r8
561	shat(i,k) = s(i,k)
562	qhat(i,k) = q(i,k)
563	end do
564	end do
565
566	do i = 1,ncol
567	capeg(i) = 0._r8
568	lclg(i) = 1
569	lelg(i) = pver
570	maxg(i) = 1
571	tlg(i) = 400._r8
572	dsubcld(i) = 0._r8
573	end do
574
575	! Evaluate Tparcel, qsat(Tparcel), buoyancy and CAPE,
576	! lcl, lel, parcel launch level at index maxi()=hmax
577
578	call buoyan_dilute(lchnk ,ncol , &
579	q ,t ,p ,z ,pf , &
580	tp ,qstp ,tl ,rl ,cape , &
581	pblt ,lcl ,lel ,lon ,maxi , &
582	rgas ,grav ,cpres ,msg , &
583	tpert )
584
585	!
586	! determine whether grid points will undergo some deep convection
587	! (ideep=1) or not (ideep=0), based on values of cape,lcl,lel
588	! (require cape.gt. 0 and lel<lcl as minimum conditions).
589	!
590	lengath = 0
591	do i=1,ncol
592	if (cape(i) > capelmt) then
593	lengath = lengath + 1
594	index(lengath) = i
595	end if
596	end do
597
598	if (lengath.eq.0) return
599	do ii=1,lengath
600	i=index(ii)
601	ideep(ii)=i
602	end do
603	!
604	! obtain gathered arrays necessary for ensuing calculations.
605	!
606	do k = 1,pver
607	do i = 1,lengath
608	dp(i,k) = 0.01_r8*dpp(ideep(i),k)
609	qg(i,k) = q(ideep(i),k)
610	tg(i,k) = t(ideep(i),k)
611	pg(i,k) = p(ideep(i),k)
612	zg(i,k) = z(ideep(i),k)
613	sg(i,k) = s(ideep(i),k)
614	tpg(i,k) = tp(ideep(i),k)
615	zfg(i,k) = zf(ideep(i),k)
616	qstpg(i,k) = qstp(ideep(i),k)
617	ug(i,k) = 0._r8
618	vg(i,k) = 0._r8
619	end do
620	end do
621	!
622	do i = 1,lengath
623	zfg(i,pver+1) = zf(ideep(i),pver+1)
624	end do
625	do i = 1,lengath
626	capeg(i) = cape(ideep(i))
627	lclg(i) = lcl(ideep(i))
628	lelg(i) = lel(ideep(i))
629	maxg(i) = maxi(ideep(i))
630	tlg(i) = tl(ideep(i))
631	landfracg(i) = landfrac(ideep(i))
632	end do
633	!
634	! calculate sub-cloud layer pressure "thickness" for use in
635	! closure and tendency routines.
636	!
637	do k = msg + 1,pver
638	do i = 1,lengath
639	if (k >= maxg(i)) then
640	dsubcld(i) = dsubcld(i) + dp(i,k)
641	end if
642	end do
643	end do
644	!
645	! define array of factors (alpha) which defines interfacial
646	! values, as well as interfacial values for (q,s) used in
647	! subsequent routines.
648	!
649	do k = msg + 2,pver
650	do i = 1,lengath
651	! alpha(i,k) = 0.5
652	sdifr = 0._r8
653	qdifr = 0._r8
654	if (sg(i,k) > 0._r8 .or. sg(i,k-1) > 0._r8) &
655	sdifr = abs((sg(i,k)-sg(i,k-1))/max(sg(i,k-1),sg(i,k)))
656	if (qg(i,k) > 0._r8 .or. qg(i,k-1) > 0._r8) &
657	qdifr = abs((qg(i,k)-qg(i,k-1))/max(qg(i,k-1),qg(i,k)))
658	if (sdifr > 1.E-6_r8) then
659	shat(i,k) = log(sg(i,k-1)/sg(i,k))sg(i,k-1)sg(i,k)/(sg(i,k-1)-sg(i,k))
660	else
661	shat(i,k) = 0.5_r8* (sg(i,k)+sg(i,k-1))
662	end if
663	if (qdifr > 1.E-6_r8) then
664	qhat(i,k) = log(qg(i,k-1)/qg(i,k))qg(i,k-1)qg(i,k)/(qg(i,k-1)-qg(i,k))
665	else
666	qhat(i,k) = 0.5_r8* (qg(i,k)+qg(i,k-1))
667	end if
668	end do
669	end do
670	!
671	! obtain cloud properties.
672	!
673
674	call cldprp(lchnk , &
675	qg ,tg ,ug ,vg ,pg , &
676	zg ,sg ,mu ,eu ,du , &
677	md ,ed ,sd ,qd ,mc , &
678	qu ,su ,zfg ,qs ,hmn , &
679	hsat ,shat ,qlg , &
680	cmeg ,maxg ,lelg ,jt ,jlcl , &
681	maxg ,j0 ,jd ,rl ,lengath , &
682	rgas ,grav ,cpres ,msg , &
683	pflxg ,evpg ,cug ,rprdg ,limcnv ,landfracg)
684	!
685	! convert detrainment from units of "1/m" to "1/mb".
686	!
687	do k = msg + 1,pver
688	do i = 1,lengath
689	du (i,k) = du (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
690	eu (i,k) = eu (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
691	ed (i,k) = ed (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
692	cug (i,k) = cug (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
693	cmeg (i,k) = cmeg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
694	rprdg(i,k) = rprdg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
695	evpg (i,k) = evpg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
696	end do
697	end do
698
699	call closure(lchnk , &
700	qg ,tg ,pg ,zg ,sg , &
701	tpg ,qs ,qu ,su ,mc , &
702	du ,mu ,md ,qd ,sd , &
703	qhat ,shat ,dp ,qstpg ,zfg , &
704	qlg ,dsubcld ,mb ,capeg ,tlg , &
705	lclg ,lelg ,jt ,maxg ,1 , &
706	lengath ,rgas ,grav ,cpres ,rl , &
707	msg ,capelmt )
708	!
709	! limit cloud base mass flux to theoretical upper bound.
710	!
711	do i=1,lengath
712	mumax(i) = 0
713	end do
714	do k=msg + 2,pver
715	do i=1,lengath
716	mumax(i) = max(mumax(i), mu(i,k)/dp(i,k))
717	end do
718	end do
719
720	do i=1,lengath
721	if (mumax(i) > 0._r8) then
722	mb(i) = min(mb(i),0.5_r8/(delt*mumax(i)))
723	else
724	mb(i) = 0._r8
725	endif
726	end do
727	! If no_deep_pbl = .true., don't allow convection entirely
728	! within PBL (suggestion of Bjorn Stevens, 8-2000)
729
730	if (no_deep_pbl) then
731	do i=1,lengath
732	if (zm(ideep(i),jt(i)) < pblh(ideep(i))) mb(i) = 0
733	end do
734	end if
735
736
737	do k=msg+1,pver
738	do i=1,lengath
739	mu (i,k) = mu (i,k)*mb(i)
740	md (i,k) = md (i,k)*mb(i)
741	mc (i,k) = mc (i,k)*mb(i)
742	du (i,k) = du (i,k)*mb(i)
743	eu (i,k) = eu (i,k)*mb(i)
744	ed (i,k) = ed (i,k)*mb(i)
745	cmeg (i,k) = cmeg (i,k)*mb(i)
746	rprdg(i,k) = rprdg(i,k)*mb(i)
747	cug (i,k) = cug (i,k)*mb(i)
748	evpg (i,k) = evpg (i,k)*mb(i)
749	pflxg(i,k+1)= pflxg(i,k+1)mb(i)100._r8/grav
750	end do
751	end do
752	!
753	! compute temperature and moisture changes due to convection.
754	!
755	call q1q2_pjr(lchnk , &
756	dqdt ,dsdt ,qg ,qs ,qu , &
757	su ,du ,qhat ,shat ,dp , &
758	mu ,md ,sd ,qd ,qlg , &
759	dsubcld ,jt ,maxg ,1 ,lengath , &
760	cpres ,rl ,msg , &
761	dlg ,evpg ,cug )
762	!
763	! gather back temperature and mixing ratio.
764	!
765	do k = msg + 1,pver
766	!DIR$ CONCURRENT
767	do i = 1,lengath
768	!
769	! q is updated to compute net precip.
770	!
771	q(ideep(i),k) = qh(ideep(i),k,1) + 2._r8deltdqdt(i,k)
772	qtnd(ideep(i),k) = dqdt (i,k)
773	cme (ideep(i),k) = cmeg (i,k)
774	rprd(ideep(i),k) = rprdg(i,k)
775	zdu (ideep(i),k) = du (i,k)
776	mcon(ideep(i),k) = mc (i,k)
777	heat(ideep(i),k) = dsdt (i,k)*cpres
778	dlf (ideep(i),k) = dlg (i,k)
779	pflx(ideep(i),k) = pflxg(i,k)
780	ql (ideep(i),k) = qlg (i,k)
781	end do
782	end do
783	!
784	!DIR$ CONCURRENT
785	do i = 1,lengath
786	jctop(ideep(i)) = jt(i)
787	!++bee
788	jcbot(ideep(i)) = maxg(i)
789	!--bee
790	pflx(ideep(i),pverp) = pflxg(i,pverp)
791	end do
792
793	! Compute precip by integrating change in water vapor minus detrained cloud water
794	do k = pver,msg + 1,-1
795	do i = 1,ncol
796	prec(i) = prec(i) - dpp(i,k)* (q(i,k)-qh(i,k,1)) - dpp(i,k)dlf(i,k)2*delt
797	end do
798	end do
799
800	! obtain final precipitation rate in m/s.
801	do i = 1,ncol
802	prec(i) = rgravmax(prec(i),0._r8)/ (2._r8delt)/1000._r8
803	end do
804
805	! Compute reserved liquid (not yet in cldliq) for energy integrals.
806	! Treat rliq as flux out bottom, to be added back later.
807	do k = 1, pver
808	do i = 1, ncol
809	rliq(i) = rliq(i) + dlf(i,k)*dpp(i,k)/gravit
810	end do
811	end do
812	rliq(:ncol) = rliq(:ncol) /1000._r8
813
814	return
815	end subroutine zm_convr
816
817	!===============================================================================
818	subroutine zm_conv_evap(ncol,lchnk, &
819	t,pmid,pdel,q, &
820	tend_s, tend_s_snwprd, tend_s_snwevmlt, tend_q, &
821	prdprec, cldfrc, deltat, &
822	prec, snow, ntprprd, ntsnprd, flxprec, flxsnow )
823
824	!-----------------------------------------------------------------------
825	! Compute tendencies due to evaporation of rain from ZM scheme
826	!--
827	! Compute the total precipitation and snow fluxes at the surface.
828	! Add in the latent heat of fusion for snow formation and melt, since it not dealt with
829	! in the Zhang-MacFarlane parameterization.
830	! Evaporate some of the precip directly into the environment using a Sundqvist type algorithm
831	!-----------------------------------------------------------------------
832
833	use wv_saturation, only: aqsat
834	#ifndef WRF_PORT
835	use phys_grid, only: get_rlat_all_p
836	#endif
837
838	!------------------------------Arguments--------------------------------
839	integer,intent(in) :: ncol, lchnk ! number of columns and chunk index
840	real(r8),intent(in), dimension(pcols,pver) :: t ! temperature (K)
841	real(r8),intent(in), dimension(pcols,pver) :: pmid ! midpoint pressure (Pa)
842	real(r8),intent(in), dimension(pcols,pver) :: pdel ! layer thickness (Pa)
843	real(r8),intent(in), dimension(pcols,pver) :: q ! water vapor (kg/kg)
844	real(r8),intent(inout), dimension(pcols,pver) :: tend_s ! heating rate (J/kg/s)
845	real(r8),intent(inout), dimension(pcols,pver) :: tend_q ! water vapor tendency (kg/kg/s)
846	real(r8),intent(out ), dimension(pcols,pver) :: tend_s_snwprd ! Heating rate of snow production
847	real(r8),intent(out ), dimension(pcols,pver) :: tend_s_snwevmlt ! Heating rate of evap/melting of snow
848
849
850
851	real(r8), intent(in ) :: prdprec(pcols,pver)! precipitation production (kg/ks/s)
852	real(r8), intent(in ) :: cldfrc(pcols,pver) ! cloud fraction
853	real(r8), intent(in ) :: deltat ! time step
854
855	real(r8), intent(inout) :: prec(pcols) ! Convective-scale preciptn rate
856	real(r8), intent(out) :: snow(pcols) ! Convective-scale snowfall rate
857	!
858	!---------------------------Local storage-------------------------------
859
860	real(r8) :: est (pcols,pver) ! Saturation vapor pressure
861	real(r8) :: fice (pcols,pver) ! ice fraction in precip production
862	real(r8) :: fsnow_conv(pcols,pver) ! snow fraction in precip production
863	real(r8) :: qsat (pcols,pver) ! saturation specific humidity
864	real(r8),intent(out) :: flxprec(pcols,pverp) ! Convective-scale flux of precip at interfaces (kg/m2/s)
865	real(r8),intent(out) :: flxsnow(pcols,pverp) ! Convective-scale flux of snow at interfaces (kg/m2/s)
866	real(r8),intent(out) :: ntprprd(pcols,pver) ! net precip production in layer
867	real(r8),intent(out) :: ntsnprd(pcols,pver) ! net snow production in layer
868	real(r8) :: work1 ! temp variable (pjr)
869	real(r8) :: work2 ! temp variable (pjr)
870
871	real(r8) :: evpvint(pcols) ! vertical integral of evaporation
872	real(r8) :: evpprec(pcols) ! evaporation of precipitation (kg/kg/s)
873	real(r8) :: evpsnow(pcols) ! evaporation of snowfall (kg/kg/s)
874	real(r8) :: snowmlt(pcols) ! snow melt tendency in layer
875	real(r8) :: flxsntm(pcols) ! flux of snow into layer, after melting
876
877	real(r8) :: evplimit ! temp variable for evaporation limits
878	real(r8) :: rlat(pcols)
879
880	integer :: i,k ! longitude,level indices
881
882
883	!-----------------------------------------------------------------------
884
885	! convert input precip to kg/m2/s
886	prec(:ncol) = prec(:ncol)*1000._r8
887
888	! determine saturation vapor pressure
889	call aqsat (t ,pmid ,est ,qsat ,pcols , &
890	ncol ,pver ,1 ,pver )
891
892	! determine ice fraction in rain production (use cloud water parameterization fraction at present)
893	call cldwat_fice(ncol, t, fice, fsnow_conv)
894
895	! zero the flux integrals on the top boundary
896	flxprec(:ncol,1) = 0._r8
897	flxsnow(:ncol,1) = 0._r8
898	evpvint(:ncol) = 0._r8
899
900	do k = 1, pver
901	do i = 1, ncol
902
903	! Melt snow falling into layer, if necessary.
904	if (t(i,k) > tmelt) then
905	flxsntm(i) = 0._r8
906	snowmlt(i) = flxsnow(i,k) * gravit/ pdel(i,k)
907	else
908	flxsntm(i) = flxsnow(i,k)
909	snowmlt(i) = 0._r8
910	end if
911
912	! relative humidity depression must be > 0 for evaporation
913	evplimit = max(1._r8 - q(i,k)/qsat(i,k), 0._r8)
914
915	! total evaporation depends on flux in the top of the layer
916	! flux prec is the net production above layer minus evaporation into environmet
917	evpprec(i) = ke * (1._r8 - cldfrc(i,k)) * evplimit * sqrt(flxprec(i,k))
918	!**********************************************************
919	!! evpprec(i) = 0. ! turn off evaporation for now
920	!**********************************************************
921
922	! Don't let evaporation supersaturate layer (approx). Layer may already be saturated.
923	! Currently does not include heating/cooling change to qsat
924	evplimit = max(0._r8, (qsat(i,k)-q(i,k)) / deltat)
925
926	! Don't evaporate more than is falling into the layer - do not evaporate rain formed
927	! in this layer but if precip production is negative, remove from the available precip
928	! Negative precip production occurs because of evaporation in downdrafts.
929	!!$ evplimit = flxprec(i,k) * gravit / pdel(i,k) + min(prdprec(i,k), 0.)
930	evplimit = min(evplimit, flxprec(i,k) * gravit / pdel(i,k))
931
932	! Total evaporation cannot exceed input precipitation
933	evplimit = min(evplimit, (prec(i) - evpvint(i)) * gravit / pdel(i,k))
934
935	evpprec(i) = min(evplimit, evpprec(i))
936
937	! evaporation of snow depends on snow fraction of total precipitation in the top after melting
938	if (flxprec(i,k) > 0._r8) then
939	! evpsnow(i) = evpprec(i) * flxsntm(i) / flxprec(i,k)
940	! prevent roundoff problems
941	work1 = min(max(0._r8,flxsntm(i)/flxprec(i,k)),1._r8)
942	evpsnow(i) = evpprec(i) * work1
943	else
944	evpsnow(i) = 0._r8
945	end if
946
947	! vertically integrated evaporation
948	evpvint(i) = evpvint(i) + evpprec(i) * pdel(i,k)/gravit
949
950	! net precip production is production - evaporation
951	ntprprd(i,k) = prdprec(i,k) - evpprec(i)
952	! net snow production is precip production * ice fraction - evaporation - melting
953	!pjrworks ntsnprd(i,k) = prdprec(i,k)*fice(i,k) - evpsnow(i) - snowmlt(i)
954	!pjrwrks2 ntsnprd(i,k) = prdprec(i,k)*fsnow_conv(i,k) - evpsnow(i) - snowmlt(i)
955	! the small amount added to flxprec in the work1 expression has been increased from
956	! 1e-36 to 8.64e-11 (1e-5 mm/day). This causes the temperature based partitioning
957	! scheme to be used for small flxprec amounts. This is to address error growth problems.
958	#ifdef PERGRO
959	work1 = min(max(0._r8,flxsnow(i,k)/(flxprec(i,k)+8.64e-11_r8)),1._r8)
960	#else
961	if (flxprec(i,k).gt.0._r8) then
962	work1 = min(max(0._r8,flxsnow(i,k)/flxprec(i,k)),1._r8)
963	else
964	work1 = 0._r8
965	endif
966	#endif
967	work2 = max(fsnow_conv(i,k), work1)
968	if (snowmlt(i).gt.0._r8) work2 = 0._r8
969	! work2 = fsnow_conv(i,k)
970	ntsnprd(i,k) = prdprec(i,k)*work2 - evpsnow(i) - snowmlt(i)
971	tend_s_snwprd (i,k) = prdprec(i,k)work2latice
972	tend_s_snwevmlt(i,k) = - ( evpsnow(i) + snowmlt(i) )*latice
973
974	! precipitation fluxes
975	flxprec(i,k+1) = flxprec(i,k) + ntprprd(i,k) * pdel(i,k)/gravit
976	flxsnow(i,k+1) = flxsnow(i,k) + ntsnprd(i,k) * pdel(i,k)/gravit
977
978	! protect against rounding error
979	flxprec(i,k+1) = max(flxprec(i,k+1), 0._r8)
980	flxsnow(i,k+1) = max(flxsnow(i,k+1), 0._r8)
981	! more protection (pjr)
982	! flxsnow(i,k+1) = min(flxsnow(i,k+1), flxprec(i,k+1))
983
984	! heating (cooling) and moistening due to evaporation
985	! - latent heat of vaporization for precip production has already been accounted for
986	! - snow is contained in prec
987	tend_s(i,k) =-evpprec(i)latvap + ntsnprd(i,k)latice
988	tend_q(i,k) = evpprec(i)
989	end do
990	end do
991
992	! set output precipitation rates (m/s)
993	prec(:ncol) = flxprec(:ncol,pver+1) / 1000._r8
994	snow(:ncol) = flxsnow(:ncol,pver+1) / 1000._r8
995
996	!**********************************************************
997	!!$ tend_s(:ncol,:) = 0. ! turn heating off
998	!**********************************************************
999
1000	end subroutine zm_conv_evap
1001
1002
1003
1004	subroutine convtran(lchnk , &
1005	doconvtran,q ,ncnst ,mu ,md , &
1006	du ,eu ,ed ,dp ,dsubcld , &
1007	jt ,mx ,ideep ,il1g ,il2g , &
1008	nstep ,fracis ,dqdt ,dpdry )
1009	!-----------------------------------------------------------------------
1010	!
1011	! Purpose:
1012	! Convective transport of trace species
1013	!
1014	! Mixing ratios may be with respect to either dry or moist air
1015	!
1016	! Method:
1017	! <Describe the algorithm(s) used in the routine.>
1018	! <Also include any applicable external references.>
1019	!
1020	! Author: P. Rasch
1021	!
1022	!-----------------------------------------------------------------------
1023	use shr_kind_mod, only: r8 => shr_kind_r8
1024	use constituents, only: cnst_get_type_byind
1025	#ifndef WRF_PORT
1026	use ppgrid
1027	use abortutils, only: endrun
1028	!~~~#else
1029	! use module_cam_support, only: cnst_get_type_byind
1030	#endif
1031
1032	implicit none
1033	!-----------------------------------------------------------------------
1034	!
1035	! Input arguments
1036	!
1037	integer, intent(in) :: lchnk ! chunk identifier
1038	integer, intent(in) :: ncnst ! number of tracers to transport
1039	logical, intent(in) :: doconvtran(ncnst) ! flag for doing convective transport
1040	real(r8), intent(in) :: q(pcols,pver,ncnst) ! Tracer array including moisture
1041	real(r8), intent(in) :: mu(pcols,pver) ! Mass flux up
1042	real(r8), intent(in) :: md(pcols,pver) ! Mass flux down
1043	real(r8), intent(in) :: du(pcols,pver) ! Mass detraining from updraft
1044	real(r8), intent(in) :: eu(pcols,pver) ! Mass entraining from updraft
1045	real(r8), intent(in) :: ed(pcols,pver) ! Mass entraining from downdraft
1046	real(r8), intent(in) :: dp(pcols,pver) ! Delta pressure between interfaces
1047	real(r8), intent(in) :: dsubcld(pcols) ! Delta pressure from cloud base to sfc
1048	real(r8), intent(in) :: fracis(pcols,pver,ncnst) ! fraction of tracer that is insoluble
1049
1050	integer, intent(in) :: jt(pcols) ! Index of cloud top for each column
1051	integer, intent(in) :: mx(pcols) ! Index of cloud top for each column
1052	integer, intent(in) :: ideep(pcols) ! Gathering array
1053	integer, intent(in) :: il1g ! Gathered min lon indices over which to operate
1054	integer, intent(in) :: il2g ! Gathered max lon indices over which to operate
1055	integer, intent(in) :: nstep ! Time step index
1056
1057	real(r8), intent(in) :: dpdry(pcols,pver) ! Delta pressure between interfaces
1058
1059
1060	! input/output
1061
1062	real(r8), intent(out) :: dqdt(pcols,pver,ncnst) ! Tracer tendency array
1063
1064	!--------------------------Local Variables------------------------------
1065
1066	integer i ! Work index
1067	integer k ! Work index
1068	integer kbm ! Highest altitude index of cloud base
1069	integer kk ! Work index
1070	integer kkp1 ! Work index
1071	integer km1 ! Work index
1072	integer kp1 ! Work index
1073	integer ktm ! Highest altitude index of cloud top
1074	integer m ! Work index
1075
1076	real(r8) cabv ! Mix ratio of constituent above
1077	real(r8) cbel ! Mix ratio of constituent below
1078	real(r8) cdifr ! Normalized diff between cabv and cbel
1079	real(r8) chat(pcols,pver) ! Mix ratio in env at interfaces
1080	real(r8) cond(pcols,pver) ! Mix ratio in downdraft at interfaces
1081	real(r8) const(pcols,pver) ! Gathered tracer array
1082	real(r8) fisg(pcols,pver) ! gathered insoluble fraction of tracer
1083	real(r8) conu(pcols,pver) ! Mix ratio in updraft at interfaces
1084	real(r8) dcondt(pcols,pver) ! Gathered tend array
1085	real(r8) small ! A small number
1086	real(r8) mbsth ! Threshold for mass fluxes
1087	real(r8) mupdudp ! A work variable
1088	real(r8) minc ! A work variable
1089	real(r8) maxc ! A work variable
1090	real(r8) fluxin ! A work variable
1091	real(r8) fluxout ! A work variable
1092	real(r8) netflux ! A work variable
1093
1094	real(r8) dutmp(pcols,pver) ! Mass detraining from updraft
1095	real(r8) eutmp(pcols,pver) ! Mass entraining from updraft
1096	real(r8) edtmp(pcols,pver) ! Mass entraining from downdraft
1097	real(r8) dptmp(pcols,pver) ! Delta pressure between interfaces
1098	!-----------------------------------------------------------------------
1099	!
1100	small = 1.e-36_r8
1101
1102	! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
1103	mbsth = 1.e-15_r8
1104
1105	! Find the highest level top and bottom levels of convection
1106	ktm = pver
1107	kbm = pver
1108	do i = il1g, il2g
1109	ktm = min(ktm,jt(i))
1110	kbm = min(kbm,mx(i))
1111	end do
1112
1113	! Loop ever each constituent
1114	do m = 2, ncnst
1115	if (doconvtran(m)) then
1116
1117	if (cnst_get_type_byind(m).eq.'dry') then
1118	do k = 1,pver
1119	do i =il1g,il2g
1120	dptmp(i,k) = dpdry(i,k)
1121	dutmp(i,k) = du(i,k)*dp(i,k)/dpdry(i,k)
1122	eutmp(i,k) = eu(i,k)*dp(i,k)/dpdry(i,k)
1123	edtmp(i,k) = ed(i,k)*dp(i,k)/dpdry(i,k)
1124	end do
1125	end do
1126	else
1127	do k = 1,pver
1128	do i =il1g,il2g
1129	dptmp(i,k) = dp(i,k)
1130	dutmp(i,k) = du(i,k)
1131	eutmp(i,k) = eu(i,k)
1132	edtmp(i,k) = ed(i,k)
1133	end do
1134	end do
1135	endif
1136	! dptmp = dp
1137
1138	! Gather up the constituent and set tend to zero
1139	do k = 1,pver
1140	do i =il1g,il2g
1141	const(i,k) = q(ideep(i),k,m)
1142	fisg(i,k) = fracis(ideep(i),k,m)
1143	end do
1144	end do
1145
1146	! From now on work only with gathered data
1147
1148	! Interpolate environment tracer values to interfaces
1149	do k = 1,pver
1150	km1 = max(1,k-1)
1151	do i = il1g, il2g
1152	minc = min(const(i,km1),const(i,k))
1153	maxc = max(const(i,km1),const(i,k))
1154	if (minc < 0) then
1155	cdifr = 0._r8
1156	else
1157	cdifr = abs(const(i,k)-const(i,km1))/max(maxc,small)
1158	endif
1159
1160	! If the two layers differ significantly use a geometric averaging
1161	! procedure
1162	if (cdifr > 1.E-6_r8) then
1163	cabv = max(const(i,km1),maxc*1.e-12_r8)
1164	cbel = max(const(i,k),maxc*1.e-12_r8)
1165	chat(i,k) = log(cabv/cbel)/(cabv-cbel)cabvcbel
1166
1167	else ! Small diff, so just arithmetic mean
1168	chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1))
1169	end if
1170
1171	! Provisional up and down draft values
1172	conu(i,k) = chat(i,k)
1173	cond(i,k) = chat(i,k)
1174
1175	! provisional tends
1176	dcondt(i,k) = 0._r8
1177
1178	end do
1179	end do
1180
1181	! Do levels adjacent to top and bottom
1182	k = 2
1183	km1 = 1
1184	kk = pver
1185	do i = il1g,il2g
1186	mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
1187	if (mupdudp > mbsth) then
1188	conu(i,kk) = (+eutmp(i,kk)fisg(i,kk)const(i,kk)*dptmp(i,kk))/mupdudp
1189	endif
1190	if (md(i,k) < -mbsth) then
1191	cond(i,k) = (-edtmp(i,km1)fisg(i,km1)const(i,km1)*dptmp(i,km1))/md(i,k)
1192	endif
1193	end do
1194
1195	! Updraft from bottom to top
1196	do kk = pver-1,1,-1
1197	kkp1 = min(pver,kk+1)
1198	do i = il1g,il2g
1199	mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk)
1200	if (mupdudp > mbsth) then
1201	conu(i,kk) = ( mu(i,kkp1)conu(i,kkp1)+eutmp(i,kk)fisg(i,kk)* &
1202	const(i,kk)*dptmp(i,kk) )/mupdudp
1203	endif
1204	end do
1205	end do
1206
1207	! Downdraft from top to bottom
1208	do k = 3,pver
1209	km1 = max(1,k-1)
1210	do i = il1g,il2g
1211	if (md(i,k) < -mbsth) then
1212	cond(i,k) = ( md(i,km1)cond(i,km1)-edtmp(i,km1)fisg(i,km1)*const(i,km1) &
1213	*dptmp(i,km1) )/md(i,k)
1214	endif
1215	end do
1216	end do
1217
1218
1219	do k = ktm,pver
1220	km1 = max(1,k-1)
1221	kp1 = min(pver,k+1)
1222	do i = il1g,il2g
1223
1224	! version 1 hard to check for roundoff errors
1225	! dcondt(i,k) =
1226	! $ +(+mu(i,kp1)* (conu(i,kp1)-chat(i,kp1))
1227	! $ -mu(i,k)* (conu(i,k)-chat(i,k))
1228	! $ +md(i,kp1)* (cond(i,kp1)-chat(i,kp1))
1229	! $ -md(i,k)* (cond(i,k)-chat(i,k))
1230	! $ )/dp(i,k)
1231
1232	! version 2 hard to limit fluxes
1233	! fluxin = mu(i,kp1)conu(i,kp1) + mu(i,k)chat(i,k)
1234	! $ -(md(i,k) cond(i,k) + md(i,kp1)chat(i,kp1))
1235	! fluxout = mu(i,k)conu(i,k) + mu(i,kp1)chat(i,kp1)
1236	! $ -(md(i,kp1)cond(i,kp1) + md(i,k)chat(i,k))
1237
1238	! version 3 limit fluxes outside convection to mass in appropriate layer
1239	! these limiters are probably only safe for positive definite quantitities
1240	! it assumes that mu and md already satify a courant number limit of 1
1241	fluxin = mu(i,kp1)conu(i,kp1)+ mu(i,k)min(chat(i,k),const(i,km1)) &
1242	-(md(i,k) cond(i,k) + md(i,kp1)min(chat(i,kp1),const(i,kp1)))
1243	fluxout = mu(i,k)conu(i,k) + mu(i,kp1)min(chat(i,kp1),const(i,k)) &
1244	-(md(i,kp1)cond(i,kp1) + md(i,k)min(chat(i,k),const(i,k)))
1245
1246	netflux = fluxin - fluxout
1247	if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then
1248	netflux = 0._r8
1249	endif
1250	dcondt(i,k) = netflux/dptmp(i,k)
1251	end do
1252	end do
1253	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1254	!
1255	!DIR$ NOINTERCHANGE
1256	do k = kbm,pver
1257	km1 = max(1,k-1)
1258	do i = il1g,il2g
1259	if (k == mx(i)) then
1260
1261	! version 1
1262	! dcondt(i,k) = (1./dsubcld(i))*
1263	! $ (-mu(i,k)*(conu(i,k)-chat(i,k))
1264	! $ -md(i,k)*(cond(i,k)-chat(i,k))
1265	! $ )
1266
1267	! version 2
1268	! fluxin = mu(i,k)chat(i,k) - md(i,k)cond(i,k)
1269	! fluxout = mu(i,k)conu(i,k) - md(i,k)chat(i,k)
1270	! version 3
1271	fluxin = mu(i,k)min(chat(i,k),const(i,km1)) - md(i,k)cond(i,k)
1272	fluxout = mu(i,k)conu(i,k) - md(i,k)min(chat(i,k),const(i,k))
1273
1274	netflux = fluxin - fluxout
1275	if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then
1276	netflux = 0._r8
1277	endif
1278	! dcondt(i,k) = netflux/dsubcld(i)
1279	dcondt(i,k) = netflux/dptmp(i,k)
1280	else if (k > mx(i)) then
1281	! dcondt(i,k) = dcondt(i,k-1)
1282	dcondt(i,k) = 0._r8
1283	end if
1284	end do
1285	end do
1286
1287	! Initialize to zero everywhere, then scatter tendency back to full array
1288	dqdt(:,:,m) = 0._r8
1289	do k = 1,pver
1290	kp1 = min(pver,k+1)
1291	!DIR$ CONCURRENT
1292	do i = il1g,il2g
1293	dqdt(ideep(i),k,m) = dcondt(i,k)
1294	end do
1295	end do
1296
1297	end if ! for doconvtran
1298
1299	end do
1300
1301	return
1302	end subroutine convtran
1303
1304	!=========================================================================================
1305
1306	subroutine momtran(lchnk, ncol, &
1307	domomtran,q ,ncnst ,mu ,md , &
1308	du ,eu ,ed ,dp ,dsubcld , &
1309	jt ,mx ,ideep ,il1g ,il2g , &
1310	nstep ,dqdt ,pguall ,pgdall, icwu, icwd, dt, seten )
1311	!-----------------------------------------------------------------------
1312	!
1313	! Purpose:
1314	! Convective transport of momentum
1315	!
1316	! Mixing ratios may be with respect to either dry or moist air
1317	!
1318	! Method:
1319	! Based on the convtran subroutine by P. Rasch
1320	! <Also include any applicable external references.>
1321	!
1322	! Author: J. Richter and P. Rasch
1323	!
1324	!-----------------------------------------------------------------------
1325	use shr_kind_mod, only: r8 => shr_kind_r8
1326	#ifndef WRF_PORT
1327	use constituents, only: cnst_get_type_byind
1328	use ppgrid
1329	use abortutils, only: endrun
1330	#endif
1331	implicit none
1332	!-----------------------------------------------------------------------
1333	!
1334	! Input arguments
1335	!
1336	integer, intent(in) :: lchnk ! chunk identifier
1337	integer, intent(in) :: ncol ! number of atmospheric columns
1338	integer, intent(in) :: ncnst ! number of tracers to transport
1339	logical, intent(in) :: domomtran(ncnst) ! flag for doing convective transport
1340	real(r8), intent(in) :: q(pcols,pver,ncnst) ! Wind array
1341	real(r8), intent(in) :: mu(pcols,pver) ! Mass flux up
1342	real(r8), intent(in) :: md(pcols,pver) ! Mass flux down
1343	real(r8), intent(in) :: du(pcols,pver) ! Mass detraining from updraft
1344	real(r8), intent(in) :: eu(pcols,pver) ! Mass entraining from updraft
1345	real(r8), intent(in) :: ed(pcols,pver) ! Mass entraining from downdraft
1346	real(r8), intent(in) :: dp(pcols,pver) ! Delta pressure between interfaces
1347	real(r8), intent(in) :: dsubcld(pcols) ! Delta pressure from cloud base to sfc
1348	real(r8), intent(in) :: dt ! time step in seconds : 2*delta_t
1349
1350	integer, intent(in) :: jt(pcols) ! Index of cloud top for each column
1351	integer, intent(in) :: mx(pcols) ! Index of cloud top for each column
1352	integer, intent(in) :: ideep(pcols) ! Gathering array
1353	integer, intent(in) :: il1g ! Gathered min lon indices over which to operate
1354	integer, intent(in) :: il2g ! Gathered max lon indices over which to operate
1355	integer, intent(in) :: nstep ! Time step index
1356
1357
1358
1359	! input/output
1360
1361	real(r8), intent(out) :: dqdt(pcols,pver,ncnst) ! Tracer tendency array
1362
1363	!--------------------------Local Variables------------------------------
1364
1365	integer i ! Work index
1366	integer k ! Work index
1367	integer kbm ! Highest altitude index of cloud base
1368	integer kk ! Work index
1369	integer kkp1 ! Work index
1370	integer kkm1 ! Work index
1371	integer km1 ! Work index
1372	integer kp1 ! Work index
1373	integer ktm ! Highest altitude index of cloud top
1374	integer m ! Work index
1375	integer ii ! Work index
1376
1377	real(r8) cabv ! Mix ratio of constituent above
1378	real(r8) cbel ! Mix ratio of constituent below
1379	real(r8) cdifr ! Normalized diff between cabv and cbel
1380	real(r8) chat(pcols,pver) ! Mix ratio in env at interfaces
1381	real(r8) cond(pcols,pver) ! Mix ratio in downdraft at interfaces
1382	real(r8) const(pcols,pver) ! Gathered wind array
1383	real(r8) conu(pcols,pver) ! Mix ratio in updraft at interfaces
1384	real(r8) dcondt(pcols,pver) ! Gathered tend array
1385	real(r8) small ! A small number
1386	real(r8) mbsth ! Threshold for mass fluxes
1387	real(r8) mupdudp ! A work variable
1388	real(r8) minc ! A work variable
1389	real(r8) maxc ! A work variable
1390	real(r8) fluxin ! A work variable
1391	real(r8) fluxout ! A work variable
1392	real(r8) netflux ! A work variable
1393
1394	real(r8) momcu ! constant for updraft pressure gradient term
1395	real(r8) momcd ! constant for downdraft pressure gradient term
1396	real(r8) sum ! sum
1397	real(r8) sum2 ! sum2
1398
1399	real(r8) mududp(pcols,pver) ! working variable
1400	real(r8) mddudp(pcols,pver) ! working variable
1401
1402	real(r8) pgu(pcols,pver) ! Pressure gradient term for updraft
1403	real(r8) pgd(pcols,pver) ! Pressure gradient term for downdraft
1404
1405	real(r8),intent(out) :: pguall(pcols,pver,ncnst) ! Apparent force from updraft PG
1406	real(r8),intent(out) :: pgdall(pcols,pver,ncnst) ! Apparent force from downdraft PG
1407
1408	real(r8),intent(out) :: icwu(pcols,pver,ncnst) ! In-cloud winds in updraft
1409	real(r8),intent(out) :: icwd(pcols,pver,ncnst) ! In-cloud winds in downdraft
1410
1411	real(r8),intent(out) :: seten(pcols,pver) ! Dry static energy tendency
1412	real(r8) gseten(pcols,pver) ! Gathered dry static energy tendency
1413
1414	real(r8) mflux(pcols,pverp,ncnst) ! Gathered momentum flux
1415
1416	real(r8) wind0(pcols,pver,ncnst) ! gathered wind before time step
1417	real(r8) windf(pcols,pver,ncnst) ! gathered wind after time step
1418	real(r8) fkeb, fket, ketend_cons, ketend, utop, ubot, vtop, vbot, gset2
1419
1420
1421	!-----------------------------------------------------------------------
1422	!
1423
1424	! Initialize outgoing fields
1425	pguall(:,:,:) = 0.0_r8
1426	pgdall(:,:,:) = 0.0_r8
1427	! Initialize in-cloud winds to environmental wind
1428	icwu(:ncol,:,:) = q(:ncol,:,:)
1429	icwd(:ncol,:,:) = q(:ncol,:,:)
1430
1431	! Initialize momentum flux and final winds
1432	mflux(:,:,:) = 0.0_r8
1433	wind0(:,:,:) = 0.0_r8
1434	windf(:,:,:) = 0.0_r8
1435
1436	! Initialize dry static energy
1437
1438	seten(:,:) = 0.0_r8
1439	gseten(:,:) = 0.0_r8
1440
1441	! Define constants for parameterization
1442
1443	momcu = 0.4_r8
1444	momcd = 0.4_r8
1445
1446	small = 1.e-36_r8
1447	! mbsth is the threshold below which we treat the mass fluxes as zero (in mb/s)
1448	mbsth = 1.e-15_r8
1449
1450	! Find the highest level top and bottom levels of convection
1451	ktm = pver
1452	kbm = pver
1453	do i = il1g, il2g
1454	ktm = min(ktm,jt(i))
1455	kbm = min(kbm,mx(i))
1456	end do
1457
1458	! Loop ever each wind component
1459	do m = 1, ncnst !start at m = 1 to transport momentum
1460	if (domomtran(m)) then
1461
1462	! Gather up the winds and set tend to zero
1463	do k = 1,pver
1464	do i =il1g,il2g
1465	const(i,k) = q(ideep(i),k,m)
1466	wind0(i,k,m) = const(i,k)
1467	end do
1468	end do
1469
1470
1471	! From now on work only with gathered data
1472
1473	! Interpolate winds to interfaces
1474
1475	do k = 1,pver
1476	km1 = max(1,k-1)
1477	do i = il1g, il2g
1478
1479	! use arithmetic mean
1480	chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1))
1481
1482	! Provisional up and down draft values
1483	conu(i,k) = chat(i,k)
1484	cond(i,k) = chat(i,k)
1485
1486	! provisional tends
1487	dcondt(i,k) = 0._r8
1488
1489	end do
1490	end do
1491
1492
1493	!
1494	! Pressure Perturbation Term
1495	!
1496
1497	!Top boundary: assume mu is zero
1498
1499	k=1
1500	pgu(:il2g,k) = 0.0_r8
1501	pgd(:il2g,k) = 0.0_r8
1502
1503	do k=2,pver-1
1504	km1 = max(1,k-1)
1505	kp1 = min(pver,k+1)
1506	do i = il1g,il2g
1507
1508	!interior points
1509
1510	mududp(i,k) = ( mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
1511	+ mu(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
1512
1513	pgu(i,k) = - momcu * 0.5_r8 * mududp(i,k)
1514
1515
1516	mddudp(i,k) = ( md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) &
1517	+ md(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k))
1518
1519	pgd(i,k) = - momcd * 0.5_r8 * mddudp(i,k)
1520
1521
1522	end do
1523	end do
1524
1525	! bottom boundary
1526	k = pver
1527	km1 = max(1,k-1)
1528	do i=il1g,il2g
1529
1530	mududp(i,k) = mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
1531	pgu(i,k) = - momcu * mududp(i,k)
1532
1533	mddudp(i,k) = md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1)
1534
1535	pgd(i,k) = - momcd * mddudp(i,k)
1536
1537	end do
1538
1539
1540	!
1541	! In-cloud velocity calculations
1542	!
1543
1544	! Do levels adjacent to top and bottom
1545	k = 2
1546	km1 = 1
1547	kk = pver
1548	kkm1 = max(1,kk-1)
1549	do i = il1g,il2g
1550	mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
1551	if (mupdudp > mbsth) then
1552
1553	conu(i,kk) = (+eu(i,kk)const(i,kk)dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp
1554	endif
1555	if (md(i,k) < -mbsth) then
1556	cond(i,k) = (-ed(i,km1)const(i,km1)dp(i,km1))-pgd(i,km1)*dp(i,km1)/md(i,k)
1557	endif
1558
1559
1560	end do
1561
1562
1563
1564	! Updraft from bottom to top
1565	do kk = pver-1,1,-1
1566	kkm1 = max(1,kk-1)
1567	kkp1 = min(pver,kk+1)
1568	do i = il1g,il2g
1569	mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk)
1570	if (mupdudp > mbsth) then
1571
1572	conu(i,kk) = ( mu(i,kkp1)conu(i,kkp1)+eu(i,kk) &
1573	const(i,kk)dp(i,kk)+pgu(i,kk)dp(i,kk))/mupdudp
1574	endif
1575	end do
1576
1577	end do
1578
1579
1580	! Downdraft from top to bottom
1581	do k = 3,pver
1582	km1 = max(1,k-1)
1583	do i = il1g,il2g
1584	if (md(i,k) < -mbsth) then
1585
1586	cond(i,k) = ( md(i,km1)cond(i,km1)-ed(i,km1)const(i,km1) &
1587	dp(i,km1)-pgd(i,km1)dp(i,km1) )/md(i,k)
1588
1589	endif
1590	end do
1591	end do
1592
1593
1594	sum = 0._r8
1595	sum2 = 0._r8
1596
1597
1598	do k = ktm,pver
1599	km1 = max(1,k-1)
1600	kp1 = min(pver,k+1)
1601	do i = il1g,il2g
1602	ii = ideep(i)
1603
1604	! version 1 hard to check for roundoff errors
1605	dcondt(i,k) = &
1606	+(mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) &
1607	-mu(i,k)* (conu(i,k)-chat(i,k)) &
1608	+md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) &
1609	-md(i,k)* (cond(i,k)-chat(i,k)) &
1610	)/dp(i,k)
1611
1612	end do
1613	end do
1614
1615	! dcont for bottom layer
1616	!
1617	!DIR$ NOINTERCHANGE
1618	do k = kbm,pver
1619	km1 = max(1,k-1)
1620	do i = il1g,il2g
1621	if (k == mx(i)) then
1622
1623	! version 1
1624	dcondt(i,k) = (1./dp(i,k))* &
1625	(-mu(i,k)*(conu(i,k)-chat(i,k)) &
1626	-md(i,k)*(cond(i,k)-chat(i,k)) &
1627	)
1628	end if
1629	end do
1630	end do
1631
1632	! Initialize to zero everywhere, then scatter tendency back to full array
1633	dqdt(:,:,m) = 0._r8
1634
1635	do k = 1,pver
1636	do i = il1g,il2g
1637	ii = ideep(i)
1638	dqdt(ii,k,m) = dcondt(i,k)
1639	! Output apparent force on the mean flow from pressure gradient
1640	pguall(ii,k,m) = -pgu(i,k)
1641	pgdall(ii,k,m) = -pgd(i,k)
1642	icwu(ii,k,m) = conu(i,k)
1643	icwd(ii,k,m) = cond(i,k)
1644	end do
1645	end do
1646
1647	! Calculate momentum flux in units of mb*m/s2
1648
1649	do k = ktm,pver
1650	do i = il1g,il2g
1651	ii = ideep(i)
1652	mflux(i,k,m) = &
1653	-mu(i,k)* (conu(i,k)-chat(i,k)) &
1654	-md(i,k)* (cond(i,k)-chat(i,k))
1655	end do
1656	end do
1657
1658
1659	! Calculate winds at the end of the time step
1660
1661	do k = ktm,pver
1662	do i = il1g,il2g
1663	ii = ideep(i)
1664	km1 = max(1,k-1)
1665	kp1 = k+1
1666	windf(i,k,m) = const(i,k) - (mflux(i,kp1,m) - mflux(i,k,m)) * dt /dp(i,k)
1667
1668	end do
1669	end do
1670
1671	end if ! for domomtran
1672	end do
1673
1674	! Need to add an energy fix to account for the dissipation of kinetic energy
1675	! Formulation follows from Boville and Bretherton (2003)
1676	! formulation by PJR
1677
1678	do k = ktm,pver
1679	km1 = max(1,k-1)
1680	kp1 = min(pver,k+1)
1681	do i = il1g,il2g
1682
1683	ii = ideep(i)
1684
1685	! calculate the KE fluxes at top and bot of layer
1686	! based on a discrete approximation to b&b eq(35) F_KE = uF_u + vF_v at interface
1687	utop = (wind0(i,k,1)+wind0(i,km1,1))/2.
1688	vtop = (wind0(i,k,2)+wind0(i,km1,2))/2.
1689	ubot = (wind0(i,kp1,1)+wind0(i,k,1))/2.
1690	vbot = (wind0(i,kp1,2)+wind0(i,k,2))/2.
1691	fket = utopmflux(i,k,1) + vtopmflux(i,k,2) ! top of layer
1692	fkeb = ubotmflux(i,k+1,1) + vbotmflux(i,k+1,2) ! bot of layer
1693
1694	! divergence of these fluxes should give a conservative redistribution of KE
1695	ketend_cons = (fket-fkeb)/dp(i,k)
1696
1697	! tendency in kinetic energy resulting from the momentum transport
1698	ketend = ((windf(i,k,1)2 + windf(i,k,2)2) - (wind0(i,k,1)2 + wind0(i,k,2)2))*0.5/dt
1699
1700	! the difference should be the dissipation
1701	gset2 = ketend_cons - ketend
1702	gseten(i,k) = gset2
1703
1704	end do
1705
1706	end do
1707
1708	! Scatter dry static energy to full array
1709	do k = 1,pver
1710	do i = il1g,il2g
1711	ii = ideep(i)
1712	seten(ii,k) = gseten(i,k)
1713
1714	end do
1715	end do
1716
1717	return
1718	end subroutine momtran
1719
1720	#ifndef WRF_PORT
1721	!=========================================================================================
1722
1723	subroutine buoyan(lchnk ,ncol , &
1724	q ,t ,p ,z ,pf , &
1725	tp ,qstp ,tl ,rl ,cape , &
1726	pblt ,lcl ,lel ,lon ,mx , &
1727	rd ,grav ,cp ,msg , &
1728	tpert )
1729	!-----------------------------------------------------------------------
1730	!
1731	! Purpose:
1732	! <Say what the routine does>
1733	!
1734	! Method:
1735	! <Describe the algorithm(s) used in the routine.>
1736	! <Also include any applicable external references.>
1737	!
1738	! Author:
1739	! This is contributed code not fully standardized by the CCM core group.
1740	! The documentation has been enhanced to the degree that we are able.
1741	! Reviewed: P. Rasch, April 1996
1742	!
1743	!-----------------------------------------------------------------------
1744	implicit none
1745	!-----------------------------------------------------------------------
1746	!
1747	! input arguments
1748	!
1749	integer, intent(in) :: lchnk ! chunk identifier
1750	integer, intent(in) :: ncol ! number of atmospheric columns
1751
1752	real(r8), intent(in) :: q(pcols,pver) ! spec. humidity
1753	real(r8), intent(in) :: t(pcols,pver) ! temperature
1754	real(r8), intent(in) :: p(pcols,pver) ! pressure
1755	real(r8), intent(in) :: z(pcols,pver) ! height
1756	real(r8), intent(in) :: pf(pcols,pver+1) ! pressure at interfaces
1757	real(r8), intent(in) :: pblt(pcols) ! index of pbl depth
1758	real(r8), intent(in) :: tpert(pcols) ! perturbation temperature by pbl processes
1759
1760	!
1761	! output arguments
1762	!
1763	real(r8), intent(out) :: tp(pcols,pver) ! parcel temperature
1764	real(r8), intent(out) :: qstp(pcols,pver) ! saturation mixing ratio of parcel
1765	real(r8), intent(out) :: tl(pcols) ! parcel temperature at lcl
1766	real(r8), intent(out) :: cape(pcols) ! convective aval. pot. energy.
1767	integer lcl(pcols) !
1768	integer lel(pcols) !
1769	integer lon(pcols) ! level of onset of deep convection
1770	integer mx(pcols) ! level of max moist static energy
1771	!
1772	!--------------------------Local Variables------------------------------
1773	!
1774	real(r8) capeten(pcols,5) ! provisional value of cape
1775	real(r8) tv(pcols,pver) !
1776	real(r8) tpv(pcols,pver) !
1777	real(r8) buoy(pcols,pver)
1778
1779	real(r8) a1(pcols)
1780	real(r8) a2(pcols)
1781	real(r8) estp(pcols)
1782	real(r8) pl(pcols)
1783	real(r8) plexp(pcols)
1784	real(r8) hmax(pcols)
1785	real(r8) hmn(pcols)
1786	real(r8) y(pcols)
1787
1788	logical plge600(pcols)
1789	integer knt(pcols)
1790	integer lelten(pcols,5)
1791
1792	real(r8) cp
1793	real(r8) e
1794	real(r8) grav
1795
1796	integer i
1797	integer k
1798	integer msg
1799	integer n
1800
1801	real(r8) rd
1802	real(r8) rl
1803	#ifdef PERGRO
1804	real(r8) rhd
1805	#endif
1806	!
1807	!-----------------------------------------------------------------------
1808	!
1809	do n = 1,5
1810	do i = 1,ncol
1811	lelten(i,n) = pver
1812	capeten(i,n) = 0._r8
1813	end do
1814	end do
1815	!
1816	do i = 1,ncol
1817	lon(i) = pver
1818	knt(i) = 0
1819	lel(i) = pver
1820	mx(i) = lon(i)
1821	cape(i) = 0._r8
1822	hmax(i) = 0._r8
1823	end do
1824
1825	tp(:ncol,:) = t(:ncol,:)
1826	qstp(:ncol,:) = q(:ncol,:)
1827
1828	!!! RBN - Initialize tv and buoy for output.
1829	!!! tv=tv : tpv=tpv : qstp=q : buoy=0.
1830	tv(:ncol,:) = t(:ncol,:) (1._r8+1.608q(:ncol,:))/ (1._r8+q(:ncol,:))
1831	tpv(:ncol,:) = tv(:ncol,:)
1832	buoy(:ncol,:) = 0._r8
1833
1834	!
1835	! set "launching" level(mx) to be at maximum moist static energy.
1836	! search for this level stops at planetary boundary layer top.
1837	!
1838	#ifdef PERGRO
1839	do k = pver,msg + 1,-1
1840	do i = 1,ncol
1841	hmn(i) = cpt(i,k) + gravz(i,k) + rl*q(i,k)
1842	!
1843	! Reset max moist static energy level when relative difference exceeds 1.e-4
1844	!
1845	rhd = (hmn(i) - hmax(i))/(hmn(i) + hmax(i))
1846	if (k >= nint(pblt(i)) .and. k <= lon(i) .and. rhd > -1.e-4_r8) then
1847	hmax(i) = hmn(i)
1848	mx(i) = k
1849	end if
1850	end do
1851	end do
1852	#else
1853	do k = pver,msg + 1,-1
1854	do i = 1,ncol
1855	hmn(i) = cpt(i,k) + gravz(i,k) + rl*q(i,k)
1856	if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then
1857	hmax(i) = hmn(i)
1858	mx(i) = k
1859	end if
1860	end do
1861	end do
1862	#endif
1863	!
1864	do i = 1,ncol
1865	lcl(i) = mx(i)
1866	e = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
1867	tl(i) = 2840._r8/ (3.5_r8*log(t(i,mx(i)))-log(e)-4.805_r8) + 55._r8
1868	if (tl(i) < t(i,mx(i))) then
1869	plexp(i) = (1._r8/ (0.2854_r8* (1._r8-0.28_r8*q(i,mx(i)))))
1870	pl(i) = p(i,mx(i))* (tl(i)/t(i,mx(i)))**plexp(i)
1871	else
1872	tl(i) = t(i,mx(i))
1873	pl(i) = p(i,mx(i))
1874	end if
1875	end do
1876
1877	!
1878	! calculate lifting condensation level (lcl).
1879	!
1880	do k = pver,msg + 2,-1
1881	do i = 1,ncol
1882	if (k <= mx(i) .and. (p(i,k) > pl(i) .and. p(i,k-1) <= pl(i))) then
1883	lcl(i) = k - 1
1884	end if
1885	end do
1886	end do
1887	!
1888	! if lcl is above the nominal level of non-divergence (600 mbs),
1889	! no deep convection is permitted (ensuing calculations
1890	! skipped and cape retains initialized value of zero).
1891	!
1892	do i = 1,ncol
1893	plge600(i) = pl(i).ge.600._r8
1894	end do
1895	!
1896	! initialize parcel properties in sub-cloud layer below lcl.
1897	!
1898	do k = pver,msg + 1,-1
1899	do i=1,ncol
1900	if (k > lcl(i) .and. k <= mx(i) .and. plge600(i)) then
1901	tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
1902	qstp(i,k) = q(i,mx(i))
1903	tp(i,k) = t(i,mx(i))* (p(i,k)/p(i,mx(i)))*(0.2854_r8 (1._r8-0.28_r8*q(i,mx(i))))
1904	!
1905	! buoyancy is increased by 0.5 k as in tiedtke
1906	!
1907	!-jjh tpv (i,k)=tp(i,k)(1.+1.608q(i,mx(i)))/
1908	!-jjh 1 (1.+q(i,mx(i)))
1909	tpv(i,k) = (tp(i,k)+tpert(i))(1._r8+1.608_r8q(i,mx(i)))/ (1._r8+q(i,mx(i)))
1910	buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
1911	end if
1912	end do
1913	end do
1914
1915	!
1916	! define parcel properties at lcl (i.e. level immediately above pl).
1917	!
1918	do k = pver,msg + 1,-1
1919	do i=1,ncol
1920	if (k == lcl(i) .and. plge600(i)) then
1921	tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
1922	qstp(i,k) = q(i,mx(i))
1923	tp(i,k) = tl(i)* (p(i,k)/pl(i))*(0.2854_r8 (1._r8-0.28_r8*qstp(i,k)))
1924	! estp(i) =exp(a-b/tp(i,k))
1925	! use of different formulas for est has about 1 g/kg difference
1926	! in qs at t= 300k, and 0.02 g/kg at t=263k, with the formula
1927	! above giving larger qs.
1928	!
1929	estp(i) = c1exp((c2 (tp(i,k)-tfreez))/((tp(i,k)-tfreez)+c3))
1930
1931	qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i))
1932	a1(i) = cp / rl + qstp(i,k) * (1._r8+ qstp(i,k) / eps1) * rl * eps1 / &
1933	(rd * tp(i,k) ** 2)
1934	a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1qstp(i,k)) &
1935	(1._r8+qstp(i,k)/eps1)eps12rl*rl/ &
1936	(rd*2tp(i,k)*4)-qstp(i,k) &
1937	(1._r8+qstp(i,k)/eps1)2._r8eps1*rl/ &
1938	(rdtp(i,k)*3))
1939	a1(i) = 1._r8/a1(i)
1940	a2(i) = -a2(i)a1(i)*3
1941	y(i) = q(i,mx(i)) - qstp(i,k)
1942	tp(i,k) = tp(i,k) + a1(i)y(i) + a2(i)y(i)**2
1943	! estp(i) =exp(a-b/tp(i,k))
1944	estp(i) = c1exp((c2 (tp(i,k)-tfreez))/ ((tp(i,k)-tfreez)+c3))
1945
1946	qstp(i,k) = eps1*estp(i) / (p(i,k)-estp(i))
1947	!
1948	! buoyancy is increased by 0.5 k in cape calculation.
1949	! dec. 9, 1994
1950	!-jjh tpv(i,k) =tp(i,k)(1.+1.608qstp(i,k))/(1.+q(i,mx(i)))
1951	!
1952	tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k)) / (1._r8+q(i,mx(i)))
1953	buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
1954	end if
1955	end do
1956	end do
1957	!
1958	! main buoyancy calculation.
1959	!
1960	do k = pver - 1,msg + 1,-1
1961	do i=1,ncol
1962	if (k < lcl(i) .and. plge600(i)) then
1963	tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
1964	qstp(i,k) = qstp(i,k+1)
1965	tp(i,k) = tp(i,k+1)* (p(i,k)/p(i,k+1))*(0.2854_r8 (1._r8-0.28_r8*qstp(i,k)))
1966	! estp(i) = exp(a-b/tp(i,k))
1967	estp(i) = c1exp((c2 (tp(i,k)-tfreez))/((tp(i,k)-tfreez)+c3))
1968
1969	qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i))
1970	a1(i) = cp/rl + qstp(i,k)* (1._r8+qstp(i,k)/eps1)rleps1/ (rdtp(i,k)*2)
1971	a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1qstp(i,k)) &
1972	(1._r8+qstp(i,k)/eps1)eps12rl*rl/ &
1973	(rd*2tp(i,k)*4)-qstp(i,k) &
1974	(1._r8+qstp(i,k)/eps1)2._r8eps1*rl/ &
1975	(rdtp(i,k)*3))
1976	a1(i) = 1._r8/a1(i)
1977	a2(i) = -a2(i)a1(i)*3
1978	y(i) = qstp(i,k+1) - qstp(i,k)
1979	tp(i,k) = tp(i,k) + a1(i)y(i) + a2(i)y(i)**2
1980	! estp(i) =exp(a-b/tp(i,k))
1981	estp(i) = c1exp((c2 (tp(i,k)-tfreez))/ ((tp(i,k)-tfreez)+c3))
1982
1983	qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i))
1984	!-jjh tpv(i,k) =tp(i,k)(1.+1.608qstp(i,k))/
1985	!jt (1.+q(i,mx(i)))
1986	tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k))/(1._r8+q(i,mx(i)))
1987	buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add
1988	end if
1989	end do
1990	end do
1991
1992	!
1993	do k = msg + 2,pver
1994	do i = 1,ncol
1995	if (k < lcl(i) .and. plge600(i)) then
1996	if (buoy(i,k+1) > 0._r8 .and. buoy(i,k) <= 0._r8) then
1997	knt(i) = min(5,knt(i) + 1)
1998	lelten(i,knt(i)) = k
1999	end if
2000	end if
2001	end do
2002	end do
2003	!
2004	! calculate convective available potential energy (cape).
2005	!
2006	do n = 1,5
2007	do k = msg + 1,pver
2008	do i = 1,ncol
2009	if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then
2010	capeten(i,n) = capeten(i,n) + rdbuoy(i,k)log(pf(i,k+1)/pf(i,k))
2011	end if
2012	end do
2013	end do
2014	end do
2015	!
2016	! find maximum cape from all possible tentative capes from
2017	! one sounding,
2018	! and use it as the final cape, april 26, 1995
2019	!
2020	do n = 1,5
2021	do i = 1,ncol
2022	if (capeten(i,n) > cape(i)) then
2023	cape(i) = capeten(i,n)
2024	lel(i) = lelten(i,n)
2025	end if
2026	end do
2027	end do
2028	!
2029	! put lower bound on cape for diagnostic purposes.
2030	!
2031	do i = 1,ncol
2032	cape(i) = max(cape(i), 0._r8)
2033	end do
2034	!
2035	return
2036	end subroutine buoyan
2037
2038	#endif
2039
2040	subroutine cldprp(lchnk , &
2041	q ,t ,u ,v ,p , &
2042	z ,s ,mu ,eu ,du , &
2043	md ,ed ,sd ,qd ,mc , &
2044	qu ,su ,zf ,qst ,hmn , &
2045	hsat ,shat ,ql , &
2046	cmeg ,jb ,lel ,jt ,jlcl , &
2047	mx ,j0 ,jd ,rl ,il2g , &
2048	rd ,grav ,cp ,msg , &
2049	pflx ,evp ,cu ,rprd ,limcnv ,landfrac)
2050	!-----------------------------------------------------------------------
2051	!
2052	! Purpose:
2053	! <Say what the routine does>
2054	!
2055	! Method:
2056	! may 09/91 - guang jun zhang, m.lazare, n.mcfarlane.
2057	! original version cldprop.
2058	!
2059	! Author: See above, modified by P. Rasch
2060	! This is contributed code not fully standardized by the CCM core group.
2061	!
2062	! this code is very much rougher than virtually anything else in the CCM
2063	! there are debug statements left strewn about and code segments disabled
2064	! these are to facilitate future development. We expect to release a
2065	! cleaner code in a future release
2066	!
2067	! the documentation has been enhanced to the degree that we are able
2068	!
2069	!-----------------------------------------------------------------------
2070
2071	implicit none
2072
2073	!------------------------------------------------------------------------------
2074	!
2075	! Input arguments
2076	!
2077	integer, intent(in) :: lchnk ! chunk identifier
2078
2079	real(r8), intent(in) :: q(pcols,pver) ! spec. humidity of env
2080	real(r8), intent(in) :: t(pcols,pver) ! temp of env
2081	real(r8), intent(in) :: p(pcols,pver) ! pressure of env
2082	real(r8), intent(in) :: z(pcols,pver) ! height of env
2083	real(r8), intent(in) :: s(pcols,pver) ! normalized dry static energy of env
2084	real(r8), intent(in) :: zf(pcols,pverp) ! height of interfaces
2085	real(r8), intent(in) :: u(pcols,pver) ! zonal velocity of env
2086	real(r8), intent(in) :: v(pcols,pver) ! merid. velocity of env
2087	real(r8), intent(in) :: landfrac(pcols) ! RBN Landfrac
2088
2089	integer, intent(in) :: jb(pcols) ! updraft base level
2090	integer, intent(in) :: lel(pcols) ! updraft launch level
2091	integer, intent(out) :: jt(pcols) ! updraft plume top
2092	integer, intent(out) :: jlcl(pcols) ! updraft lifting cond level
2093	integer, intent(in) :: mx(pcols) ! updraft base level (same is jb)
2094	integer, intent(out) :: j0(pcols) ! level where updraft begins detraining
2095	integer, intent(out) :: jd(pcols) ! level of downdraft
2096	integer, intent(in) :: limcnv ! convection limiting level
2097	integer, intent(in) :: il2g !CORE GROUP REMOVE
2098	integer, intent(in) :: msg ! missing moisture vals (always 0)
2099	real(r8), intent(in) :: rl ! latent heat of vap
2100	real(r8), intent(in) :: shat(pcols,pver) ! interface values of dry stat energy
2101	!
2102	! output
2103	!
2104	real(r8), intent(out) :: rprd(pcols,pver) ! rate of production of precip at that layer
2105	real(r8), intent(out) :: du(pcols,pver) ! detrainement rate of updraft
2106	real(r8), intent(out) :: ed(pcols,pver) ! entrainment rate of downdraft
2107	real(r8), intent(out) :: eu(pcols,pver) ! entrainment rate of updraft
2108	real(r8), intent(out) :: hmn(pcols,pver) ! moist stat energy of env
2109	real(r8), intent(out) :: hsat(pcols,pver) ! sat moist stat energy of env
2110	real(r8), intent(out) :: mc(pcols,pver) ! net mass flux
2111	real(r8), intent(out) :: md(pcols,pver) ! downdraft mass flux
2112	real(r8), intent(out) :: mu(pcols,pver) ! updraft mass flux
2113	real(r8), intent(out) :: pflx(pcols,pverp) ! precipitation flux thru layer
2114	real(r8), intent(out) :: qd(pcols,pver) ! spec humidity of downdraft
2115	real(r8), intent(out) :: ql(pcols,pver) ! liq water of updraft
2116	real(r8), intent(out) :: qst(pcols,pver) ! saturation spec humidity of env.
2117	real(r8), intent(out) :: qu(pcols,pver) ! spec hum of updraft
2118	real(r8), intent(out) :: sd(pcols,pver) ! normalized dry stat energy of downdraft
2119	real(r8), intent(out) :: su(pcols,pver) ! normalized dry stat energy of updraft
2120
2121
2122	real(r8) rd ! gas constant for dry air
2123	real(r8) grav ! gravity
2124	real(r8) cp ! heat capacity of dry air
2125
2126	!
2127	! Local workspace
2128	!
2129	real(r8) gamma(pcols,pver)
2130	real(r8) dz(pcols,pver)
2131	real(r8) iprm(pcols,pver)
2132	real(r8) hu(pcols,pver)
2133	real(r8) hd(pcols,pver)
2134	real(r8) eps(pcols,pver)
2135	real(r8) f(pcols,pver)
2136	real(r8) k1(pcols,pver)
2137	real(r8) i2(pcols,pver)
2138	real(r8) ihat(pcols,pver)
2139	real(r8) i3(pcols,pver)
2140	real(r8) idag(pcols,pver)
2141	real(r8) i4(pcols,pver)
2142	real(r8) qsthat(pcols,pver)
2143	real(r8) hsthat(pcols,pver)
2144	real(r8) gamhat(pcols,pver)
2145	real(r8) cu(pcols,pver)
2146	real(r8) evp(pcols,pver)
2147	real(r8) cmeg(pcols,pver)
2148	real(r8) qds(pcols,pver)
2149	! RBN For c0mask
2150	real(r8) c0mask(pcols)
2151	real(r8) hmin(pcols)
2152	real(r8) expdif(pcols)
2153	real(r8) expnum(pcols)
2154	real(r8) ftemp(pcols)
2155	real(r8) eps0(pcols)
2156	real(r8) rmue(pcols)
2157	real(r8) zuef(pcols)
2158	real(r8) zdef(pcols)
2159	real(r8) epsm(pcols)
2160	real(r8) ratmjb(pcols)
2161	real(r8) est(pcols)
2162	real(r8) totpcp(pcols)
2163	real(r8) totevp(pcols)
2164	real(r8) alfa(pcols)
2165	real(r8) ql1
2166	real(r8) tu
2167	real(r8) estu
2168	real(r8) qstu
2169
2170	real(r8) small
2171	real(r8) mdt
2172
2173	integer khighest
2174	integer klowest
2175	integer kount
2176	integer i,k
2177
2178	logical doit(pcols)
2179	logical done(pcols)
2180	!
2181	!------------------------------------------------------------------------------
2182	!
2183	do i = 1,il2g
2184	ftemp(i) = 0._r8
2185	expnum(i) = 0._r8
2186	expdif(i) = 0._r8
2187	c0mask(i) = c0_ocn * (1._r8-landfrac(i)) + c0_lnd * landfrac(i)
2188	end do
2189	!
2190	!jr Change from msg+1 to 1 to prevent blowup
2191	!
2192	do k = 1,pver
2193	do i = 1,il2g
2194	dz(i,k) = zf(i,k) - zf(i,k+1)
2195	end do
2196	end do
2197
2198	!
2199	! initialize many output and work variables to zero
2200	!
2201	pflx(:il2g,1) = 0
2202
2203	do k = 1,pver
2204	do i = 1,il2g
2205	k1(i,k) = 0._r8
2206	i2(i,k) = 0._r8
2207	i3(i,k) = 0._r8
2208	i4(i,k) = 0._r8
2209	mu(i,k) = 0._r8
2210	f(i,k) = 0._r8
2211	eps(i,k) = 0._r8
2212	eu(i,k) = 0._r8
2213	du(i,k) = 0._r8
2214	ql(i,k) = 0._r8
2215	cu(i,k) = 0._r8
2216	evp(i,k) = 0._r8
2217	cmeg(i,k) = 0._r8
2218	qds(i,k) = q(i,k)
2219	md(i,k) = 0._r8
2220	ed(i,k) = 0._r8
2221	sd(i,k) = s(i,k)
2222	qd(i,k) = q(i,k)
2223	mc(i,k) = 0._r8
2224	qu(i,k) = q(i,k)
2225	su(i,k) = s(i,k)
2226	! est(i)=exp(a-b/t(i,k))
2227	est(i) = c1exp((c2 (t(i,k)-tfreez))/((t(i,k)-tfreez)+c3))
2228	!++bee
2229	if ( p(i,k)-est(i) > 0._r8 ) then
2230	qst(i,k) = eps1*est(i)/ (p(i,k)-est(i))
2231	else
2232	qst(i,k) = 1.0_r8
2233	end if
2234	!--bee
2235	gamma(i,k) = qst(i,k)(1._r8 + qst(i,k)/eps1)eps1rl/(rdt(i,k)*2)rl/cp
2236	hmn(i,k) = cpt(i,k) + gravz(i,k) + rl*q(i,k)
2237	hsat(i,k) = cpt(i,k) + gravz(i,k) + rl*qst(i,k)
2238	hu(i,k) = hmn(i,k)
2239	hd(i,k) = hmn(i,k)
2240	rprd(i,k) = 0._r8
2241	end do
2242	end do
2243	!
2244	!jr Set to zero things which make this routine blow up
2245	!
2246	do k=1,msg
2247	do i=1,il2g
2248	rprd(i,k) = 0._r8
2249	end do
2250	end do
2251	!
2252	! interpolate the layer values of qst, hsat and gamma to
2253	! layer interfaces
2254	!
2255	do i = 1,il2g
2256	hsthat(i,msg+1) = hsat(i,msg+1)
2257	qsthat(i,msg+1) = qst(i,msg+1)
2258	gamhat(i,msg+1) = gamma(i,msg+1)
2259	totpcp(i) = 0._r8
2260	totevp(i) = 0._r8
2261	end do
2262	do k = msg + 2,pver
2263	do i = 1,il2g
2264	if (abs(qst(i,k-1)-qst(i,k)) > 1.E-6_r8) then
2265	qsthat(i,k) = log(qst(i,k-1)/qst(i,k))qst(i,k-1)qst(i,k)/ (qst(i,k-1)-qst(i,k))
2266	else
2267	qsthat(i,k) = qst(i,k)
2268	end if
2269	hsthat(i,k) = cpshat(i,k) + rlqsthat(i,k)
2270	if (abs(gamma(i,k-1)-gamma(i,k)) > 1.E-6_r8) then
2271	gamhat(i,k) = log(gamma(i,k-1)/gamma(i,k))gamma(i,k-1)gamma(i,k)/ &
2272	(gamma(i,k-1)-gamma(i,k))
2273	else
2274	gamhat(i,k) = gamma(i,k)
2275	end if
2276	end do
2277	end do
2278	!
2279	! initialize cloud top to highest plume top.
2280	!jr changed hard-wired 4 to limcnv+1 (not to exceed pver)
2281	!
2282	jt(:) = pver
2283	do i = 1,il2g
2284	jt(i) = max(lel(i),limcnv+1)
2285	jt(i) = min(jt(i),pver)
2286	jd(i) = pver
2287	jlcl(i) = lel(i)
2288	hmin(i) = 1.E6_r8
2289	end do
2290	!
2291	! find the level of minimum hsat, where detrainment starts
2292	!
2293	do k = msg + 1,pver
2294	do i = 1,il2g
2295	if (hsat(i,k) <= hmin(i) .and. k >= jt(i) .and. k <= jb(i)) then
2296	hmin(i) = hsat(i,k)
2297	j0(i) = k
2298	end if
2299	end do
2300	end do
2301	do i = 1,il2g
2302	j0(i) = min(j0(i),jb(i)-2)
2303	j0(i) = max(j0(i),jt(i)+2)
2304	!
2305	! Fix from Guang Zhang to address out of bounds array reference
2306	!
2307	j0(i) = min(j0(i),pver)
2308	end do
2309	!
2310	! Initialize certain arrays inside cloud
2311	!
2312	do k = msg + 1,pver
2313	do i = 1,il2g
2314	if (k >= jt(i) .and. k <= jb(i)) then
2315	hu(i,k) = hmn(i,mx(i)) + cp*tiedke_add
2316	su(i,k) = s(i,mx(i)) + tiedke_add
2317	end if
2318	end do
2319	end do
2320	!
2321	! *********************************************************
2322	! compute taylor series for approximate eps(z) below
2323	! *********************************************************
2324	!
2325	do k = pver - 1,msg + 1,-1
2326	do i = 1,il2g
2327	if (k < jb(i) .and. k >= jt(i)) then
2328	k1(i,k) = k1(i,k+1) + (hmn(i,mx(i))-hmn(i,k))*dz(i,k)
2329	ihat(i,k) = 0.5_r8* (k1(i,k+1)+k1(i,k))
2330	i2(i,k) = i2(i,k+1) + ihat(i,k)*dz(i,k)
2331	idag(i,k) = 0.5_r8* (i2(i,k+1)+i2(i,k))
2332	i3(i,k) = i3(i,k+1) + idag(i,k)*dz(i,k)
2333	iprm(i,k) = 0.5_r8* (i3(i,k+1)+i3(i,k))
2334	i4(i,k) = i4(i,k+1) + iprm(i,k)*dz(i,k)
2335	end if
2336	end do
2337	end do
2338	!
2339	! re-initialize hmin array for ensuing calculation.
2340	!
2341	do i = 1,il2g
2342	hmin(i) = 1.E6_r8
2343	end do
2344	do k = msg + 1,pver
2345	do i = 1,il2g
2346	if (k >= j0(i) .and. k <= jb(i) .and. hmn(i,k) <= hmin(i)) then
2347	hmin(i) = hmn(i,k)
2348	expdif(i) = hmn(i,mx(i)) - hmin(i)
2349	end if
2350	end do
2351	end do
2352	!
2353	! *********************************************************
2354	! compute approximate eps(z) using above taylor series
2355	! *********************************************************
2356	!
2357	do k = msg + 2,pver
2358	do i = 1,il2g
2359	expnum(i) = 0._r8
2360	ftemp(i) = 0._r8
2361	if (k < jt(i) .or. k >= jb(i)) then
2362	k1(i,k) = 0._r8
2363	expnum(i) = 0._r8
2364	else
2365	expnum(i) = hmn(i,mx(i)) - (hsat(i,k-1)*(zf(i,k)-z(i,k)) + &
2366	hsat(i,k)* (z(i,k-1)-zf(i,k)))/(z(i,k-1)-z(i,k))
2367	end if
2368	if ((expdif(i) > 100._r8 .and. expnum(i) > 0._r8) .and. &
2369	k1(i,k) > expnum(i)*dz(i,k)) then
2370	ftemp(i) = expnum(i)/k1(i,k)
2371	f(i,k) = ftemp(i) + i2(i,k)/k1(i,k)ftemp(i)*2 + &
2372	(2._r8i2(i,k)2-k1(i,k)i3(i,k))/k1(i,k)*2 &
2373	ftemp(i)*3 + (-5._r8k1(i,k)i2(i,k)i3(i,k)+ &
2374	5._r8i2(i,k)3+k1(i,k)2i4(i,k))/ &
2375	k1(i,k)*3ftemp(i)**4
2376	f(i,k) = max(f(i,k),0._r8)
2377	f(i,k) = min(f(i,k),0.0002_r8)
2378	end if
2379	end do
2380	end do
2381	do i = 1,il2g
2382	if (j0(i) < jb(i)) then
2383	if (f(i,j0(i)) < 1.E-6_r8 .and. f(i,j0(i)+1) > f(i,j0(i))) j0(i) = j0(i) + 1
2384	end if
2385	end do
2386	do k = msg + 2,pver
2387	do i = 1,il2g
2388	if (k >= jt(i) .and. k <= j0(i)) then
2389	f(i,k) = max(f(i,k),f(i,k-1))
2390	end if
2391	end do
2392	end do
2393	do i = 1,il2g
2394	eps0(i) = f(i,j0(i))
2395	eps(i,jb(i)) = eps0(i)
2396	end do
2397	!
2398	! This is set to match the Rasch and Kristjansson paper
2399	!
2400	do k = pver,msg + 1,-1
2401	do i = 1,il2g
2402	if (k >= j0(i) .and. k <= jb(i)) then
2403	eps(i,k) = f(i,j0(i))
2404	end if
2405	end do
2406	end do
2407	do k = pver,msg + 1,-1
2408	do i = 1,il2g
2409	if (k < j0(i) .and. k >= jt(i)) eps(i,k) = f(i,k)
2410	end do
2411	end do
2412	!
2413	! specify the updraft mass flux mu, entrainment eu, detrainment du
2414	! and moist static energy hu.
2415	! here and below mu, eu,du, md and ed are all normalized by mb
2416	!
2417	do i = 1,il2g
2418	if (eps0(i) > 0._r8) then
2419	mu(i,jb(i)) = 1._r8
2420	eu(i,jb(i)) = mu(i,jb(i))/dz(i,jb(i))
2421	end if
2422	end do
2423	do k = pver,msg + 1,-1
2424	do i = 1,il2g
2425	if (eps0(i) > 0._r8 .and. (k >= jt(i) .and. k < jb(i))) then
2426	zuef(i) = zf(i,k) - zf(i,jb(i))
2427	rmue(i) = (1._r8/eps0(i))* (exp(eps(i,k+1)*zuef(i))-1._r8)/zuef(i)
2428	mu(i,k) = (1._r8/eps0(i))* (exp(eps(i,k )*zuef(i))-1._r8)/zuef(i)
2429	eu(i,k) = (rmue(i)-mu(i,k+1))/dz(i,k)
2430	du(i,k) = (rmue(i)-mu(i,k))/dz(i,k)
2431	end if
2432	end do
2433	end do
2434	!
2435	khighest = pverp
2436	klowest = 1
2437	do i=1,il2g
2438	khighest = min(khighest,lel(i))
2439	klowest = max(klowest,jb(i))
2440	end do
2441	do k = klowest-1,khighest,-1
2442	!cdir$ ivdep
2443	do i = 1,il2g
2444	if (k <= jb(i)-1 .and. k >= lel(i) .and. eps0(i) > 0._r8) then
2445	if (mu(i,k) < 0.01_r8) then
2446	hu(i,k) = hu(i,jb(i))
2447	mu(i,k) = 0._r8
2448	eu(i,k) = 0._r8
2449	du(i,k) = mu(i,k+1)/dz(i,k)
2450	else
2451	hu(i,k) = mu(i,k+1)/mu(i,k)*hu(i,k+1) + &
2452	dz(i,k)/mu(i,k)* (eu(i,k)hmn(i,k)- du(i,k)hsat(i,k))
2453	end if
2454	end if
2455	end do
2456	end do
2457	!
2458	! reset cloud top index beginning from two layers above the
2459	! cloud base (i.e. if cloud is only one layer thick, top is not reset
2460	!
2461	do i=1,il2g
2462	doit(i) = .true.
2463	end do
2464	do k=klowest-2,khighest-1,-1
2465	do i=1,il2g
2466	if (doit(i) .and. k <= jb(i)-2 .and. k >= lel(i)-1) then
2467	if (hu(i,k) <= hsthat(i,k) .and. hu(i,k+1) > hsthat(i,k+1) &
2468	.and. mu(i,k) >= 0.02_r8) then
2469	if (hu(i,k)-hsthat(i,k) < -2000._r8) then
2470	jt(i) = k + 1
2471	doit(i) = .false.
2472	else
2473	jt(i) = k
2474	if (eps0(i) <= 0._r8) doit(i) = .false.
2475	end if
2476	else if (hu(i,k) > hu(i,jb(i)) .or. mu(i,k) < 0.01_r8) then
2477	jt(i) = k + 1
2478	doit(i) = .false.
2479	end if
2480	end if
2481	end do
2482	end do
2483	do k = pver,msg + 1,-1
2484	!cdir$ ivdep
2485	do i = 1,il2g
2486	if (k >= lel(i) .and. k <= jt(i) .and. eps0(i) > 0._r8) then
2487	mu(i,k) = 0._r8
2488	eu(i,k) = 0._r8
2489	du(i,k) = 0._r8
2490	hu(i,k) = hu(i,jb(i))
2491	end if
2492	if (k == jt(i) .and. eps0(i) > 0._r8) then
2493	du(i,k) = mu(i,k+1)/dz(i,k)
2494	eu(i,k) = 0._r8
2495	mu(i,k) = 0._r8
2496	end if
2497	end do
2498	end do
2499	!
2500	! specify downdraft properties (no downdrafts if jd.ge.jb).
2501	! scale down downward mass flux profile so that net flux
2502	! (up-down) at cloud base in not negative.
2503	!
2504	do i = 1,il2g
2505	!
2506	! in normal downdraft strength run alfa=0.2. In test4 alfa=0.1
2507	!
2508	alfa(i) = 0.1_r8
2509	jt(i) = min(jt(i),jb(i)-1)
2510	jd(i) = max(j0(i),jt(i)+1)
2511	jd(i) = min(jd(i),jb(i))
2512	hd(i,jd(i)) = hmn(i,jd(i)-1)
2513	if (jd(i) < jb(i) .and. eps0(i) > 0._r8) then
2514	epsm(i) = eps0(i)
2515	md(i,jd(i)) = -alfa(i)*epsm(i)/eps0(i)
2516	end if
2517	end do
2518	do k = msg + 1,pver
2519	do i = 1,il2g
2520	if ((k > jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8) then
2521	zdef(i) = zf(i,jd(i)) - zf(i,k)
2522	md(i,k) = -alfa(i)/ (2._r8eps0(i))(exp(2._r8epsm(i)zdef(i))-1._r8)/zdef(i)
2523	end if
2524	end do
2525	end do
2526	do k = msg + 1,pver
2527	do i = 1,il2g
2528	if ((k >= jt(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8 .and. jd(i) < jb(i)) then
2529	ratmjb(i) = min(abs(mu(i,jb(i))/md(i,jb(i))),1._r8)
2530	md(i,k) = md(i,k)*ratmjb(i)
2531	end if
2532	end do
2533	end do
2534
2535	small = 1.e-20_r8
2536	do k = msg + 1,pver
2537	do i = 1,il2g
2538	if ((k >= jt(i) .and. k <= pver) .and. eps0(i) > 0._r8) then
2539	ed(i,k-1) = (md(i,k-1)-md(i,k))/dz(i,k-1)
2540	mdt = min(md(i,k),-small)
2541	hd(i,k) = (md(i,k-1)hd(i,k-1) - dz(i,k-1)ed(i,k-1)*hmn(i,k-1))/mdt
2542	end if
2543	end do
2544	end do
2545	!
2546	! calculate updraft and downdraft properties.
2547	!
2548	do k = msg + 2,pver
2549	do i = 1,il2g
2550	if ((k >= jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8 .and. jd(i) < jb(i)) then
2551	qds(i,k) = qsthat(i,k) + gamhat(i,k)*(hd(i,k)-hsthat(i,k))/ &
2552	(rl*(1._r8 + gamhat(i,k)))
2553	end if
2554	end do
2555	end do
2556	!
2557	do i = 1,il2g
2558	done(i) = .false.
2559	end do
2560	kount = 0
2561	do k = pver,msg + 2,-1
2562	do i = 1,il2g
2563	if (( .not. done(i) .and. k > jt(i) .and. k < jb(i)) .and. eps0(i) > 0._r8) then
2564	su(i,k) = mu(i,k+1)/mu(i,k)*su(i,k+1) + &
2565	dz(i,k)/mu(i,k)* (eu(i,k)-du(i,k))*s(i,k)
2566	qu(i,k) = mu(i,k+1)/mu(i,k)qu(i,k+1) + dz(i,k)/mu(i,k) (eu(i,k)*q(i,k)- &
2567	du(i,k)*qst(i,k))
2568	tu = su(i,k) - grav/cp*zf(i,k)
2569	estu = c1exp((c2 (tu-tfreez))/ ((tu-tfreez)+c3))
2570	qstu = eps1*estu/ ((p(i,k)+p(i,k-1))/2._r8-estu)
2571	if (qu(i,k) >= qstu) then
2572	jlcl(i) = k
2573	kount = kount + 1
2574	done(i) = .true.
2575	end if
2576	end if
2577	end do
2578	if (kount >= il2g) goto 690
2579	end do
2580	690 continue
2581	do k = msg + 2,pver
2582	do i = 1,il2g
2583	if (k == jb(i) .and. eps0(i) > 0._r8) then
2584	qu(i,k) = q(i,mx(i))
2585	su(i,k) = (hu(i,k)-rl*qu(i,k))/cp
2586	end if
2587	if ((k > jt(i) .and. k <= jlcl(i)) .and. eps0(i) > 0._r8) then
2588	su(i,k) = shat(i,k) + (hu(i,k)-hsthat(i,k))/(cp* (1._r8+gamhat(i,k)))
2589	qu(i,k) = qsthat(i,k) + gamhat(i,k)*(hu(i,k)-hsthat(i,k))/ &
2590	(rl* (1._r8+gamhat(i,k)))
2591	end if
2592	end do
2593	end do
2594
2595	! compute condensation in updraft
2596	do k = pver,msg + 2,-1
2597	do i = 1,il2g
2598	if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._r8) then
2599	cu(i,k) = ((mu(i,k)su(i,k)-mu(i,k+1)su(i,k+1))/ &
2600	dz(i,k)- (eu(i,k)-du(i,k))*s(i,k))/(rl/cp)
2601	if (k == jt(i)) cu(i,k) = 0._r8
2602	cu(i,k) = max(0._r8,cu(i,k))
2603	end if
2604	end do
2605	end do
2606
2607	! compute condensed liquid, rain production rate
2608	! accumulate total precipitation (condensation - detrainment of liquid)
2609	! Note ql1 = ql(k) + rprd(k)*dz(k)/mu(k)
2610	! The differencing is somewhat strange (e.g. du(i,k)*ql(i,k+1)) but is
2611	! consistently applied.
2612	! mu, ql are interface quantities
2613	! cu, du, eu, rprd are midpoint quantites
2614	do k = pver,msg + 2,-1
2615	do i = 1,il2g
2616	rprd(i,k) = 0._r8
2617	if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._r8 .and. mu(i,k) >= 0.0_r8) then
2618	if (mu(i,k) > 0._r8) then
2619	ql1 = 1._r8/mu(i,k)* (mu(i,k+1)*ql(i,k+1)- &
2620	dz(i,k)du(i,k)ql(i,k+1)+dz(i,k)*cu(i,k))
2621	ql(i,k) = ql1/ (1._r8+dz(i,k)*c0mask(i))
2622	else
2623	ql(i,k) = 0._r8
2624	end if
2625	totpcp(i) = totpcp(i) + dz(i,k)(cu(i,k)-du(i,k)ql(i,k+1))
2626	rprd(i,k) = c0mask(i)mu(i,k)ql(i,k)
2627	end if
2628	end do
2629	end do
2630	!
2631	do i = 1,il2g
2632	qd(i,jd(i)) = qds(i,jd(i))
2633	sd(i,jd(i)) = (hd(i,jd(i)) - rl*qd(i,jd(i)))/cp
2634	end do
2635	!
2636	do k = msg + 2,pver
2637	do i = 1,il2g
2638	if (k >= jd(i) .and. k < jb(i) .and. eps0(i) > 0._r8) then
2639	qd(i,k+1) = qds(i,k+1)
2640	evp(i,k) = -ed(i,k)q(i,k) + (md(i,k)qd(i,k)-md(i,k+1)*qd(i,k+1))/dz(i,k)
2641	evp(i,k) = max(evp(i,k),0._r8)
2642	mdt = min(md(i,k+1),-small)
2643	sd(i,k+1) = ((rl/cpevp(i,k)-ed(i,k)s(i,k))dz(i,k) + md(i,k)sd(i,k))/mdt
2644	totevp(i) = totevp(i) - dz(i,k)ed(i,k)q(i,k)
2645	end if
2646	end do
2647	end do
2648	do i = 1,il2g
2649	!guang totevp(i) = totevp(i) + md(i,jd(i))q(i,jd(i)-1) -
2650	totevp(i) = totevp(i) + md(i,jd(i))qd(i,jd(i)) - md(i,jb(i))qd(i,jb(i))
2651	end do
2652	!!$ if (.true.) then
2653	if (.false.) then
2654	do i = 1,il2g
2655	k = jb(i)
2656	if (eps0(i) > 0._r8) then
2657	evp(i,k) = -ed(i,k)q(i,k) + (md(i,k)qd(i,k))/dz(i,k)
2658	evp(i,k) = max(evp(i,k),0._r8)
2659	totevp(i) = totevp(i) - dz(i,k)ed(i,k)q(i,k)
2660	end if
2661	end do
2662	endif
2663
2664	do i = 1,il2g
2665	totpcp(i) = max(totpcp(i),0._r8)
2666	totevp(i) = max(totevp(i),0._r8)
2667	end do
2668	!
2669	do k = msg + 2,pver
2670	do i = 1,il2g
2671	if (totevp(i) > 0._r8 .and. totpcp(i) > 0._r8) then
2672	md(i,k) = md (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i)))
2673	ed(i,k) = ed (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i)))
2674	evp(i,k) = evp(i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i)))
2675	else
2676	md(i,k) = 0._r8
2677	ed(i,k) = 0._r8
2678	evp(i,k) = 0._r8
2679	end if
2680	! cmeg is the cloud water condensed - rain water evaporated
2681	! rprd is the cloud water converted to rain - (rain evaporated)
2682	cmeg(i,k) = cu(i,k) - evp(i,k)
2683	rprd(i,k) = rprd(i,k)-evp(i,k)
2684	end do
2685	end do
2686
2687	! compute the net precipitation flux across interfaces
2688	pflx(:il2g,1) = 0._r8
2689	do k = 2,pverp
2690	do i = 1,il2g
2691	pflx(i,k) = pflx(i,k-1) + rprd(i,k-1)*dz(i,k-1)
2692	end do
2693	end do
2694	!
2695	do k = msg + 1,pver
2696	do i = 1,il2g
2697	mc(i,k) = mu(i,k) + md(i,k)
2698	end do
2699	end do
2700	!
2701	return
2702	end subroutine cldprp
2703
2704	subroutine closure(lchnk , &
2705	q ,t ,p ,z ,s , &
2706	tp ,qs ,qu ,su ,mc , &
2707	du ,mu ,md ,qd ,sd , &
2708	qhat ,shat ,dp ,qstp ,zf , &
2709	ql ,dsubcld ,mb ,cape ,tl , &
2710	lcl ,lel ,jt ,mx ,il1g , &
2711	il2g ,rd ,grav ,cp ,rl , &
2712	msg ,capelmt )
2713	!-----------------------------------------------------------------------
2714	!
2715	! Purpose:
2716	! <Say what the routine does>
2717	!
2718	! Method:
2719	! <Describe the algorithm(s) used in the routine.>
2720	! <Also include any applicable external references.>
2721	!
2722	! Author: G. Zhang and collaborators. CCM contact:P. Rasch
2723	! This is contributed code not fully standardized by the CCM core group.
2724	!
2725	! this code is very much rougher than virtually anything else in the CCM
2726	! We expect to release cleaner code in a future release
2727	!
2728	! the documentation has been enhanced to the degree that we are able
2729	!
2730	!-----------------------------------------------------------------------
2731	#ifndef WRF_PORT
2732	use dycore, only: dycore_is, get_resolution
2733	#endif
2734
2735	implicit none
2736
2737	!
2738	!-----------------------------Arguments---------------------------------
2739	!
2740	integer, intent(in) :: lchnk ! chunk identifier
2741
2742	real(r8), intent(inout) :: q(pcols,pver) ! spec humidity
2743	real(r8), intent(inout) :: t(pcols,pver) ! temperature
2744	real(r8), intent(inout) :: p(pcols,pver) ! pressure (mb)
2745	real(r8), intent(inout) :: mb(pcols) ! cloud base mass flux
2746	real(r8), intent(in) :: z(pcols,pver) ! height (m)
2747	real(r8), intent(in) :: s(pcols,pver) ! normalized dry static energy
2748	real(r8), intent(in) :: tp(pcols,pver) ! parcel temp
2749	real(r8), intent(in) :: qs(pcols,pver) ! sat spec humidity
2750	real(r8), intent(in) :: qu(pcols,pver) ! updraft spec. humidity
2751	real(r8), intent(in) :: su(pcols,pver) ! normalized dry stat energy of updraft
2752	real(r8), intent(in) :: mc(pcols,pver) ! net convective mass flux
2753	real(r8), intent(in) :: du(pcols,pver) ! detrainment from updraft
2754	real(r8), intent(in) :: mu(pcols,pver) ! mass flux of updraft
2755	real(r8), intent(in) :: md(pcols,pver) ! mass flux of downdraft
2756	real(r8), intent(in) :: qd(pcols,pver) ! spec. humidity of downdraft
2757	real(r8), intent(in) :: sd(pcols,pver) ! dry static energy of downdraft
2758	real(r8), intent(in) :: qhat(pcols,pver) ! environment spec humidity at interfaces
2759	real(r8), intent(in) :: shat(pcols,pver) ! env. normalized dry static energy at intrfcs
2760	real(r8), intent(in) :: dp(pcols,pver) ! pressure thickness of layers
2761	real(r8), intent(in) :: qstp(pcols,pver) ! spec humidity of parcel
2762	real(r8), intent(in) :: zf(pcols,pver+1) ! height of interface levels
2763	real(r8), intent(in) :: ql(pcols,pver) ! liquid water mixing ratio
2764
2765	real(r8), intent(in) :: cape(pcols) ! available pot. energy of column
2766	real(r8), intent(in) :: tl(pcols)
2767	real(r8), intent(in) :: dsubcld(pcols) ! thickness of subcloud layer
2768
2769	integer, intent(in) :: lcl(pcols) ! index of lcl
2770	integer, intent(in) :: lel(pcols) ! index of launch leve
2771	integer, intent(in) :: jt(pcols) ! top of updraft
2772	integer, intent(in) :: mx(pcols) ! base of updraft
2773	!
2774	!--------------------------Local variables------------------------------
2775	!
2776	real(r8) dtpdt(pcols,pver)
2777	real(r8) dqsdtp(pcols,pver)
2778	real(r8) dtmdt(pcols,pver)
2779	real(r8) dqmdt(pcols,pver)
2780	real(r8) dboydt(pcols,pver)
2781	real(r8) thetavp(pcols,pver)
2782	real(r8) thetavm(pcols,pver)
2783
2784	real(r8) dtbdt(pcols),dqbdt(pcols),dtldt(pcols)
2785	real(r8) beta
2786	real(r8) capelmt
2787	real(r8) cp
2788	real(r8) dadt(pcols)
2789	real(r8) debdt
2790	real(r8) dltaa
2791	real(r8) eb
2792	real(r8) grav
2793
2794	integer i
2795	integer il1g
2796	integer il2g
2797	integer k, kmin, kmax
2798	integer msg
2799
2800	real(r8) rd
2801	real(r8) rl
2802	! change of subcloud layer properties due to convection is
2803	! related to cumulus updrafts and downdrafts.
2804	! mc(z)=f(z)mb, mub=betaumb, mdb=betad*mb are used
2805	! to define betau, betad and f(z).
2806	! note that this implies all time derivatives are in effect
2807	! time derivatives per unit cloud-base mass flux, i.e. they
2808	! have units of 1/mb instead of 1/sec.
2809	!
2810	do i = il1g,il2g
2811	mb(i) = 0._r8
2812	eb = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
2813	dtbdt(i) = (1._r8/dsubcld(i))* (mu(i,mx(i))*(shat(i,mx(i))-su(i,mx(i)))+ &
2814	md(i,mx(i))* (shat(i,mx(i))-sd(i,mx(i))))
2815	dqbdt(i) = (1._r8/dsubcld(i))* (mu(i,mx(i))*(qhat(i,mx(i))-qu(i,mx(i)))+ &
2816	md(i,mx(i))* (qhat(i,mx(i))-qd(i,mx(i))))
2817	debdt = eps1p(i,mx(i))/ (eps1+q(i,mx(i)))2dqbdt(i)
2818	dtldt(i) = -2840._r8* (3.5_r8/t(i,mx(i))*dtbdt(i)-debdt/eb)/ &
2819	(3.5_r8log(t(i,mx(i)))-log(eb)-4.805_r8)*2
2820	end do
2821	!
2822	! dtmdt and dqmdt are cumulus heating and drying.
2823	!
2824	do k = msg + 1,pver
2825	do i = il1g,il2g
2826	dtmdt(i,k) = 0._r8
2827	dqmdt(i,k) = 0._r8
2828	end do
2829	end do
2830	!
2831	do k = msg + 1,pver - 1
2832	do i = il1g,il2g
2833	if (k == jt(i)) then
2834	dtmdt(i,k) = (1._r8/dp(i,k))(mu(i,k+1) (su(i,k+1)-shat(i,k+1)- &
2835	rl/cpql(i,k+1))+md(i,k+1) (sd(i,k+1)-shat(i,k+1)))
2836	dqmdt(i,k) = (1._r8/dp(i,k))(mu(i,k+1) (qu(i,k+1)- &
2837	qhat(i,k+1)+ql(i,k+1))+md(i,k+1)*(qd(i,k+1)-qhat(i,k+1)))
2838	end if
2839	end do
2840	end do
2841	!
2842	beta = 0._r8
2843	do k = msg + 1,pver - 1
2844	do i = il1g,il2g
2845	if (k > jt(i) .and. k < mx(i)) then
2846	dtmdt(i,k) = (mc(i,k)* (shat(i,k)-s(i,k))+mc(i,k+1)* (s(i,k)-shat(i,k+1)))/ &
2847	dp(i,k) - rl/cpdu(i,k)(betaql(i,k)+ (1-beta)ql(i,k+1))
2848	! dqmdt(i,k)=(mc(i,k)*(qhat(i,k)-q(i,k))
2849	! 1 +mc(i,k+1)*(q(i,k)-qhat(i,k+1)))/dp(i,k)
2850	! 2 +du(i,k)*(qs(i,k)-q(i,k))
2851	! 3 +du(i,k)(betaql(i,k)+(1-beta)*ql(i,k+1))
2852
2853	dqmdt(i,k) = (mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)+cp/rl* (su(i,k+1)-s(i,k)))- &
2854	mu(i,k)* (qu(i,k)-qhat(i,k)+cp/rl(su(i,k)-s(i,k)))+md(i,k+1) &
2855	(qd(i,k+1)-qhat(i,k+1)+cp/rl(sd(i,k+1)-s(i,k)))-md(i,k) &
2856	(qd(i,k)-qhat(i,k)+cp/rl*(sd(i,k)-s(i,k))))/dp(i,k) + &
2857	du(i,k)* (betaql(i,k)+(1-beta)ql(i,k+1))
2858	end if
2859	end do
2860	end do
2861	!
2862	do k = msg + 1,pver
2863	do i = il1g,il2g
2864	if (k >= lel(i) .and. k <= lcl(i)) then
2865	thetavp(i,k) = tp(i,k)* (1000._r8/p(i,k))** (rd/cp)(1._r8+1.608_r8qstp(i,k)-q(i,mx(i)))
2866	thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)(1._r8+0.608_r8q(i,k))
2867	dqsdtp(i,k) = qstp(i,k)* (1._r8+qstp(i,k)/eps1)eps1rl/(rdtp(i,k)*2)
2868	!
2869	! dtpdt is the parcel temperature change due to change of
2870	! subcloud layer properties during convection.
2871	!
2872	dtpdt(i,k) = tp(i,k)/ (1._r8+rl/cp* (dqsdtp(i,k)-qstp(i,k)/tp(i,k)))* &
2873	(dtbdt(i)/t(i,mx(i))+rl/cp* (dqbdt(i)/tl(i)-q(i,mx(i))/ &
2874	tl(i)*2dtldt(i)))
2875	!
2876	! dboydt is the integrand of cape change.
2877	!
2878	dboydt(i,k) = ((dtpdt(i,k)/tp(i,k)+1._r8/(1._r8+1.608_r8qstp(i,k)-q(i,mx(i))) &
2879	(1.608_r8 * dqsdtp(i,k) * dtpdt(i,k) -dqbdt(i))) - (dtmdt(i,k)/t(i,k)+0.608_r8/ &
2880	(1._r8+0.608_r8q(i,k))dqmdt(i,k)))gravthetavp(i,k)/thetavm(i,k)
2881	end if
2882	end do
2883	end do
2884	!
2885	do k = msg + 1,pver
2886	do i = il1g,il2g
2887	if (k > lcl(i) .and. k < mx(i)) then
2888	thetavp(i,k) = tp(i,k)* (1000._r8/p(i,k))** (rd/cp)(1._r8+0.608_r8q(i,mx(i)))
2889	thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)(1._r8+0.608_r8q(i,k))
2890	!
2891	! dboydt is the integrand of cape change.
2892	!
2893	dboydt(i,k) = (dtbdt(i)/t(i,mx(i))+0.608_r8/ (1._r8+0.608_r8q(i,mx(i)))dqbdt(i)- &
2894	dtmdt(i,k)/t(i,k)-0.608_r8/ (1._r8+0.608_r8q(i,k))dqmdt(i,k))* &
2895	grav*thetavp(i,k)/thetavm(i,k)
2896	end if
2897	end do
2898	end do
2899
2900	!
2901	! buoyant energy change is set to 2/3*excess cape per 3 hours
2902	!
2903	dadt(il1g:il2g) = 0._r8
2904	kmin = minval(lel(il1g:il2g))
2905	kmax = maxval(mx(il1g:il2g)) - 1
2906	do k = kmin, kmax
2907	do i = il1g,il2g
2908	if ( k >= lel(i) .and. k <= mx(i) - 1) then
2909	dadt(i) = dadt(i) + dboydt(i,k)* (zf(i,k)-zf(i,k+1))
2910	endif
2911	end do
2912	end do
2913	do i = il1g,il2g
2914	dltaa = -1._r8* (cape(i)-capelmt)
2915	if (dadt(i) /= 0._r8) mb(i) = max(dltaa/tau/dadt(i),0._r8)
2916	end do
2917	!
2918	return
2919	end subroutine closure
2920
2921	subroutine q1q2_pjr(lchnk , &
2922	dqdt ,dsdt ,q ,qs ,qu , &
2923	su ,du ,qhat ,shat ,dp , &
2924	mu ,md ,sd ,qd ,ql , &
2925	dsubcld ,jt ,mx ,il1g ,il2g , &
2926	cp ,rl ,msg , &
2927	dl ,evp ,cu )
2928
2929
2930	implicit none
2931
2932	!-----------------------------------------------------------------------
2933	!
2934	! Purpose:
2935	! <Say what the routine does>
2936	!
2937	! Method:
2938	! <Describe the algorithm(s) used in the routine.>
2939	! <Also include any applicable external references.>
2940	!
2941	! Author: phil rasch dec 19 1995
2942	!
2943	!-----------------------------------------------------------------------
2944
2945
2946	real(r8), intent(in) :: cp
2947
2948	integer, intent(in) :: lchnk ! chunk identifier
2949	integer, intent(in) :: il1g
2950	integer, intent(in) :: il2g
2951	integer, intent(in) :: msg
2952
2953	real(r8), intent(in) :: q(pcols,pver)
2954	real(r8), intent(in) :: qs(pcols,pver)
2955	real(r8), intent(in) :: qu(pcols,pver)
2956	real(r8), intent(in) :: su(pcols,pver)
2957	real(r8), intent(in) :: du(pcols,pver)
2958	real(r8), intent(in) :: qhat(pcols,pver)
2959	real(r8), intent(in) :: shat(pcols,pver)
2960	real(r8), intent(in) :: dp(pcols,pver)
2961	real(r8), intent(in) :: mu(pcols,pver)
2962	real(r8), intent(in) :: md(pcols,pver)
2963	real(r8), intent(in) :: sd(pcols,pver)
2964	real(r8), intent(in) :: qd(pcols,pver)
2965	real(r8), intent(in) :: ql(pcols,pver)
2966	real(r8), intent(in) :: evp(pcols,pver)
2967	real(r8), intent(in) :: cu(pcols,pver)
2968	real(r8), intent(in) :: dsubcld(pcols)
2969
2970	real(r8),intent(out) :: dqdt(pcols,pver),dsdt(pcols,pver)
2971	real(r8),intent(out) :: dl(pcols,pver)
2972	integer kbm
2973	integer ktm
2974	integer jt(pcols)
2975	integer mx(pcols)
2976	!
2977	! work fields:
2978	!
2979	integer i
2980	integer k
2981
2982	real(r8) emc
2983	real(r8) rl
2984	!-------------------------------------------------------------------
2985	do k = msg + 1,pver
2986	do i = il1g,il2g
2987	dsdt(i,k) = 0._r8
2988	dqdt(i,k) = 0._r8
2989	dl(i,k) = 0._r8
2990	end do
2991	end do
2992	!
2993	! find the highest level top and bottom levels of convection
2994	!
2995	ktm = pver
2996	kbm = pver
2997	do i = il1g, il2g
2998	ktm = min(ktm,jt(i))
2999	kbm = min(kbm,mx(i))
3000	end do
3001
3002	do k = ktm,pver-1
3003	do i = il1g,il2g
3004	emc = -cu (i,k) & ! condensation in updraft
3005	+evp(i,k) ! evaporating rain in downdraft
3006
3007	dsdt(i,k) = -rl/cp*emc &
3008	+ (+mu(i,k+1)* (su(i,k+1)-shat(i,k+1)) &
3009	-mu(i,k)* (su(i,k)-shat(i,k)) &
3010	+md(i,k+1)* (sd(i,k+1)-shat(i,k+1)) &
3011	-md(i,k)* (sd(i,k)-shat(i,k)) &
3012	)/dp(i,k)
3013
3014	dqdt(i,k) = emc + &
3015	(+mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)) &
3016	-mu(i,k)* (qu(i,k)-qhat(i,k)) &
3017	+md(i,k+1)* (qd(i,k+1)-qhat(i,k+1)) &
3018	-md(i,k)* (qd(i,k)-qhat(i,k)) &
3019	)/dp(i,k)
3020
3021	dl(i,k) = du(i,k)*ql(i,k+1)
3022
3023	end do
3024	end do
3025
3026	!
3027	!DIR$ NOINTERCHANGE!
3028	do k = kbm,pver
3029	do i = il1g,il2g
3030	if (k == mx(i)) then
3031	dsdt(i,k) = (1._r8/dsubcld(i))* &
3032	(-mu(i,k)* (su(i,k)-shat(i,k)) &
3033	-md(i,k)* (sd(i,k)-shat(i,k)) &
3034	)
3035	dqdt(i,k) = (1._r8/dsubcld(i))* &
3036	(-mu(i,k)*(qu(i,k)-qhat(i,k)) &
3037	-md(i,k)*(qd(i,k)-qhat(i,k)) &
3038	)
3039	else if (k > mx(i)) then
3040	dsdt(i,k) = dsdt(i,k-1)
3041	dqdt(i,k) = dqdt(i,k-1)
3042	end if
3043	end do
3044	end do
3045	!
3046	return
3047	end subroutine q1q2_pjr
3048
3049	subroutine buoyan_dilute(lchnk ,ncol , &
3050	q ,t ,p ,z ,pf , &
3051	tp ,qstp ,tl ,rl ,cape , &
3052	pblt ,lcl ,lel ,lon ,mx , &
3053	rd ,grav ,cp ,msg , &
3054	tpert )
3055	!-----------------------------------------------------------------------
3056	!
3057	! Purpose:
3058	! Calculates CAPE the lifting condensation level and the convective top
3059	! where buoyancy is first -ve.
3060	!
3061	! Method: Calculates the parcel temperature based on a simple constant
3062	! entraining plume model. CAPE is integrated from buoyancy.
3063	! 09/09/04 - Simplest approach using an assumed entrainment rate for
3064	! testing (dmpdp).
3065	! 08/04/05 - Swap to convert dmpdz to dmpdp
3066	!
3067	! SCAM Logical Switches - DILUTE:RBN - Now Disabled
3068	! ---------------------
3069	! switch(1) = .T. - Uses the dilute parcel calculation to obtain tendencies.
3070	! switch(2) = .T. - Includes entropy/q changes due to condensate loss and freezing.
3071	! switch(3) = .T. - Adds the PBL Tpert for the parcel temperature at all levels.
3072	!
3073	! References:
3074	! Raymond and Blythe (1992) JAS
3075	!
3076	! Author:
3077	! Richard Neale - September 2004
3078	!
3079	!-----------------------------------------------------------------------
3080	implicit none
3081	!-----------------------------------------------------------------------
3082	!
3083	! input arguments
3084	!
3085	integer, intent(in) :: lchnk ! chunk identifier
3086	integer, intent(in) :: ncol ! number of atmospheric columns
3087
3088	real(r8), intent(in) :: q(pcols,pver) ! spec. humidity
3089	real(r8), intent(in) :: t(pcols,pver) ! temperature
3090	real(r8), intent(in) :: p(pcols,pver) ! pressure
3091	real(r8), intent(in) :: z(pcols,pver) ! height
3092	real(r8), intent(in) :: pf(pcols,pver+1) ! pressure at interfaces
3093	real(r8), intent(in) :: pblt(pcols) ! index of pbl depth
3094	real(r8), intent(in) :: tpert(pcols) ! perturbation temperature by pbl processes
3095
3096	!
3097	! output arguments
3098	!
3099	real(r8), intent(out) :: tp(pcols,pver) ! parcel temperature
3100	real(r8), intent(out) :: qstp(pcols,pver) ! saturation mixing ratio of parcel (only above lcl, just q below).
3101	real(r8), intent(out) :: tl(pcols) ! parcel temperature at lcl
3102	real(r8), intent(out) :: cape(pcols) ! convective aval. pot. energy.
3103	integer lcl(pcols) !
3104	integer lel(pcols) !
3105	integer lon(pcols) ! level of onset of deep convection
3106	integer mx(pcols) ! level of max moist static energy
3107	!
3108	!--------------------------Local Variables------------------------------
3109	!
3110	real(r8) capeten(pcols,5) ! provisional value of cape
3111	real(r8) tv(pcols,pver) !
3112	real(r8) tpv(pcols,pver) !
3113	real(r8) buoy(pcols,pver)
3114
3115	real(r8) a1(pcols)
3116	real(r8) a2(pcols)
3117	real(r8) estp(pcols)
3118	real(r8) pl(pcols)
3119	real(r8) plexp(pcols)
3120	real(r8) hmax(pcols)
3121	real(r8) hmn(pcols)
3122	real(r8) y(pcols)
3123
3124	logical plge600(pcols)
3125	integer knt(pcols)
3126	integer lelten(pcols,5)
3127
3128	real(r8) cp
3129	real(r8) e
3130	real(r8) grav
3131
3132	integer i
3133	integer k
3134	integer msg
3135	integer n
3136
3137	real(r8) rd
3138	real(r8) rl
3139	#ifdef PERGRO
3140	real(r8) rhd
3141	#endif
3142	!
3143	!-----------------------------------------------------------------------
3144	!
3145	do n = 1,5
3146	do i = 1,ncol
3147	lelten(i,n) = pver
3148	capeten(i,n) = 0._r8
3149	end do
3150	end do
3151	!
3152	do i = 1,ncol
3153	lon(i) = pver
3154	knt(i) = 0
3155	lel(i) = pver
3156	mx(i) = lon(i)
3157	cape(i) = 0._r8
3158	hmax(i) = 0._r8
3159	end do
3160
3161	tp(:ncol,:) = t(:ncol,:)
3162	qstp(:ncol,:) = q(:ncol,:)
3163
3164	!!! RBN - Initialize tv and buoy for output.
3165	!!! tv=tv : tpv=tpv : qstp=q : buoy=0.
3166	tv(:ncol,:) = t(:ncol,:) (1._r8+1.608_r8q(:ncol,:))/ (1._r8+q(:ncol,:))
3167	tpv(:ncol,:) = tv(:ncol,:)
3168	buoy(:ncol,:) = 0._r8
3169
3170	!
3171	! set "launching" level(mx) to be at maximum moist static energy.
3172	! search for this level stops at planetary boundary layer top.
3173	!
3174	#ifdef PERGRO
3175	do k = pver,msg + 1,-1
3176	do i = 1,ncol
3177	hmn(i) = cpt(i,k) + gravz(i,k) + rl*q(i,k)
3178	!
3179	! Reset max moist static energy level when relative difference exceeds 1.e-4
3180	!
3181	rhd = (hmn(i) - hmax(i))/(hmn(i) + hmax(i))
3182	if (k >= nint(pblt(i)) .and. k <= lon(i) .and. rhd > -1.e-4_r8) then
3183	hmax(i) = hmn(i)
3184	mx(i) = k
3185	end if
3186	end do
3187	end do
3188	#else
3189	do k = pver,msg + 1,-1
3190	do i = 1,ncol
3191	hmn(i) = cpt(i,k) + gravz(i,k) + rl*q(i,k)
3192	if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then
3193	hmax(i) = hmn(i)
3194	mx(i) = k
3195	end if
3196	end do
3197	end do
3198	#endif
3199
3200	! LCL dilute calculation - initialize to mx(i)
3201	! Determine lcl in parcel_dilute and get pl,tl after parcel_dilute
3202	! Original code actually sets LCL as level above wher condensate forms.
3203	! Therefore in parcel_dilute lcl(i) will be at first level where qsmix < qtmix.
3204
3205	do i = 1,ncol ! Initialise LCL variables.
3206	lcl(i) = mx(i)
3207	tl(i) = t(i,mx(i))
3208	pl(i) = p(i,mx(i))
3209	end do
3210
3211	!
3212	! main buoyancy calculation.
3213	!
3214	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3215	!!! DILUTE PLUME CALCULATION USING ENTRAINING PLUME !!!
3216	!!! RBN 9/9/04 !!!
3217
3218	call parcel_dilute(lchnk, ncol, msg, mx, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl)
3219
3220
3221	! If lcl is above the nominal level of non-divergence (600 mbs),
3222	! no deep convection is permitted (ensuing calculations
3223	! skipped and cape retains initialized value of zero).
3224	!
3225	do i = 1,ncol
3226	plge600(i) = pl(i).ge.600._r8 ! Just change to always allow buoy calculation.
3227	end do
3228
3229	!
3230	! Main buoyancy calculation.
3231	!
3232	do k = pver,msg + 1,-1
3233	do i=1,ncol
3234	if (k <= mx(i) .and. plge600(i)) then ! Define buoy from launch level to cloud top.
3235	tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k))
3236	buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add ! +0.5K or not?
3237	else
3238	qstp(i,k) = q(i,k)
3239	tp(i,k) = t(i,k)
3240	tpv(i,k) = tv(i,k)
3241	endif
3242	end do
3243	end do
3244
3245
3246
3247	!-------------------------------------------------------------------------------
3248
3249	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3250
3251
3252
3253	!
3254	do k = msg + 2,pver
3255	do i = 1,ncol
3256	if (k < lcl(i) .and. plge600(i)) then
3257	if (buoy(i,k+1) > 0. .and. buoy(i,k) <= 0._r8) then
3258	knt(i) = min(5,knt(i) + 1)
3259	lelten(i,knt(i)) = k
3260	end if
3261	end if
3262	end do
3263	end do
3264	!
3265	! calculate convective available potential energy (cape).
3266	!
3267	do n = 1,5
3268	do k = msg + 1,pver
3269	do i = 1,ncol
3270	if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then
3271	capeten(i,n) = capeten(i,n) + rdbuoy(i,k)log(pf(i,k+1)/pf(i,k))
3272	end if
3273	end do
3274	end do
3275	end do
3276	!
3277	! find maximum cape from all possible tentative capes from
3278	! one sounding,
3279	! and use it as the final cape, april 26, 1995
3280	!
3281	do n = 1,5
3282	do i = 1,ncol
3283	if (capeten(i,n) > cape(i)) then
3284	cape(i) = capeten(i,n)
3285	lel(i) = lelten(i,n)
3286	end if
3287	end do
3288	end do
3289	!
3290	! put lower bound on cape for diagnostic purposes.
3291	!
3292	do i = 1,ncol
3293	cape(i) = max(cape(i), 0._r8)
3294	end do
3295	!
3296	return
3297	end subroutine buoyan_dilute
3298
3299	subroutine parcel_dilute (lchnk, ncol, msg, klaunch, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl)
3300
3301	! Routine to determine
3302	! 1. Tp - Parcel temperature
3303	! 2. qstp - Saturated mixing ratio at the parcel temperature.
3304
3305	!--------------------
3306	implicit none
3307	!--------------------
3308
3309	integer, intent(in) :: lchnk
3310	integer, intent(in) :: ncol
3311	integer, intent(in) :: msg
3312
3313	integer, intent(in), dimension(pcols) :: klaunch(pcols)
3314
3315	real(r8), intent(in), dimension(pcols,pver) :: p
3316	real(r8), intent(in), dimension(pcols,pver) :: t
3317	real(r8), intent(in), dimension(pcols,pver) :: q
3318	real(r8), intent(in), dimension(pcols) :: tpert ! PBL temperature perturbation.
3319
3320	real(r8), intent(inout), dimension(pcols,pver) :: tp ! Parcel temp.
3321	real(r8), intent(inout), dimension(pcols,pver) :: qstp ! Parcel water vapour (sat value above lcl).
3322	real(r8), intent(inout), dimension(pcols) :: tl ! Actual temp of LCL.
3323	real(r8), intent(inout), dimension(pcols) :: pl ! Actual pressure of LCL.
3324
3325	integer, intent(inout), dimension(pcols) :: lcl ! Lifting condesation level (first model level with saturation).
3326
3327	real(r8), intent(out), dimension(pcols,pver) :: tpv ! Define tpv within this routine.
3328
3329	!--------------------
3330
3331	! Have to be careful as s is also dry static energy.
3332
3333
3334	! If we are to retain the fact that CAM loops over grid-points in the internal
3335	! loop then we need to dimension sp,atp,mp,xsh2o with ncol.
3336
3337
3338	real(r8) tmix(pcols,pver) ! Tempertaure of the entraining parcel.
3339	real(r8) qtmix(pcols,pver) ! Total water of the entraining parcel.
3340	real(r8) qsmix(pcols,pver) ! Saturated mixing ratio at the tmix.
3341	real(r8) smix(pcols,pver) ! Entropy of the entraining parcel.
3342	real(r8) xsh2o(pcols,pver) ! Precipitate lost from parcel.
3343	real(r8) ds_xsh2o(pcols,pver) ! Entropy change due to loss of condensate.
3344	real(r8) ds_freeze(pcols,pver) ! Entropy change sue to freezing of precip.
3345
3346	real(r8) mp(pcols) ! Parcel mass flux.
3347	real(r8) qtp(pcols) ! Parcel total water.
3348	real(r8) sp(pcols) ! Parcel entropy.
3349
3350	real(r8) sp0(pcols) ! Parcel launch entropy.
3351	real(r8) qtp0(pcols) ! Parcel launch total water.
3352	real(r8) mp0(pcols) ! Parcel launch relative mass flux.
3353
3354	real(r8) lwmax ! Maximum condesate that can be held in cloud before rainout.
3355	real(r8) dmpdp ! Parcel fractional mass entrainment rate (/mb).
3356	!real(r8) dmpdpc ! In cloud parcel mass entrainment rate (/mb).
3357	real(r8) dmpdz ! Parcel fractional mass entrainment rate (/m)
3358	real(r8) dpdz,dzdp ! Hydrstatic relation and inverse of.
3359	real(r8) senv ! Environmental entropy at each grid point.
3360	real(r8) qtenv ! Environmental total water " " ".
3361	real(r8) penv ! Environmental total pressure " " ".
3362	real(r8) tenv ! Environmental total temperature " " ".
3363	real(r8) new_s ! Hold value for entropy after condensation/freezing adjustments.
3364	real(r8) new_q ! Hold value for total water after condensation/freezing adjustments.
3365	real(r8) dp ! Layer thickness (center to center)
3366	real(r8) tfguess ! First guess for entropy inversion - crucial for efficiency!
3367	real(r8) tscool ! Super cooled temperature offset (in degC) (eg -35).
3368
3369	real(r8) qxsk, qxskp1 ! LCL excess water (k, k+1)
3370	real(r8) dsdp, dqtdp, dqxsdp ! LCL s, qt, p gradients (k, k+1)
3371	real(r8) slcl,qtlcl,qslcl ! LCL s, qt, qs values.
3372
3373	integer rcall ! Number of ientropy call for errors recording
3374	integer nit_lheat ! Number of iterations for condensation/freezing loop.
3375	integer i,k,ii ! Loop counters.
3376
3377	!======================================================================
3378	! SUMMARY
3379	!
3380	! 9/9/04 - Assumes parcel is initiated from level of maxh (klaunch)
3381	! and entrains at each level with a specified entrainment rate.
3382	!
3383	! 15/9/04 - Calculates lcl(i) based on k where qsmix is first < qtmix.
3384	!
3385	!======================================================================
3386	!
3387	! Set some values that may be changed frequently.
3388	!
3389
3390	nit_lheat = 2 ! iterations for ds,dq changes from condensation freezing.
3391	dmpdz=-1.e-3_r8 ! Entrainment rate. (-ve for /m)
3392	!dmpdpc = 3.e-2_r8 ! In cloud entrainment rate (/mb).
3393	lwmax = 1.e-3_r8 ! Need to put formula in for this.
3394	tscool = 0.0_r8 ! Temp at which water loading freezes in the cloud.
3395
3396	qtmix=0._r8
3397	smix=0._r8
3398
3399	qtenv = 0._r8
3400	senv = 0._r8
3401	tenv = 0._r8
3402	penv = 0._r8
3403
3404	qtp0 = 0._r8
3405	sp0 = 0._r8
3406	mp0 = 0._r8
3407
3408	qtp = 0._r8
3409	sp = 0._r8
3410	mp = 0._r8
3411
3412	new_q = 0._r8
3413	new_s = 0._r8
3414
3415	! ** Begin loops **
3416
3417	do k = pver, msg+1, -1
3418	do i=1,ncol
3419
3420	! Initialize parcel values at launch level.
3421
3422	if (k == klaunch(i)) then
3423	qtp0(i) = q(i,k) ! Parcel launch total water (assuming subsaturated) - OK????.
3424	sp0(i) = entropy(t(i,k),p(i,k),qtp0(i)) ! Parcel launch entropy.
3425	mp0(i) = 1._r8 ! Parcel launch relative mass (i.e. 1 parcel stays 1 parcel for dmpdp=0, undilute).
3426	smix(i,k) = sp0(i)
3427	qtmix(i,k) = qtp0(i)
3428	tfguess = t(i,k)
3429	rcall = 1
3430	call ientropy (rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess)
3431	end if
3432
3433	! Entraining levels
3434
3435	if (k < klaunch(i)) then
3436
3437	! Set environmental values for this level.
3438
3439	dp = (p(i,k)-p(i,k+1)) ! In -ve mb as p decreasing with height - difference between center of layers.
3440	qtenv = 0.5_r8*(q(i,k)+q(i,k+1)) ! Total water of environment.
3441	tenv = 0.5_r8*(t(i,k)+t(i,k+1))
3442	penv = 0.5_r8*(p(i,k)+p(i,k+1))
3443
3444	senv = entropy(tenv,penv,qtenv) ! Entropy of environment.
3445
3446	! Determine fractional entrainment rate /pa given value /m.
3447
3448	dpdz = -(penvgrav)/(rgastenv) ! in mb/m since p in mb.
3449	dzdp = 1._r8/dpdz ! in m/mb
3450	dmpdp = dmpdz*dzdp ! /mb Fractional entrainment
3451
3452	! Sum entrainment to current level
3453	! entrains q,s out of intervening dp layers, in which linear variation is assumed
3454	! so really it entrains the mean of the 2 stored values.
3455
3456	sp(i) = sp(i) - dmpdpdpsenv
3457	qtp(i) = qtp(i) - dmpdpdpqtenv
3458	mp(i) = mp(i) - dmpdp*dp
3459
3460	! Entrain s and qt to next level.
3461
3462	smix(i,k) = (sp0(i) + sp(i)) / (mp0(i) + mp(i))
3463	qtmix(i,k) = (qtp0(i) + qtp(i)) / (mp0(i) + mp(i))
3464
3465	! Invert entropy from s and q to determine T and saturation-capped q of mixture.
3466	! t(i,k) used as a first guess so that it converges faster.
3467
3468	tfguess = tmix(i,k+1)
3469	rcall = 2
3470	call ientropy(rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess)
3471
3472	!
3473	! Determine if this is lcl of this column if qsmix <= qtmix.
3474	! FIRST LEVEL where this happens on ascending.
3475
3476	if (qsmix(i,k) <= qtmix(i,k) .and. qsmix(i,k+1) > qtmix(i,k+1)) then
3477	lcl(i) = k
3478	qxsk = qtmix(i,k) - qsmix(i,k)
3479	qxskp1 = qtmix(i,k+1) - qsmix(i,k+1)
3480	dqxsdp = (qxsk - qxskp1)/dp
3481	pl(i) = p(i,k+1) - qxskp1/dqxsdp ! pressure level of actual lcl.
3482	dsdp = (smix(i,k) - smix(i,k+1))/dp
3483	dqtdp = (qtmix(i,k) - qtmix(i,k+1))/dp
3484	slcl = smix(i,k+1) + dsdp* (pl(i)-p(i,k+1))
3485	qtlcl = qtmix(i,k+1) + dqtdp*(pl(i)-p(i,k+1))
3486
3487	tfguess = tmix(i,k)
3488	rcall = 3
3489	call ientropy (rcall,i,lchnk,slcl,pl(i),qtlcl,tl(i),qslcl,tfguess)
3490
3491	! write(iulog,*)' '
3492	! write(iulog,*)' p',p(i,k+1),pl(i),p(i,lcl(i))
3493	! write(iulog,*)' t',tmix(i,k+1),tl(i),tmix(i,lcl(i))
3494	! write(iulog,*)' s',smix(i,k+1),slcl,smix(i,lcl(i))
3495	! write(iulog,*)'qt',qtmix(i,k+1),qtlcl,qtmix(i,lcl(i))
3496	! write(iulog,*)'qs',qsmix(i,k+1),qslcl,qsmix(i,lcl(i))
3497
3498	endif
3499	!
3500	end if ! k < klaunch
3501
3502
3503	end do ! Levels loop
3504	end do ! Columns loop
3505
3506	!!!!!!!!!!!!!!!!!!!!!!!!!!END ENTRAINMENT LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3507
3508	!! Could stop now and test with this as it will provide some estimate of buoyancy
3509	!! without the effects of freezing/condensation taken into account for tmix.
3510
3511	!! So we now have a profile of entropy and total water of the entraining parcel
3512	!! Varying with height from the launch level klaunch parcel=environment. To the
3513	!! top allowed level for the existence of convection.
3514
3515	!! Now we have to adjust these values such that the water held in vaopor is < or
3516	!! = to qsmix. Therefore, we assume that the cloud holds a certain amount of
3517	!! condensate (lwmax) and the rest is rained out (xsh2o). This, obviously
3518	!! provides latent heating to the mixed parcel and so this has to be added back
3519	!! to it. But does this also increase qsmix as well? Also freezing processes
3520
3521
3522	xsh2o = 0._r8
3523	ds_xsh2o = 0._r8
3524	ds_freeze = 0._r8
3525
3526	!!!!!!!!!!!!!!!!!!!!!!!!!PRECIPITATION/FREEZING LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!
3527	!! Iterate solution twice for accuracy
3528
3529
3530
3531	do k = pver, msg+1, -1
3532	do i=1,ncol
3533
3534	! Initialize variables at k=klaunch
3535
3536	if (k == klaunch(i)) then
3537
3538	! Set parcel values at launch level assume no liquid water.
3539
3540	tp(i,k) = tmix(i,k)
3541	qstp(i,k) = q(i,k)
3542	tpv(i,k) = (tp(i,k) + tpert(i)) * (1._r8+1.608_r8*qstp(i,k)) / (1._r8+qstp(i,k))
3543
3544	end if
3545
3546	if (k < klaunch(i)) then
3547
3548	! Initiaite loop if switch(2) = .T. - RBN:DILUTE - TAKEN OUT BUT COULD BE RETURNED LATER.
3549
3550	! Iterate nit_lheat times for s,qt changes.
3551
3552	do ii=0,nit_lheat-1
3553
3554	! Rain (xsh2o) is excess condensate, bar LWMAX (Accumulated loss from qtmix).
3555
3556	xsh2o(i,k) = max (0._r8, qtmix(i,k) - qsmix(i,k) - lwmax)
3557
3558	! Contribution to ds from precip loss of condensate (Accumulated change from smix).(-ve)
3559
3560	ds_xsh2o(i,k) = ds_xsh2o(i,k+1) - cpliq * log (tmix(i,k)/tfreez) * max(0._r8,(xsh2o(i,k)-xsh2o(i,k+1)))
3561	!
3562	! Entropy of freezing: latice times amount of water involved divided by T.
3563	!
3564
3565	if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) == 0._r8) then ! One off freezing of condensate.
3566	ds_freeze(i,k) = (latice/tmix(i,k)) * max(0._r8,qtmix(i,k)-qsmix(i,k)-xsh2o(i,k)) ! Gain of LH
3567	end if
3568
3569	if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) /= 0._r8) then ! Continual freezing of additional condensate.
3570	ds_freeze(i,k) = ds_freeze(i,k+1)+(latice/tmix(i,k)) * max(0._r8,(qsmix(i,k+1)-qsmix(i,k)))
3571	end if
3572
3573	! Adjust entropy and accordingly to sum of ds (be careful of signs).
3574
3575	new_s = smix(i,k) + ds_xsh2o(i,k) + ds_freeze(i,k)
3576
3577	! Adjust liquid water and accordingly to xsh2o.
3578
3579	new_q = qtmix(i,k) - xsh2o(i,k)
3580
3581	! Invert entropy to get updated Tmix and qsmix of parcel.
3582
3583	tfguess = tmix(i,k)
3584	rcall =4
3585	call ientropy (rcall,i,lchnk,new_s, p(i,k), new_q, tmix(i,k), qsmix(i,k), tfguess)
3586
3587	end do ! Iteration loop for freezing processes.
3588
3589	! tp - Parcel temp is temp of mixture.
3590	! tpv - Parcel v. temp should be density temp with new_q total water.
3591
3592	tp(i,k) = tmix(i,k)
3593
3594	! tpv = tprho in the presence of condensate (i.e. when new_q > qsmix)
3595
3596	if (new_q > qsmix(i,k)) then ! Super-saturated so condensate present - reduces buoyancy.
3597	qstp(i,k) = qsmix(i,k)
3598	else ! Just saturated/sub-saturated - no condensate virtual effects.
3599	qstp(i,k) = new_q
3600	end if
3601
3602	tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k)) / (1._r8+ new_q)
3603
3604	end if ! k < klaunch
3605
3606	end do ! Loop for columns
3607
3608	end do ! Loop for vertical levels.
3609
3610
3611	return
3612	end subroutine parcel_dilute
3613
3614	!-----------------------------------------------------------------------------------------
3615	real(r8) function entropy(TK,p,qtot)
3616	!-----------------------------------------------------------------------------------------
3617	!
3618	! TK(K),p(mb),qtot(kg/kg)
3619	! from Raymond and Blyth 1992
3620	!
3621	real(r8), intent(in) :: p,qtot,TK
3622	real(r8) :: qv,qsat,e,esat,L,eref,pref
3623
3624	pref = 1000.0_r8 ! mb
3625	eref = 6.106_r8 ! sat p at tfreez (mb)
3626
3627	L = rl - (cpliq - cpwv)*(TK-tfreez) ! T IN CENTIGRADE
3628
3629	! Replace call to satmixutils.
3630
3631	esat = c1exp(c2(TK-tfreez)/(c3+TK-tfreez)) ! esat(T) in mb
3632	qsat=eps1*esat/(p-esat) ! Sat. mixing ratio (in kg/kg).
3633
3634	qv = min(qtot,qsat) ! Partition qtot into vapor part only.
3635	e = qv*p / (eps1 +qv)
3636
3637	entropy = (cpres + qtotcpliq)log( TK/tfreez) - rgas*log( (p-e)/pref ) + &
3638	Lqv/TK - qvrh2o*log(qv/qsat)
3639	!
3640	return
3641	end FUNCTION entropy
3642
3643	!
3644	!-----------------------------------------------------------------------------------------
3645	SUBROUTINE ientropy (rcall,icol,lchnk,s,p,qt,T,qsat,Tfg)
3646	!-----------------------------------------------------------------------------------------
3647	!
3648	! p(mb), Tfg/T(K), qt/qv(kg/kg), s(J/kg).
3649	! Inverts entropy, pressure and total water qt
3650	! for T and saturated vapor mixing ratio
3651	!
3652	#ifndef WRF_PORT
3653	use phys_grid, only: get_rlon_p, get_rlat_p
3654	#endif
3655
3656	integer, intent(in) :: icol, lchnk, rcall
3657	real(r8), intent(in) :: s, p, Tfg, qt
3658	real(r8), intent(out) :: qsat, T
3659	real(r8) :: qv,Ts,dTs,fs1,fs2,esat
3660	real(r8) :: pref,eref,L,e
3661	real(r8) :: this_lat,this_lon
3662	integer :: LOOPMAX,i
3663
3664	LOOPMAX = 100 !* max number of iteration loops
3665
3666	! Values for entropy
3667	pref = 1000.0_r8 ! mb ref pressure.
3668	eref = 6.106_r8 ! sat p at tfreez (mb)
3669
3670	! Invert the entropy equation -- use Newton's method
3671
3672	Ts = Tfg ! Better first guess based on Tprofile from conv.
3673
3674	converge: do i=0, LOOPMAX
3675
3676	L = rl - (cpliq - cpwv)*(Ts-tfreez)
3677
3678	esat = c1exp(c2(Ts-tfreez)/(c3+Ts-tfreez)) ! Bolton (eq. 10)
3679	qsat = eps1*esat/(p-esat)
3680	qv = min(qt,qsat)
3681	e = qv*p / (eps1 +qv) ! Bolton (eq. 16)
3682	fs1 = (cpres + qtcpliq)log( Ts/tfreez ) - rgas*log( (p-e)/pref ) + &
3683	Lqv/Ts - qvrh2o*log(qv/qsat) - s
3684
3685	L = rl - (cpliq - cpwv)*(Ts-1._r8-tfreez)
3686
3687	esat = c1exp(c2(Ts-1._r8-tfreez)/(c3+Ts-1._r8-tfreez))
3688	qsat = eps1*esat/(p-esat)
3689	qv = min(qt,qsat)
3690	e = qv*p / (eps1 +qv)
3691	fs2 = (cpres + qtcpliq)log( (Ts-1._r8)/tfreez ) - rgas*log( (p-e)/pref ) + &
3692	Lqv/(Ts-1._r8) - qvrh2o*log(qv/qsat) - s
3693
3694	dTs = fs1/(fs2 - fs1)
3695	Ts = Ts+dTs
3696	if (abs(dTs).lt.0.001_r8) exit converge
3697	if (i .eq. LOOPMAX - 1) then
3698	#ifndef WRF_PORT
3699	this_lat = get_rlat_p(lchnk, icol)*57.296_r8
3700	this_lon = get_rlon_p(lchnk, icol)*57.296_r8
3701	#else
3702	!Do not worry about the specific lat/lon in WRF
3703	this_lat = 0.
3704	this_lon = 0.
3705	#endif
3706	write(iulog,) ' ZM_CONV: IENTROPY: Failed and about to exit, info follows **'
3707	#ifdef WRF_PORT
3708	call wrf_message(iulog)
3709	#endif
3710	write(iulog,100) 'ZM_CONV: IENTROPY. Details: call#,lchnk,icol= ',rcall,lchnk,icol, &
3711	' lat: ',this_lat,' lon: ',this_lon, &
3712	' P(mb)= ', p, ' Tfg(K)= ', Tfg, ' qt(g/kg) = ', 1000._r8*qt, &
3713	' qsat(g/kg) = ', 1000._r8*qsat,', s(J/kg) = ',s
3714	#ifdef WRF_PORT
3715	call wrf_message(iulog)
3716	#endif
3717	call endrun('** ZM_CONV IENTROPY: Tmix did not converge **')
3718	end if
3719	enddo converge
3720
3721	! Replace call to satmixutils.
3722
3723	esat = c1exp(c2(Ts-tfreez)/(c3+Ts-tfreez))
3724	qsat=eps1*esat/(p-esat)
3725
3726	qv = min(qt,qsat) ! /* check for saturation */
3727	T = Ts
3728
3729	100 format (A,I1,I4,I4,7(A,F6.2))
3730
3731	return
3732	end SUBROUTINE ientropy
3733
3734	end module module_cu_camzm

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

Download in other formats:

Original Format