Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_camzm_driver.F @ 1

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 38.6 KB

Line
1	MODULE module_cu_camzm_driver
2
3	!Roughly based on zm_conv_intr.F90 from CAM
4
5	USE module_cam_support, only: pcnst, pcols, pver, pverp
6
7	IMPLICIT NONE
8
9	PRIVATE !Default to private
10	PUBLIC :: & !Public entities
11	camzm_driver, &
12	zm_conv_init
13
14	CONTAINS
15
16	!------------------------------------------------------------------------
17	SUBROUTINE camzm_driver( &
18	ids,ide, jds,jde, kds,kde &
19	,ims,ime, jms,jme, kms,kme &
20	,its,ite, jts,jte, kts,kte &
21	,itimestep, bl_pbl_physics, sf_sfclay_physics &
22	,th, t_phy, tsk, tke_pbl, ust, qv, qc, qi &
23	,mavail, kpbl, pblh, xland, z &
24	,z_at_w, dz8w, ht, p, p8w, pi_phy, psfc &
25	,u_phy, v_phy, hfx, qfx, cldfra &
26	,tpert_camuwpbl &
27	,dx, dt, stepcu, cudt, curr_secs, adapt_step_flag&
28	,cape_out, mu_out, md_out, zmdt, zmdq &
29	,rliq_out, dlf_out &
30	,pconvb, pconvt, cubot, cutop, raincv, pratec &
31	,rucuten, rvcuten &
32	,rthcuten, rqvcuten, rqccuten, rqicuten &
33	,evaptzm, fzsntzm, evsntzm, evapqzm, zmflxprc &
34	,zmflxsnw, zmntprpd, zmntsnpd, zmeiheat &
35	,cmfmc, cmfmcdzm, preccdzm, precz &
36	,zmmtu, zmmtv, zmupgu, zmupgd, zmvpgu, zmvpgd &
37	,zmicuu, zmicud, zmicvu, zmicvd &
38	,zmdice, zmdliq &
39	)
40	! This routine is based on zm_conv_tend in CAM. It handles the mapping of
41	! variables from the WRF to the CAM framework for the Zhang-McFarlane
42	! convective parameterization.
43	!
44	! Author: William.Gustafson@pnl.gov, Nov. 2009
45	! Last modified: William.Gustafson@pnl.gov, Nov. 2010
46	!------------------------------------------------------------------------
47	USE shr_kind_mod, only: r8 => shr_kind_r8
48	USE physconst, only: cpair, gravit
49	USE module_cu_camzm, only: convtran, momtran, zm_conv_evap, zm_convr
50
51	! Subroutine arguments...
52	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
53	ims,ime, jms,jme, kms,kme, &
54	its,ite, jts,jte, kts,kte, &
55	bl_pbl_physics, & !pbl scheme
56	itimestep, & !time step index
57	sf_sfclay_physics, & !surface layer scheme
58	stepcu !number of time steps between Cu calls
59
60	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: &
61	cldfra, & !cloud fraction
62	dz8w, & !height between interfaces (m)
63	p, & !pressure at mid-level (Pa)
64	p8w, & !pressure at level interface (Pa)
65	pi_phy, & !exner function, (p0/p)^(R/cpair) (none)
66	qv, & !water vapor mixing ratio (kg/kg-dry air)
67	th, & !potential temperature (K)
68	tke_pbl, & !turbulent kinetic energy from PBL (m2/s2)
69	t_phy, & !temperature (K)
70	u_phy, & !zonal wind component on T points (m/s)
71	v_phy, & !meridional wind component on T points (m/s)
72	z, & !height above sea level at mid-level (m)
73	z_at_w !height above sea level at interface (m)
74
75	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN), OPTIONAL :: &
76	qc, & !cloud droplet mixing ratio (kg/kg-dry air)
77	qi !cloud ice crystal mixing ratio (kg/kg-dry air)
78
79	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: &
80	dlf_out, & !detraining cloud water tendendcy
81	evaptzm, & !temp. tendency - evap/snow prod from ZM (K/s)
82	fzsntzm, & !temp. tendency - rain to snow conversion from ZM (K/s)
83	evsntzm, & !temp. tendency - snow to rain conversion from ZM (K/s)
84	evapqzm, & !spec. humidity tend. - evaporation from ZM (kg/kg/s)
85	zmflxprc, & !flux of precipitation from ZM (kg/m2/s)
86	zmflxsnw, & !flux of snow from ZM (kg/m2/s)
87	zmntprpd, & !net precipitation production from ZM (kg/kg/s)
88	zmntsnpd, & !net snow production from ZM (kg/kg/s)
89	zmeiheat, & !heating by ice and evaporation from ZM (W/kg)
90	cmfmc, & !convective mass flux--m sub c, deep here but ultimately deep+shallow (kg/m2/s)
91	cmfmcdzm, & !convection mass flux from ZM deep (kg/m2/s)
92	md_out, & !output convective downdraft mass flux (kg/m2/s)
93	mu_out, & !output convective updraft mass flux (kg/m2/s)
94	rucuten, & !UNcoupled zonal wind tendency due to Cu scheme (m/s2)
95	rvcuten, & !UNcoupled meridional wind tendency due to Cu scheme (m/s2)
96	rthcuten, & !UNcoupled potential temperature tendendcy due to cu scheme (K/s)
97	rqvcuten, & !UNcoupled water vapor mixing ratio tendency due to Cu scheme (kg/kg/s)
98	zmdt, & !temp. tendency from moist convection (K/s)
99	zmdq, & !spec. humidity tendency from moist convection (kg/kg/s)
100	zmmtu, & !U tendency from ZM convective momentum transport (m/s2)
101	zmmtv, & !V tendency from ZM convective momentum transport (m/s2)
102	zmupgu, & !zonal force from ZM updraft pressure gradient term (m/s2)
103	zmupgd, & !zonal force from ZM downdraft pressure gradient term (m/s2)
104	zmvpgu, & !meridional force from ZM updraft pressure gradient term (m/s2)
105	zmvpgd, & !meridional force from ZM downdraft pressure gradient term (m/s2)
106	zmicuu, & !ZM in-cloud U updrafts (m/s)
107	zmicud, & !ZM in-cloud U downdrafts (m/s)
108	zmicvu, & !ZM in-cloud V updrafts (m/s)
109	zmicvd, & !ZM in-cloud V downdrafts (m/s)
110	zmdice, & !ZM cloud ice tendency (kg/kg/s)
111	zmdliq !ZM cloud liquid tendency (kg/kg/s)
112
113	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT), OPTIONAL :: &
114	rqccuten, & !UNcoupled cloud droplet mixing ratio tendency due to Cu scheme (kg/kg/s)
115	rqicuten !UNcoupled ice crystal mixing ratio tendency due to Cu scheme (kg/kg/s)
116
117	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: &
118	hfx, & !upward heat flux at sfc (W/m2)
119	ht, & !terrain height (m)
120	xland, & !land/water mask, 1=land, 2=water
121	mavail, & !soil moisture availability
122	pblh, & !planetary boundary layer height (m)
123	psfc, & !surface pressure (Pa)
124	qfx, & !upward moisture flux at sfc (kg/m2/s)
125	tpert_camuwpbl, & !temperature perturbation from UW CAM PBL
126	tsk, & !skin temperature (K)
127	ust !u* in similarity theory (m/s)
128
129	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
130	cape_out, & !convective available potential energy (J/kg)
131	cubot, & !level number of base of convection
132	cutop, & !level number of top of convection
133	pconvb, & !pressure of base of convection (Pa)
134	pconvt, & !pressure of top of convection (Pa)
135	pratec, & !rain rate returned to WRF for accumulation (mm/s)
136	preccdzm, & !convection precipitation rate from ZM deep (m/s)
137	precz, & !total precipitation rate from ZM (m/s)
138	raincv, & !time-step convective rain amount (mm)
139	rliq_out !vertical integral of reserved cloud water
140
141	REAL, INTENT(IN) :: &
142	cudt, & !cumulus time step (min)
143	curr_secs, & !current forecast time (s)
144	dt, & !domain time step (s)
145	dx !grid spacing (m)
146
147	INTEGER, DIMENSION( ims:ime, jms:jme), INTENT(IN) :: &
148	kpbl !index of PBL level
149
150	LOGICAL, INTENT(IN) :: &
151	adapt_step_flag !using adaptive time steps?
152
153	! Local variables...
154	!Variables dimensioned for input to ZM routines
155	REAL(r8), DIMENSION(pcols, kte+1) :: &
156	mcon, & !convective mass flux--m sub c (sub arg out in CAM)
157	pflx, & !scattered precip flux at each level (sub arg out in CAM)
158	pint8, & !pressure at layer interface (Pa)
159	zi8 !height above sea level at mid-level (m)
160
161	REAL(r8), DIMENSION(pcols, kte, pcnst) :: &
162	qh8 !specific humidity (kg/kg-moist air)
163
164	REAL(r8), DIMENSION(pcols, kte, 3) :: &
165	cloud, & !holder for cloud water and ice (q in CAM)
166	cloudtnd, & !holder for cloud tendencies (ptend_loc%q in CAM)
167	fracis !fraction of cloud species that are insoluble
168
169	REAL(r8), DIMENSION(pcols, kte, 2) :: &
170	icwu, & !in-cloud winds in upraft (m/s)
171	icwd, & !in-cloud winds in downdraft (m/s)
172	pguall, & !apparent force from updraft pres. gradient (~units?)
173	pgdall, & !apparent force from downdraft pres. gradient (~units?)
174	winds, & !wind components (m/s)
175	wind_tends !wind component tendencies (m/s2)
176
177	REAL(r8), DIMENSION(pcols, kte) :: &
178	cld8, & !cloud fraction
179	cme, & !cmf condensation - evaporation (~units?, sub arg out in CAM)
180	dlf, & !scattered version of the detraining cld h2o tendency (~units?)
181	fake_dpdry, & !place holder for dpdry, delta pres. of dry atmos.
182	flxprec, & !evap outfld: convective-scale flux of precip at interfaces (kg/m2/s)
183	flxsnow, & !evap outfld: convective-scale flux of snow at interfaces (kg/m2/s)
184	ntprprd, & !evap outfld: net precip production in layer (kg/kg/s)
185	ntsnprd, & !evap outfld: net snow production in layer (kg/kg/s)
186	pdel8, & !pressure thickness of layer (between interfaces, Pa)
187	pmid8, & !pressure at layer middle (Pa)
188	ql8, & !cloud liquid water (~units?)
189	qi8, & !cloud ice (~units?)
190	t8, & !temperature (K)
191	zm8, & !height above ground at mid-level (m)
192	qtnd, & !specific humidity tendency (kg/kg/s)
193	rprd, & !rain production rate (kg/kg/s, comes from pbuf array in CAM, add to restart?~)
194	stnd, & !heating rate (dry static energy tendency, W/kg)
195	tend_s_snwprd, & !heating rate of snow production (~units?)
196	tend_s_snwevmlt, & !heating rate of evap/melting of snow (~units?)
197	zdu !detraining mass flux (~units? sub arg out in CAM)
198
199	REAL(r8), DIMENSION(pcols) :: &
200	cape, & !convective available potential energy (J/kg)
201	jctop, & !row of top-of-deep-convection indices passed out (sub arg out in CAM)
202	jcbot, & !row of base of cloud indices passed out (sub arg out in CAM)
203	landfrac, & !land fraction
204	pblh8, & !planetary boundary layer height (m)
205	phis, & !geopotential at terrain height (m2/s2)
206	prec, & !convective-scale precipitation rate from ZM (m/s)
207	rliq, & !reserved liquid (not yet in cldliq) for energy integrals (units? sub arg out in CAM)
208	snow, & !convective-scale snowfall rate from ZM (m/s)
209	tpert !thermal (convective) temperature excess (K)
210
211	!Variables that were declared at the module level of zm_conv_intr in
212	!CAM. In that context they needed to be held between calls to
213	!zm_conv_tend and zm_conv_tend2 at the chunk level. In WRF, these vars
214	!are setup to hold the whole "memory" dimension, but as a 1-D vector
215	!instead of a 2-D array as is typically done in WRF. This allows us to
216	!leave the CAM routines in tact. For now, this forces us to immediately
217	!call zm_conv_tend2 before leaving this module, but it allows us to
218	!follow the WRF rules. We can deal with generalizing this for handling
219	!tracer convective transport of aerosols later.~
220	REAL(r8), DIMENSION(pcols, kte, (ime-ims+1)*(jme-jms+1)) :: &
221	dp, & !layer pres. thickness between interfaces (mb)
222	du, &
223	ed, &
224	eu, &
225	md, &
226	mu
227
228	REAL(r8), DIMENSION(pcols, (ime-ims+1)*(jme-jms+1)) :: &
229	dsubcld ! layer pres. thickness between LCL and maxi (mb)
230
231	INTEGER, DIMENSION(pcols, (ime-ims+1)*(jme-jms+1)) :: &
232	ideep, & ! holds position of gathered points
233	jt, & ! top-level index of deep cumulus convection
234	maxg ! gathered values of maxi
235
236	INTEGER, DIMENSION((ime-ims+1)*(jme-jms+1)) :: &
237	lengath ! number of gathered points
238
239	!Other local vars
240	INTEGER :: i, j, k, kflip, n, ncnst
241	INTEGER :: lchnk !chunk identifier, used to map 2-D to 1-D arrays in WRF
242	INTEGER :: ncol !number of atmospheric columns in chunk
243	LOGICAL, DIMENSION(3) :: l_qt !logical switches for constituent tendency presence
244	LOGICAL, DIMENSION(2) :: l_windt !logical switches for wind tendency presence
245	LOGICAL :: run_param !flag for handling alternate cumulus call freq.
246	!
247	! Check to see if this is a convection timestep...
248	!
249	if (adapt_step_flag) then
250	if ( (itimestep==0) .or. (cudt<=0) .or. &
251	( curr_secs+dt >= ( int(curr_secs/( cudt60 )) + 1 )cudt*60 ) ) then
252	run_param = .TRUE.
253	else
254	run_param = .FALSE.
255	endif
256
257	else
258	if (mod(itimestep,stepcu)==0 .or. itimestep==0) then
259	run_param = .TRUE.
260	else
261	run_param = .FALSE.
262	endif
263	endif
264
265	!Leave the subroutine if it is not yet time to call the cumulus param
266	if( .not. run_param ) return
267	!
268	! Initialize...
269	!
270	ncol = 1 !chunk size in WRF is 1 since we loop over all columns in a tile
271
272	cape_out(its:ite, jts:jte) = 0.
273	mu_out(its:ite, kts:kte, jts:jte) = 0.
274	md_out(its:ite, kts:kte, jts:jte) = 0.
275	zmdt(its:ite, kts:kte, jts:jte) = 0.
276	!
277	! Map variables to inputs for zm_convr and call it...
278	! Loop over the points in the tile and treat them each as a CAM chunk.
279	!
280	do j = jts,jte
281	do i = its,ite
282	lchnk = (j-jts)*(ite-its+1) + (i-its+1) !1-D index location from the 2-D tile
283
284	!Flip variables on the layer middles
285	do k = kts,kte
286	kflip = kte-k+1
287
288	cld8(1,kflip) = cldfra(i,k,j)
289	pdel8(1,kflip) = p8w(i,k,j) - p8w(i,k+1,j)
290	pmid8(1,kflip) = p(i,k,j)
291	qh8(1,kflip,1) = max( qv(i,k,j)/(1.+qv(i,k,j)), 1e-30 ) !values of 0 cause a crash in entropy
292	if( present(qc) ) then
293	ql8(1,kflip) = qc(i,k,j)/(1.+qv(i,k,j)) !Convert to moist mix ratio
294	else
295	ql8(1,kflip) = 0.
296	end if
297	if( present(qi) ) then
298	qi8(1,kflip) = qi(i,k,j)/(1.+qv(i,k,j)) !Used in convtran, ditto for conversion
299	else
300	qi8(1,kflip) = 0.
301	end if
302	t8(1,kflip) = t_phy(i,k,j)
303	zm8(1,kflip) = z(i,k,j) - ht(i,j) !Height above the ground at midlevels
304	end do
305
306	!Flip variables on the layer interfaces
307	do k = kts,kte+1
308	kflip = kte-k+2
309
310	pint8(1,kflip) = p8w(i,k,j)
311	zi8(1,kflip) = z_at_w(i,k,j) -ht(i,j) !Height above the ground at interfaces
312	end do
313
314	!Other necessary conversions for input to ZM
315	if( xland(i,j)==2 ) then
316	landfrac(1) = 1. !land, WRF is all or nothing
317	else
318	landfrac(1) = 0. !water
319	end if
320	pblh8(1) = pblh(i,j)
321	phis(1) = ht(i,j)*gravit
322
323	call get_tpert(bl_pbl_physics, sf_sfclay_physics, dx, &
324	mavail(i,j), kpbl(i,j), pblh(i,j), &
325	dz8w(i,1,j), psfc(i,j), qv(i,1,j), t_phy(i,1,j), &
326	th(i,1,j), tsk(i,j), tke_pbl(i,:,j), ust(i,j), &
327	u_phy(i,1,j), v_phy(i,1,j), hfx(i,j), qfx(i,j), &
328	tpert_camuwpbl(i,j), kte, &
329	tpert(1))
330
331	!Call the Zhang-McFarlane (1996) convection parameterization
332	!NOTE: The 0.5*dt is correct and is a nuance of CAM typically
333	! using 2*dt for physics tendencies. Everywhere in zm_convr
334	! that dt is used, the dt is multiplied by 2 to get back to
335	! the 2*dt. Everywhere else in the CAM ZM interface the full
336	! 2dt is passed into the subroutines. In WRF we use 1dt
337	! in place of CAM's 2*dt since the adjustment is made
338	! elsewhere.
339	call zm_convr( lchnk, ncol, &
340	t8, qh8, prec, jctop, jcbot, &
341	pblh8, zm8, phis, zi8, qtnd, &
342	stnd, pmid8, pint8, pdel8, &
343	0.5_r8*real(dt,r8), mcon, cme, cape, &
344	tpert, dlf, pflx, zdu, rprd, &
345	mu(:,:,lchnk), md(:,:,lchnk),du(:,:,lchnk),eu(:,:,lchnk),ed(:,:,lchnk), &
346	dp(:,:,lchnk), dsubcld(:,lchnk), jt(:,lchnk), maxg(:,lchnk), ideep(:,lchnk), &
347	lengath(lchnk), ql8, rliq, landfrac )
348
349	!Start mapping CAM output to WRF output variables. We follow the
350	!same order as in CAM's zm_conv_tend to ease maintenance...
351	do k=kts,kte
352	kflip = kte-k+1
353	dlf_out(i,k,j) = dlf(1,kflip)
354	end do
355	cape_out(i,j) = cape(1)
356	rliq_out(i,j) = rliq(1)
357
358	!Convert mass flux from reported mb/s to kg/m2/s
359	mcon(:ncol,:pver) = mcon(:ncol,:pver) * 100._r8/gravit
360
361	! Store upward and downward mass fluxes in un-gathered arrays
362	! + convert from mb/s to kg/m2/s
363	do n=1,lengath(lchnk)
364	do k=kts,kte
365	kflip = kte-k+1
366	! ii = ideep(n,lchnk) <--in WRF this is always 1 because chunk size=1
367	mu_out(i,k,j) = mu(n,kflip,lchnk) * 100._r8/gravit
368	md_out(i,k,j) = md(n,kflip,lchnk) * 100._r8/gravit
369	end do
370	end do
371
372	do k=kts,kte
373	kflip = kte-k+1
374	zmdt(i,k,j) = stnd(1,kflip)/cpair
375	zmdq(i,k,j) = qtnd(1,kflip)
376	end do
377
378	!Top and bottom pressure of convection
379	pconvt(i,j) = p8w(i,1,j)
380	pconvb(i,j) = p8w(i,1,j)
381	do n = 1,lengath(lchnk)
382	if (maxg(n,lchnk).gt.jt(n,lchnk)) then
383	pconvt(i,j) = pmid8(ideep(n,lchnk),jt(n,lchnk)) ! gathered array (or jctop ungathered)
384	pconvb(i,j) = pmid8(ideep(n,lchnk),maxg(n,lchnk))! gathered array
385	endif
386	end do
387	cutop(i,j) = jctop(1)
388	cubot(i,j) = jcbot(1)
389
390	!Add tendency from this process to tendencies arrays. Also,
391	!increment the local state arrays to reflect additional tendency
392	!from zm_convr, i.e. the physics update call in CAM. Note that
393	!we are not readjusting the pressure levels to hydrostatic
394	!balance for the new virtual temperature like is done in CAM. We
395	!will let WRF take care of such things later for the total
396	!tendency.
397	do k=kts,kte
398	kflip = kte-k+1
399
400	!Convert temperature to potential temperature and
401	!specific humidity to water vapor mixing ratio
402	rthcuten(i,k,j) = zmdt(i,k,j)/pi_phy(i,k,j)
403	rqvcuten(i,k,j) = zmdq(i,k,j)/(1._r8 - zmdq(i,k,j))
404
405	t8(1,kflip) = t8(1,kflip) + zmdt(i,k,j)*real(dt,r8)
406	qh8(1,kflip,1) = qh8(1,kflip,1) + zmdq(i,k,j)*real(dt,r8)
407	end do
408
409	! Determine the phase of the precipitation produced and add latent heat of fusion
410	! Evaporate some of the precip directly into the environment (Sundqvist)
411	! Allow this to use the updated state1 (t8 and qh8 in WRF) and the fresh ptend_loc type
412	! heating and specific humidity tendencies produced
413	qtnd = 0._r8 !re-initialize tendencies (i.e. physics_ptend_init(ptend_loc))
414	stnd = 0._r8
415	call zm_conv_evap(ncol, lchnk, &
416	t8, pmid8, pdel8, qh8(:,:,1), &
417	stnd, tend_s_snwprd, tend_s_snwevmlt, qtnd, &
418	rprd, cld8, real(dt,r8), &
419	prec, snow, ntprprd, ntsnprd , flxprec, flxsnow)
420
421	! Parse output variables from zm_conv_evap
422	do k=kts,kte
423	kflip = kte-k+1
424
425	evaptzm(i,k,j) = stnd(1,kflip)/cpair
426	fzsntzm(i,k,j) = tend_s_snwprd(1,kflip)/cpair
427	evsntzm(i,k,j) = tend_s_snwevmlt(1,kflip)/cpair
428	evapqzm(i,k,j) = qtnd(1,kflip)
429	zmflxprc(i,k,j) = flxprec(1,kflip)
430	zmflxsnw(i,k,j) = flxsnow(1,kflip)
431	zmntprpd(i,k,j) = ntprprd(1,kflip)
432	zmntsnpd(i,k,j) = ntsnprd(1,kflip)
433	zmeiheat(i,k,j) = stnd(1,kflip) !Do we really need this and evaptzm?
434	cmfmc(i,k,j) = mcon(1,kflip) !Set to deep value here, shallow added in UW scheme
435	cmfmcdzm(i,k,j) = mcon(1,kflip)
436	preccdzm(i,j) = prec(1) !Rain rate from just deep
437	precz(i,j) = prec(1) !Rain rate for total convection (just deep right now)
438	pratec(i,j) = prec(1)*1e3 !Rain rate used in WRF for accumulation (mm/s)
439	raincv(i,j) = pratec(i,j)*dt !Rain amount for time step returned back to WRF
440	end do
441
442	!Add tendency from zm_conv_evap to tendencies arrays. Also,
443	!increment the local state arrays to reflect additional tendency
444	!Note that we are not readjusting the pressure levels to hydrostatic
445	!balance for the new virtual temperature like is done in CAM. We
446	!will let WRF take care of such things later for the total
447	!tendency.
448	do k=kts,kte
449	kflip = kte-k+1
450
451	!Convert temperature to potential temperature and
452	!specific humidity to water vapor mixing ratio
453	rthcuten(i,k,j) = rthcuten(i,k,j) + &
454	evaptzm(i,k,j)/pi_phy(i,k,j)
455	rqvcuten(i,k,j) = rqvcuten(i,k,j) + &
456	evapqzm(i,k,j)/(1. - qv(i,k,j))
457
458	t8(1,kflip) = t8(1,kflip) + evaptzm(i,k,j)*real(dt,r8)
459	qh8(1,kflip,1) = qh8(1,kflip,1) + evapqzm(i,k,j)*real(dt,r8)
460	end do
461
462	! Momentum transport
463	stnd = 0._r8 !Zero relevant tendencies in preparation
464	wind_tends = 0._r8
465	do k=kts,kte
466	kflip = kte-k+1
467	winds(1,k,1) = u_phy(i,kflip,j)
468	winds(1,k,2) = v_phy(i,kflip,j)
469	end do
470	l_windt(1:2) = .true.
471
472	call momtran (lchnk, ncol, &
473	l_windt, winds, 2, mu(:,:,lchnk), md(:,:,lchnk), &
474	du(:,:,lchnk), eu(:,:,lchnk), ed(:,:,lchnk), dp(:,:,lchnk), dsubcld(:,lchnk), &
475	jt(:,lchnk),maxg(:,lchnk), ideep(:,lchnk), 1, lengath(lchnk), &
476	itimestep, wind_tends, pguall, pgdall, icwu, icwd, real(dt,r8), stnd )
477
478	!Add tendency from momtran to tendencies arrays. Also,
479	!increment the local state arrays to reflect additional tendency
480	!Note that we are not readjusting the pressure levels to hydrostatic
481	!balance for the new virtual temperature like is done in CAM. We
482	!will let WRF take care of such things later for the total
483	!tendency.
484	do k=kts,kte
485	kflip = kte-k+1
486
487	!Convert temperature to potential temperature and
488	!specific humidity to water vapor mixing ratio
489	rucuten(i,k,j) = wind_tends(1,kflip,1)
490	rvcuten(i,k,j) = wind_tends(1,kflip,2)
491	rthcuten(i,k,j) = rthcuten(i,k,j) + &
492	stnd(1,kflip)/cpair/pi_phy(i,k,j)
493
494	t8(1,kflip) = t8(1,kflip) + stnd(1,kflip)/cpair*real(dt,r8)
495	!winds is not used again so do not bother updating them
496	end do
497
498	!Parse output arrays for momtran
499	do k=kts,kte
500	kflip = kte-k+1
501
502	zmmtu(i,k,j) = wind_tends(1,kflip,1) !wind tendencies...
503	zmmtv(i,k,j) = wind_tends(1,kflip,2)
504
505	zmupgu(i,k,j) = pguall(1,kflip,1) !apparent force pres. grads...
506	zmupgd(i,k,j) = pgdall(1,kflip,1)
507	zmvpgu(i,k,j) = pguall(1,kflip,2)
508	zmvpgd(i,k,j) = pgdall(1,kflip,2)
509
510	zmicuu(i,k,j) = icwu(1,kflip,1) !in-cloud winds...
511	zmicud(i,k,j) = icwd(1,kflip,1)
512	zmicvu(i,k,j) = icwu(1,kflip,2)
513	zmicvd(i,k,j) = icwd(1,kflip,2)
514	end do
515
516	!Setup for convective transport of cloud water and ice
517	!~We should set this up to use an integer pointer instead of
518	! hard-coded numbers for each species so that we can easily
519	! handle other species. Something to come back to later...
520	l_qt(1) = .false. !do not mix water vapor
521	l_qt(2:3) = .true. !do mix cloud water and ice
522	cloudtnd = 0._r8
523	fracis(1,:,1:3) = 0._r8 !all cloud liquid & ice is soluble
524	ncnst = 3 !number of constituents in cloud array (including vapor)
525	fake_dpdry = 0._r8 !delta pres. for dry atmos. (0 since assuming moist mix ratios)
526	do k=kts,kte
527	kflip = kte-k+1
528
529	cloud(1,kflip,1) = qh8(1,kflip,1) !We can either use moist mix ratios, as is
530	cloud(1,kflip,2) = ql8(1,kflip) !done here, or else use dry mix ratios, send
531	cloud(1,kflip,3) = qi8(1,kflip) !in appropriate dpdry values, and return the
532	!approp. response from cnst_get_type_byind
533	end do
534
535	call convtran (lchnk, &
536	l_qt, cloud, ncnst, mu(:,:,lchnk), md(:,:,lchnk), &
537	du(:,:,lchnk), eu(:,:,lchnk), ed(:,:,lchnk), dp(:,:,lchnk), dsubcld(:,lchnk), &
538	jt(:,lchnk), maxg(:,lchnk), ideep(:,lchnk), 1, lengath(lchnk), &
539	itimestep, fracis, cloudtnd, fake_dpdry)
540
541	!Parse output for convtran and increment tendencies
542	do k=kts,kte
543	kflip = kte-k+1
544
545	zmdice(i,k,j) = cloudtnd(1,kflip,3)
546	zmdliq(i,k,j) = cloudtnd(1,kflip,2)
547
548	!Convert cloud tendencies from wet to dry mix ratios
549	if( present(rqccuten) ) then
550	rqccuten(i,k,j) = cloudtnd(1,kflip,2)/(1. - qv(i,k,j))
551	end if
552	if( present(rqicuten) ) then
553	rqicuten(i,k,j) = cloudtnd(1,kflip,3)/(1. - qv(i,k,j))
554	end if
555	end do
556
557	end do !i-loop
558	end do !j-loop
559	END SUBROUTINE camzm_driver
560
561
562	!------------------------------------------------------------------------
563	SUBROUTINE zm_conv_init(rucuten, rvcuten, rthcuten, rqvcuten, &
564	rqccuten, rqicuten, &
565	p_qc, p_qi, param_first_scalar, &
566	restart, &
567	ids, ide, jds, jde, kds, kde, &
568	ims, ime, jms, jme, kms, kme, &
569	its, ite, jts, jte, kts, kte )
570	! Initialization routine for Zhang-McFarlane cumulus parameterization
571	! from CAM. The routine with this name in CAM is revamped here to give
572	! the needed functionality within WRF.
573	!
574	! Author: William.Gustafson@pnl.gov, Nov. 2009
575	!------------------------------------------------------------------------
576	USE module_cam_esinti, only: esinti
577	USE physconst, only: epsilo, latvap, latice, rh2o, cpair, tmelt
578	USE module_bl_camuwpbl_driver, only: vd_register
579	USE module_cu_camzm, only: zm_convi, zmconv_readnl
580
581	LOGICAL , INTENT(IN) :: restart
582	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
583	ims, ime, jms, jme, kms, kme, &
584	its, ite, jts, jte, kts, kte, &
585	p_qc, p_qi, param_first_scalar
586
587	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
588	rucuten, &
589	rvcuten, &
590	rthcuten, &
591	rqvcuten, &
592	rqccuten, &
593	rqicuten
594
595	integer :: i, itf, j, jtf, k, ktf
596	integer :: limcnv
597
598	jtf = min(jte,jde-1)
599	ktf = min(kte,kde-1)
600	itf = min(ite,ide-1)
601	!
602	! Initialize module_cam_support variables...
603	! This could be moved to a master "cam-init" subroutine once multiple CAM
604	! parameterizations are implemented in WRF. For now, it doesn't hurt to
605	! have these possibly initialized more than once since they are
606	! essentially constant.
607	!
608	pver = ktf-kts+1
609	pverp = pver+1
610	!
611	! Initialize the saturation vapor pressure look-up table...
612	! This could be moved to a master cam-init subroutine too if it is needed
613	! by more than one CAM parameterization. In CAM this is called from
614	! phys_init.
615	!
616	call esinti(epsilo, latvap, latice, rh2o, cpair, tmelt)
617
618	!~Need code here to set limcnv to max convective level of 40 mb... To
619	! properly set an average value for the whole domain would involve doing
620	! a collective operation across the whole domain to get pressure averages
621	! by level, which is difficult within the WRF framework. So, we will just
622	! use the model top for now.
623	!
624	! Technically, limcnv should be calculated for each grid point at each
625	! time because the levels in WRF do not have a constant pressure, as
626	! opposed to CAM. But, that would involve making this variable local
627	! instead of at the module level as is currently done in CAM. For now,
628	! we will not worry about this because WRF rarely has domains higher than
629	! 40 mb. There is also little variability between grid points near the
630	! model top due to the underlying topography so a constant value would
631	! be OK across the comain. Also, remember that CAM levels are reversed in
632	! the vertical from WRF. So, setting limcnv to 1 means the top of the
633	! domain. Note that because of a bug in the parcel_dilute routine, limcnv
634	! must be >=2 instead of 1. Otherwise, an array out-of-bounds occurs.
635	limcnv = 2
636	!
637	! Get the ZM namelist variables (hard-wired for now)...
638	!
639	call zmconv_readnl("hard-wired")
640	!
641	!~need to determine if convection should happen in PBL and set
642	! no_deep_pbl_in accordingly
643	call zm_convi(limcnv, NO_DEEP_PBL_IN=.false.)
644
645	!
646	! Set initial values for arrays dependent on restart conditions
647	!
648	if(.not.restart)then
649	do j=jts,jtf
650	do k=kts,ktf
651	do i=its,itf
652	rucuten(i,k,j) = 0.
653	rvcuten(i,k,j) = 0.
654	rthcuten(i,k,j) = 0.
655	rqvcuten(i,k,j) = 0.
656	if( p_qc > param_first_scalar ) rqccuten(i,k,j) = 0.
657	if( p_qi > param_first_scalar ) rqicuten(i,k,j) = 0.
658	enddo
659	enddo
660	enddo
661	end if
662
663	END SUBROUTINE zm_conv_init
664
665
666	!------------------------------------------------------------------------
667	SUBROUTINE get_tpert(bl_pbl_physics, sf_sfclay_physics, dx, &
668	mavail, kpbl, pblh, dzlowest, &
669	psfc, qvlowest, t_phylowest, thlowest, tsk, tke_pbl, ust, &
670	u_phylowest, v_phylowest, hfx, qfx, tpert_camuwpbl, kte, tpert)
671	! Calculates the temperature perturbation used to trigger convection.
672	! This should only be used if tpert is not already provided by CAM's PBL
673	! scheme. Right now, this only works with the Mellor-Yamada boundary
674	! layer scheme in combination with the MYJ surface scheme. This is due to
675	! the need of TKE for the vertical velocity perturbation. This scheme has
676	! not been generalized to handle all schemes that produce TKE at this
677	! time, and the non-TKE schemes, e.g. YSU, could probably have an
678	! appropriate tpert deduced but we have not taken the time yet to do it.
679	!
680	! Author: William.Gustafson@pnl.gov, Nov. 2009
681	! Last updated: William.Gustafson@pnl.gov, Nov. 2010
682	!------------------------------------------------------------------------
683	USE shr_kind_mod, only: r8 => shr_kind_r8
684	USE module_model_constants, only: cp, ep_1, ep_2, g, r_d, rcp, &
685	svp1, svp2, svp3, svpt0, xlv
686	USE module_state_description, ONLY : SFCLAYSCHEME &
687	,MYJSFCSCHEME &
688	,GFSSFCSCHEME &
689	,SLABSCHEME &
690	,LSMSCHEME &
691	,RUCLSMSCHEME &
692	,MYJPBLSCHEME &
693	,CAMUWPBLSCHEME
694	!
695	! Subroutine arguments...
696	!
697	real, dimension(:), intent(in) :: tke_pbl
698	real, intent(in) :: dx, dzlowest, hfx, mavail, pblh, psfc, qvlowest, &
699	tpert_camuwpbl, tsk, t_phylowest, thlowest, ust, u_phylowest, &
700	v_phylowest
701	integer, intent(in) :: bl_pbl_physics, kpbl, kte, sf_sfclay_physics
702	real(r8),intent(inout) :: tpert
703	!
704	! Local vars...
705	!
706	real, parameter :: fak = 8.5 !Constant in surface temperature excess
707	real, parameter :: tfac = 1. !Ratio of 'tpert' to (w't')/wpert
708	real, parameter :: wfac = 1. !Ratio of 'wpert' to sqrt(tke)
709	real, parameter :: wpertmin = 1.e-6 !Min PBL eddy vertical velocity perturbation
710	real :: ebrk !In CAM, net mean TKE of CL including
711	!entrainment effect (m2/s2). In WRF,
712	!average TKE within the PBL
713	real :: br2, dthvdz, e1, flux, govrth, psfccmb, qfx, qsfc, rhox, thgb, &
714	thv, tskv, tvcon, vconv, vsgd, wpert, wspd, za
715	integer :: k
716	character(len=250) :: msg
717	logical :: UnstableOrNeutral
718	!
719	! We can get the perturbation values directly from CAMUWPBLSCHEME if
720	! available. Otherwise, we have to calculate them.
721	!
722	if( bl_pbl_physics==CAMUWPBLSCHEME ) then
723	tpert = tpert_camuwpbl
724	!
725	!...non-CAMUWPBL. Need MYJ SFC & PBL for now until other schemes
726	! get coded to give perturbations too.
727	! First, we need to determine if the conditions are stable or unstable.
728	! We will do this by mimicing the bulk Richardson calculation from
729	! module_sf_sfclay.F because the MYJ scheme does not provide this info.
730	! The logic borrowed from the CuP shallow cumulus scheme. Non-MYJ sfc
731	! scheme code is commented out for now since we also require MYJ PBL
732	! scheme.
733	!
734	elseif( bl_pbl_physics==MYJPBLSCHEME ) then
735
736	UnstableOrNeutral = .false.
737	sfclay_case: SELECT CASE (sf_sfclay_physics)
738	CASE (MYJSFCSCHEME)
739	! The MYJ sfc scheme does not already provide a stability criteria.
740	! So, we will mimic the bulk Richardson calculation from
741	! module_sf_sfclay.F.
742
743	if( pblh <= 0. ) call wrf_error_fatal( &
744	"CAMZMSCHEME needs a PBL height from a PBL scheme.")
745
746	za = 0.5*dzlowest
747
748	e1 = svp1exp(svp2(tsk-svpt0)/(tsk-svp3))
749	psfccmb=psfc/1000. !converts from Pa to cmb
750	qsfc = ep_2*e1/(psfccmb-e1)
751	thgb = tsk(100./psfccmb)*rcp
752	tskv = thgb(1.+ep_1qsfc*mavail)
753	tvcon = 1.+ep_1*qvlowest
754	thv = thlowest*tvcon
755	dthvdz = (thv-tskv)
756
757	govrth = g/thlowest
758
759	rhox = psfc/(r_dt_phylowesttvcon)
760	flux = max(hfx/rhox/cp + ep_1tskvqfx/rhox,0.)
761	vconv = (g/tskpblhflux)**.33
762	vsgd = 0.32 * (max(dx/5000.-1.,0.))**.33
763	wspd = sqrt(u_phylowestu_phylowest+v_phylowestv_phylowest)
764	wspd = sqrt(wspdwspd+vconvvconv+vsgd*vsgd)
765	wspd = max(wspd,0.1)
766
767	!And finally, the bulk Richardson number...
768	br2 = govrthzadthvdz/(wspd*wspd)
769
770	if( br2 <= 0. ) UnstableOrNeutral = .true.
771
772	CASE DEFAULT
773	call wrf_error_fatal("CAMZMSCHEME requires MYJSFCSCHEME or else CAMUWPBLSCHEME.")
774
775	END SELECT sfclay_case
776	!
777	! The perturbation temperature for unstable conditions...
778	! The calculation follows the one in caleddy inside eddy_diff.F90 from
779	! CAM.
780	!
781	!Check that we are using the MJY BL scheme since we are hard-wired to
782	!use TKE and u* from this scheme. At some point this dependency should
783	!be broken and a way needs to be found for other schemes to provide
784	!reasonable tpert values too.
785	if( bl_pbl_physics /= MYJPBLSCHEME ) &
786	call wrf_error_fatal("CAMZMSCHEME requires MYJPBLSCHEME or CAMUWPBLSCHEME")
787
788
789	!Recalculate rhox in case a non-MYJ scheme is used to get
790	!stability and rhox is not calculated above. Right now, this is
791	!technically reduncant, but cheap.
792	tvcon = 1.+ep_1*qvlowest
793	rhox = psfc/(r_dt_phylowesttvcon)
794
795	if( UnstableOrNeutral ) then
796	!1st, get avg TKE w/n the PBL as a proxy for ebrk variable in CAM
797	ebrk = 0.
798	do k=1,min(kpbl,kte)
799	ebrk = ebrk + tke_pbl(k)
800	end do
801	ebrk = ebrk/real(kpbl)
802
803	wpert = max( wfac*sqrt(ebrk), wpertmin )
804	tpert = max( abs(hfx/rhox/cp)*tfac/wpert, 0. )
805	!
806	! Or, the perturbation temperature for stable conditions...
807	!
808	else !Stable conditions
809	tpert = max( hfx/rhox/cp*fak/ust, 0. )
810	end if
811
812	else
813	call wrf_error_fatal("CAMZMSCHEME requires MYJPBLSCHEME or CAMUWPBLSCHEME")
814
815	end if !PBL choice
816
817	END SUBROUTINE get_tpert
818
819	END MODULE module_cu_camzm_driver

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

Download in other formats:

Original Format