Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_mynn.F @ 1314

Last change on this file since 1314 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: 53.7 KB

Line
1	! translated from NN f77 to F90 and put into WRF by Mariusz Pagowski
2	! NOAA/GSD & CIRA/CSU, Feb 2008
3	! changes to original code:
4	! 1. code is 1d (in z)
5	! 2. no advection of TKE, covariances and variances
6	! 3. Cranck-Nicholson replaced with the implicit scheme
7	! 4. removed terrain dependent grid since input in WRF in actual
8	! distances in z[m]
9	! 5. cosmetic changes to adhere to WRF standard (remove common blocks,
10	! intent etc)
11
12	MODULE module_bl_mynn
13
14	USE module_model_constants, only: &
15	&karman, g, p1000mb, &
16	&cp, r_d, rcp, xlv, &
17	&svp1, svp2, svp3, svpt0, ep_1, ep_2
18
19	IMPLICIT NONE
20
21	! The parameters below depend on stability functions of module_sf_mynn.
22	REAL, PARAMETER :: cphm_st=5.0, cphm_unst=16.0, &
23	cphh_st=5.0, cphh_unst=16.0
24
25	REAL, PARAMETER :: xlvcp=xlv/cp, ev=xlv, rd=r_d, rk=cp/rd, &
26	&svp11=svp1*1.e3, p608=ep_1, ep_3=1.-ep_2
27
28	REAL, PARAMETER :: tref=300.0 ! reference temperature (K)
29	REAL, PARAMETER :: tv0=p608tref, tv1=(1.+p608)tref, gtr=g/tref
30
31	! Closure constants
32	REAL, PARAMETER :: &
33	&vk = karman, &
34	&pr = 0.74, &
35	&g1 = 0.235, &
36	&b1 = 24.0, &
37	&b2 = 15.0, &
38	&c2 = 0.75, &
39	&c3 = 0.352, &
40	&c4 = 0.0, &
41	&c5 = 0.2, &
42	&a1 = b1( 1.0-3.0g1 )/6.0, &
43	! &c1 = g1 -1.0/( 3.0a1b1**(1.0/3.0) ), &
44	&c1 = g1 -1.0/( 3.0a12.88449914061481660), &
45	&a2 = a1( g1-c1 )/( g1pr ), &
46	&g2 = b2/b1( 1.0-c3 ) +2.0a1/b1( 3.0-2.0c2 )
47
48	REAL, PARAMETER :: &
49	&cc2 = 1.0-c2, &
50	&cc3 = 1.0-c3, &
51	&e1c = 3.0a2b2*cc3, &
52	&e2c = 9.0a1a2*cc2, &
53	&e3c = 9.0a2a2cc2( 1.0-c5 ), &
54	&e4c = 12.0a1a2*cc2, &
55	&e5c = 6.0a1a1
56
57	! Constants for length scale
58	REAL, PARAMETER :: qmin=0.0, zmax=1.0, cns=2.7, &
59	&alp1=0.23, alp2=1.0, alp3=5.0, alp4=100.0
60
61	! Constants for gravitational settling
62	! REAL, PARAMETER :: gno=1.e6/(1.e8)**(2./3.), gpw=5./3., qcgmin=1.e-8
63	REAL, PARAMETER :: gno=4.64158883361278196
64	REAL, PARAMETER :: gpw=5./3., qcgmin=1.e-8,qkemin=1.e-12
65	! REAL, PARAMETER :: pblh_ref=1500.
66
67	REAL, PARAMETER :: rr2=0.7071068, rrp=0.3989423
68
69	INTEGER :: mynn_level
70
71	CONTAINS
72
73	! **********************************************************************
74	! * An improved Mellor-Yamada turbulence closure model *
75	! * *
76	! * Aug/2005 M. Nakanishi (N.D.A) *
77	! * Modified: Dec/2005 M. Nakanishi (N.D.A) *
78	! * naka@nda.ac.jp *
79	! * *
80	! * Contents: *
81	! * 1. mym_initialize (to be called once initially) *
82	! * gives the closure constants and initializes the turbulent *
83	! * quantities. *
84	! * (2) mym_level2 (called in the other subroutines) *
85	! * calculates the stability functions at Level 2. *
86	! * (3) mym_length (called in the other subroutines) *
87	! * calculates the master length scale. *
88	! * 4. mym_turbulence *
89	! * calculates the vertical diffusivity coefficients and the *
90	! * production terms for the turbulent quantities. *
91	! * 5. mym_predict *
92	! * predicts the turbulent quantities at the next step. *
93	! * 6. mym_condensation *
94	! * determines the liquid water content and the cloud fraction *
95	! * diagnostically. *
96	! * *
97	! * call mym_initialize *
98	! * \| *
99	! * \|<----------------+ *
100	! * \| \| *
101	! * call mym_condensation \| *
102	! * call mym_turbulence \| *
103	! * call mym_predict \| *
104	! * \| \| *
105	! * \|-----------------+ *
106	! * \| *
107	! * end *
108	! * *
109	! * Variables worthy of special mention: *
110	! * tref : Reference temperature *
111	! * thl : Liquid water potential temperature *
112	! * qw : Total water (water vapor+liquid water) content *
113	! * ql : Liquid water content *
114	! * vt, vq : Functions for computing the buoyancy flux *
115	! * *
116	! * If the water contents are unnecessary, e.g., in the case of *
117	! * ocean models, thl is the potential temperature and qw, ql, vt *
118	! * and vq are all zero. *
119	! * *
120	! * Grid arrangement: *
121	! * k+1 +---------+ *
122	! * \| \| i = 1 - nx *
123	! * (k) \| * \| j = 1 - ny *
124	! * \| \| k = 1 - nz *
125	! * k +---------+ *
126	! * i (i) i+1 *
127	! * *
128	! * All the predicted variables are defined at the center () of
129	! * the grid boxes. The diffusivity coefficients are, however, *
130	! * defined on the walls of the grid boxes. *
131	! * # Upper boundary values are given at k=nz. *
132	! * *
133	! * References: *
134	! * 1. Nakanishi, M., 2001: *
135	! * Boundary-Layer Meteor., 99, 349-378. *
136	! * 2. Nakanishi, M. and H. Niino, 2004: *
137	! * Boundary-Layer Meteor., 112, 1-31. *
138	! * 3. Nakanishi, M. and H. Niino, 2006: *
139	! * Boundary-Layer Meteor., (in press). *
140	! **********************************************************************
141	!
142	!
143	! SUBROUTINE mym_initialize:
144	!
145	! Input variables:
146	! iniflag : <>0; turbulent quantities will be initialized
147	! = 0; turbulent quantities have been already
148	! given, i.e., they will not be initialized
149	! mx, my : Maximum numbers of grid boxes
150	! in the x and y directions, respectively
151	! nx, ny, nz : Numbers of the actual grid boxes
152	! in the x, y and z directions, respectively
153	! tref : Reference temperature (K)
154	! dz(nz) : Vertical grid spacings (m)
155	! # dz(nz)=dz(nz-1)
156	! zw(nz+1) : Heights of the walls of the grid boxes (m)
157	! # zw(1)=0.0 and zw(k)=zw(k-1)+dz(k-1)
158	! h(mx,ny) : G^(1/2) in the terrain-following coordinate
159	! # h=1-zg/zt, where zg is the height of the
160	! terrain and zt the top of the model domain
161	! pi0(mx,my,nz) : Exner function at zw*h+zg (J/kg K)
162	! defined by c_p*( p_basic/1000hPa )^kappa
163	! This is usually computed by integrating
164	! d(pi0)/dz = -h*g/tref.
165	! rmo(mx,ny) : Inverse of the Obukhov length (m^(-1))
166	! flt, flq(mx,ny) : Turbulent fluxes of sensible and latent heat,
167	! respectively, e.g., flt=-u_Theta_ (K m/s)
168	!! flt - liquid water potential temperature surface flux
169	!! flq - total water flux surface flux
170	! ust(mx,ny) : Friction velocity (m/s)
171	! pmz(mx,ny) : phi_m-zeta at z1h+z0, where z1 (=0.5dz(1))
172	! is the first grid point above the surafce, z0
173	! the roughness length and zeta=(z1h+z0)rmo
174	! phh(mx,ny) : phi_h at z1*h+z0
175	! u, v(mx,my,nz): Components of the horizontal wind (m/s)
176	! thl(mx,my,nz) : Liquid water potential temperature
177	! (K)
178	! qw(mx,my,nz) : Total water content Q_w (kg/kg)
179	!
180	! Output variables:
181	! ql(mx,my,nz) : Liquid water content (kg/kg)
182	! v?(mx,my,nz) : Functions for computing the buoyancy flux
183	! qke(mx,my,nz) : Twice the turbulent kineti! energy q^2
184	! (m^2/s^2)
185	! tsq(mx,my,nz) : Variance of Theta_l (K^2)
186	! qsq(mx,my,nz) : Variance of Q_w
187	! cov(mx,my,nz) : Covariance of Theta_l and Q_w (K)
188	! el(mx,my,nz) : Master length scale L (m)
189	! defined on the walls of the grid boxes
190	! bsh : no longer used
191	! via common : Closure constants
192	!
193	! Work arrays: see subroutine mym_level2
194	! pd?(mx,my,nz) : Half of the production terms at Level 2
195	! defined on the walls of the grid boxes
196	! qkw(mx,my,nz) : q on the walls of the grid boxes (m/s)
197	!
198	! # As to dtl, ...gh, see subroutine mym_turbulence.
199	!
200
201	SUBROUTINE mym_initialize ( kts,kte,&
202	& dz, zw, &
203	& u, v, thl, qw, &
204	! & ust, rmo, pmz, phh, flt, flq,&
205	& ust, rmo, &
206	& Qke, Tsq, Qsq, Cov)
207
208	!
209
210	INTEGER, INTENT(IN) :: kts,kte
211	! REAL, INTENT(IN) :: ust, rmo, pmz, phh, flt, flq
212	REAL, INTENT(IN) :: ust, rmo
213	REAL, DIMENSION(kts:kte), INTENT(in) :: dz
214	REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
215	REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw
216
217	REAL, DIMENSION(kts:kte), INTENT(out) :: qke,tsq,qsq,cov
218
219
220	REAL, DIMENSION(kts:kte) :: &
221	&ql,el,pdk,pdt,pdq,pdc,dtl,dqw,dtv,&
222	&gm,gh,sm,sh,qkw,vt,vq
223	INTEGER :: k,l,lmax
224	REAL :: phm,vkz,elq,elv,b1l,b2l,pmz=1.,phh=1.,flt=0.,flq=0.,tmpq
225
226
227	! At first ql, vt and vq are set to zero.
228	DO k = kts,kte
229	ql(k) = 0.0
230	vt(k) = 0.0
231	vq(k) = 0.0
232	END DO
233	!
234	CALL mym_level2 ( kts,kte,&
235	& dz, &
236	& u, v, thl, qw, &
237	& ql, vt, vq, &
238	& dtl, dqw, dtv, gm, gh, sm, sh )
239	!
240	! Preliminary setting
241
242	el (kts) = 0.0
243	qke(kts) = ust*2 ( b1pmz )*(2.0/3.0)
244	!
245	phm = phhb2 / ( b1pmz )**(1.0/3.0)
246	tsq(kts) = phm( flt/ust )*2
247	qsq(kts) = phm( flq/ust )*2
248	cov(kts) = phm( flt/ust )( flq/ust )
249	!
250	DO k = kts+1,kte
251	vkz = vk*zw(k)
252	el (k) = vkz/( 1.0 + vkz/100.0 )
253	qke(k) = 0.0
254	!
255	tsq(k) = 0.0
256	qsq(k) = 0.0
257	cov(k) = 0.0
258	END DO
259	!
260	! Initialization with an iterative manner
261	! lmax is the iteration count. This is arbitrary.
262	lmax = 5 !!kte +3
263	!
264	DO l = 1,lmax
265	!
266	CALL mym_length ( kts,kte,&
267	& dz, zw, &
268	& rmo, flt, flq, &
269	& vt, vq, &
270	& qke, &
271	& dtv, &
272	& el, &
273	& qkw)
274	!
275	DO k = kts+1,kte
276	elq = el(k)*qkw(k)
277	pdk(k) = elq( sm(k)gm (k)+&
278	&sh(k)*gh (k) )
279	pdt(k) = elq* sh(k)dtl(k)*2
280	pdq(k) = elq* sh(k)dqw(k)*2
281	pdc(k) = elq* sh(k)dtl(k)dqw(k)
282	END DO
283	!
284	! ** Strictly, vkzh(i,j) -> vk( 0.5dz(1)h(i,j)+z0 ) **
285	vkz = vk0.5dz(kts)
286	!
287	elv = 0.5*( el(kts+1)+el(kts) ) / vkz
288	qke(kts) = ust*2 ( b1pmzelv )**(2.0/3.0)
289	!
290	phm = phhb2 / ( b1pmz/elv2 )(1.0/3.0)
291	tsq(kts) = phm( flt/ust )*2
292	qsq(kts) = phm( flq/ust )*2
293	cov(kts) = phm( flt/ust )( flq/ust )
294	!
295	DO k = kts+1,kte-1
296	b1l = b10.25( el(k+1)+el(k) )
297	tmpq=MAX(b1l*( pdk(k+1)+pdk(k) ),qkemin)
298	! PRINT *,'tmpqqqqq',tmpq,pdk(k+1),pdk(k)
299	qke(k) = tmpq**(2.0/3.0)
300
301	!
302	IF ( qke(k) .LE. 0.0 ) THEN
303	b2l = 0.0
304	ELSE
305	b2l = b2*( b1l/b1 ) / SQRT( qke(k) )
306	END IF
307	!
308	tsq(k) = b2l*( pdt(k+1)+pdt(k) )
309	qsq(k) = b2l*( pdq(k+1)+pdq(k) )
310	cov(k) = b2l*( pdc(k+1)+pdc(k) )
311	END DO
312
313	!
314	END DO
315
316	!! qke(kts)=qke(kts+1)
317	!! tsq(kts)=tsq(kts+1)
318	!! qsq(kts)=qsq(kts+1)
319	!! cov(kts)=cov(kts+1)
320
321	qke(kte)=qke(kte-1)
322	tsq(kte)=tsq(kte-1)
323	qsq(kte)=qsq(kte-1)
324	cov(kte)=cov(kte-1)
325
326	!
327	! RETURN
328
329	END SUBROUTINE mym_initialize
330
331	!
332	! SUBROUTINE mym_level2:
333	!
334	! Input variables: see subroutine mym_initialize
335	!
336	! Output variables:
337	! dtl(mx,my,nz) : Vertical gradient of Theta_l (K/m)
338	! dqw(mx,my,nz) : Vertical gradient of Q_w
339	! dtv(mx,my,nz) : Vertical gradient of Theta_V (K/m)
340	! gm (mx,my,nz) : G_M divided by L^2/q^2 (s^(-2))
341	! gh (mx,my,nz) : G_H divided by L^2/q^2 (s^(-2))
342	! sm (mx,my,nz) : Stability function for momentum, at Level 2
343	! sh (mx,my,nz) : Stability function for heat, at Level 2
344	!
345	! These are defined on the walls of the grid boxes.
346	!
347	SUBROUTINE mym_level2 (kts,kte,&
348	& dz, &
349	& u, v, thl, qw, &
350	& ql, vt, vq, &
351	& dtl, dqw, dtv, gm, gh, sm, sh )
352	!
353
354	INTEGER, INTENT(IN) :: kts,kte
355	REAL, DIMENSION(kts:kte), INTENT(in) :: dz
356	REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,ql,vt,vq
357
358
359	REAL, DIMENSION(kts:kte), INTENT(out) :: &
360	&dtl,dqw,dtv,gm,gh,sm,sh
361
362	INTEGER :: k
363
364	REAL :: rfc,f1,f2,rf1,rf2,smc,shc,&
365	&ri1,ri2,ri3,ri4,duz,dtz,dqz,vtt,vqq,dtq,dzk,afk,abk,ri,rf
366
367	! ev = 2.5e6
368	! tv0 = 0.61*tref
369	! tv1 = 1.61*tref
370	! gtr = 9.81/tref
371	!
372	rfc = g1/( g1+g2 )
373	f1 = b1( g1-c1 ) +3.0a2( 1.0 -c2 )( 1.0-c5 ) &
374	& +2.0a1( 3.0-2.0*c2 )
375	f2 = b1( g1+g2 ) -3.0a1*( 1.0 -c2 )
376	rf1 = b1*( g1-c1 )/f1
377	rf2 = b1* g1 /f2
378	smc = a1 /a2* f1/f2
379	shc = 3.0a2( g1+g2 )
380	!
381	ri1 = 0.5/smc
382	ri2 = rf1*smc
383	ri3 = 4.0rf2smc -2.0*ri2
384	ri4 = ri2**2
385	!
386	DO k = kts+1,kte
387	dzk = 0.5 *( dz(k)+dz(k-1) )
388	afk = dz(k)/( dz(k)+dz(k-1) )
389	abk = 1.0 -afk
390	duz = ( u(k)-u(k-1) )2 +( v(k)-v(k-1) )2
391	duz = duz /dzk**2
392	dtz = ( thl(k)-thl(k-1) )/( dzk )
393	dqz = ( qw(k)-qw(k-1) )/( dzk )
394	!
395	vtt = 1.0 +vt(k)abk +vt(k-1)afk
396	vqq = tv0 +vq(k)abk +vq(k-1)afk
397	dtq = vttdtz +vqqdqz
398	!
399	dtl(k) = dtz
400	dqw(k) = dqz
401	dtv(k) = dtq
402	!? dtv(i,j,k) = dtz +tv0*dqz
403	!? : +( ev/pi0(i,j,k)-tv1 )
404	!? : ( ql(i,j,k)-ql(i,j,k-1) )/( dzkh(i,j) )
405	!
406	gm (k) = duz
407	gh (k) = -dtq*gtr
408	!
409	! Gradient Richardson number
410	ri = -gh(k)/MAX( duz, 1.0e-10 )
411	! Flux Richardson number
412	rf = MIN( ri1( ri+ri2-SQRT(ri2-ri3ri+ri4) ), rfc )
413	!
414	sh (k) = shc*( rfc-rf )/( 1.0-rf )
415	sm (k) = smc( rf1-rf )/( rf2-rf ) sh(k)
416	END DO
417	!
418	RETURN
419
420	END SUBROUTINE mym_level2
421
422	!
423	! SUBROUTINE mym_length:
424	!
425	! Input variables: see subroutine mym_initialize
426	!
427	! Output variables: see subroutine mym_initialize
428	!
429	! Work arrays:
430	! elt(mx,ny) : Length scale depending on the PBL depth (m)
431	! vsc(mx,ny) : Velocity scale q_c (m/s)
432	! at first, used for computing elt
433	!
434	SUBROUTINE mym_length ( kts,kte,&
435	& dz, zw, &
436	& rmo, flt, flq, &
437	& vt, vq, &
438	& qke, &
439	& dtv, &
440	& el, &
441	& qkw)
442
443	INTEGER, INTENT(IN) :: kts,kte
444	REAL, DIMENSION(kts:kte), INTENT(in) :: dz
445	REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
446	REAL, INTENT(in) :: rmo,flt,flq
447	REAL, DIMENSION(kts:kte), INTENT(IN) :: qke,vt,vq
448
449
450	REAL, DIMENSION(kts:kte), INTENT(out) :: &
451	&qkw, el
452	REAL, DIMENSION(kts:kte), INTENT(in) :: dtv
453
454	REAL :: elt,vsc
455
456	INTEGER :: i,j,k
457	REAL :: afk,abk,zwk,dzk,qdz,vflx,bv,elb,els,elf
458
459	! tv0 = 0.61*tref
460	! gtr = 9.81/tref
461	!
462	DO k = kts+1,kte
463	afk = dz(k)/( dz(k)+dz(k-1) )
464	abk = 1.0 -afk
465	qkw(k) = SQRT(MAX(qke(k)abk+qke(k-1)&
466	&afk,1.0e-10))
467	END DO
468	!
469	elt = 1.0e-5
470	vsc = 1.0e-5
471	!
472	! ** Strictly, zwkh(i,j) -> ( zwkh(i,j)+z0 ) **
473	DO k = kts+1,kte-1
474	zwk = zw(k)
475	dzk = 0.5*( dz(k)+dz(k-1) )
476	qdz = MAX( qkw(k)-qmin, 0.0 )*dzk
477	elt = elt +qdz*zwk
478	vsc = vsc +qdz
479	END DO
480	!
481	vflx = ( vt(kts)+1.0 )flt +( vq(kts)+tv0 )flq
482	elt = alp1*elt/vsc
483	vsc = ( gtreltMAX( vflx, 0.0 ) )**(1.0/3.0)
484	!
485	! Strictly, el(i,j,1) is not zero.
486	el(kts) = 0.0
487	!
488	DO k = kts+1,kte
489	zwk = zw(k)
490	! Length scale limited by the buoyancy effect
491	IF ( dtv(k) .GT. 0.0 ) THEN
492	bv = SQRT( gtr*dtv(k) )
493	elb = alp2*qkw(k) / bv &
494	& ( 1.0 + alp3/alp2&
495	&SQRT( vsc/( bv*elt ) ) )
496	elf=elb
497	elf = alp2 * qkw(k)/bv
498	ELSE
499	elb = 1.0e10
500	elf =elb
501	END IF
502	!
503	! Length scale in the surface layer
504	IF ( rmo .GT. 0.0 ) THEN
505	els = vk*zwk &
506	& /( 1.0 + cnsMIN( zwkrmo, zmax ) )
507	ELSE
508	els = vk*zwk &
509	& ( 1.0 - alp4 zwkrmo )*0.2
510	END IF
511	!
512	el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf)
513	!! el(k) = elb/( elb/elt+elb/els+1.0 )
514
515	END DO
516	!
517	RETURN
518
519	END SUBROUTINE mym_length
520
521	!
522	! SUBROUTINE mym_turbulence:
523	!
524	! Input variables: see subroutine mym_initialize
525	! levflag : <>3; Level 2.5
526	! = 3; Level 3
527	!
528	! # ql, vt, vq, qke, tsq, qsq and cov are changed to input variables.
529	!
530	! Output variables: see subroutine mym_initialize
531	! dfm(mx,my,nz) : Diffusivity coefficient for momentum,
532	! divided by dz (not dz*h(i,j)) (m/s)
533	! dfh(mx,my,nz) : Diffusivity coefficient for heat,
534	! divided by dz (not dz*h(i,j)) (m/s)
535	! dfq(mx,my,nz) : Diffusivity coefficient for q^2,
536	! divided by dz (not dz*h(i,j)) (m/s)
537	! tcd(mx,my,nz) : Countergradient diffusion term for Theta_l
538	! (K/s)
539	! qcd(mx,my,nz) : Countergradient diffusion term for Q_w
540	! (kg/kg s)
541	! pd?(mx,my,nz) : Half of the production terms
542	!
543	! Only tcd and qcd are defined at the center of the grid boxes
544	!
545	! # DO NOT forget that tcd and qcd are added on the right-hand side
546	! of the equations for Theta_l and Q_w, respectively.
547	!
548	! Work arrays: see subroutine mym_initialize and level2
549	!
550	! # dtl, dqw, dtv, gm and gh are allowed to share storage units with
551	! dfm, dfh, dfq, tcd and qcd, respectively, for saving memory.
552	!
553	SUBROUTINE mym_turbulence ( kts,kte,&
554	& levflag, &
555	& dz, zw, &
556	& u, v, thl, ql, qw, &
557	& qke, tsq, qsq, cov, &
558	& vt, vq,&
559	& rmo, flt, flq, &
560	& El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc)
561
562	!
563	INTEGER, INTENT(IN) :: kts,kte
564	INTEGER, INTENT(IN) :: levflag
565	REAL, DIMENSION(kts:kte), INTENT(in) :: dz
566	REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
567	REAL, INTENT(in) :: rmo,flt,flq
568	REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,&
569	&ql,vt,vq,qke,tsq,qsq,cov
570
571	REAL, DIMENSION(kts:kte), INTENT(out) :: dfm,dfh,dfq,&
572	&pdk,pdt,pdq,pdc,tcd,qcd,el
573
574	REAL, DIMENSION(kts:kte) :: qkw,dtl,dqw,dtv,gm,gh,sm,sh
575
576	INTEGER :: k
577	! REAL :: cc2,cc3,e1c,e2c,e3c,e4c,e5c
578	REAL :: e6c,dzk,afk,abk,vtt,vqq,&
579	&cw25,clow,cupp,gamt,gamq,smd,gamv,elq,elh
580
581	DOUBLE PRECISION q2sq, t2sq, r2sq, c2sq, elsq, gmel, ghel
582	DOUBLE PRECISION q3sq, t3sq, r3sq, c3sq, dlsq, qdiv
583	DOUBLE PRECISION e1, e2, e3, e4, enum, eden, wden
584	!
585	! tv0 = 0.61*tref
586	! gtr = 9.81/tref
587	!
588	! cc2 = 1.0-c2
589	! cc3 = 1.0-c3
590	! e1c = 3.0a2b2*cc3
591	! e2c = 9.0a1a2*cc2
592	! e3c = 9.0a2a2cc2( 1.0-c5 )
593	! e4c = 12.0a1a2*cc2
594	! e5c = 6.0a1a1
595	!
596
597	CALL mym_level2 (kts,kte,&
598	& dz, &
599	& u, v, thl, qw, &
600	& ql, vt, vq, &
601	& dtl, dqw, dtv, gm, gh, sm, sh )
602	!
603	CALL mym_length (kts,kte, &
604	& dz, zw, &
605	& rmo, flt, flq, &
606	& vt, vq, &
607	& qke, &
608	& dtv, &
609	& el, &
610	& qkw)
611	!
612	DO k = kts+1,kte
613	dzk = 0.5 *( dz(k)+dz(k-1) )
614	afk = dz(k)/( dz(k)+dz(k-1) )
615	abk = 1.0 -afk
616	elsq = el (k)**2
617	q2sq = b1elsq( sm(k)gm(k)+sh(k)gh(k) )
618	q3sq = qkw(k)**2
619	!
620	! Modified: Dec/22/2005, from here, (dlsq -> elsq)
621	gmel = gm (k)*elsq
622	ghel = gh (k)*elsq
623	! Modified: Dec/22/2005, up to here
624	!
625	! Since qkw is set to more than 0.0, q3sq > 0.0.
626	IF ( q3sq .LT. q2sq ) THEN
627	qdiv = SQRT( q3sq/q2sq )
628	sm(k) = sm(k) * qdiv
629	sh(k) = sh(k) * qdiv
630	!
631	e1 = q3sq - e1cghel qdiv**2
632	e2 = q3sq - e2cghel qdiv**2
633	e3 = e1 + e3cghel qdiv**2
634	e4 = e1 - e4cghel qdiv**2
635	eden = e2e4 + e3e5cgmel qdiv**2
636	eden = MAX( eden, 1.0d-20 )
637	ELSE
638	e1 = q3sq - e1c*ghel
639	e2 = q3sq - e2c*ghel
640	e3 = e1 + e3c*ghel
641	e4 = e1 - e4c*ghel
642	eden = e2e4 + e3e5c*gmel
643	eden = MAX( eden, 1.0d-20 )
644	!
645	qdiv = 1.0
646	sm(k) = q3sqa1( e3-3.0c1e4 )/eden
647	sh(k) = q3sqa2( e2+3.0c1e5c*gmel )/eden
648	END IF
649	!
650	! Level 3 : start
651	IF ( levflag .EQ. 3 ) THEN
652	t2sq = qdivb2elsqsh(k)dtl(k)**2
653	r2sq = qdivb2elsqsh(k)dqw(k)**2
654	c2sq = qdivb2elsqsh(k)dtl(k)*dqw(k)
655	t3sq = MAX( tsq(k)abk+tsq(k-1)afk, 0.0 )
656	r3sq = MAX( qsq(k)abk+qsq(k-1)afk, 0.0 )
657	c3sq = cov(k)abk+cov(k-1)afk
658	!
659	! Modified: Dec/22/2005, from here
660	c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq )
661	!
662	vtt = 1.0 +vt(k)abk +vt(k-1)afk
663	vqq = tv0 +vq(k)abk +vq(k-1)afk
664	t2sq = vttt2sq +vqqc2sq
665	r2sq = vttc2sq +vqqr2sq
666	c2sq = MAX( vttt2sq+vqqr2sq, 0.0d0 )
667	t3sq = vttt3sq +vqqc3sq
668	r3sq = vttc3sq +vqqr3sq
669	c3sq = MAX( vttt3sq+vqqr3sq, 0.0d0 )
670	!
671	cw25 = e1( e2 + 3.0c1e5cgmelqdiv2 )/( 3.0eden )
672	!
673	! Limitation on q, instead of L/q
674	dlsq = elsq
675	IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k)
676	!
677	! Limitation on c3sq (0.12 =< cw =< 0.76)
678	e2 = q3sq - e2cghel qdiv**2
679	e3 = q3sq + e3cghel qdiv**2
680	e4 = q3sq - e4cghel qdiv**2
681	eden = e2e4 + e3 e5cgmel qdiv**2
682	!
683	wden = cc3gtr2 dlsq*2/elsq qdiv**2 &
684	& ( e2e4c - e3ce5cgmel * qdiv**2 )
685	!
686	IF ( wden .NE. 0.0 ) THEN
687	clow = q3sq( 0.12-cw25 )eden/wden
688	cupp = q3sq( 0.76-cw25 )eden/wden
689	!
690	IF ( wden .GT. 0.0 ) THEN
691	c3sq = MIN( MAX( c3sq, c2sq+clow ), c2sq+cupp )
692	ELSE
693	c3sq = MAX( MIN( c3sq, c2sq+clow ), c2sq+cupp )
694	END IF
695	END IF
696	!
697	e1 = e2 + e5cgmel qdiv**2
698	eden = MAX( eden, 1.0d-20 )
699	! Modified: Dec/22/2005, up to here
700	!
701	e6c = 3.0a2cc3gtr dlsq/elsq
702	!
703	! for Gamma_theta
704	!! enum = qdive6c( t3sq-t2sq )
705	IF ( t2sq .GE. 0.0 ) THEN
706	enum = MAX( qdive6c( t3sq-t2sq ), 0.0d0 )
707	ELSE
708	enum = MIN( qdive6c( t3sq-t2sq ), 0.0d0 )
709	ENDIF
710
711	gamt =-e1 *enum /eden
712	!
713	! for Gamma_q
714	!! enum = qdive6c( r3sq-r2sq )
715	IF ( r2sq .GE. 0.0 ) THEN
716	enum = MAX( qdive6c( r3sq-r2sq ), 0.0d0 )
717	ELSE
718	enum = MIN( qdive6c( r3sq-r2sq ), 0.0d0 )
719	ENDIF
720
721	gamq =-e1 *enum /eden
722	!
723	! for Sm' and Sh'd(Theta_V)/dz
724	!! enum = qdive6c( c3sq-c2sq )
725	enum = MAX( qdive6c( c3sq-c2sq ), 0.0d0)
726
727	smd = dlsqenumgtr/eden * qdiv*2 (e3c+e4c)*a1/a2
728	gamv = e1 enumgtr/eden
729	!
730
731	sm(k) = sm(k) +smd
732	!
733	! For elh (see below), qdiv at Level 3 is reset to 1.0.
734	qdiv = 1.0
735	! Level 3 : end
736	!
737	ELSE
738	! At Level 2.5, qdiv is not reset.
739	gamt = 0.0
740	gamq = 0.0
741	gamv = 0.0
742	END IF
743	!
744	elq = el(k)*qkw(k)
745	elh = elq*qdiv
746	!
747	pdk(k) = elq( sm(k)gm (k) &
748	& +sh(k)*gh (k)+gamv )
749	pdt(k) = elh( sh(k)dtl(k)+gamt )*dtl(k)
750	pdq(k) = elh( sh(k)dqw(k)+gamq )*dqw(k)
751	pdc(k) = elh( sh(k)dtl(k)+gamt )&
752	&dqw(k)0.5 &
753	&+elh( sh(k)dqw(k)+gamq )dtl(k)0.5
754	!
755	tcd(k) = elq*gamt
756	qcd(k) = elq*gamq
757	!
758	dfm(k) = elq*sm (k) / dzk
759	dfh(k) = elq*sh (k) / dzk
760	! Modified: Dec/22/2005, from here
761	! In sub.mym_predict, dfq for the TKE and scalar variance
762	! ** are set to 3.0dfm and 1.0dfm, respectively. **
763	dfq(k) = dfm(k)
764	! Modified: Dec/22/2005, up to here
765	END DO
766	!
767	dfm(kts) = 0.0
768	dfh(kts) = 0.0
769	dfq(kts) = 0.0
770	tcd(kts) = 0.0
771	qcd(kts) = 0.0
772
773	tcd(kte) = 0.0
774	qcd(kte) = 0.0
775
776	!
777	DO k = kts,kte-1
778	dzk = dz(k)
779	tcd(k) = ( tcd(k+1)-tcd(k) )/( dzk )
780	qcd(k) = ( qcd(k+1)-qcd(k) )/( dzk )
781	END DO
782	!
783	RETURN
784
785	END SUBROUTINE mym_turbulence
786
787	!
788	! SUBROUTINE mym_predict:
789	!
790	! Input variables: see subroutine mym_initialize and turbulence
791	! qke(mx,my,nz) : qke at (n)th time level
792	! tsq, ...cov : ditto
793	!
794	! Output variables:
795	! qke(mx,my,nz) : qke at (n+1)th time level
796	! tsq, ...cov : ditto
797	!
798	! Work arrays:
799	! qkw(mx,my,nz) : q at the center of the grid boxes (m/s)
800	! bp (mx,my,nz) : = 1/2*F, see below
801	! rp (mx,my,nz) : = P-1/2FQ, see below
802	!
803	! # The equation for a turbulent quantity Q can be expressed as
804	! dQ/dt + Ah + Av = Dh + Dv + P - F*Q, (1)
805	! where A is the advection, D the diffusion, P the production,
806	! F*Q the dissipation and h and v denote horizontal and vertical,
807	! respectively. If Q is q^2, F is 2q/B_1L.
808	! Using the Crank-Nicholson scheme for Av, Dv and F*Q, a finite
809	! difference equation is written as
810	! Q{n+1} - Q{n} = dt *( Dh{n} - Ah{n} + P{n} )
811	! + dt/2( Dv{n} - Av{n} - FQ{n} )
812	! + dt/2( Dv{n+1} - Av{n+1} - FQ{n+1} ), (2)
813	! where n denotes the time level.
814	! When the advection and diffusion terms are discretized as
815	! dt/2*( Dv - Av ) = a(k)Q(k+1) - b(k)Q(k) + c(k)Q(k-1), (3)
816	! Eq.(2) can be rewritten as
817	! - a(k)Q(k+1) + [ 1 + b(k) + dt/2*F ]Q(k) - c(k)Q(k-1)
818	! = Q{n} + dt *( Dh{n} - Ah{n} + P{n} )
819	! + dt/2( Dv{n} - Av{n} - FQ{n} ), (4)
820	! where Q on the left-hand side is at (n+1)th time level.
821	!
822	! In this subroutine, a(k), b(k) and c(k) are obtained from
823	! subprogram coefvu and are passed to subprogram tinteg via
824	! common. 1/2F and P-1/2F*Q are stored in bp and rp,
825	! respectively. Subprogram tinteg solves Eq.(4).
826	!
827	! Modify this subroutine according to your numerical integration
828	! scheme (program).
829	!
830
831
832	SUBROUTINE mym_predict (kts,kte,&
833	& levflag, &
834	& delt,&
835	& dz, &
836	& ust, flt, flq, pmz, phh, &
837	& el, dfq, &
838	& pdk, pdt, pdq, pdc,&
839	& qke, tsq, qsq, cov)
840
841	INTEGER, INTENT(IN) :: kts,kte
842	INTEGER, INTENT(IN) :: levflag
843	REAL, INTENT(IN) :: delt
844	REAL, DIMENSION(kts:kte), INTENT(IN) :: dz, dfq,el
845	REAL, DIMENSION(kts:kte), INTENT(INOUT) :: pdk, pdt, pdq, pdc
846	REAL, INTENT(IN) :: flt, flq, ust, pmz, phh
847	REAL, DIMENSION(kts:kte), INTENT(INOUT) :: qke,tsq, qsq, cov
848
849	INTEGER :: k,nz
850	REAL, DIMENSION(kts:kte) :: qkw, bp, rp, df3q
851	REAL :: vkz,pdk1,phm,pdt1,pdq1,pdc1,b1l,b2l
852	REAL, DIMENSION(kts:kte) :: dtz
853	REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d
854
855	nz=kte-kts+1
856
857	! ** Strictly, vkzh(i,j) -> vk( 0.5dz(1)h(i,j)+z0 ) **
858	vkz = vk0.5dz(kts)
859	!
860	! Modified: Dec/22/2005, from here
861	! ** dfq for the TKE is 3.0dfm. *
862	! CALL coefvu ( dfq, 3.0 ) ! make change here
863	! Modified: Dec/22/2005, up to here
864	!
865	DO k = kts,kte
866	!! qke(k) = MAX(qke(k), 0.0)
867	qkw(k) = SQRT( MAX( qke(k), 0.0 ) )
868	df3q(k)=3.*dfq(k)
869	dtz(k)=delt/dz(k)
870	END DO
871	!
872	pdk1 = 2.0ust3pmz/( vkz )
873	phm = 2.0/ust *phh/( vkz )
874	pdt1 = phmflt*2
875	pdq1 = phmflq*2
876	pdc1 = phmfltflq
877	!
878	! pdk(i,j,1)+pdk(i,j,2) corresponds to pdk1.
879	pdk(kts) = pdk1 -pdk(kts+1)
880
881	!! pdt(kts) = pdt1 -pdt(kts+1)
882	!! pdq(kts) = pdq1 -pdq(kts+1)
883	!! pdc(kts) = pdc1 -pdc(kts+1)
884	pdt(kts) = pdt(kts+1)
885	pdq(kts) = pdq(kts+1)
886	pdc(kts) = pdc(kts+1)
887	!
888	! Prediction of twice the turbulent kinetic energy
889	!! DO k = kts+1,kte-1
890	DO k = kts,kte-1
891	b1l = b10.5( el(k+1)+el(k) )
892	bp(k) = 2.*qkw(k) / b1l
893	rp(k) = pdk(k+1) + pdk(k)
894	END DO
895
896	!! a(1)=0.
897	!! b(1)=1.
898	!! c(1)=-1.
899	!! d(1)=0.
900
901	! Since df3q(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)df3q(k+1)+bp(k)delt.
902	DO k=kts,kte-1
903	a(k-kts+1)=-dtz(k)*df3q(k)
904	b(k-kts+1)=1.+dtz(k)(df3q(k)+df3q(k+1))+bp(k)delt
905	c(k-kts+1)=-dtz(k)*df3q(k+1)
906	d(k-kts+1)=rp(k)*delt + qke(k)
907	ENDDO
908
909	!! DO k=kts+1,kte-1
910	!! a(k-kts+1)=-dtz(k)*df3q(k)
911	!! b(k-kts+1)=1.+dtz(k)*(df3q(k)+df3q(k+1))
912	!! c(k-kts+1)=-dtz(k)*df3q(k+1)
913	!! d(k-kts+1)=rp(k)delt + qke(k) - qke(k)bp(k)*delt
914	!! ENDDO
915
916	a(nz)=-1. !0.
917	b(nz)=1.
918	c(nz)=0.
919	d(nz)=0.
920
921	CALL tridiag(nz,a,b,c,d)
922
923	DO k=kts,kte
924	qke(k)=d(k-kts+1)
925	ENDDO
926
927
928	IF ( levflag .EQ. 3 ) THEN
929	!
930	! Modified: Dec/22/2005, from here
931	! ** dfq for the scalar variance is 1.0dfm. *
932	! CALL coefvu ( dfq, 1.0 ) make change here
933	! Modified: Dec/22/2005, up to here
934	!
935	! Prediction of the temperature variance
936	!! DO k = kts+1,kte-1
937	DO k = kts,kte-1
938	b2l = b20.5( el(k+1)+el(k) )
939	bp(k) = 2.*qkw(k) / b2l
940	rp(k) = pdt(k+1) + pdt(k)
941	END DO
942
943	!zero gradient for tsq at bottom and top
944
945	!! a(1)=0.
946	!! b(1)=1.
947	!! c(1)=-1.
948	!! d(1)=0.
949
950	! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)dfq(k+1)+bp(k)delt.
951	DO k=kts,kte-1
952	a(k-kts+1)=-dtz(k)*dfq(k)
953	b(k-kts+1)=1.+dtz(k)(dfq(k)+dfq(k+1))+bp(k)delt
954	c(k-kts+1)=-dtz(k)*dfq(k+1)
955	d(k-kts+1)=rp(k)*delt + tsq(k)
956	ENDDO
957
958	!! DO k=kts+1,kte-1
959	!! a(k-kts+1)=-dtz(k)*dfq(k)
960	!! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
961	!! c(k-kts+1)=-dtz(k)*dfq(k+1)
962	!! d(k-kts+1)=rp(k)delt + tsq(k) - tsq(k)bp(k)*delt
963	!! ENDDO
964
965	a(nz)=-1. !0.
966	b(nz)=1.
967	c(nz)=0.
968	d(nz)=0.
969
970	CALL tridiag(nz,a,b,c,d)
971
972	DO k=kts,kte
973	tsq(k)=d(k-kts+1)
974	ENDDO
975
976	! Prediction of the moisture variance
977	!! DO k = kts+1,kte-1
978	DO k = kts,kte-1
979	b2l = b20.5( el(k+1)+el(k) )
980	bp(k) = 2.*qkw(k) / b2l
981	rp(k) = pdq(k+1) +pdq(k)
982	END DO
983
984	!zero gradient for qsq at bottom and top
985
986	!! a(1)=0.
987	!! b(1)=1.
988	!! c(1)=-1.
989	!! d(1)=0.
990
991	! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)dfq(k+1)+bp(k)delt.
992	DO k=kts,kte-1
993	a(k-kts+1)=-dtz(k)*dfq(k)
994	b(k-kts+1)=1.+dtz(k)(dfq(k)+dfq(k+1))+bp(k)delt
995	c(k-kts+1)=-dtz(k)*dfq(k+1)
996	d(k-kts+1)=rp(k)*delt + qsq(k)
997	ENDDO
998
999	!! DO k=kts+1,kte-1
1000	!! a(k-kts+1)=-dtz(k)*dfq(k)
1001	!! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
1002	!! c(k-kts+1)=-dtz(k)*dfq(k+1)
1003	!! d(k-kts+1)=rp(k)delt + qsq(k) -qsq(k)bp(k)*delt
1004	!! ENDDO
1005
1006	a(nz)=-1. !0.
1007	b(nz)=1.
1008	c(nz)=0.
1009	d(nz)=0.
1010
1011	CALL tridiag(nz,a,b,c,d)
1012
1013	DO k=kts,kte
1014	qsq(k)=d(k-kts+1)
1015	ENDDO
1016
1017	! Prediction of the temperature-moisture covariance
1018	!! DO k = kts+1,kte-1
1019	DO k = kts,kte-1
1020	b2l = b20.5( el(k+1)+el(k) )
1021	bp(k) = 2.*qkw(k) / b2l
1022	rp(k) = pdc(k+1) + pdc(k)
1023	END DO
1024
1025	!zero gradient for tqcov at bottom and top
1026
1027	!! a(1)=0.
1028	!! b(1)=1.
1029	!! c(1)=-1.
1030	!! d(1)=0.
1031
1032	! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)dfq(k+1)+bp(k)delt.
1033	DO k=kts,kte-1
1034	a(k-kts+1)=-dtz(k)*dfq(k)
1035	b(k-kts+1)=1.+dtz(k)(dfq(k)+dfq(k+1))+bp(k)delt
1036	c(k-kts+1)=-dtz(k)*dfq(k+1)
1037	d(k-kts+1)=rp(k)*delt + cov(k)
1038	ENDDO
1039
1040	!! DO k=kts+1,kte-1
1041	!! a(k-kts+1)=-dtz(k)*dfq(k)
1042	!! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
1043	!! c(k-kts+1)=-dtz(k)*dfq(k+1)
1044	!! d(k-kts+1)=rp(k)delt + cov(k) - cov(k)bp(k)*delt
1045	!! ENDDO
1046
1047	a(nz)=-1. !0.
1048	b(nz)=1.
1049	c(nz)=0.
1050	d(nz)=0.
1051
1052	CALL tridiag(nz,a,b,c,d)
1053
1054	DO k=kts,kte
1055	cov(k)=d(k-kts+1)
1056	ENDDO
1057
1058	ELSE
1059	!! DO k = kts+1,kte-1
1060	DO k = kts,kte-1
1061	IF ( qkw(k) .LE. 0.0 ) THEN
1062	b2l = 0.0
1063	ELSE
1064	b2l = b20.25( el(k+1)+el(k) )/qkw(k)
1065	END IF
1066	!
1067	tsq(k) = b2l*( pdt(k+1)+pdt(k) )
1068	qsq(k) = b2l*( pdq(k+1)+pdq(k) )
1069	cov(k) = b2l*( pdc(k+1)+pdc(k) )
1070	END DO
1071
1072	!! tsq(kts)=tsq(kts+1)
1073	!! qsq(kts)=qsq(kts+1)
1074	!! cov(kts)=cov(kts+1)
1075
1076	tsq(kte)=tsq(kte-1)
1077	qsq(kte)=qsq(kte-1)
1078	cov(kte)=cov(kte-1)
1079
1080	END IF
1081
1082	END SUBROUTINE mym_predict
1083
1084	! SUBROUTINE mym_condensation:
1085	!
1086	! Input variables: see subroutine mym_initialize and turbulence
1087	! pi (mx,my,nz) : Perturbation of the Exner function (J/kg K)
1088	! defined on the walls of the grid boxes
1089	! This is usually computed by integrating
1090	! d(pi)/dz = hgtv/tref**2
1091	! from the upper boundary, where tv is the
1092	! virtual potential temperature minus tref.
1093	!
1094	! Output variables: see subroutine mym_initialize
1095	! cld(mx,my,nz) : Cloud fraction
1096	!
1097	! Work arrays:
1098	! qmq(mx,my,nz) : Q_w-Q_{sl}, where Q_{sl} is the saturation
1099	! specific humidity at T=Tl
1100	! alp(mx,my,nz) : Functions in the condensation process
1101	! bet(mx,my,nz) : ditto
1102	! sgm(mx,my,nz) : Combined standard deviation sigma_s
1103	! multiplied by 2/alp
1104	!
1105	! # qmq, alp, bet and sgm are allowed to share storage units with
1106	! any four of other work arrays for saving memory.
1107	!
1108	! # Results are sensitive particularly to values of cp and rd.
1109	! Set these values to those adopted by you.
1110	!
1111
1112
1113
1114	SUBROUTINE mym_condensation (kts,kte, &
1115	& dz, &
1116	& thl, qw, &
1117	& p,exner, &
1118	& tsq, qsq, cov, &
1119	& Vt, Vq)
1120
1121	INTEGER, INTENT(IN) :: kts,kte
1122
1123	REAL, DIMENSION(kts:kte), INTENT(IN) :: dz
1124	REAL, DIMENSION(kts:kte), INTENT(IN) :: p,exner, thl, qw, &
1125	&tsq, qsq, cov
1126
1127	REAL, DIMENSION(kts:kte), INTENT(OUT) :: vt,vq
1128
1129
1130	REAL, DIMENSION(kts:kte) :: qmq,alp,bet,sgm,ql,cld
1131
1132	DOUBLE PRECISION :: t3sq, r3sq, c3sq
1133	!
1134
1135	REAL :: p2a,t,esl,qsl,dqsl,q1,cld0,eq1,qll,&
1136	&q2p,pt,rac,qt
1137	INTEGER :: i,j,k
1138
1139	REAL :: erf
1140
1141	! Note: kte needs to be larger than kts, i.e., kte >= kts+1.
1142
1143	DO k = kts,kte-1
1144	p2a = exner(k)
1145	t = thl(k)*p2a
1146
1147	!x if ( ct .gt. 0.0 ) then
1148	! a = 17.27
1149	! b = 237.3
1150	!x else
1151	!x a = 21.87
1152	!x b = 265.5
1153	!x end if
1154	!
1155	! ** 3.8 = 0.6226.11 (hPa) *
1156	esl=svp11EXP(svp2(t-svpt0)/(t-svp3))
1157	qsl=ep_2esl/(p(k)-ep_3esl)
1158	! qsl = 3.8EXP( act/( b+ct ) ) / ( 1000.0p2a*rk )
1159	dqsl = qslep_2ev/( rdt*2 )
1160	!
1161	qmq(k) = qw(k) -qsl
1162
1163	alp(k) = 1.0/( 1.0+dqsl*xlvcp )
1164	bet(k) = dqsl*p2a
1165	!
1166	t3sq = MAX( tsq(k), 0.0 )
1167	r3sq = MAX( qsq(k), 0.0 )
1168	c3sq = cov(k)
1169	c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq )
1170	!
1171	r3sq = r3sq +bet(k)*2t3sq -2.0bet(k)c3sq
1172	sgm(k) = SQRT( MAX( r3sq, 1.0d-10 ) )
1173	END DO
1174	!
1175	DO k = kts,kte-1
1176	q1 = qmq(k) / sgm(k)
1177	cld0 = 0.5( 1.0+erf( q1rr2 ) )
1178	! q1=0.
1179	! cld0=0.
1180
1181	eq1 = rrpEXP( -0.5q1*q1 )
1182	qll = MAX( cld0*q1 + eq1, 0.0 )
1183
1184	cld(k) = cld0
1185	ql (k) = alp(k)sgm(k)qll
1186	!
1187	q2p = xlvcp/exner( k )
1188	pt = thl(k) +q2p*ql(k)
1189	qt = 1.0 +p608qw(k) -(1.+p608)ql(k)
1190	rac = alp(k)( cld0-qlleq1 )( q2pqt-(1.+p608)*pt )
1191	!
1192	vt (k) = qt-1.0 -rac*bet(k)
1193	vq (k) = p608*pt-tv0 +rac
1194	END DO
1195	!
1196
1197	cld(kte) = cld(kte-1)
1198	ql(kte) = ql(kte-1)
1199	vt(kte) = vt(kte-1)
1200	vq(kte) = vq(kte-1)
1201
1202	RETURN
1203
1204	END SUBROUTINE mym_condensation
1205
1206	SUBROUTINE mynn_tendencies(kts,kte,&
1207	&levflag,grav_settling,&
1208	&delt,&
1209	&dz,&
1210	&u,v,th,qv,qc,p,exner,&
1211	&thl,sqv,sqc,sqw,&
1212	&ust,flt,flq,wspd,qcg,&
1213	&tsq,qsq,cov,&
1214	&tcd,qcd,&
1215	&dfm,dfh,dfq,&
1216	&Du,Dv,Dth,Dqv,Dqc)
1217
1218	INTEGER, INTENT(in) :: kts,kte
1219	INTEGER, INTENT(in) :: grav_settling,levflag
1220
1221	!! grav_settling = 1 for gravitational settling of droplets
1222	!! grav_settling = 0 otherwise
1223	! thl - liquid water potential temperature
1224	! qw - total water
1225	! dfm,dfh,dfq - as above
1226	! flt - surface flux of thl
1227	! flq - surface flux of qw
1228
1229
1230	REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,th,qv,qc,p,exner,&
1231	&dfm,dfh,dfq,dz,tsq,qsq,cov,tcd,qcd
1232	REAL, DIMENSION(kts:kte), INTENT(inout) :: thl,sqw,sqv,sqc
1233	REAL, DIMENSION(kts:kte), INTENT(out) :: du,dv,dth,dqv,dqc
1234	REAL, INTENT(IN) :: delt,ust,flt,flq,wspd,qcg
1235
1236	! REAL, INTENT(IN) :: delt,ust,flt,flq,qcg,&
1237	! &gradu_top,gradv_top,gradth_top,gradqv_top
1238
1239
1240	!local vars
1241
1242	REAL, DIMENSION(kts:kte) :: dtz,vt,vq
1243
1244	REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d
1245
1246	REAL :: rhs,gfluxm,gfluxp,dztop
1247	INTEGER :: k,kk,nz
1248
1249	nz=kte-kts+1
1250
1251	dztop=.5*(dz(kte)+dz(kte-1))
1252
1253	DO k=kts,kte
1254	dtz(k)=delt/dz(k)
1255	ENDDO
1256
1257	!! u
1258
1259	k=kts
1260
1261	a(1)=0.
1262	b(1)=1.+dtz(k)(dfm(k+1)+ust*2/wspd)
1263	c(1)=-dtz(k)*dfm(k+1)
1264	d(1)=u(k)
1265
1266	!! a(1)=0.
1267	!! b(1)=1.+dtz(k)*dfm(k+1)
1268	!! c(1)=-dtz(k)*dfm(k+1)
1269	!! d(1)=u(k)(1.-ust2/wspddtz(k))
1270
1271	DO k=kts+1,kte-1
1272	kk=k-kts+1
1273	a(kk)=-dtz(k)*dfm(k)
1274	b(kk)=1.+dtz(k)*(dfm(k)+dfm(k+1))
1275	c(kk)=-dtz(k)*dfm(k+1)
1276	d(kk)=u(k)
1277	ENDDO
1278
1279	!! no flux at the top
1280
1281	! a(nz)=-1.
1282	! b(nz)=1.
1283	! c(nz)=0.
1284	! d(nz)=0.
1285
1286	!! specified gradient at the top
1287
1288	! a(nz)=-1.
1289	! b(nz)=1.
1290	! c(nz)=0.
1291	! d(nz)=gradu_top*dztop
1292
1293	!! prescribed value
1294
1295	a(nz)=0
1296	b(nz)=1.
1297	c(nz)=0.
1298	d(nz)=u(kte)
1299
1300	CALL tridiag(nz,a,b,c,d)
1301
1302	DO k=kts,kte
1303	du(k)=(d(k-kts+1)-u(k))/delt
1304	ENDDO
1305
1306	!! v
1307
1308	k=kts
1309
1310	a(1)=0.
1311	b(1)=1.+dtz(k)(dfm(k+1)+ust*2/wspd)
1312	c(1)=-dtz(k)*dfm(k+1)
1313	d(1)=v(k)
1314
1315	!! a(1)=0.
1316	!! b(1)=1.+dtz(k)*dfm(k+1)
1317	!! c(1)=-dtz(k)*dfm(k+1)
1318	!! d(1)=v(k)(1.-ust2/wspddtz(k))
1319
1320	DO k=kts+1,kte-1
1321	kk=k-kts+1
1322	a(kk)=-dtz(k)*dfm(k)
1323	b(kk)=1.+dtz(k)*(dfm(k)+dfm(k+1))
1324	c(kk)=-dtz(k)*dfm(k+1)
1325	d(kk)=v(k)
1326	ENDDO
1327
1328	!! no flux at the top
1329
1330	! a(nz)=-1.
1331	! b(nz)=1.
1332	! c(nz)=0.
1333	! d(nz)=0.
1334
1335
1336	!! specified gradient at the top
1337
1338	! a(nz)=-1.
1339	! b(nz)=1.
1340	! c(nz)=0.
1341	! d(nz)=gradv_top*dztop
1342
1343	!! prescribed value
1344
1345	a(nz)=0
1346	b(nz)=1.
1347	c(nz)=0.
1348	d(nz)=v(kte)
1349
1350	CALL tridiag(nz,a,b,c,d)
1351
1352	DO k=kts,kte
1353	dv(k)=(d(k-kts+1)-v(k))/delt
1354	ENDDO
1355
1356	!! thl
1357
1358	k=kts
1359
1360	a(1)=0.
1361	b(1)=1.+dtz(k)*dfh(k+1)
1362	c(1)=-dtz(k)*dfh(k+1)
1363
1364	! if qcg not used then assume constant flux in the surface layer
1365
1366	IF (qcg < qcgmin) THEN
1367	IF (sqc(k) > qcgmin) THEN
1368	gfluxm=grav_settlinggnosqc(k)**gpw
1369	ELSE
1370	gfluxm=0.
1371	ENDIF
1372	ELSE
1373	gfluxm=grav_settlinggno(qcg/(1.+qcg))**gpw
1374	ENDIF
1375
1376	IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
1377	gfluxp=grav_settlinggno(.5(sqc(k+1)+sqc(k)))*gpw
1378	ELSE
1379	gfluxp=0.
1380	ENDIF
1381
1382	rhs=-xlvcp/exner(k)&
1383	&*( &
1384	&(gfluxp - gfluxm)/dz(k)&
1385	& ) + tcd(k)
1386
1387	d(1)=thl(k)+dtz(k)flt+rhsdelt
1388
1389	DO k=kts+1,kte-1
1390	kk=k-kts+1
1391	a(kk)=-dtz(k)*dfh(k)
1392	b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1))
1393	c(kk)=-dtz(k)*dfh(k+1)
1394
1395	IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
1396	gfluxp=grav_settlinggno(.5(sqc(k+1)+sqc(k)))*gpw
1397	ELSE
1398	gfluxp=0.
1399	ENDIF
1400
1401	IF (.5*(sqc(k-1)+sqc(k)) > qcgmin) THEN
1402	gfluxm=grav_settlinggno(.5(sqc(k-1)+sqc(k)))*gpw
1403	ELSE
1404	gfluxm=0.
1405	ENDIF
1406
1407	rhs=-xlvcp/exner(k)&
1408	&*( &
1409	&(gfluxp - gfluxm)/dz(k)&
1410	& ) + tcd(k)
1411	d(kk)=thl(k)+rhs*delt
1412	ENDDO
1413
1414	!! no flux at the top
1415
1416	! a(nz)=-1.
1417	! b(nz)=1.
1418	! c(nz)=0.
1419	! d(nz)=0.
1420
1421	!! specified gradient at the top
1422
1423	!assume gradthl_top=gradth_top
1424
1425	! a(nz)=-1.
1426	! b(nz)=1.
1427	! c(nz)=0.
1428	! d(nz)=gradth_top*dztop
1429
1430	!! prescribed value
1431
1432	a(nz)=0.
1433	b(nz)=1.
1434	c(nz)=0.
1435	d(nz)=thl(kte)
1436
1437	CALL tridiag(nz,a,b,c,d)
1438
1439	DO k=kts,kte
1440	thl(k)=d(k-kts+1)
1441	ENDDO
1442
1443	!! total water
1444
1445	k=kts
1446
1447	a(1)=0.
1448	b(1)=1.+dtz(k)*dfh(k+1)
1449	c(1)=-dtz(k)*dfh(k+1)
1450
1451	IF (qcg < qcgmin) THEN
1452	IF (sqc(k) > qcgmin) THEN
1453	gfluxm=grav_settlinggnosqc(k)**gpw
1454	ELSE
1455	gfluxm=0.
1456	ENDIF
1457	ELSE
1458	gfluxm=grav_settlinggno(qcg/(1.+qcg))**gpw
1459	ENDIF
1460
1461	IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
1462	gfluxp=grav_settlinggno(.5(sqc(k+1)+sqc(k)))*gpw
1463	ELSE
1464	gfluxp=0.
1465	ENDIF
1466
1467	rhs=&
1468	&( &
1469	&(gfluxp - gfluxm)/dz(k)&
1470	& ) + qcd(k)
1471
1472	d(1)=sqw(k)+dtz(k)flq+rhsdelt
1473
1474	DO k=kts+1,kte-1
1475	kk=k-kts+1
1476	a(kk)=-dtz(k)*dfh(k)
1477	b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1))
1478	c(kk)=-dtz(k)*dfh(k+1)
1479
1480	IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
1481	gfluxp=grav_settlinggno(.5(sqc(k+1)+sqc(k)))*gpw
1482	ELSE
1483	gfluxp=0.
1484	ENDIF
1485
1486	IF (.5*(sqc(k-1)+sqc(k)) > qcgmin) THEN
1487	gfluxm=grav_settlinggno(.5(sqc(k-1)+sqc(k)))*gpw
1488	ELSE
1489	gfluxm=0.
1490	ENDIF
1491
1492	rhs=&
1493	&( &
1494	&(gfluxp - gfluxm)/dz(k)&
1495	& ) + qcd(k)
1496	d(kk)=sqw(k) + rhs*delt
1497	ENDDO
1498
1499
1500	!! no flux at the top
1501
1502	! a(nz)=-1.
1503	! b(nz)=1.
1504	! c(nz)=0.
1505	! d(nz)=0.
1506
1507	!! specified gradient at the top
1508	!assume gradqw_top=gradqv_top
1509
1510	! a(nz)=-1.
1511	! b(nz)=1.
1512	! c(nz)=0.
1513	! d(nz)=gradqv_top*dztop
1514
1515	!! prescribed value
1516
1517	a(nz)=0.
1518	b(nz)=1.
1519	c(nz)=0.
1520	d(nz)=sqw(kte)
1521
1522	CALL tridiag(nz,a,b,c,d)
1523
1524	!convert to mixing ratios for wrf
1525
1526	DO k=kts,kte
1527	sqw(k)=d(k-kts+1)
1528	sqv(k)=sqw(k)-sqc(k)
1529	Dqv(k)=(sqv(k)/(1.-sqv(k))-qv(k))/delt
1530	! Dqc(k)=(sqc(k)/(1.-sqc(k))-qc(k))/delt
1531	Dqc(k)=0.
1532	Dth(k)=(thl(k)+xlvcp/exner(k)*sqc(k)-th(k))/delt
1533	ENDDO
1534
1535	END SUBROUTINE mynn_tendencies
1536
1537	SUBROUTINE retrieve_exchange_coeffs(kts,kte,&
1538	&dfm,dfh,dfq,dz,&
1539	&K_m,K_h,K_q)
1540
1541	INTEGER , INTENT(in) :: kts,kte
1542
1543	REAL, DIMENSION(KtS:KtE), INTENT(in) :: dz,dfm,dfh,dfq
1544
1545	REAL, DIMENSION(KtS:KtE), INTENT(out) :: &
1546	&K_m, K_h, K_q
1547
1548
1549	INTEGER :: k
1550	REAL :: dzk
1551
1552	K_m(kts)=0.
1553	K_h(kts)=0.
1554	K_q(kts)=0.
1555
1556	DO k=kts+1,kte
1557	dzk = 0.5 *( dz(k)+dz(k-1) )
1558	K_m(k)=dfm(k)*dzk
1559	K_h(k)=dfh(k)*dzk
1560	K_q(k)=dfq(k)*dzk
1561	ENDDO
1562
1563	END SUBROUTINE retrieve_exchange_coeffs
1564
1565
1566	SUBROUTINE tridiag(n,a,b,c,d)
1567
1568	!! to solve system of linear eqs on tridiagonal matrix n times n
1569	!! after Peaceman and Rachford, 1955
1570	!! a,b,c,d - are vectors of order n
1571	!! a,b,c - are coefficients on the LHS
1572	!! d - is initially RHS on the output becomes a solution vector
1573
1574	INTEGER, INTENT(in):: n
1575	REAL, DIMENSION(n), INTENT(in) :: a,b
1576	REAL, DIMENSION(n), INTENT(inout) :: c,d
1577
1578	INTEGER :: i
1579	REAL :: p
1580	REAL, DIMENSION(n) :: q
1581
1582	c(n)=0.
1583	q(1)=-c(1)/b(1)
1584	d(1)=d(1)/b(1)
1585
1586	DO i=2,n
1587	p=1./(b(i)+a(i)*q(i-1))
1588	q(i)=-c(i)*p
1589	d(i)=(d(i)-a(i)d(i-1))p
1590	ENDDO
1591
1592	DO i=n-1,1,-1
1593	d(i)=d(i)+q(i)*d(i+1)
1594	ENDDO
1595
1596	END SUBROUTINE tridiag
1597
1598	SUBROUTINE mynn_bl_driver(&
1599	&initflag,&
1600	&grav_settling,&
1601	&delt,&
1602	&dz,&
1603	&u,v,th,qv,qc,&
1604	&p,exner,rho,&
1605	&xland,ts,qsfc,qcg,ps,&
1606	&ust,ch,hfx,qfx,rmol,wspd,&
1607	&Qke,Tsq,Qsq,Cov,&
1608	&Du,Dv,Dth,&
1609	&Dqv,Dqc,&
1610	! &K_m,K_h,K_q&
1611	&K_h,k_m,&
1612	&Pblh&
1613	&,IDS,IDE,JDS,JDE,KDS,KDE &
1614	&,IMS,IME,JMS,JME,KMS,KME &
1615	&,ITS,ITE,JTS,JTE,KTS,KTE)
1616
1617	INTEGER, INTENT(in) :: initflag
1618	INTEGER, INTENT(in) :: grav_settling
1619
1620	INTEGER,INTENT(IN) :: &
1621	& IDS,IDE,JDS,JDE,KDS,KDE &
1622	&,IMS,IME,JMS,JME,KMS,KME &
1623	&,ITS,ITE,JTS,JTE,KTS,KTE
1624
1625
1626	! initflag > 0 for TRUE
1627	! else for FALSE
1628	! levflag : <>3; Level 2.5
1629	! = 3; Level 3
1630	! grav_settling = 1 when gravitational settling accounted for
1631	! grav_settling = 0 when gravitational settling NOT accounted for
1632
1633	REAL, INTENT(in) :: delt
1634	REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(in) :: dz,&
1635	&u,v,th,qv,qc,p,exner,rho
1636	REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(in) :: xland,ust,&
1637	&ch,rmol,ts,qsfc,qcg,ps,hfx,qfx, wspd
1638	!
1639	REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: &
1640	&Qke,Tsq,Qsq,Cov
1641	REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: &
1642	&Du,Dv,Dth,Dqv,Dqc
1643
1644	REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: &
1645	&K_h,K_m
1646
1647	REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(inout) :: &
1648	&Pblh
1649
1650	!local vars
1651	INTEGER :: ITF,JTF,KTF
1652	INTEGER :: i,j,k
1653	REAL, DIMENSION(KMS:KME) :: thl,sqv,sqc,sqw,&
1654	&El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, Vt, Vq
1655
1656	REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME) :: K_q
1657
1658	REAL, DIMENSION(KMS:KME+1) :: zw
1659
1660	REAL :: cpm,sqcg,flt,flq,pmz,phh,exnerg,zet
1661
1662	REAL, DIMENSION(KMS:KME) :: thetav
1663	REAL :: thsfc
1664
1665
1666	INTEGER, SAVE :: levflag
1667
1668	JTF=MIN0(JTE,JDE-1)
1669	ITF=MIN0(ITE,IDE-1)
1670	KTF=MIN0(KTE,KDE-1)
1671
1672	levflag=mynn_level
1673
1674	IF (initflag > 0) THEN
1675
1676	DO j=JTS,JTF
1677	DO i=ITS,ITF
1678	DO k=KTS,KTF
1679	sqv(k)=qv(i,k,j)/(1.+qv(i,k,j))
1680	thl(k)=th(i,k,j)
1681	!JOE-PBLH test
1682	thetav(k)=th(i,k,j)(1.+0.61qv(i,k,j))
1683	!JOE-end
1684	IF (k==kts) THEN
1685	zw(k)=0.
1686	ELSE
1687	zw(k)=zw(k-1)+dz(i,k-1,j)
1688	ENDIF
1689
1690	k_m(i,k,j)=0.
1691	k_h(i,k,j)=0.
1692	k_q(i,k,j)=0.
1693
1694	ENDDO
1695
1696	zw(ktf+1)=zw(ktf)+dz(i,ktf,j)
1697
1698	CALL mym_initialize ( kts,kte,&
1699	&dz(i,kts:kte,j), zw, &
1700	&u(i,kts:kte,j), v(i,kts:kte,j), &
1701	&thl, sqv,&
1702	&ust(i,j), rmol(i,j),&
1703	&Qke(i,kts:kte,j), Tsq(i,kts:kte,j), &
1704	&Qsq(i,kts:kte,j), Cov(i,kts:kte,j))
1705
1706	! pblh(i,j)=pblh_ref
1707
1708	!JOE-PBLH test begin
1709	PBLH(i,j)=0.
1710	thsfc = thetav(1)
1711	DO k = kts+1,kte
1712	IF (PBLH(i,j)==0. .AND. thetav(k)>(thsfc+0.5))THEN
1713	PBLH(i,j)=zw(k) - ((thetav(k)-(thsfc+0.5)) * &
1714	&(dz(i,k-1,j)/(MAX(thetav(k)-thetav(k-1),0.01))))
1715	ENDIF
1716	ENDDO
1717	!JOE--PBLH test end
1718
1719	ENDDO
1720	ENDDO
1721
1722	ENDIF
1723
1724	DO j=JTS,JTF
1725	DO i=ITS,ITF
1726	DO k=KTS,KTF
1727	sqv(k)=qv(i,k,j)/(1.+qv(i,k,j))
1728	sqc(k)=qc(i,k,j)/(1.+qc(i,k,j))
1729	sqw(k)=sqv(k)+sqc(k)
1730	thl(k)=th(i,k,j)-xlvcp/exner(i,k,j)*sqc(k)
1731	!JOE-PBLH test
1732	thetav(k)=th(i,k,j)(1.+0.61qv(i,k,j))
1733	!JOE-end
1734
1735	IF (k==kts) THEN
1736	zw(k)=0.
1737	ELSE
1738	zw(k)=zw(k-1)+dz(i,k-1,j)
1739	ENDIF
1740	ENDDO
1741
1742	zw(ktf+1)=zw(ktf)+dz(i,ktf,j)
1743
1744	sqcg=qcg(i,j)/(1.+qcg(i,j))
1745	cpm=cp(1.+0.8qv(i,kts,j))
1746
1747	! The exchange coefficient for cloud water is assumed to be the same as
1748	! that for heat. CH is multiplied by WSPD. See module_sf_mynn.F
1749	exnerg=(ps(i,j)/p1000mb)**rcp
1750	flt = hfx(i,j)/( rho(i,kts,j)*cpm ) &
1751	+xlvcpch(i,j)(sqc(kts)/exner(i,kts,j)-sqcg/exnerg)
1752	flq = qfx(i,j)/ rho(i,kts,j) &
1753	-ch(i,j)*(sqc(kts) -sqcg )
1754
1755	!!!!!
1756	zet = 0.5dz(i,kts,j)rmol(i,j)
1757	if ( zet >= 0.0 ) then
1758	pmz = 1.0 + (cphm_st-1.0) * zet
1759	phh = 1.0 + cphh_st * zet
1760	else
1761	pmz = 1.0/ (1.0-cphm_unstzet)*0.25 - zet
1762	phh = 1.0/SQRT(1.0-cphh_unst*zet)
1763	end if
1764	!!!!!
1765
1766	CALL mym_condensation ( kts,kte,&
1767	&dz(i,kts:kte,j), &
1768	&thl, sqw, &
1769	&p(i,kts:kte,j),exner(i,kts:kte,j), &
1770	&tsq(i,kts:kte,j), qsq(i,kts:kte,j), cov(i,kts:kte,j), &
1771	&Vt, Vq)
1772
1773	CALL mym_turbulence ( kts,kte,&
1774	&levflag, &
1775	&dz(i,kts:kte,j), zw, &
1776	&u(i,kts:kte,j), v(i,kts:kte,j), thl, sqc, sqw, &
1777	&qke(i,kts:kte,j), tsq(i,kts:kte,j), &
1778	&qsq(i,kts:kte,j), cov(i,kts:kte,j), &
1779	&vt, vq,&
1780	&rmol(i,j), flt, flq, &
1781	&El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc)
1782
1783
1784	CALL mym_predict (kts,kte,&
1785	&levflag, &
1786	&delt,&
1787	&dz(i,kts:kte,j), &
1788	&ust(i,j), flt, flq, pmz, phh, &
1789	&el, dfq, &
1790	&pdk, pdt, pdq, pdc,&
1791	&Qke(i,kts:kte,j), Tsq(i,kts:kte,j), &
1792	&Qsq(i,kts:kte,j), Cov(i,kts:kte,j))
1793
1794	CALL mynn_tendencies(kts,kte,&
1795	&levflag,grav_settling,&
1796	&delt,&
1797	&dz(i,kts:kte,j),&
1798	&u(i,kts:kte,j),v(i,kts:kte,j),&
1799	&th(i,kts:kte,j),qv(i,kts:kte,j),qc(i,kts:kte,j),&
1800	&p(i,kts:kte,j),exner(i,kts:kte,j),&
1801	&thl,sqv,sqc,sqw,&
1802	&ust(i,j),flt,flq,wspd(i,j),qcg(i,j),&
1803	&tsq(i,kts:kte,j),qsq(i,kts:kte,j),cov(i,kts:kte,j),&
1804	&tcd,qcd,&
1805	&dfm,dfh,dfq,&
1806	&Du(i,kts:kte,j),Dv(i,kts:kte,j),Dth(i,kts:kte,j),&
1807	&Dqv(i,kts:kte,j),Dqc(i,kts:kte,j))
1808
1809	CALL retrieve_exchange_coeffs(kts,kte,&
1810	&dfm,dfh,dfq,dz(i,kts:kte,j),&
1811	&K_m(i,kts:kte,j),K_h(i,kts:kte,j),K_q(i,kts:kte,j))
1812
1813
1814	!JOE-PBLH test begin
1815	PBLH(i,j)=0.
1816	thsfc = thetav(1)
1817	DO k = kts+1,kte
1818	IF (PBLH(i,j)==0. .AND. thetav(k)>(thsfc+0.5))THEN
1819	PBLH(i,j)=zw(k) - ((thetav(k)-(thsfc+0.5)) * &
1820	&(dz(i,k-1,j)/(MAX(thetav(k)-thetav(k-1),0.01))))
1821	ENDIF
1822	ENDDO
1823	!JOE--PBLH test end
1824
1825	! pblh(i,j)=pblh_ref
1826
1827	ENDDO
1828	ENDDO
1829
1830	END SUBROUTINE mynn_bl_driver
1831
1832	SUBROUTINE mynn_bl_init_driver(&
1833	&Du,Dv,Dth,&
1834	&Dqv,Dqc&
1835	&,RESTART,ALLOWED_TO_READ,LEVEL&
1836	&,IDS,IDE,JDS,JDE,KDS,KDE &
1837	&,IMS,IME,JMS,JME,KMS,KME &
1838	&,ITS,ITE,JTS,JTE,KTS,KTE)
1839
1840	LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
1841	INTEGER,INTENT(IN) :: LEVEL
1842
1843	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
1844	& IMS,IME,JMS,JME,KMS,KME, &
1845	& ITS,ITE,JTS,JTE,KTS,KTE
1846
1847
1848	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: &
1849	&Du,Dv,Dth,Dqv,Dqc
1850
1851	INTEGER :: I,J,K,ITF,JTF,KTF
1852
1853	JTF=MIN0(JTE,JDE-1)
1854	KTF=MIN0(KTE,KDE-1)
1855	ITF=MIN0(ITE,IDE-1)
1856
1857	IF(.NOT.RESTART)THEN
1858	DO J=JTS,JTF
1859	DO K=KTS,KTF
1860	DO I=ITS,ITF
1861	Du(i,k,j)=0.
1862	Dv(i,k,j)=0.
1863	Dth(i,k,j)=0.
1864	Dqv(i,k,j)=0.
1865	Dqc(i,k,j)=0.
1866	ENDDO
1867	ENDDO
1868	ENDDO
1869	ENDIF
1870
1871	mynn_level=level
1872
1873	END SUBROUTINE mynn_bl_init_driver
1874
1875	END MODULE module_bl_mynn

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

Download in other formats:

Original Format