Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_diffusion_em.F @ 14

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

Line
1	! WRF:MODEL_LAYER:PHYSICS
2
3	MODULE module_diffusion_em
4
5	USE module_configure
6	USE module_bc
7	USE module_state_description
8	USE module_big_step_utilities_em
9	USE module_model_constants
10	USE module_wrf_error
11
12	CONTAINS
13
14	!=======================================================================
15	!=======================================================================
16
17	SUBROUTINE cal_deform_and_div( config_flags, u, v, w, div, &
18	defor11, defor22, defor33, &
19	defor12, defor13, defor23, &
20	nba_rij, n_nba_rij, & !JDM
21	u_base, v_base, msfux, msfuy, &
22	msfvx, msfvy, msftx, msfty, &
23	rdx, rdy, dn, dnw, rdz, rdzw, &
24	fnm, fnp, cf1, cf2, cf3, zx, zy, &
25	ids, ide, jds, jde, kds, kde, &
26	ims, ime, jms, jme, kms, kme, &
27	its, ite, jts, jte, kts, kte )
28
29	! History: Sep 2003 Changes by Jason Knievel and George Bryan, NCAR
30	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
31	! ... ...
32
33	! Purpose: This routine calculates deformation and 3-d divergence.
34
35	! References: Klemp and Wilhelmson (JAS 1978)
36	! Chen and Dudhia (NCAR WRF physics report 2000)
37
38	!-----------------------------------------------------------------------
39	! Comments 10-MAR-05
40	! Equations 13a-f, Chen and Dudhia 2000, Appendix A:
41	! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
42	! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
43	! Eqn 13c: D33=defor33= 2 * partial dw/dz [SIMPLER FORM]
44	! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
45	! partial dpsi/dx * partial dv^/dpsi +
46	! partial dpsi/dy * partial du^/dpsi)
47	! Eqn 13e: D13=defor13= m^2 * (partial dw^/dX + partial dpsi/dx * partial dw^/dpsi)
48	! + partial du/dz [SIMPLER FORM]
49	! Eqn 13f: D23=defor23= m^2 * (partial dw^/dY + partial dpsi/dy * partial dw^/dpsi)
50	! + partial dv/dz [SIMPLER FORM]
51	!-----------------------------------------------------------------------
52	! Begin declarations.
53
54	IMPLICIT NONE
55
56	TYPE( grid_config_rec_type ), INTENT( IN ) &
57	:: config_flags
58
59	INTEGER, INTENT( IN ) &
60	:: ids, ide, jds, jde, kds, kde, &
61	ims, ime, jms, jme, kms, kme, &
62	its, ite, jts, jte, kts, kte
63
64	REAL, INTENT( IN ) &
65	:: rdx, rdy, cf1, cf2, cf3
66
67	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
68	:: fnm, fnp, dn, dnw, u_base, v_base
69
70	REAL, DIMENSION( ims:ime , jms:jme ), INTENT( IN ) &
71	:: msfux, msfuy, msfvx, msfvy, msftx, msfty
72
73	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
74	:: u, v, w, zx, zy, rdz, rdzw
75
76	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
77	:: defor11, defor22, defor33, defor12, defor13, defor23, div
78
79	INTEGER, INTENT( IN ) :: n_nba_rij !JDM
80
81	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_rij), INTENT(INOUT) & !JDM
82	:: nba_rij
83
84
85	! Local variables.
86
87	INTEGER &
88	:: i, j, k, ktf, ktes1, ktes2, i_start, i_end, j_start, j_end
89
90	REAL &
91	:: tmp, tmpzx, tmpzy, tmpzeta_z, cft1, cft2
92
93	REAL, DIMENSION( its:ite, jts:jte ) &
94	:: mm, zzavg, zeta_zd12
95
96	REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) &
97	:: tmp1, hat, hatavg
98
99	! End declarations.
100	!-----------------------------------------------------------------------
101
102	! Comments 10-MAR-2005
103	! Treat all differentials as 'div-style' [or 'curl-style'],
104	! i.e., du/dX becomes (in map coordinate space) mxmy d(u/my)/dx,
105	! NB - all equations referred to here are from Chen and Dudhia 2002, from the
106	! WRF physics documents web pages:
107	! http://www.mmm.ucar.edu/wrf/users/docs/wrf-doc-physics.pdf
108
109	!=======================================================================
110	! In the following section, calculate 3-d divergence and the first three
111	! (defor11, defor22, defor33) of six deformation terms.
112
113	ktes1 = kte-1
114	ktes2 = kte-2
115
116	cft2 = - 0.5 * dnw(ktes1) / dn(ktes1)
117	cft1 = 1.0 - cft2
118
119	ktf = MIN( kte, kde-1 )
120
121	i_start = its
122	i_end = MIN( ite, ide-1 )
123	j_start = jts
124	j_end = MIN( jte, jde-1 )
125
126	! Square the map scale factor.
127
128	DO j = j_start, j_end
129	DO i = i_start, i_end
130	mm(i,j) = msftx(i,j) * msfty(i,j)
131	END DO
132	END DO
133
134	!-----------------------------------------------------------------------
135	! Calculate du/dx.
136
137	! Apply a coordinate transformation to zonal velocity, u.
138
139	DO j = j_start, j_end
140	DO k = kts, ktf
141	DO i = i_start, i_end+1
142	hat(i,k,j) = u(i,k,j) / msfuy(i,j)
143	END DO
144	END DO
145	END DO
146
147	! Average in x and z.
148
149	DO j=j_start,j_end
150	DO k=kts+1,ktf
151	DO i=i_start,i_end
152	hatavg(i,k,j) = 0.5 * &
153	( fnm(k) * ( hat(i,k ,j) + hat(i+1, k,j) ) + &
154	fnp(k) * ( hat(i,k-1,j) + hat(i+1,k-1,j) ) )
155	END DO
156	END DO
157	END DO
158
159	! Extrapolate to top and bottom of domain (to w levels).
160
161	DO j = j_start, j_end
162	DO i = i_start, i_end
163	hatavg(i,1,j) = 0.5 * ( &
164	cf1 * hat(i ,1,j) + &
165	cf2 * hat(i ,2,j) + &
166	cf3 * hat(i ,3,j) + &
167	cf1 * hat(i+1,1,j) + &
168	cf2 * hat(i+1,2,j) + &
169	cf3 * hat(i+1,3,j) )
170	hatavg(i,kte,j) = 0.5 * ( &
171	cft1 * ( hat(i,ktes1,j) + hat(i+1,ktes1,j) ) + &
172	cft2 * ( hat(i,ktes2,j) + hat(i+1,ktes2,j) ) )
173	END DO
174	END DO
175
176	! Comments 10-MAR-05
177	! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
178	! Below, D11 is set = 2*tmp1
179	! => tmp1 = m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
180	! tmpzx = averaged value of dpsi/dx (=zx)
181
182	DO j = j_start, j_end
183	DO k = kts, ktf
184	DO i = i_start, i_end
185	tmpzx = 0.25 * ( &
186	zx(i,k ,j) + zx(i+1,k ,j) + &
187	zx(i,k+1,j) + zx(i+1,k+1,j) )
188	tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) tmpzx rdzw(i,k,j)
189	! tmp1 to here = partial dpsi/dx * partial du^/dpsi:
190	END DO
191	END DO
192	END DO
193
194	DO j = j_start, j_end
195	DO k = kts, ktf
196	DO i = i_start, i_end
197	tmp1(i,k,j) = mm(i,j) * ( rdx * ( hat(i+1,k,j) - hat(i,k,j) ) - &
198	tmp1(i,k,j))
199	END DO
200	END DO
201	END DO
202
203	! End calculation of du/dx.
204	!-----------------------------------------------------------------------
205
206	!-----------------------------------------------------------------------
207	! Calculate defor11 (2*du/dx).
208	! Comments 10-MAR-05
209	! Eqn 13a: D11=defor11= 2 m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi)
210	! = 2*tmp1
211
212	DO j = j_start, j_end
213	DO k = kts, ktf
214	DO i = i_start, i_end
215	defor11(i,k,j) = 2.0 * tmp1(i,k,j)
216	END DO
217	END DO
218	END DO
219
220	! End calculation of defor11.
221	!-----------------------------------------------------------------------
222
223	!-----------------------------------------------------------------------
224	! Calculate zonal divergence (du/dx) and add it to the divergence array.
225
226	DO j = j_start, j_end
227	DO k = kts, ktf
228	DO i = i_start, i_end
229	div(i,k,j) = tmp1(i,k,j)
230	END DO
231	END DO
232	END DO
233
234	! End calculation of zonal divergence.
235	!-----------------------------------------------------------------------
236
237	!-----------------------------------------------------------------------
238	! Calculate dv/dy.
239
240	! Apply a coordinate transformation to meridional velocity, v.
241
242	DO j = j_start, j_end+1
243	DO k = kts, ktf
244	DO i = i_start, i_end
245	! Because msfvx at the poles will be undefined (1./0.), we will have
246	! trouble. But we are OK since v at the poles is 0., and that takes
247	! precedence in this case.
248	IF ((config_flags%polar) .AND. ((j == jds) .OR. (j == jde))) THEN
249	hat(i,k,j) = 0.
250	ELSE ! normal code
251	hat(i,k,j) = v(i,k,j) / msfvx(i,j)
252	ENDIF
253	END DO
254	END DO
255	END DO
256
257	! Account for the slope in y of eta surfaces.
258
259	DO j=j_start,j_end
260	DO k=kts+1,ktf
261	DO i=i_start,i_end
262	hatavg(i,k,j) = 0.5 * ( &
263	fnm(k) * ( hat(i,k ,j) + hat(i,k ,j+1) ) + &
264	fnp(k) * ( hat(i,k-1,j) + hat(i,k-1,j+1) ) )
265	END DO
266	END DO
267	END DO
268
269	! Extrapolate to top and bottom of domain (to w levels).
270
271	DO j = j_start, j_end
272	DO i = i_start, i_end
273	hatavg(i,1,j) = 0.5 * ( &
274	cf1 * hat(i,1,j ) + &
275	cf2 * hat(i,2,j ) + &
276	cf3 * hat(i,3,j ) + &
277	cf1 * hat(i,1,j+1) + &
278	cf2 * hat(i,2,j+1) + &
279	cf3 * hat(i,3,j+1) )
280	hatavg(i,kte,j) = 0.5 * ( &
281	cft1 * ( hat(i,ktes1,j) + hat(i,ktes1,j+1) ) + &
282	cft2 * ( hat(i,ktes2,j) + hat(i,ktes2,j+1) ) )
283	END DO
284	END DO
285
286	! Comments 10-MAR-05
287	! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
288	! Below, D22 is set = 2*tmp1
289	! => tmp1 = m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
290	! tmpzy = averaged value of dpsi/dy (=zy)
291
292	DO j = j_start, j_end
293	DO k = kts, ktf
294	DO i = i_start, i_end
295	tmpzy = 0.25 * ( &
296	zy(i,k ,j) + zy(i,k ,j+1) + &
297	zy(i,k+1,j) + zy(i,k+1,j+1) )
298	tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * tmpzy * rdzw(i,k,j)
299	! tmp1 to here = partial dpsi/dy * partial dv^/dpsi:
300	END DO
301	END DO
302	END DO
303
304	DO j = j_start, j_end
305	DO k = kts, ktf
306	DO i = i_start, i_end
307	tmp1(i,k,j) = mm(i,j) * ( &
308	rdy * ( hat(i,k,j+1) - hat(i,k,j) ) - tmp1(i,k,j) )
309	END DO
310	END DO
311	END DO
312
313	! End calculation of dv/dy.
314	!-----------------------------------------------------------------------
315
316	!-----------------------------------------------------------------------
317	! Calculate defor22 (2*dv/dy).
318	! Comments 10-MAR-05
319	! Eqn 13b: D22=defor22= 2 m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi)
320	! = 2*tmp1
321
322	DO j = j_start, j_end
323	DO k = kts, ktf
324	DO i = i_start, i_end
325	defor22(i,k,j) = 2.0 * tmp1(i,k,j)
326	END DO
327	END DO
328	END DO
329
330	! End calculation of defor22.
331	!-----------------------------------------------------------------------
332
333	!-----------------------------------------------------------------------
334	! Calculate meridional divergence (dv/dy) and add it to the divergence
335	! array.
336
337	DO j = j_start, j_end
338	DO k = kts, ktf
339	DO i = i_start, i_end
340	div(i,k,j) = div(i,k,j) + tmp1(i,k,j)
341	END DO
342	END DO
343	END DO
344
345	! End calculation of meridional divergence.
346	!-----------------------------------------------------------------------
347
348	!-----------------------------------------------------------------------
349	! Comments 10-MAR-05
350	! Eqn 13c: D33=defor33= 2 * partial dw/dz
351	! Below, D33 is set = 2*tmp1
352	! => tmp1 = partial dw/dz
353
354	! Calculate dw/dz.
355
356	DO j = j_start, j_end
357	DO k = kts, ktf
358	DO i = i_start, i_end
359	tmp1(i,k,j) = ( w(i,k+1,j) - w(i,k,j) ) * rdzw(i,k,j)
360	END DO
361	END DO
362	END DO
363
364	! End calculation of dw/dz.
365	!-----------------------------------------------------------------------
366
367	!-----------------------------------------------------------------------
368	! Calculate defor33 (2*dw/dz).
369
370	DO j = j_start, j_end
371	DO k = kts, ktf
372	DO i = i_start, i_end
373	defor33(i,k,j) = 2.0 * tmp1(i,k,j)
374	END DO
375	END DO
376	END DO
377
378	! End calculation of defor33.
379	!-----------------------------------------------------------------------
380
381	!-----------------------------------------------------------------------
382	! Calculate vertical divergence (dw/dz) and add it to the divergence
383	! array.
384
385	DO j = j_start, j_end
386	DO k = kts, ktf
387	DO i = i_start, i_end
388	div(i,k,j) = div(i,k,j) + tmp1(i,k,j)
389	END DO
390	END DO
391	END DO
392
393	! End calculation of vertical divergence.
394	!-----------------------------------------------------------------------
395
396	! Three-dimensional divergence is now finished and values are in array
397	! "div." Also, the first three (defor11, defor22, defor33) of six
398	! deformation terms are now calculated at pressure points.
399	!=======================================================================
400
401	! Comments 10-MAR-2005
402	! Treat all differentials as 'div-style' [or 'curl-style'],
403	! i.e., du/dY becomes (in map coordinate space) mxmy d(u/mx)/dy,
404	! dv/dX becomes (in map coordinate space) mxmy d(v/my)/dx,
405	! (see e.g. Haltiner and Williams p. 441)
406
407	!=======================================================================
408	! Calculate the final three deformations (defor12, defor13, defor23) at
409	! vorticity points.
410
411	i_start = its
412	i_end = ite
413	j_start = jts
414	j_end = jte
415
416	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
417	config_flags%nested) i_start = MAX( ids+1, its )
418	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
419	config_flags%nested) i_end = MIN( ide-1, ite )
420	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
421	config_flags%nested) j_start = MAX( jds+1, jts )
422	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
423	config_flags%nested) j_end = MIN( jde-1, jte )
424	IF ( config_flags%periodic_x ) i_start = its
425	IF ( config_flags%periodic_x ) i_end = ite
426
427
428	!-----------------------------------------------------------------------
429	! Calculate du/dy.
430
431	! First, calculate an average mapscale factor.
432
433	! Comments 10-MAR-05
434	! du/dy => need u map scale factor in x (which is defined at u points)
435	! averaged over j and j-1
436	! dv/dx => need v map scale factor in y (which is defined at v points)
437	! averaged over i and i-1
438
439	DO j = j_start, j_end
440	DO i = i_start, i_end
441	mm(i,j) = 0.25 * ( msfux(i,j-1) + msfux(i,j) ) * ( msfvy(i-1,j) + msfvy(i,j) )
442	END DO
443	END DO
444
445	! Apply a coordinate transformation to zonal velocity, u.
446
447	DO j =j_start-1, j_end
448	DO k =kts, ktf
449	DO i =i_start, i_end
450	! Fixes to set_physical_bc2/3d for polar boundary conditions
451	! remove issues with loop over j
452	hat(i,k,j) = u(i,k,j) / msfux(i,j)
453	END DO
454	END DO
455	END DO
456
457	! Average in y and z.
458
459	DO j=j_start,j_end
460	DO k=kts+1,ktf
461	DO i=i_start,i_end
462	hatavg(i,k,j) = 0.5 * ( &
463	fnm(k) * ( hat(i,k ,j-1) + hat(i,k ,j) ) + &
464	fnp(k) * ( hat(i,k-1,j-1) + hat(i,k-1,j) ) )
465	END DO
466	END DO
467	END DO
468
469	! Extrapolate to top and bottom of domain (to w levels).
470
471	DO j = j_start, j_end
472	DO i = i_start, i_end
473	hatavg(i,1,j) = 0.5 * ( &
474	cf1 * hat(i,1,j-1) + &
475	cf2 * hat(i,2,j-1) + &
476	cf3 * hat(i,3,j-1) + &
477	cf1 * hat(i,1,j ) + &
478	cf2 * hat(i,2,j ) + &
479	cf3 * hat(i,3,j ) )
480	hatavg(i,kte,j) = 0.5 * ( &
481	cft1 * ( hat(i,ktes1,j-1) + hat(i,ktes1,j) ) + &
482	cft2 * ( hat(i,ktes2,j-1) + hat(i,ktes2,j) ) )
483	END DO
484	END DO
485
486	! tmpzy = averaged value of dpsi/dy (=zy) on vorticity grid
487	! tmp1 = partial dpsi/dy * partial du^/dpsi
488	DO j = j_start, j_end
489	DO k = kts, ktf
490	DO i = i_start, i_end
491	tmpzy = 0.25 * ( &
492	zy(i-1,k ,j) + zy(i,k ,j) + &
493	zy(i-1,k+1,j) + zy(i,k+1,j) )
494	tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * &
495	0.25 * tmpzy * ( rdzw(i,k,j) + rdzw(i-1,k,j) + &
496	rdzw(i-1,k,j-1) + rdzw(i,k,j-1) )
497	END DO
498	END DO
499	END DO
500
501	! End calculation of du/dy.
502	!----------------------------------------------------------------------
503
504	!-----------------------------------------------------------------------
505	! Add the first term to defor12 (du/dy+dv/dx) at vorticity points.
506
507	! Comments 10-MAR-05
508	! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
509	! partial dpsi/dx * partial dv^/dpsi +
510	! partial dpsi/dy * partial du^/dpsi)
511	! Here deal with m^2 * (partial du^/dY + partial dpsi/dy * partial du^/dpsi)
512	! Still need to add v^ terms:
513	! m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi)
514
515	DO j = j_start, j_end
516	DO k = kts, ktf
517	DO i = i_start, i_end
518	defor12(i,k,j) = mm(i,j) * ( &
519	rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) )
520	END DO
521	END DO
522	END DO
523
524	! End addition of the first term to defor12.
525	!-----------------------------------------------------------------------
526
527	!-----------------------------------------------------------------------
528	! Calculate dv/dx.
529
530	! Apply a coordinate transformation to meridional velocity, v.
531
532	DO j = j_start, j_end
533	DO k = kts, ktf
534	DO i = i_start-1, i_end
535	hat(i,k,j) = v(i,k,j) / msfvy(i,j)
536	END DO
537	END DO
538	END DO
539
540	! Account for the slope in x of eta surfaces.
541
542	DO j = j_start, j_end
543	DO k = kts+1, ktf
544	DO i = i_start, i_end
545	hatavg(i,k,j) = 0.5 * ( &
546	fnm(k) * ( hat(i-1,k ,j) + hat(i,k ,j) ) + &
547	fnp(k) * ( hat(i-1,k-1,j) + hat(i,k-1,j) ) )
548	END DO
549	END DO
550	END DO
551
552	! Extrapolate to top and bottom of domain (to w levels).
553
554	DO j = j_start, j_end
555	DO i = i_start, i_end
556	hatavg(i,1,j) = 0.5 * ( &
557	cf1 * hat(i-1,1,j) + &
558	cf2 * hat(i-1,2,j) + &
559	cf3 * hat(i-1,3,j) + &
560	cf1 * hat(i ,1,j) + &
561	cf2 * hat(i ,2,j) + &
562	cf3 * hat(i ,3,j) )
563	hatavg(i,kte,j) = 0.5 * ( &
564	cft1 * ( hat(i,ktes1,j) + hat(i-1,ktes1,j) ) + &
565	cft2 * ( hat(i,ktes2,j) + hat(i-1,ktes2,j) ) )
566	END DO
567	END DO
568
569	! Fixes to set_physical_bc2/3d have made any check for polar B.C.'s
570	! unnecessary in this place. zx, rdzw, and hatavg are all defined
571	! in places they need to be and the values at the poles are replications
572	! of the values one grid point in, so the averaging over j and j-1 works
573	! to act as just using the value at j or j-1 (with out extra code).
574	!
575	! tmpzx = averaged value of dpsi/dx (=zx) on vorticity grid
576	! tmp1 = partial dpsi/dx * partial dv^/dpsi
577	DO j = j_start, j_end
578	DO k = kts, ktf
579	DO i = i_start, i_end
580	tmpzx = 0.25 * ( &
581	zx(i,k ,j-1) + zx(i,k ,j) + &
582	zx(i,k+1,j-1) + zx(i,k+1,j) )
583	tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * &
584	0.25 * tmpzx * ( rdzw(i,k,j) + rdzw(i,k,j-1) + &
585	rdzw(i-1,k,j-1) + rdzw(i-1,k,j) )
586	END DO
587	END DO
588	END DO
589
590	! End calculation of dv/dx.
591	!-----------------------------------------------------------------------
592
593	!-----------------------------------------------------------------------
594	! Add the second term to defor12 (du/dy+dv/dx) at vorticity points.
595
596	! Comments 10-MAR-05
597	! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY +
598	! partial dpsi/dx * partial dv^/dpsi +
599	! partial dpsi/dy * partial du^/dpsi)
600	! Here adding v^ terms:
601	! m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi)
602
603	IF ( config_flags%sfs_opt .GT. 0 ) THEN ! NBA--
604
605	!JDM____________________________________________________________________
606	!
607	! s12 = du/dy + dv/dx
608	! = (du/dy - dz/dydu/dz) + (dv/dx - dz/dxdv/dz)
609	! ______defor12______ ___tmp1___
610	!
611	! r12 = du/dy - dv/dx
612	! = (du/dy - dz/dydu/dz) - (dv/dx - dz/dxdv/dz)
613	! ______defor12______ ___tmp1___
614	!_______________________________________________________________________
615
616
617	DO j = j_start, j_end
618	DO k = kts, ktf
619	DO i = i_start, i_end
620
621	nba_rij(i,k,j,P_r12) = defor12(i,k,j) - &
622	mm(i,j) * ( &
623	rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
624
625	defor12(i,k,j) = defor12(i,k,j) + &
626	mm(i,j) * ( &
627	rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
628	END DO
629	END DO
630	END DO
631
632	! End addition of the second term to defor12.
633	!-----------------------------------------------------------------------
634
635	!-----------------------------------------------------------------------
636	! Update the boundary for defor12 (might need to change later).
637
638	IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1 ) THEN
639	DO j = jts, jte
640	DO k = kts, kte
641	defor12(ids,k,j) = defor12(ids+1,k,j)
642	nba_rij(ids,k,j,P_r12) = nba_rij(ids+1,k,j,P_r12)
643	END DO
644	END DO
645	END IF
646
647	IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
648	DO k = kts, kte
649	DO i = its, ite
650	defor12(i,k,jds) = defor12(i,k,jds+1)
651	nba_rij(i,k,jds,P_r12) = nba_rij(i,k,jds+1,P_r12)
652	END DO
653	END DO
654	END IF
655
656	IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
657	DO j = jts, jte
658	DO k = kts, kte
659	defor12(ide,k,j) = defor12(ide-1,k,j)
660	nba_rij(ide,k,j,P_r12) = nba_rij(ide-1,k,j,P_r12)
661	END DO
662	END DO
663	END IF
664
665	IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
666	DO k = kts, kte
667	DO i = its, ite
668	defor12(i,k,jde) = defor12(i,k,jde-1)
669	nba_rij(i,k,jde,P_r12) = nba_rij(i,k,jde-1,P_r12)
670	END DO
671	END DO
672	END IF
673
674	ELSE ! NOT NBA--------------------------------------------------------
675
676	DO j = j_start, j_end
677	DO k = kts, ktf
678	DO i = i_start, i_end
679	defor12(i,k,j) = defor12(i,k,j) + &
680	mm(i,j) * ( &
681	rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
682	END DO
683	END DO
684	END DO
685
686	! End addition of the second term to defor12.
687	!-----------------------------------------------------------------------
688
689	!-----------------------------------------------------------------------
690	! Update the boundary for defor12 (might need to change later).
691
692	IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1 ) THEN
693	DO j = jts, jte
694	DO k = kts, kte
695	defor12(ids,k,j) = defor12(ids+1,k,j)
696	END DO
697	END DO
698	END IF
699
700	IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
701	DO k = kts, kte
702	DO i = its, ite
703	defor12(i,k,jds) = defor12(i,k,jds+1)
704	END DO
705	END DO
706	END IF
707
708	IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
709	DO j = jts, jte
710	DO k = kts, kte
711	defor12(ide,k,j) = defor12(ide-1,k,j)
712	END DO
713	END DO
714	END IF
715
716	IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
717	DO k = kts, kte
718	DO i = its, ite
719	defor12(i,k,jde) = defor12(i,k,jde-1)
720	END DO
721	END DO
722	END IF
723
724	ENDIF ! NBA-----------------------------------------------------------
725
726	! End update of boundary for defor12.
727	!-----------------------------------------------------------------------
728
729	! Comments 10-MAR-05
730	! Further deformation terms not needed for 2-dimensional Smagorinsky diffusion,
731	! so those terms have not been dealt with yet.
732	! A "y" has simply been added to all map scale factors to allow the model to
733	! compile without errors.
734
735	!-----------------------------------------------------------------------
736	! Calculate dw/dx.
737
738	i_start = its
739	i_end = MIN( ite, ide-1 )
740	j_start = jts
741	j_end = MIN( jte, jde-1 )
742
743	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
744	config_flags%nested) i_start = MAX( ids+1, its )
745	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
746	config_flags%nested) j_start = MAX( jds+1, jts )
747
748	IF ( config_flags%periodic_x ) i_start = its
749	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide )
750	IF ( config_flags%periodic_y ) j_end = MIN( jte, jde )
751
752	! Square the mapscale factor.
753
754	DO j = jts, jte
755	DO i = its, ite
756	mm(i,j) = msfux(i,j) * msfuy(i,j)
757	END DO
758	END DO
759
760	! Apply a coordinate transformation to vertical velocity, w. This is for both
761	! defor13 and defor23.
762
763	DO j = j_start, j_end
764	DO k = kts, kte
765	DO i = i_start, i_end
766	hat(i,k,j) = w(i,k,j) / msfty(i,j)
767	END DO
768	END DO
769	END DO
770
771	i = i_start-1
772	DO j = j_start, MIN( jte, jde-1 )
773	DO k = kts, kte
774	hat(i,k,j) = w(i,k,j) / msfty(i,j)
775	END DO
776	END DO
777
778	j = j_start-1
779	DO k = kts, kte
780	DO i = i_start, MIN( ite, ide-1 )
781	hat(i,k,j) = w(i,k,j) / msfty(i,j)
782	END DO
783	END DO
784
785	! QUESTION: What is this for?
786
787	DO j = j_start, j_end
788	DO k = kts, ktf
789	DO i = i_start, i_end
790	hatavg(i,k,j) = 0.25 * ( &
791	hat(i ,k ,j) + &
792	hat(i ,k+1,j) + &
793	hat(i-1,k ,j) + &
794	hat(i-1,k+1,j) )
795	END DO
796	END DO
797	END DO
798
799	! Calculate dw/dx.
800
801	DO j = j_start, j_end
802	DO k = kts+1, ktf
803	DO i = i_start, i_end
804	tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zx(i,k,j) * &
805	0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
806	END DO
807	END DO
808	END DO
809
810	! End calculation of dw/dx.
811	!-----------------------------------------------------------------------
812
813	!-----------------------------------------------------------------------
814	! Add the first term (dw/dx) to defor13 (dw/dx+du/dz) at vorticity
815	! points.
816
817	DO j = j_start, j_end
818	DO k = kts+1, ktf
819	DO i = i_start, i_end
820	defor13(i,k,j) = mm(i,j) * ( &
821	rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
822	END DO
823	END DO
824	END DO
825
826	DO j = j_start, j_end
827	DO i = i_start, i_end
828	defor13(i,kts,j ) = 0.0
829	defor13(i,ktf+1,j) = 0.0
830	END DO
831	END DO
832
833	! End addition of the first term to defor13.
834	!-----------------------------------------------------------------------
835
836	!-----------------------------------------------------------------------
837	! Calculate du/dz.
838
839	IF ( config_flags%mix_full_fields ) THEN
840
841	DO j = j_start, j_end
842	DO k = kts+1, ktf
843	DO i = i_start, i_end
844	tmp1(i,k,j) = ( u(i,k,j) - u(i,k-1,j) ) * &
845	0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
846	END DO
847	END DO
848	END DO
849
850	ELSE
851
852	DO j = j_start, j_end
853	DO k = kts+1, ktf
854	DO i = i_start, i_end
855	tmp1(i,k,j) = ( u(i,k,j) - u_base(k) - u(i,k-1,j) + u_base(k-1) ) * &
856	0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) )
857	END DO
858	END DO
859	END DO
860
861	END IF
862
863	!-----------------------------------------------------------------------
864	! Add the second term (du/dz) to defor13 (dw/dx+du/dz) at vorticity
865	! points.
866
867
868	IF ( config_flags%sfs_opt .GT. 0 ) THEN ! NBA--
869
870	!JDM____________________________________________________________________
871	!
872	! s13 = du/dz + dw/dx
873	! = du/dz + (dw/dx - dz/dx*dw/dz)
874	! = tmp1 + ______defor13______
875	!
876	! r13 = du/dz - dw/dx
877	! = du/dz - (dw/dx - dz/dx*dw/dz)
878	! = tmp1 - ______defor13______
879	!_______________________________________________________________________
880
881	DO j = j_start, j_end
882	DO k = kts+1, ktf
883	DO i = i_start, i_end
884	nba_rij(i,k,j,P_r13) = tmp1(i,k,j) - defor13(i,k,j)
885	defor13(i,k,j) = defor13(i,k,j) + tmp1(i,k,j)
886	END DO
887	END DO
888	END DO
889
890	DO j = j_start, j_end !change for different surface B. C.
891	DO i = i_start, i_end
892	nba_rij(i,kts ,j,P_r13) = 0.0
893	nba_rij(i,ktf+1,j,P_r13) = 0.0
894	END DO
895	END DO
896
897	ELSE ! NOT NBA--------------------------------------------------------
898
899	DO j = j_start, j_end
900	DO k = kts+1, ktf
901	DO i = i_start, i_end
902	defor13(i,k,j) = defor13(i,k,j) + tmp1(i,k,j)
903	END DO
904	END DO
905	END DO
906
907	ENDIF ! NBA-----------------------------------------------------------
908
909	! End addition of the second term to defor13.
910	!-----------------------------------------------------------------------
911
912	!-----------------------------------------------------------------------
913	! Calculate dw/dy.
914
915	i_start = its
916	i_end = MIN( ite, ide-1 )
917	j_start = jts
918	j_end = MIN( jte, jde-1 )
919
920	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
921	config_flags%nested) i_start = MAX( ids+1, its )
922	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
923	config_flags%nested) j_start = MAX( jds+1, jts )
924	IF ( config_flags%periodic_y ) j_end = MIN( jte, jde )
925	IF ( config_flags%periodic_x ) i_start = its
926	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
927
928	! Square mapscale factor.
929
930	DO j = jts, jte
931	DO i = its, ite
932	mm(i,j) = msfvx(i,j) * msfvy(i,j)
933	END DO
934	END DO
935
936	! Apply a coordinate transformation to vertical velocity, w. Added by CW 7/19/07
937
938	DO j = j_start, j_end
939	DO k = kts, kte
940	DO i = i_start, i_end
941	hat(i,k,j) = w(i,k,j) / msftx(i,j)
942	END DO
943	END DO
944	END DO
945
946	i = i_start-1
947	DO j = j_start, MIN( jte, jde-1 )
948	DO k = kts, kte
949	hat(i,k,j) = w(i,k,j) / msftx(i,j)
950	END DO
951	END DO
952
953	j = j_start-1
954	DO k = kts, kte
955	DO i = i_start, MIN( ite, ide-1 )
956	hat(i,k,j) = w(i,k,j) / msftx(i,j)
957	END DO
958	END DO
959
960	! QUESTION: What is this for?
961
962	DO j = j_start, j_end
963	DO k = kts, ktf
964	DO i = i_start, i_end
965	hatavg(i,k,j) = 0.25 * ( &
966	hat(i,k ,j ) + &
967	hat(i,k+1,j ) + &
968	hat(i,k ,j-1) + &
969	hat(i,k+1,j-1) )
970	END DO
971	END DO
972	END DO
973
974	! Calculate dw/dy and store in tmp1.
975
976	DO j = j_start, j_end
977	DO k = kts+1, ktf
978	DO i = i_start, i_end
979	tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zy(i,k,j) * &
980	0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
981	END DO
982	END DO
983	END DO
984
985	! End calculation of dw/dy.
986	!-----------------------------------------------------------------------
987
988	!-----------------------------------------------------------------------
989	! Add the first term (dw/dy) to defor23 (dw/dy+dv/dz) at vorticity
990	! points.
991
992	DO j = j_start, j_end
993	DO k = kts+1, ktf
994	DO i = i_start, i_end
995	defor23(i,k,j) = mm(i,j) * ( &
996	rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) )
997	END DO
998	END DO
999	END DO
1000
1001	DO j = j_start, j_end
1002	DO i = i_start, i_end
1003	defor23(i,kts,j ) = 0.0
1004	defor23(i,ktf+1,j) = 0.0
1005	END DO
1006	END DO
1007
1008	! End addition of the first term to defor23.
1009	!-----------------------------------------------------------------------
1010
1011	!-----------------------------------------------------------------------
1012	! Calculate dv/dz.
1013
1014	IF ( config_flags%mix_full_fields ) THEN
1015
1016	DO j = j_start, j_end
1017	DO k = kts+1, ktf
1018	DO i = i_start, i_end
1019	tmp1(i,k,j) = ( v(i,k,j) - v(i,k-1,j) ) * &
1020	0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
1021	END DO
1022	END DO
1023	END DO
1024
1025	ELSE
1026
1027	DO j = j_start, j_end
1028	DO k = kts+1, ktf
1029	DO i = i_start, i_end
1030	tmp1(i,k,j) = ( v(i,k,j) - v_base(k) - v(i,k-1,j) + v_base(k-1) ) * &
1031	0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) )
1032	END DO
1033	END DO
1034	END DO
1035
1036	END IF
1037
1038	! End calculation of dv/dz.
1039	!-----------------------------------------------------------------------
1040
1041	!-----------------------------------------------------------------------
1042	! Add the second term (dv/dz) to defor23 (dw/dy+dv/dz) at vorticity
1043	! points.
1044
1045	IF ( config_flags%sfs_opt .GT. 0 ) THEN ! NBA--
1046
1047	!JDM___________________________________________________________________
1048	!
1049	! s23 = dv/dz + dw/dy
1050	! = dv/dz + (dw/dy - dz/dy*dw/dz)
1051	! tmp1 + ______defor23______
1052	!
1053	! r23 = dv/dz - dw/dy
1054	! = dv/dz - (dw/dy - dz/dy*dw/dz)
1055	! = tmp1 - ______defor23______
1056
1057	! Add tmp1 to defor23.
1058
1059	DO j = j_start, j_end
1060	DO k = kts+1, ktf
1061	DO i = i_start, i_end
1062	nba_rij(i,k,j,P_r23) = tmp1(i,k,j) - defor23(i,k,j)
1063	defor23(i,k,j) = defor23(i,k,j) + tmp1(i,k,j)
1064	END DO
1065	END DO
1066	END DO
1067
1068	DO j = j_start, j_end
1069	DO i = i_start, i_end
1070	nba_rij(i,kts ,j,P_r23) = 0.0
1071	nba_rij(i,ktf+1,j,P_r23) = 0.0
1072	END DO
1073	END DO
1074
1075	! End addition of the second term to defor23.
1076	!-----------------------------------------------------------------------
1077
1078	!-----------------------------------------------------------------------
1079	! Update the boundary for defor13 and defor23 (might need to change
1080	! later).
1081
1082	IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1) THEN
1083	DO j = jts, jte
1084	DO k = kts, kte
1085	defor13(ids,k,j) = defor13(ids+1,k,j)
1086	defor23(ids,k,j) = defor23(ids+1,k,j)
1087	nba_rij(ids,k,j,P_r13) = nba_rij(ids+1,k,j,P_r13)
1088	nba_rij(ids,k,j,P_r23) = nba_rij(ids+1,k,j,P_r23)
1089	END DO
1090	END DO
1091	END IF
1092
1093	IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
1094	DO k = kts, kte
1095	DO i = its, ite
1096	defor13(i,k,jds) = defor13(i,k,jds+1)
1097	defor23(i,k,jds) = defor23(i,k,jds+1)
1098	nba_rij(i,k,jds,P_r13) = nba_rij(i,k,jds+1,P_r13)
1099	nba_rij(i,k,jds,P_r23) = nba_rij(i,k,jds+1,P_r23)
1100	END DO
1101	END DO
1102	END IF
1103
1104	IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
1105	DO j = jts, jte
1106	DO k = kts, kte
1107	defor13(ide,k,j) = defor13(ide-1,k,j)
1108	defor23(ide,k,j) = defor23(ide-1,k,j)
1109	nba_rij(ide,k,j,P_r13) = nba_rij(ide-1,k,j,P_r13)
1110	nba_rij(ide,k,j,P_r23) = nba_rij(ide-1,k,j,P_r23)
1111	END DO
1112	END DO
1113	END IF
1114
1115	IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
1116	DO k = kts, kte
1117	DO i = its, ite
1118	defor13(i,k,jde) = defor13(i,k,jde-1)
1119	defor23(i,k,jde) = defor23(i,k,jde-1)
1120	nba_rij(i,k,jde,P_r13) = nba_rij(i,k,jde-1,P_r13)
1121	nba_rij(i,k,jde,P_r23) = nba_rij(i,k,jde-1,P_r23)
1122	END DO
1123	END DO
1124	END IF
1125
1126	ELSE ! NOT NBA--------------------------------------------------------
1127
1128	! Add tmp1 to defor23.
1129
1130	DO j = j_start, j_end
1131	DO k = kts+1, ktf
1132	DO i = i_start, i_end
1133	defor23(i,k,j) = defor23(i,k,j) + tmp1(i,k,j)
1134	END DO
1135	END DO
1136	END DO
1137
1138	! End addition of the second term to defor23.
1139	!-----------------------------------------------------------------------
1140
1141	!-----------------------------------------------------------------------
1142	! Update the boundary for defor13 and defor23 (might need to change
1143	! later).
1144
1145	IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1) THEN
1146	DO j = jts, jte
1147	DO k = kts, kte
1148	defor13(ids,k,j) = defor13(ids+1,k,j)
1149	defor23(ids,k,j) = defor23(ids+1,k,j)
1150	END DO
1151	END DO
1152	END IF
1153
1154	IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN
1155	DO k = kts, kte
1156	DO i = its, ite
1157	defor13(i,k,jds) = defor13(i,k,jds+1)
1158	defor23(i,k,jds) = defor23(i,k,jds+1)
1159	END DO
1160	END DO
1161	END IF
1162
1163	IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN
1164	DO j = jts, jte
1165	DO k = kts, kte
1166	defor13(ide,k,j) = defor13(ide-1,k,j)
1167	defor23(ide,k,j) = defor23(ide-1,k,j)
1168	END DO
1169	END DO
1170	END IF
1171
1172	IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN
1173	DO k = kts, kte
1174	DO i = its, ite
1175	defor13(i,k,jde) = defor13(i,k,jde-1)
1176	defor23(i,k,jde) = defor23(i,k,jde-1)
1177	END DO
1178	END DO
1179	END IF
1180
1181	ENDIF ! NBA-----------------------------------------------------------
1182
1183	! End update of boundary for defor13 and defor23.
1184	!-----------------------------------------------------------------------
1185
1186	! The second three (defor12, defor13, defor23) of six deformation terms
1187	! are now calculated at vorticity points.
1188	!=======================================================================
1189
1190	END SUBROUTINE cal_deform_and_div
1191
1192	!=======================================================================
1193	!=======================================================================
1194
1195	SUBROUTINE calculate_km_kh( config_flags, dt, &
1196	dampcoef, zdamp, damp_opt, &
1197	xkmh, xkmv, xkhh, xkhv, &
1198	BN2, khdif, kvdif, div, &
1199	defor11, defor22, defor33, &
1200	defor12, defor13, defor23, &
1201	tke, p8w, t8w, theta, t, p, moist, &
1202	dn, dnw, dx, dy, rdz, rdzw, isotropic, &
1203	n_moist, cf1, cf2, cf3, warm_rain, &
1204	mix_upper_bound, &
1205	msftx, msfty, &
1206	ids, ide, jds, jde, kds, kde, &
1207	ims, ime, jms, jme, kms, kme, &
1208	its, ite, jts, jte, kts, kte )
1209
1210	! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
1211	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
1212	! ... ...
1213
1214	! Purpose: This routine calculates exchange coefficients for the TKE
1215	! scheme.
1216
1217	! References: Klemp and Wilhelmson (JAS 1978)
1218	! Deardorff (B-L Meteor 1980)
1219	! Chen and Dudhia (NCAR WRF physics report 2000)
1220
1221	!-----------------------------------------------------------------------
1222	! Begin declarations.
1223
1224	IMPLICIT NONE
1225
1226	TYPE( grid_config_rec_type ), INTENT( IN ) &
1227	:: config_flags
1228
1229	INTEGER, INTENT( IN ) &
1230	:: n_moist, damp_opt, isotropic, &
1231	ids, ide, jds, jde, kds, kde, &
1232	ims, ime, jms, jme, kms, kme, &
1233	its, ite, jts, jte, kts, kte
1234
1235	LOGICAL, INTENT( IN ) &
1236	:: warm_rain
1237
1238	REAL, INTENT( IN ) &
1239	:: dx, dy, zdamp, dt, dampcoef, cf1, cf2, cf3, khdif, kvdif
1240
1241	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
1242	:: dnw, dn
1243
1244	REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist ), INTENT( INOUT ) &
1245	:: moist
1246
1247	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
1248	:: xkmv, xkmh, xkhv, xkhh, BN2
1249
1250	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT( IN ) &
1251	:: defor11, defor22, defor33, defor12, defor13, defor23, &
1252	div, rdz, rdzw, p8w, t8w, theta, t, p
1253
1254	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
1255	:: tke
1256
1257	REAL, INTENT( IN ) &
1258	:: mix_upper_bound
1259
1260	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
1261	:: msftx, msfty
1262
1263	! Local variables.
1264
1265	INTEGER &
1266	:: i_start, i_end, j_start, j_end, ktf, i, j, k
1267
1268	! End declarations.
1269	!-----------------------------------------------------------------------
1270
1271	ktf = MIN( kte, kde-1 )
1272	i_start = its
1273	i_end = MIN( ite, ide-1 )
1274	j_start = jts
1275	j_end = MIN( jte, jde-1 )
1276
1277	CALL calculate_N2( config_flags, BN2, moist, &
1278	theta, t, p, p8w, t8w, &
1279	dnw, dn, rdz, rdzw, &
1280	n_moist, cf1, cf2, cf3, warm_rain, &
1281	ids, ide, jds, jde, kds, kde, &
1282	ims, ime, jms, jme, kms, kme, &
1283	its, ite, jts, jte, kts, kte )
1284
1285	! Select a scheme for calculating diffusion coefficients.
1286
1287	km_coef: SELECT CASE( config_flags%km_opt )
1288
1289	CASE (1)
1290	CALL isotropic_km( config_flags, xkmh, xkmv, &
1291	xkhh, xkhv, khdif, kvdif, &
1292	ids, ide, jds, jde, kds, kde, &
1293	ims, ime, jms, jme, kms, kme, &
1294	its, ite, jts, jte, kts, kte )
1295	CASE (2)
1296	CALL tke_km( config_flags, xkmh, xkmv, &
1297	xkhh, xkhv, BN2, tke, p8w, t8w, theta, &
1298	rdz, rdzw, dx, dy, dt, isotropic, &
1299	mix_upper_bound, msftx, msfty, &
1300	ids, ide, jds, jde, kds, kde, &
1301	ims, ime, jms, jme, kms, kme, &
1302	its, ite, jts, jte, kts, kte )
1303	CASE (3)
1304	CALL smag_km( config_flags, xkmh, xkmv, &
1305	xkhh, xkhv, BN2, div, &
1306	defor11, defor22, defor33, &
1307	defor12, defor13, defor23, &
1308	rdzw, dx, dy, dt, isotropic, &
1309	mix_upper_bound, msftx, msfty, &
1310	ids, ide, jds, jde, kds, kde, &
1311	ims, ime, jms, jme, kms, kme, &
1312	its, ite, jts, jte, kts, kte )
1313	CASE (4)
1314	CALL smag2d_km( config_flags, xkmh, xkmv, &
1315	xkhh, xkhv, defor11, defor22, defor12, &
1316	rdzw, dx, dy, msftx, msfty, &
1317	ids, ide, jds, jde, kds, kde, &
1318	ims, ime, jms, jme, kms, kme, &
1319	its, ite, jts, jte, kts, kte )
1320	CASE DEFAULT
1321	CALL wrf_error_fatal( 'Please choose diffusion coefficient scheme' )
1322
1323	END SELECT km_coef
1324
1325	IF ( damp_opt .eq. 1 ) THEN
1326	CALL cal_dampkm( config_flags, xkmh, xkhh, xkmv, xkhv, &
1327	dx, dy, dt, dampcoef, rdz, rdzw, zdamp, &
1328	msftx, msfty, &
1329	ids, ide, jds, jde, kds, kde, &
1330	ims, ime, jms, jme, kms, kme, &
1331	its, ite, jts, jte, kts, kte )
1332	END IF
1333
1334	END SUBROUTINE calculate_km_kh
1335
1336	!=======================================================================
1337
1338	SUBROUTINE cal_dampkm( config_flags,xkmh,xkhh,xkmv,xkhv, &
1339	dx,dy,dt,dampcoef, &
1340	rdz, rdzw ,zdamp, &
1341	msftx, msfty, &
1342	ids,ide, jds,jde, kds,kde, &
1343	ims,ime, jms,jme, kms,kme, &
1344	its,ite, jts,jte, kts,kte )
1345
1346	!-----------------------------------------------------------------------
1347	! Begin declarations.
1348
1349	IMPLICIT NONE
1350
1351	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
1352
1353	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1354	ims, ime, jms, jme, kms, kme, &
1355	its, ite, jts, jte, kts, kte
1356
1357	REAL , INTENT(IN ) :: zdamp,dx,dy,dt,dampcoef
1358
1359
1360	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: xkmh , &
1361	xkhh , &
1362	xkmv , &
1363	xkhv
1364
1365	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdz, &
1366	rdzw
1367
1368	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
1369	msfty
1370	! LOCAL VARS
1371
1372	INTEGER :: i_start, i_end, j_start, j_end, ktf, ktfm1, i, j, k
1373	REAL :: kmmax,kmmvmax,degrad90,dz,tmp
1374	REAL :: ds
1375	REAL , DIMENSION( its:ite ) :: deltaz
1376	REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: dampk,dampkv
1377
1378	! End declarations.
1379	!-----------------------------------------------------------------------
1380
1381	ktf = min(kte,kde-1)
1382	ktfm1 = ktf-1
1383
1384	i_start = its
1385	i_end = MIN(ite,ide-1)
1386	j_start = jts
1387	j_end = MIN(jte,jde-1)
1388
1389	! keep upper damping diffusion away from relaxation zones at boundaries if used
1390	IF(config_flags%specified .OR. config_flags%nested)THEN
1391	i_start = MAX(i_start,ids+config_flags%spec_bdy_width-1)
1392	i_end = MIN(i_end,ide-config_flags%spec_bdy_width)
1393	j_start = MAX(j_start,jds+config_flags%spec_bdy_width-1)
1394	j_end = MIN(j_end,jde-config_flags%spec_bdy_width)
1395	ENDIF
1396
1397	kmmax=dx*dx/dt
1398	degrad90=DEGRAD*90.
1399	DO j = j_start, j_end
1400
1401	k=ktf
1402	DO i = i_start, i_end
1403	! Unmodified dx used above may produce very large diffusivities
1404	! when msftx is very large. And the above formula ignores the fact
1405	! that dy may now be different from dx as well. Let's fix that by
1406	! defining a "ds" as the minimum of the "real-space" (physical
1407	! distance) values of dx and dy, and then using that smallest value
1408	! to calculate a point-by-point kmmax
1409	ds = MIN(dx/msftx(i,j),dy/msfty(i,j))
1410	kmmax=ds*ds/dt
1411
1412	! deltaz(i)=0.5*dnw(k)/zeta_z(i,j)
1413	! dz=dnw(k)/zeta_z(i,j)
1414	dz = 1./rdzw(i,k,j)
1415	deltaz(i) = 0.5*dz
1416
1417	kmmvmax=dz*dz/dt
1418	tmp=min(deltaz(i)/zdamp,1.)
1419	dampk(i,k,j)=cos(degrad90tmp)cos(degrad90tmp)kmmax*dampcoef
1420	dampkv(i,k,j)=cos(degrad90tmp)cos(degrad90tmp)kmmvmax*dampcoef
1421	! set upper limit on vertical K (based on horizontal K)
1422	dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j))
1423
1424	ENDDO
1425
1426	DO k = ktfm1,kts,-1
1427	DO i = i_start, i_end
1428	! Unmodified dx used above may produce very large diffusivities
1429	! when msftx is very large. And the above formula ignores the fact
1430	! that dy may now be different from dx as well. Let's fix that by
1431	! defining a "ds" as the minimum of the "real-space" (physical
1432	! distance) values of dx and dy, and then using that smallest value
1433	! to calculate a point-by-point kmmax
1434	ds = MIN(dx/msftx(i,j),dy/msfty(i,j))
1435	kmmax=ds*ds/dt
1436
1437	! deltaz(i)=deltaz(i)+dn(k)/zeta_z(i,j)
1438	! dz=dnw(k)/zeta_z(i,j)
1439	dz = 1./rdz(i,k,j)
1440	deltaz(i) = deltaz(i) + dz
1441	dz = 1./rdzw(i,k,j)
1442
1443	kmmvmax=dz*dz/dt
1444	tmp=min(deltaz(i)/zdamp,1.)
1445	dampk(i,k,j)=cos(degrad90tmp)cos(degrad90tmp)kmmax*dampcoef
1446	dampkv(i,k,j)=cos(degrad90tmp)cos(degrad90tmp)kmmvmax*dampcoef
1447	! set upper limit on vertical K (based on horizontal K)
1448	dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j))
1449	ENDDO
1450	ENDDO
1451
1452	ENDDO
1453
1454	DO j = j_start, j_end
1455	DO k = kts,ktf
1456	DO i = i_start, i_end
1457	xkmh(i,k,j)=max(xkmh(i,k,j),dampk(i,k,j))
1458	xkhh(i,k,j)=max(xkhh(i,k,j),dampk(i,k,j))
1459	xkmv(i,k,j)=max(xkmv(i,k,j),dampkv(i,k,j))
1460	xkhv(i,k,j)=max(xkhv(i,k,j),dampkv(i,k,j))
1461	ENDDO
1462	ENDDO
1463	ENDDO
1464
1465	END SUBROUTINE cal_dampkm
1466
1467	!=======================================================================
1468	!=======================================================================
1469
1470	SUBROUTINE calculate_N2( config_flags, BN2, moist, &
1471	theta, t, p, p8w, t8w, &
1472	dnw, dn, rdz, rdzw, &
1473	n_moist, cf1, cf2, cf3, warm_rain, &
1474	ids, ide, jds, jde, kds, kde, &
1475	ims, ime, jms, jme, kms, kme, &
1476	its, ite, jts, jte, kts, kte )
1477
1478	!-----------------------------------------------------------------------
1479	! Begin declarations.
1480
1481	IMPLICIT NONE
1482
1483	TYPE( grid_config_rec_type ), INTENT( IN ) &
1484	:: config_flags
1485
1486	INTEGER, INTENT( IN ) &
1487	:: n_moist, &
1488	ids, ide, jds, jde, kds, kde, &
1489	ims, ime, jms, jme, kms, kme, &
1490	its, ite, jts, jte, kts, kte
1491
1492	LOGICAL, INTENT( IN ) &
1493	:: warm_rain
1494
1495	REAL, INTENT( IN ) &
1496	:: cf1, cf2, cf3
1497
1498	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
1499	:: BN2
1500
1501	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
1502	:: rdz, rdzw, theta, t, p, p8w, t8w
1503
1504	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
1505	:: dnw, dn
1506
1507	REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), INTENT( INOUT ) &
1508	:: moist
1509
1510	! Local variables.
1511
1512	INTEGER &
1513	:: i, j, k, ktf, ispe, ktes1, ktes2, &
1514	i_start, i_end, j_start, j_end
1515
1516	REAL &
1517	:: coefa, thetaep1, thetaem1, qc_cr, es, tc, qlpqi, qsw, qsi, &
1518	tmpdz, xlvqv, thetaesfc, thetasfc, qvtop, qvsfc, thetatop, thetaetop
1519
1520	REAL, DIMENSION( its:ite, jts:jte ) &
1521	:: tmp1sfc, tmp1top
1522
1523	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
1524	:: tmp1, qvs, qctmp
1525
1526	! End declarations.
1527	!-----------------------------------------------------------------------
1528
1529	qc_cr = 0.00001 ! in Kg/Kg
1530
1531	ktf = MIN( kte, kde-1 )
1532	ktes1 = kte-1
1533	ktes2 = kte-2
1534
1535	i_start = its
1536	i_end = MIN( ite, ide-1 )
1537	j_start = jts
1538	j_end = MIN( jte, jde-1 )
1539
1540	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
1541	config_flags%nested) i_start = MAX( ids+1, its )
1542	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
1543	config_flags%nested) i_end = MIN( ide-2, ite )
1544	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
1545	config_flags%nested) j_start = MAX( jds+1, jts )
1546	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
1547	config_flags%nested) j_end = MIN( jde-2 ,jte )
1548	IF ( config_flags%periodic_x ) i_start = its
1549	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
1550
1551	IF ( P_QC .GT. PARAM_FIRST_SCALAR) THEN
1552	DO j = j_start, j_end
1553	DO k = kts, ktf
1554	DO i = i_start, i_end
1555	qctmp(i,k,j) = moist(i,k,j,P_QC)
1556	END DO
1557	END DO
1558	END DO
1559	ELSE
1560	DO j = j_start, j_end
1561	DO k = kts, ktf
1562	DO i = i_start, i_end
1563	qctmp(i,k,j) = 0.0
1564	END DO
1565	END DO
1566	END DO
1567	END IF
1568
1569	DO j = jts, jte
1570	DO k = kts, kte
1571	DO i = its, ite
1572	tmp1(i,k,j) = 0.0
1573	END DO
1574	END DO
1575	END DO
1576
1577	DO j = jts,jte
1578	DO i = its,ite
1579	tmp1sfc(i,j) = 0.0
1580	tmp1top(i,j) = 0.0
1581	END DO
1582	END DO
1583
1584	DO ispe = PARAM_FIRST_SCALAR, n_moist
1585	IF ( ispe .EQ. P_QV .OR. ispe .EQ. P_QC .OR. ispe .EQ. P_QI) THEN
1586	DO j = j_start, j_end
1587	DO k = kts, ktf
1588	DO i = i_start, i_end
1589	tmp1(i,k,j) = tmp1(i,k,j) + moist(i,k,j,ispe)
1590	END DO
1591	END DO
1592	END DO
1593
1594	DO j = j_start, j_end
1595	DO i = i_start, i_end
1596	tmp1sfc(i,j) = tmp1sfc(i,j) + &
1597	cf1 * moist(i,1,j,ispe) + &
1598	cf2 * moist(i,2,j,ispe) + &
1599	cf3 * moist(i,3,j,ispe)
1600	tmp1top(i,j) = tmp1top(i,j) + &
1601	moist(i,ktes1,j,ispe) + &
1602	( moist(i,ktes1,j,ispe) - moist(i,ktes2,j,ispe) ) * &
1603	0.5 * dnw(ktes1) / dn(ktes1)
1604	END DO
1605	END DO
1606	END IF
1607	END DO
1608
1609	! Calculate saturation mixing ratio.
1610
1611	DO j = j_start, j_end
1612	DO k = kts, ktf
1613	DO i = i_start, i_end
1614	tc = t(i,k,j) - SVPT0
1615	es = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t(i,k,j) - SVP3 ) )
1616	qvs(i,k,j) = EP_2 * es / ( p(i,k,j) - es )
1617	END DO
1618	END DO
1619	END DO
1620
1621	DO j = j_start, j_end
1622	DO k = kts+1, ktf-1
1623	DO i = i_start, i_end
1624	tmpdz = 1.0 / rdz(i,k,j) + 1.0 / rdz(i,k+1,j)
1625	IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN
1626	xlvqv = XLV * moist(i,k,j,P_QV)
1627	coefa = ( 1.0 + xlvqv / R_d / t(i,k,j) ) / &
1628	( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) / &
1629	theta(i,k,j)
1630	thetaep1 = theta(i,k+1,j) * &
1631	( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) )
1632	thetaem1 = theta(i,k-1,j) * &
1633	( 1.0 + XLV * qvs(i,k-1,j) / Cp / t(i,k-1,j) )
1634	BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaem1 ) / tmpdz - &
1635	( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz )
1636	ELSE
1637	BN2(i,k,j) = g * ( (theta(i,k+1,j) - theta(i,k-1,j) ) / &
1638	theta(i,k,j) / tmpdz + &
1639	1.61 * ( moist(i,k+1,j,P_QV) - moist(i,k-1,j,P_QV) ) / &
1640	tmpdz - &
1641	( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz )
1642	ENDIF
1643	END DO
1644	END DO
1645	END DO
1646
1647	k = kts
1648	DO j = j_start, j_end
1649	DO i = i_start, i_end
1650	tmpdz = 1.0 / rdz(i,k+1,j) + 0.5 / rdzw(i,k,j)
1651	thetasfc = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp )
1652	IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN
1653	qvsfc = cf1 * qvs(i,1,j) + &
1654	cf2 * qvs(i,2,j) + &
1655	cf3 * qvs(i,3,j)
1656	xlvqv = XLV * moist(i,k,j,P_QV)
1657	coefa = ( 1.0 + xlvqv / R_d / t(i,k,j) ) / &
1658	( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) / &
1659	theta(i,k,j)
1660	thetaep1 = theta(i,k+1,j) * &
1661	( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) )
1662	thetaesfc = thetasfc * &
1663	( 1.0 + XLV * qvsfc / Cp / t8w(i,kts,j) )
1664	BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaesfc ) / tmpdz - &
1665	( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
1666	ELSE
1667	qvsfc = cf1 * moist(i,1,j,P_QV) + &
1668	cf2 * moist(i,2,j,P_QV) + &
1669	cf3 * moist(i,3,j,P_QV)
1670	! BN2(i,k,j) = g * ( ( theta(i,k+1,j) - thetasfc ) / &
1671	! theta(i,k,j) / tmpdz + &
1672	! 1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) / &
1673	! tmpdz - &
1674	! ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
1675	!...... MARTA: change in computation of BN2 at the surface, WCS 040331
1676
1677	tmpdz= 1./rdzw(i,k,j) ! controlare come calcola rdzw
1678	BN2(i,k,j) = g * ( ( theta(i,k+1,j) - theta(i,k,j)) / &
1679	theta(i,k,j) / tmpdz + &
1680	1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) / &
1681	tmpdz - &
1682	( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz )
1683	! end of MARTA/WCS change
1684
1685	ENDIF
1686	END DO
1687	END DO
1688
1689
1690	!...... MARTA: change in computation of BN2 at the top, WCS 040331
1691	DO j = j_start, j_end
1692	DO i = i_start, i_end
1693	BN2(i,ktf,j)=BN2(i,ktf-1,j)
1694	END DO
1695	END DO
1696	! end of MARTA/WCS change
1697
1698	END SUBROUTINE calculate_N2
1699
1700	!=======================================================================
1701	!=======================================================================
1702
1703	SUBROUTINE isotropic_km( config_flags, &
1704	xkmh,xkmv,xkhh,xkhv,khdif,kvdif, &
1705	ids,ide, jds,jde, kds,kde, &
1706	ims,ime, jms,jme, kms,kme, &
1707	its,ite, jts,jte, kts,kte )
1708
1709	!-----------------------------------------------------------------------
1710	! Begin declarations.
1711
1712	IMPLICIT NONE
1713
1714	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
1715
1716	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1717	ims, ime, jms, jme, kms, kme, &
1718	its, ite, jts, jte, kts, kte
1719
1720	REAL , INTENT(IN ) :: khdif,kvdif
1721
1722	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
1723	xkmv, &
1724	xkhh, &
1725	xkhv
1726	! LOCAL VARS
1727
1728	INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
1729	REAL :: khdif3,kvdif3
1730
1731	! End declarations.
1732	!-----------------------------------------------------------------------
1733
1734	ktf = kte
1735
1736	i_start = its
1737	i_end = MIN(ite,ide-1)
1738	j_start = jts
1739	j_end = MIN(jte,jde-1)
1740
1741	! khdif3=khdif*3.
1742	! kvdif3=kvdif*3.
1743	khdif3=khdif/prandtl
1744	kvdif3=kvdif/prandtl
1745
1746	DO j = j_start, j_end
1747	DO k = kts, ktf
1748	DO i = i_start, i_end
1749	xkmh(i,k,j)=khdif
1750	xkmv(i,k,j)=kvdif
1751	xkhh(i,k,j)=khdif3
1752	xkhv(i,k,j)=kvdif3
1753	ENDDO
1754	ENDDO
1755	ENDDO
1756
1757	END SUBROUTINE isotropic_km
1758
1759	!=======================================================================
1760	!=======================================================================
1761
1762	SUBROUTINE smag_km( config_flags,xkmh,xkmv,xkhh,xkhv,BN2, &
1763	div,defor11,defor22,defor33,defor12, &
1764	defor13,defor23, &
1765	rdzw,dx,dy,dt,isotropic, &
1766	mix_upper_bound, msftx, msfty, &
1767	ids,ide, jds,jde, kds,kde, &
1768	ims,ime, jms,jme, kms,kme, &
1769	its,ite, jts,jte, kts,kte )
1770
1771	!-----------------------------------------------------------------------
1772	! Begin declarations.
1773
1774	IMPLICIT NONE
1775
1776	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
1777
1778	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1779	ims, ime, jms, jme, kms, kme, &
1780	its, ite, jts, jte, kts, kte
1781
1782	INTEGER , INTENT(IN ) :: isotropic
1783	REAL , INTENT(IN ) :: dx, dy, dt, mix_upper_bound
1784
1785
1786	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: BN2, &
1787	rdzw
1788
1789	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
1790	xkmv, &
1791	xkhh, &
1792	xkhv
1793
1794	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: &
1795	defor11, &
1796	defor22, &
1797	defor33, &
1798	defor12, &
1799	defor13, &
1800	defor23, &
1801	div
1802	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
1803	msfty
1804	! LOCAL VARS
1805
1806	INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
1807	REAL :: deltas, tmp, pr, mlen_h, mlen_v, c_s
1808
1809	REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2
1810
1811	! End declarations.
1812	!-----------------------------------------------------------------------
1813
1814	ktf = min(kte,kde-1)
1815
1816	i_start = its
1817	i_end = MIN(ite,ide-1)
1818	j_start = jts
1819	j_end = MIN(jte,jde-1)
1820
1821	IF ( config_flags%open_xs .or. config_flags%specified .or. &
1822	config_flags%nested) i_start = MAX(ids+1,its)
1823	IF ( config_flags%open_xe .or. config_flags%specified .or. &
1824	config_flags%nested) i_end = MIN(ide-2,ite)
1825	IF ( config_flags%open_ys .or. config_flags%specified .or. &
1826	config_flags%nested) j_start = MAX(jds+1,jts)
1827	IF ( config_flags%open_ye .or. config_flags%specified .or. &
1828	config_flags%nested) j_end = MIN(jde-2,jte)
1829	IF ( config_flags%periodic_x ) i_start = its
1830	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
1831
1832	pr = prandtl
1833	c_s = config_flags%c_s
1834
1835	do j=j_start,j_end
1836	do k=kts,ktf
1837	do i=i_start,i_end
1838	def2(i,k,j)=0.5(defor11(i,k,j)defor11(i,k,j) + &
1839	defor22(i,k,j)*defor22(i,k,j) + &
1840	defor33(i,k,j)*defor33(i,k,j))
1841	enddo
1842	enddo
1843	enddo
1844
1845	do j=j_start,j_end
1846	do k=kts,ktf
1847	do i=i_start,i_end
1848	tmp=0.25*(defor12(i ,k,j)+defor12(i ,k,j+1)+ &
1849	defor12(i+1,k,j)+defor12(i+1,k,j+1))
1850	def2(i,k,j)=def2(i,k,j)+tmp*tmp
1851	enddo
1852	enddo
1853	enddo
1854
1855	do j=j_start,j_end
1856	do k=kts,ktf
1857	do i=i_start,i_end
1858	tmp=0.25*(defor13(i ,k+1,j)+defor13(i ,k,j)+ &
1859	defor13(i+1,k+1,j)+defor13(i+1,k,j))
1860	def2(i,k,j)=def2(i,k,j)+tmp*tmp
1861	enddo
1862	enddo
1863	enddo
1864
1865	do j=j_start,j_end
1866	do k=kts,ktf
1867	do i=i_start,i_end
1868	tmp=0.25*(defor23(i,k+1,j )+defor23(i,k,j )+ &
1869	defor23(i,k+1,j+1)+defor23(i,k,j+1))
1870	def2(i,k,j)=def2(i,k,j)+tmp*tmp
1871	enddo
1872	enddo
1873	enddo
1874	!
1875	IF (isotropic .EQ. 0) THEN
1876	DO j = j_start, j_end
1877	DO k = kts, ktf
1878	DO i = i_start, i_end
1879	mlen_h=sqrt(dx/msftx(i,j) * dy/msfty(i,j))
1880	mlen_v= 1./rdzw(i,k,j)
1881	tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
1882	tmp=tmp**0.5
1883	xkmh(i,k,j)=max(c_sc_smlen_hmlen_htmp, 1.0E-6mlen_hmlen_h )
1884	xkmh(i,k,j)=min(xkmh(i,k,j), mix_upper_bound * mlen_h * mlen_h / dt )
1885	xkmv(i,k,j)=max(c_sc_smlen_vmlen_vtmp, 1.0E-6mlen_vmlen_v )
1886	xkmv(i,k,j)=min(xkmv(i,k,j), mix_upper_bound * mlen_v * mlen_v / dt )
1887	xkhh(i,k,j)=xkmh(i,k,j)/pr
1888	xkhh(i,k,j)=min(xkhh(i,k,j), mix_upper_bound * mlen_h * mlen_h / dt )
1889	xkhv(i,k,j)=xkmv(i,k,j)/pr
1890	xkhv(i,k,j)=min(xkhv(i,k,j), mix_upper_bound * mlen_v * mlen_v / dt )
1891	ENDDO
1892	ENDDO
1893	ENDDO
1894	ELSE
1895	DO j = j_start, j_end
1896	DO k = kts, ktf
1897	DO i = i_start, i_end
1898	deltas=(dx/msftx(i,j) * dy/msfty(i,j)/rdzw(i,k,j))**0.33333333
1899	tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
1900	tmp=tmp**0.5
1901	xkmh(i,k,j)=max(c_sc_sdeltasdeltastmp, 1.0E-6deltasdeltas )
1902	xkmh(i,k,j)=min(xkmh(i,k,j), mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt )
1903	xkmv(i,k,j)=xkmh(i,k,j)
1904	xkmv(i,k,j)=min(xkmv(i,k,j), mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt )
1905	xkhh(i,k,j)=xkmh(i,k,j)/pr
1906	xkhh(i,k,j)=min(xkhh(i,k,j), mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt )
1907	xkhv(i,k,j)=xkmv(i,k,j)/pr
1908	xkhv(i,k,j)=min(xkhv(i,k,j), mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt )
1909	ENDDO
1910	ENDDO
1911	ENDDO
1912	ENDIF
1913
1914	END SUBROUTINE smag_km
1915
1916	!=======================================================================
1917	!=======================================================================
1918
1919	SUBROUTINE smag2d_km( config_flags,xkmh,xkmv,xkhh,xkhv, &
1920	defor11,defor22,defor12, &
1921	rdzw,dx,dy,msftx, msfty, &
1922	ids,ide, jds,jde, kds,kde, &
1923	ims,ime, jms,jme, kms,kme, &
1924	its,ite, jts,jte, kts,kte )
1925
1926	!-----------------------------------------------------------------------
1927	! Begin declarations.
1928
1929	IMPLICIT NONE
1930
1931	TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
1932
1933	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1934	ims, ime, jms, jme, kms, kme, &
1935	its, ite, jts, jte, kts, kte
1936
1937	REAL , INTENT(IN ) :: dx, dy
1938
1939
1940	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rdzw
1941
1942	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
1943	xkmv, &
1944	xkhh, &
1945	xkhv
1946
1947	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: &
1948	defor11, &
1949	defor22, &
1950	defor12
1951
1952	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
1953	msfty
1954	! LOCAL VARS
1955
1956	INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
1957	REAL :: deltas, tmp, pr, mlen_h, c_s
1958
1959	REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2
1960
1961	! End declarations.
1962	!-----------------------------------------------------------------------
1963
1964	ktf = min(kte,kde-1)
1965
1966	i_start = its
1967	i_end = MIN(ite,ide-1)
1968	j_start = jts
1969	j_end = MIN(jte,jde-1)
1970
1971	IF ( config_flags%open_xs .or. config_flags%specified .or. &
1972	config_flags%nested) i_start = MAX(ids+1,its)
1973	IF ( config_flags%open_xe .or. config_flags%specified .or. &
1974	config_flags%nested) i_end = MIN(ide-2,ite)
1975	IF ( config_flags%open_ys .or. config_flags%specified .or. &
1976	config_flags%nested) j_start = MAX(jds+1,jts)
1977	IF ( config_flags%open_ye .or. config_flags%specified .or. &
1978	config_flags%nested) j_end = MIN(jde-2,jte)
1979	IF ( config_flags%periodic_x ) i_start = its
1980	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
1981
1982	pr=prandtl
1983	c_s = config_flags%c_s
1984
1985	do j=j_start,j_end
1986	do k=kts,ktf
1987	do i=i_start,i_end
1988	def2(i,k,j)=0.25((defor11(i,k,j)-defor22(i,k,j))(defor11(i,k,j)-defor22(i,k,j)))
1989	tmp=0.25*(defor12(i ,k,j)+defor12(i ,k,j+1)+ &
1990	defor12(i+1,k,j)+defor12(i+1,k,j+1))
1991	def2(i,k,j)=def2(i,k,j)+tmp*tmp
1992	enddo
1993	enddo
1994	enddo
1995	!
1996	DO j = j_start, j_end
1997	DO k = kts, ktf
1998	DO i = i_start, i_end
1999	mlen_h=sqrt(dx/msftx(i,j) * dy/msfty(i,j))
2000	tmp=sqrt(def2(i,k,j))
2001	! xkmh(i,k,j)=max(c_sc_smlen_hmlen_htmp, 1.0E-6mlen_hmlen_h )
2002	xkmh(i,k,j)=c_sc_smlen_hmlen_htmp
2003	xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen_h )
2004	xkmv(i,k,j)=0.
2005	xkhh(i,k,j)=xkmh(i,k,j)/pr
2006	xkhv(i,k,j)=0.
2007	ENDDO
2008	ENDDO
2009	ENDDO
2010
2011	END SUBROUTINE smag2d_km
2012
2013	!=======================================================================
2014	!=======================================================================
2015
2016	SUBROUTINE tke_km( config_flags, xkmh, xkmv, xkhh, xkhv, &
2017	bn2, tke, p8w, t8w, theta, &
2018	rdz, rdzw, dx,dy, dt, isotropic, &
2019	mix_upper_bound, msftx, msfty, &
2020	ids, ide, jds, jde, kds, kde, &
2021	ims, ime, jms, jme, kms, kme, &
2022	its, ite, jts, jte, kts, kte )
2023
2024	! History: Sep 2003 Changes by Jason Knievel and George Bryan, NCAR
2025	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
2026	! ... ...
2027
2028	! Purpose: This routine calculates the exchange coefficients for the
2029	! TKE turbulence parameterization.
2030
2031	! References: Klemp and Wilhelmson (JAS 1978)
2032	! Chen and Dudhia (NCAR WRF physics report 2000)
2033
2034	!-----------------------------------------------------------------------
2035	! Begin declarations.
2036
2037	IMPLICIT NONE
2038
2039	TYPE( grid_config_rec_type ), INTENT( IN ) &
2040	:: config_flags
2041
2042	INTEGER, INTENT( IN ) &
2043	:: ids, ide, jds, jde, kds, kde, &
2044	ims, ime, jms, jme, kms, kme, &
2045	its, ite, jts, jte, kts, kte
2046
2047	INTEGER, INTENT( IN ) :: isotropic
2048	REAL, INTENT( IN ) &
2049	:: dx, dy, dt
2050
2051	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
2052	:: tke, p8w, t8w, theta, rdz, rdzw, bn2
2053
2054	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
2055	:: xkmh, xkmv, xkhh, xkhv
2056
2057	REAL, INTENT( IN ) &
2058	:: mix_upper_bound
2059
2060	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
2061	msfty
2062	! Local variables.
2063
2064	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
2065	:: l_scale
2066
2067	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
2068	:: dthrdn
2069
2070	REAL &
2071	:: deltas, tmp, mlen_s, mlen_h, mlen_v, tmpdz, &
2072	thetasfc, thetatop, minkx, pr_inv, pr_inv_h, pr_inv_v, c_k
2073
2074	INTEGER &
2075	:: i_start, i_end, j_start, j_end, ktf, i, j, k
2076
2077	REAL, PARAMETER :: tke_seed_value = 1.e-06
2078	REAL :: tke_seed
2079	REAL, PARAMETER :: epsilon = 1.e-10
2080
2081	! End declarations.
2082	!-----------------------------------------------------------------------
2083
2084	ktf = MIN( kte, kde-1 )
2085	i_start = its
2086	i_end = MIN( ite, ide-1 )
2087	j_start = jts
2088	j_end = MIN( jte, jde-1 )
2089
2090	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
2091	config_flags%nested) i_start = MAX( ids+1, its )
2092	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
2093	config_flags%nested) i_end = MIN( ide-2, ite )
2094	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
2095	config_flags%nested) j_start = MAX( jds+1, jts )
2096	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
2097	config_flags%nested) j_end = MIN( jde-2, jte)
2098	IF ( config_flags%periodic_x ) i_start = its
2099	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
2100
2101	! in the absence of surface drag or a surface heat flux, there
2102	! is no way to generate tke without pre-existing tke. Use
2103	! tke_seed if the drag and flux are off.
2104
2105	c_k = config_flags%c_k
2106	tke_seed = tke_seed_value
2107	if( (config_flags%tke_drag_coefficient .gt. epsilon) .or. &
2108	(config_flags%tke_heat_flux .gt. epsilon) ) tke_seed = 0.
2109
2110	DO j = j_start, j_end
2111	DO k = kts+1, ktf-1
2112	DO i = i_start, i_end
2113	tmpdz = 1.0 / ( rdz(i,k+1,j) + rdz(i,k,j) )
2114	dthrdn(i,k,j) = ( theta(i,k+1,j) - theta(i,k-1,j) ) / tmpdz
2115	END DO
2116	END DO
2117	END DO
2118
2119	k = kts
2120	DO j = j_start, j_end
2121	DO i = i_start, i_end
2122	tmpdz = 1.0 / ( rdzw(i,k+1,j) + rdzw(i,k,j) )
2123	thetasfc = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp )
2124	dthrdn(i,k,j) = ( theta(i,k+1,j) - thetasfc ) / tmpdz
2125	END DO
2126	END DO
2127
2128	k = ktf
2129	DO j = j_start, j_end
2130	DO i = i_start, i_end
2131	tmpdz = 1.0 / rdz(i,k,j) + 0.5 / rdzw(i,k,j)
2132	thetatop = T8w(i,kde,j) / ( p8w(i,kde,j) / p1000mb )**( R_d / Cp )
2133	dthrdn(i,k,j) = ( thetatop - theta(i,k-1,j) ) / tmpdz
2134	END DO
2135	END DO
2136
2137	IF ( isotropic .EQ. 0 ) THEN
2138	DO j = j_start, j_end
2139	DO k = kts, ktf
2140	DO i = i_start, i_end
2141	mlen_h = SQRT( dx/msftx(i,j) * dy/msfty(i,j) )
2142	tmp = SQRT( MAX( tke(i,k,j), tke_seed ) )
2143	deltas = 1.0 / rdzw(i,k,j)
2144	mlen_v = deltas
2145	IF ( dthrdn(i,k,j) .GT. 0.) THEN
2146	mlen_s = 0.76 * tmp / ( ABS( g / theta(i,k,j) * dthrdn(i,k,j) ) )**0.5
2147	mlen_v = MIN( mlen_v, mlen_s )
2148	END IF
2149	xkmh(i,k,j) = MAX( c_k * tmp * mlen_h, 1.0E-6 * mlen_h * mlen_h )
2150	xkmh(i,k,j) = MIN( xkmh(i,k,j), mix_upper_bound * mlen_h *mlen_h / dt )
2151	xkmv(i,k,j) = MAX( c_k * tmp * mlen_v, 1.0E-6 * deltas * deltas )
2152	xkmv(i,k,j) = MIN( xkmv(i,k,j), mix_upper_bound * deltas *deltas / dt )
2153	pr_inv_h = 1./prandtl
2154	pr_inv_v = 1.0 + 2.0 * mlen_v / deltas
2155	xkhh(i,k,j) = xkmh(i,k,j) * pr_inv_h
2156	xkhv(i,k,j) = xkmv(i,k,j) * pr_inv_v
2157	END DO
2158	END DO
2159	END DO
2160	ELSE
2161	CALL calc_l_scale( config_flags, tke, BN2, l_scale, &
2162	i_start, i_end, ktf, j_start, j_end, &
2163	dx, dy, rdzw, msftx, msfty, &
2164	ids, ide, jds, jde, kds, kde, &
2165	ims, ime, jms, jme, kms, kme, &
2166	its, ite, jts, jte, kts, kte )
2167	DO j = j_start, j_end
2168	DO k = kts, ktf
2169	DO i = i_start, i_end
2170	tmp = SQRT( MAX( tke(i,k,j), tke_seed ) )
2171	deltas = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333
2172	xkmh(i,k,j) = c_k * tmp * l_scale(i,k,j)
2173	xkmh(i,k,j) = MIN( mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt, xkmh(i,k,j) )
2174	xkmv(i,k,j) = c_k * tmp * l_scale(i,k,j)
2175	xkmv(i,k,j) = MIN( mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt , xkmv(i,k,j) )
2176	pr_inv = 1.0 + 2.0 * l_scale(i,k,j) / deltas
2177	xkhh(i,k,j) = MIN( mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt, xkmh(i,k,j) * pr_inv )
2178	xkhv(i,k,j) = MIN( mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt, xkmv(i,k,j) * pr_inv )
2179	END DO
2180	END DO
2181	END DO
2182	END IF
2183
2184	END SUBROUTINE tke_km
2185
2186	!=======================================================================
2187	!=======================================================================
2188
2189	SUBROUTINE calc_l_scale( config_flags, tke, BN2, l_scale, &
2190	i_start, i_end, ktf, j_start, j_end, &
2191	dx, dy, rdzw, msftx, msfty, &
2192	ids, ide, jds, jde, kds, kde, &
2193	ims, ime, jms, jme, kms, kme, &
2194	its, ite, jts, jte, kts, kte )
2195
2196	! History: Sep 2003 Written by Bryan and Knievel, NCAR
2197
2198	! Purpose: This routine calculates the length scale, based on stability,
2199	! for TKE parameterization of subgrid-scale turbulence.
2200
2201	!-----------------------------------------------------------------------
2202	! Begin declarations.
2203
2204	IMPLICIT NONE
2205
2206	TYPE( grid_config_rec_type ), INTENT( IN ) &
2207	:: config_flags
2208
2209	INTEGER, INTENT( IN ) &
2210	:: i_start, i_end, ktf, j_start, j_end, &
2211	ids, ide, jds, jde, kds, kde, &
2212	ims, ime, jms, jme, kms, kme, &
2213	its, ite, jts, jte, kts, kte
2214
2215	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
2216	:: BN2, tke, rdzw
2217
2218	REAL, INTENT( IN ) &
2219	:: dx, dy
2220
2221	REAL, DIMENSION( its:ite, kts:kte, jts:jte ), INTENT( OUT ) &
2222	:: l_scale
2223
2224	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, &
2225	msfty
2226	! Local variables.
2227
2228	INTEGER &
2229	:: i, j, k
2230
2231	REAL &
2232	:: deltas, tmp
2233
2234	! End declarations.
2235	!-----------------------------------------------------------------------
2236
2237	DO j = j_start, j_end
2238	DO k = kts, ktf
2239	DO i = i_start, i_end
2240	deltas = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333
2241	l_scale(i,k,j) = deltas
2242
2243	IF ( BN2(i,k,j) .gt. 1.0e-6 ) THEN
2244	tmp = SQRT( MAX( tke(i,k,j), 1.0e-6 ) )
2245	l_scale(i,k,j) = 0.76 * tmp / SQRT( BN2(i,k,j) )
2246	l_scale(i,k,j) = MIN( l_scale(i,k,j), deltas)
2247	l_scale(i,k,j) = MAX( l_scale(i,k,j), 0.001 * deltas )
2248	END IF
2249
2250	END DO
2251	END DO
2252	END DO
2253
2254	END SUBROUTINE calc_l_scale
2255
2256	!=======================================================================
2257	!=======================================================================
2258
2259	SUBROUTINE horizontal_diffusion_2 ( rt_tendf, ru_tendf, rv_tendf, rw_tendf, &
2260	tke_tendf, &
2261	moist_tendf, n_moist, &
2262	chem_tendf, n_chem, &
2263	scalar_tendf, n_scalar, &
2264	tracer_tendf, n_tracer, &
2265	thp, theta, mu, tke, config_flags, &
2266	defor11, defor22, defor12, &
2267	defor13, defor23, &
2268	nba_mij, n_nba_mij, & !JDM
2269	div, &
2270	moist, chem, scalar,tracer, &
2271	msfux, msfuy, msfvx, msfvy, &
2272	msftx, msfty, xkmh, xkhh,km_opt, &
2273	rdx, rdy, rdz, rdzw, fnm, fnp, &
2274	cf1, cf2, cf3, zx, zy, dn, dnw, &
2275	ids, ide, jds, jde, kds, kde, &
2276	ims, ime, jms, jme, kms, kme, &
2277	its, ite, jts, jte, kts, kte )
2278
2279	!-----------------------------------------------------------------------
2280	! Begin declarations.
2281
2282	IMPLICIT NONE
2283
2284	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2285
2286	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2287	ims, ime, jms, jme, kms, kme, &
2288	its, ite, jts, jte, kts, kte
2289
2290	INTEGER , INTENT(IN ) :: n_moist,n_chem,n_scalar,n_tracer,km_opt
2291
2292	REAL , INTENT(IN ) :: cf1, cf2, cf3
2293
2294	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2295	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2296	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
2297	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
2298
2299	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux, &
2300	msfuy, &
2301	msfvx, &
2302	msfvy, &
2303	msftx, &
2304	msfty, &
2305	mu
2306
2307	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::rt_tendf,&
2308	ru_tendf,&
2309	rv_tendf,&
2310	rw_tendf,&
2311	tke_tendf
2312
2313	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
2314	INTENT(INOUT) :: moist_tendf
2315
2316	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
2317	INTENT(INOUT) :: chem_tendf
2318
2319	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar), &
2320	INTENT(INOUT) :: scalar_tendf
2321
2322	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_tracer), &
2323	INTENT(INOUT) :: tracer_tendf
2324
2325	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
2326	INTENT(IN ) :: moist
2327
2328	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
2329	INTENT(IN ) :: chem
2330
2331	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
2332	INTENT(IN ) :: scalar
2333
2334	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_tracer) , &
2335	INTENT(IN ) :: tracer
2336
2337	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
2338	defor22, &
2339	defor12, &
2340	defor13, &
2341	defor23, &
2342	div, &
2343	xkmh, &
2344	xkhh, &
2345	zx, &
2346	zy, &
2347	theta, &
2348	thp, &
2349	tke, &
2350	rdz, &
2351	rdzw
2352
2353
2354	REAL , INTENT(IN ) :: rdx, &
2355	rdy
2356	INTEGER, INTENT( IN ) :: n_nba_mij !JDM
2357
2358	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
2359	:: nba_mij
2360
2361	! LOCAL VARS
2362
2363	INTEGER :: im, ic, is
2364
2365	! REAL , DIMENSION(its-1:ite+1, kts:kte, jts-1:jte+1) :: xkhh
2366
2367	! End declarations.
2368	!-----------------------------------------------------------------------
2369
2370	!-----------------------------------------------------------------------
2371	! Call diffusion subroutines.
2372
2373	CALL horizontal_diffusion_u_2( ru_tendf, mu, config_flags, &
2374	defor11, defor12, div, &
2375	nba_mij, n_nba_mij, & !JDM
2376	tke(ims,kms,jms), &
2377	msfux, msfuy, xkmh, rdx, rdy, fnm, fnp, &
2378	zx, zy, rdzw, &
2379	ids, ide, jds, jde, kds, kde, &
2380	ims, ime, jms, jme, kms, kme, &
2381	its, ite, jts, jte, kts, kte )
2382
2383	CALL horizontal_diffusion_v_2( rv_tendf, mu, config_flags, &
2384	defor12, defor22, div, &
2385	nba_mij, n_nba_mij, & !JDM
2386	tke(ims,kms,jms), &
2387	msfvx, msfvy, xkmh, rdx, rdy, fnm, fnp, &
2388	zx, zy, rdzw, &
2389	ids, ide, jds, jde, kds, kde, &
2390	ims, ime, jms, jme, kms, kme, &
2391	its, ite, jts, jte, kts, kte )
2392
2393	CALL horizontal_diffusion_w_2( rw_tendf, mu, config_flags, &
2394	defor13, defor23, div, &
2395	nba_mij, n_nba_mij, & !JDM
2396	tke(ims,kms,jms), &
2397	msftx, msfty, xkmh, rdx, rdy, fnm, fnp, &
2398	zx, zy, rdz, &
2399	ids, ide, jds, jde, kds, kde, &
2400	ims, ime, jms, jme, kms, kme, &
2401	its, ite, jts, jte, kts, kte )
2402
2403	CALL horizontal_diffusion_s ( rt_tendf, mu, config_flags, thp, &
2404	msftx, msfty, msfux, msfuy, &
2405	msfvx, msfvy, xkhh, rdx, rdy, &
2406	fnm, fnp, cf1, cf2, cf3, &
2407	zx, zy, rdz, rdzw, dnw, dn, &
2408	.false., &
2409	ids, ide, jds, jde, kds, kde, &
2410	ims, ime, jms, jme, kms, kme, &
2411	its, ite, jts, jte, kts, kte )
2412
2413	IF (km_opt .eq. 2) &
2414	CALL horizontal_diffusion_s ( tke_tendf(ims,kms,jms), &
2415	mu, config_flags, &
2416	tke(ims,kms,jms), &
2417	msftx, msfty, msfux, msfuy, &
2418	msfvx, msfvy, xkhh, rdx, rdy, &
2419	fnm, fnp, cf1, cf2, cf3, &
2420	zx, zy, rdz, rdzw, dnw, dn, &
2421	.true., &
2422	ids, ide, jds, jde, kds, kde, &
2423	ims, ime, jms, jme, kms, kme, &
2424	its, ite, jts, jte, kts, kte )
2425
2426	IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN
2427
2428	moist_loop: do im = PARAM_FIRST_SCALAR, n_moist
2429
2430	CALL horizontal_diffusion_s( moist_tendf(ims,kms,jms,im), &
2431	mu, config_flags, &
2432	moist(ims,kms,jms,im), &
2433	msftx, msfty, msfux, msfuy, &
2434	msfvx, msfvy, xkhh, rdx, rdy, &
2435	fnm, fnp, cf1, cf2, cf3, &
2436	zx, zy, rdz, rdzw, dnw, dn, &
2437	.false., &
2438	ids, ide, jds, jde, kds, kde, &
2439	ims, ime, jms, jme, kms, kme, &
2440	its, ite, jts, jte, kts, kte )
2441
2442	ENDDO moist_loop
2443
2444	ENDIF
2445
2446	IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN
2447
2448	chem_loop: do ic = PARAM_FIRST_SCALAR, n_chem
2449
2450	CALL horizontal_diffusion_s( chem_tendf(ims,kms,jms,ic), &
2451	mu, config_flags, &
2452	chem(ims,kms,jms,ic), &
2453	msftx, msfty, msfux, msfuy, &
2454	msfvx, msfvy, xkhh, rdx, rdy, &
2455	fnm, fnp, cf1, cf2, cf3, &
2456	zx, zy, rdz, rdzw, dnw, dn, &
2457	.false., &
2458	ids, ide, jds, jde, kds, kde, &
2459	ims, ime, jms, jme, kms, kme, &
2460	its, ite, jts, jte, kts, kte )
2461
2462	ENDDO chem_loop
2463
2464	ENDIF
2465
2466	IF (n_tracer .ge. PARAM_FIRST_SCALAR) THEN
2467
2468	tracer_loop: do ic = PARAM_FIRST_SCALAR, n_tracer
2469
2470	CALL horizontal_diffusion_s( tracer_tendf(ims,kms,jms,ic), &
2471	mu, config_flags, &
2472	tracer(ims,kms,jms,ic), &
2473	msftx, msfty, msfux, msfuy, &
2474	msfvx, msfvy, xkhh, rdx, rdy, &
2475	fnm, fnp, cf1, cf2, cf3, &
2476	zx, zy, rdz, rdzw, dnw, dn, &
2477	.false., &
2478	ids, ide, jds, jde, kds, kde, &
2479	ims, ime, jms, jme, kms, kme, &
2480	its, ite, jts, jte, kts, kte )
2481
2482	ENDDO tracer_loop
2483
2484	ENDIF
2485	IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN
2486
2487	scalar_loop: do is = PARAM_FIRST_SCALAR, n_scalar
2488
2489	CALL horizontal_diffusion_s( scalar_tendf(ims,kms,jms,is), &
2490	mu, config_flags, &
2491	scalar(ims,kms,jms,is), &
2492	msftx, msfty, msfux, msfuy, &
2493	msfvx, msfvy, xkhh, rdx, rdy, &
2494	fnm, fnp, cf1, cf2, cf3, &
2495	zx, zy, rdz, rdzw, dnw, dn, &
2496	.false., &
2497	ids, ide, jds, jde, kds, kde, &
2498	ims, ime, jms, jme, kms, kme, &
2499	its, ite, jts, jte, kts, kte )
2500
2501	ENDDO scalar_loop
2502
2503	ENDIF
2504
2505	END SUBROUTINE horizontal_diffusion_2
2506
2507	!=======================================================================
2508	!=======================================================================
2509
2510	SUBROUTINE horizontal_diffusion_u_2( tendency, mu, config_flags, &
2511	defor11, defor12, div, &
2512	nba_mij, n_nba_mij, & !JDM
2513	tke, &
2514	msfux, msfuy, &
2515	xkmh, rdx, rdy, fnm, fnp, &
2516	zx, zy, rdzw, &
2517	ids, ide, jds, jde, kds, kde, &
2518	ims, ime, jms, jme, kms, kme, &
2519	its, ite, jts, jte, kts, kte )
2520
2521	!-----------------------------------------------------------------------
2522	! Begin declarations.
2523
2524	IMPLICIT NONE
2525
2526	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2527
2528	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2529	ims, ime, jms, jme, kms, kme, &
2530	its, ite, jts, jte, kts, kte
2531
2532	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2533	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2534
2535	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux, &
2536	msfuy, &
2537	mu
2538
2539	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
2540
2541	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdzw
2542
2543
2544	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
2545	defor12, &
2546	div, &
2547	tke, &
2548	xkmh, &
2549	zx, &
2550	zy
2551
2552	INTEGER, INTENT( IN ) :: n_nba_mij !JDM
2553
2554	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
2555	:: nba_mij
2556
2557	REAL , INTENT(IN ) :: rdx, &
2558	rdy
2559	! Local data
2560
2561	INTEGER :: i, j, k, ktf
2562
2563	INTEGER :: i_start, i_end, j_start, j_end
2564	INTEGER :: is_ext,ie_ext,js_ext,je_ext
2565
2566	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
2567	titau2avg, &
2568	titau1, &
2569	titau2, &
2570	xkxavg, &
2571	rravg
2572	! new
2573	! zxavg, &
2574	! zyavg
2575	REAL :: mrdx, mrdy, rcoup
2576
2577	REAL :: tmpzy, tmpzeta_z
2578
2579	REAL :: term1, term2, term3
2580
2581	! End declarations.
2582	!-----------------------------------------------------------------------
2583
2584	ktf=MIN(kte,kde-1)
2585
2586	!-----------------------------------------------------------------------
2587	! u : p (.), u(\|), w(-)
2588	!
2589	! p u p u u u
2590	!
2591	! p \| . \| . \| . \| k+1 \| . \| . \| . \| k+1
2592	!
2593	! w - 13 - - k+1 13 k+1
2594	!
2595	! p \| 11 O 11 \| . \| k \| 12 O 12 \| . \| k
2596	!
2597	! w - 13 - - k 13 k
2598	!
2599	! p \| . \| . \| . \| k-1 \| . \| . \| . \| k-1
2600	!
2601	! i-1 i i i+1 j-1 j j j+1 j+1
2602	!
2603
2604	i_start = its
2605	i_end = ite
2606	j_start = jts
2607	j_end = MIN(jte,jde-1)
2608
2609	IF ( config_flags%open_xs .or. config_flags%specified .or. &
2610	config_flags%nested) i_start = MAX(ids+1,its)
2611	IF ( config_flags%open_xe .or. config_flags%specified .or. &
2612	config_flags%nested) i_end = MIN(ide-1,ite)
2613	IF ( config_flags%open_ys .or. config_flags%specified .or. &
2614	config_flags%nested) j_start = MAX(jds+1,jts)
2615	IF ( config_flags%open_ye .or. config_flags%specified .or. &
2616	config_flags%nested) j_end = MIN(jde-2,jte)
2617	IF ( config_flags%periodic_x ) i_start = its
2618	IF ( config_flags%periodic_x ) i_end = ite
2619
2620	! titau1 = titau11
2621	is_ext=1
2622	ie_ext=0
2623	js_ext=0
2624	je_ext=0
2625	CALL cal_titau_11_22_33( config_flags, titau1, &
2626	mu, tke, xkmh, defor11, &
2627	nba_mij(ims,kms,jms,P_m11), & !JDM
2628	is_ext, ie_ext, js_ext, je_ext, &
2629	ids, ide, jds, jde, kds, kde, &
2630	ims, ime, jms, jme, kms, kme, &
2631	its, ite, jts, jte, kts, kte )
2632
2633	! titau2 = titau12
2634	is_ext=0
2635	ie_ext=0
2636	js_ext=0
2637	je_ext=1
2638	CALL cal_titau_12_21( config_flags, titau2, &
2639	mu, xkmh, defor12, &
2640	nba_mij(ims,kms,jms,P_m12), & !JDM
2641	is_ext, ie_ext, js_ext, je_ext, &
2642	ids, ide, jds, jde, kds, kde, &
2643	ims, ime, jms, jme, kms, kme, &
2644	its, ite, jts, jte, kts, kte )
2645
2646	! titau1avg = titau11avg
2647	! titau2avg = titau12avg
2648
2649	DO j = j_start, j_end
2650	DO k = kts+1,ktf
2651	DO i = i_start, i_end
2652	titau1avg(i,k,j)=0.5(fnm(k)(titau1(i-1,k ,j)+titau1(i,k ,j))+ &
2653	fnp(k)*(titau1(i-1,k-1,j)+titau1(i,k-1,j)))
2654	titau2avg(i,k,j)=0.5(fnm(k)(titau2(i,k ,j+1)+titau2(i,k ,j))+ &
2655	fnp(k)*(titau2(i,k-1,j+1)+titau2(i,k-1,j)))
2656	tmpzy = 0.25*( zy(i-1,k,j )+zy(i,k,j )+ &
2657	zy(i-1,k,j+1)+zy(i,k,j+1) )
2658	! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
2659	! titau1avg(i,k,j)=titau1avg(i,k,j)zx(i,k,j)tmpzeta_z
2660	! titau2avg(i,k,j)=titau2avg(i,k,j)tmpzy tmpzeta_z
2661
2662	titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j)
2663	titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy
2664
2665	ENDDO
2666	ENDDO
2667	ENDDO
2668	!
2669	DO j = j_start, j_end
2670	DO i = i_start, i_end
2671	titau1avg(i,kts,j)=0.
2672	titau1avg(i,ktf+1,j)=0.
2673	titau2avg(i,kts,j)=0.
2674	titau2avg(i,ktf+1,j)=0.
2675	ENDDO
2676	ENDDO
2677	!
2678	DO j = j_start, j_end
2679	DO k = kts,ktf
2680	DO i = i_start, i_end
2681
2682	mrdx=msfux(i,j)*rdx
2683	mrdy=msfuy(i,j)*rdy
2684	tendency(i,k,j)=tendency(i,k,j)- &
2685	(mrdx*(titau1(i,k,j )-titau1(i-1,k,j))+ &
2686	mrdy*(titau2(i,k,j+1)-titau2(i,k,j ))- &
2687	msfuy(i,j)rdzw(i,k,j)((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
2688	(titau2avg(i,k+1,j)-titau2avg(i,k,j)) &
2689	) )
2690	ENDDO
2691	ENDDO
2692	ENDDO
2693
2694	END SUBROUTINE horizontal_diffusion_u_2
2695
2696	!=======================================================================
2697	!=======================================================================
2698
2699	SUBROUTINE horizontal_diffusion_v_2( tendency, mu, config_flags, &
2700	defor12, defor22, div, &
2701	nba_mij, n_nba_mij, & !JDM
2702	tke, &
2703	msfvx, msfvy, &
2704	xkmh, rdx, rdy, fnm, fnp, &
2705	zx, zy, rdzw, &
2706	ids, ide, jds, jde, kds, kde, &
2707	ims, ime, jms, jme, kms, kme, &
2708	its, ite, jts, jte, kts, kte )
2709
2710	!-----------------------------------------------------------------------
2711	! Begin declarations.
2712
2713	IMPLICIT NONE
2714
2715	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2716
2717	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2718	ims, ime, jms, jme, kms, kme, &
2719	its, ite, jts, jte, kts, kte
2720
2721	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2722	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2723
2724	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvx, &
2725	msfvy, &
2726	mu
2727
2728	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
2729
2730	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor12, &
2731	defor22, &
2732	div, &
2733	tke, &
2734	xkmh, &
2735	zx, &
2736	zy, &
2737	rdzw
2738
2739	INTEGER, INTENT( IN ) :: n_nba_mij !JDM
2740
2741	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
2742	:: nba_mij
2743
2744	REAL , INTENT(IN ) :: rdx, &
2745	rdy
2746
2747	! Local data
2748
2749	INTEGER :: i, j, k, ktf
2750
2751	INTEGER :: i_start, i_end, j_start, j_end
2752	INTEGER :: is_ext,ie_ext,js_ext,je_ext
2753
2754	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
2755	titau2avg, &
2756	titau1, &
2757	titau2, &
2758	xkxavg, &
2759	rravg
2760	! new
2761	! zxavg, &
2762	! zyavg
2763
2764	REAL :: mrdx, mrdy, rcoup
2765
2766	REAL :: tmpzx, tmpzeta_z
2767
2768	! End declarations.
2769	!-----------------------------------------------------------------------
2770
2771	ktf=MIN(kte,kde-1)
2772
2773	!-----------------------------------------------------------------------
2774	! v : p (.), v(+), w(-)
2775	!
2776	! p v p v v v
2777	!
2778	! p + . + . + . + k+1 + . + . + . + k+1
2779	!
2780	! w - 23 - - k+1 23 k+1
2781	!
2782	! p + 22 O 22 + . + k + 21 O 21 + . + k
2783	!
2784	! w - 23 - - k 23 k
2785	!
2786	! p + . + . + . + k-1 + . + . + . + k-1
2787	!
2788	! j-1 j j j+1 i-1 i i i+1 i+1
2789	!
2790
2791	i_start = its
2792	i_end = MIN(ite,ide-1)
2793	j_start = jts
2794	j_end = jte
2795
2796	IF ( config_flags%open_xs .or. config_flags%specified .or. &
2797	config_flags%nested) i_start = MAX(ids+1,its)
2798	IF ( config_flags%open_xe .or. config_flags%specified .or. &
2799	config_flags%nested) i_end = MIN(ide-2,ite)
2800	IF ( config_flags%open_ys .or. config_flags%specified .or. &
2801	config_flags%nested) j_start = MAX(jds+1,jts)
2802	IF ( config_flags%open_ye .or. config_flags%specified .or. &
2803	config_flags%nested) j_end = MIN(jde-1,jte)
2804	IF ( config_flags%periodic_x ) i_start = its
2805	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2806
2807	! titau1 = titau21
2808	is_ext=0
2809	ie_ext=1
2810	js_ext=0
2811	je_ext=0
2812	CALL cal_titau_12_21( config_flags, titau1, &
2813	mu, xkmh, defor12, &
2814	nba_mij(ims,kms,jms,P_m12), & !JDM
2815	is_ext,ie_ext,js_ext,je_ext, &
2816	ids, ide, jds, jde, kds, kde, &
2817	ims, ime, jms, jme, kms, kme, &
2818	its, ite, jts, jte, kts, kte )
2819
2820	! titau2 = titau22
2821	is_ext=0
2822	ie_ext=0
2823	js_ext=1
2824	je_ext=0
2825	CALL cal_titau_11_22_33( config_flags, titau2, &
2826	mu, tke, xkmh, defor22, &
2827	nba_mij(ims,kms,jms,P_m22), & !JDM
2828	is_ext, ie_ext, js_ext, je_ext, &
2829	ids, ide, jds, jde, kds, kde, &
2830	ims, ime, jms, jme, kms, kme, &
2831	its, ite, jts, jte, kts, kte )
2832
2833	DO j = j_start, j_end
2834	DO k = kts+1,ktf
2835	DO i = i_start, i_end
2836	titau1avg(i,k,j)=0.5(fnm(k)(titau1(i+1,k ,j)+titau1(i,k ,j))+ &
2837	fnp(k)*(titau1(i+1,k-1,j)+titau1(i,k-1,j)))
2838	titau2avg(i,k,j)=0.5(fnm(k)(titau2(i,k ,j-1)+titau2(i,k ,j))+ &
2839	fnp(k)*(titau2(i,k-1,j-1)+titau2(i,k-1,j)))
2840
2841	tmpzx = 0.25*( zx(i,k,j )+zx(i+1,k,j )+ &
2842	zx(i,k,j-1)+zx(i+1,k,j-1) )
2843
2844
2845	titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
2846	titau2avg(i,k,j)=titau2avg(i,k,j)*zy(i,k,j)
2847
2848
2849	ENDDO
2850	ENDDO
2851	ENDDO
2852
2853	DO j = j_start, j_end
2854	DO i = i_start, i_end
2855	titau1avg(i,kts,j)=0.
2856	titau1avg(i,ktf+1,j)=0.
2857	titau2avg(i,kts,j)=0.
2858	titau2avg(i,ktf+1,j)=0.
2859	ENDDO
2860	ENDDO
2861	!
2862	DO j = j_start, j_end
2863	DO k = kts,ktf
2864	DO i = i_start, i_end
2865
2866	mrdx=msfvx(i,j)*rdx
2867	mrdy=msfvy(i,j)*rdy
2868	tendency(i,k,j)=tendency(i,k,j)- &
2869	(mrdy*(titau2(i ,k,j)-titau2(i,k,j-1))+ &
2870	mrdx*(titau1(i+1,k,j)-titau1(i,k,j ))- &
2871	msfvy(i,j)rdzw(i,k,j)((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
2872	(titau2avg(i,k+1,j)-titau2avg(i,k,j)) &
2873	) &
2874	)
2875
2876	ENDDO
2877	ENDDO
2878	ENDDO
2879
2880	END SUBROUTINE horizontal_diffusion_v_2
2881
2882	!=======================================================================
2883	!=======================================================================
2884
2885	SUBROUTINE horizontal_diffusion_w_2( tendency, mu, config_flags, &
2886	defor13, defor23, div, &
2887	nba_mij, n_nba_mij, & !JDM
2888	tke, &
2889	msftx, msfty, &
2890	xkmh, rdx, rdy, fnm, fnp, &
2891	zx, zy, rdz, &
2892	ids, ide, jds, jde, kds, kde, &
2893	ims, ime, jms, jme, kms, kme, &
2894	its, ite, jts, jte, kts, kte )
2895
2896	!-----------------------------------------------------------------------
2897	! Begin declarations.
2898
2899	IMPLICIT NONE
2900
2901	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2902
2903	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2904	ims, ime, jms, jme, kms, kme, &
2905	its, ite, jts, jte, kts, kte
2906
2907	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
2908	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
2909
2910	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msftx, &
2911	msfty, &
2912	mu
2913
2914	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
2915
2916	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
2917	defor23, &
2918	div, &
2919	tke, &
2920	xkmh, &
2921	zx, &
2922	zy, &
2923	rdz
2924
2925	INTEGER, INTENT( IN ) :: n_nba_mij !JDM
2926
2927	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
2928	:: nba_mij
2929
2930	REAL , INTENT(IN ) :: rdx, &
2931	rdy
2932
2933	! Local data
2934
2935	INTEGER :: i, j, k, ktf
2936
2937	INTEGER :: i_start, i_end, j_start, j_end
2938	INTEGER :: is_ext,ie_ext,js_ext,je_ext
2939
2940	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
2941	titau2avg, &
2942	titau1, &
2943	titau2, &
2944	xkxavg, &
2945	rravg
2946	! new
2947	! zxavg, &
2948	! zyavg
2949
2950	REAL :: mrdx, mrdy, rcoup
2951
2952	REAL :: tmpzx, tmpzy, tmpzeta_z
2953
2954	! End declarations.
2955	!-----------------------------------------------------------------------
2956
2957	ktf=MIN(kte,kde-1)
2958
2959	!-----------------------------------------------------------------------
2960	! w : p (.), u(\|), v(+), w(-)
2961	!
2962	! p u p u p v p v
2963	!
2964	! w - - - k+1 w - - - k+1
2965	!
2966	! p . \| 33 \| . k p . + 33 + . k
2967	!
2968	! w - 31 O 31 - k w - 32 O 32 - k
2969	!
2970	! p . \| 33 \| . k-1 p . \| 33 \| . k-1
2971	!
2972	! w - - - k-1 w - - - k-1
2973	!
2974	! i-1 i i i+1 j-1 j j j+1
2975	!
2976	i_start = its
2977	i_end = MIN(ite,ide-1)
2978	j_start = jts
2979	j_end = MIN(jte,jde-1)
2980
2981	IF ( config_flags%open_xs .or. config_flags%specified .or. &
2982	config_flags%nested) i_start = MAX(ids+1,its)
2983	IF ( config_flags%open_xe .or. config_flags%specified .or. &
2984	config_flags%nested) i_end = MIN(ide-2,ite)
2985	IF ( config_flags%open_ys .or. config_flags%specified .or. &
2986	config_flags%nested) j_start = MAX(jds+1,jts)
2987	IF ( config_flags%open_ye .or. config_flags%specified .or. &
2988	config_flags%nested) j_end = MIN(jde-2,jte)
2989	IF ( config_flags%periodic_x ) i_start = its
2990	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2991
2992	! titau1 = titau31
2993	is_ext=0
2994	ie_ext=1
2995	js_ext=0
2996	je_ext=0
2997	CALL cal_titau_13_31( config_flags, titau1, defor13, &
2998	nba_mij(ims,kms,jms,P_m13), & !JDM
2999	mu, xkmh, fnm, fnp, &
3000	is_ext, ie_ext, js_ext, je_ext, &
3001	ids, ide, jds, jde, kds, kde, &
3002	ims, ime, jms, jme, kms, kme, &
3003	its, ite, jts, jte, kts, kte )
3004
3005	! titau2 = titau32
3006	is_ext=0
3007	ie_ext=0
3008	js_ext=0
3009	je_ext=1
3010	CALL cal_titau_23_32( config_flags, titau2, defor23, &
3011	nba_mij(ims,kms,jms,P_m23), & !JDM
3012	mu, xkmh, fnm, fnp, &
3013	is_ext, ie_ext, js_ext, je_ext, &
3014	ids, ide, jds, jde, kds, kde, &
3015	ims, ime, jms, jme, kms, kme, &
3016	its, ite, jts, jte, kts, kte )
3017
3018	! titau1avg = titau31avg * zx * zeta_z = titau13avg * zx * zeta_z
3019	! titau2avg = titau32avg * zy * zeta_z = titau23avg * zy * zeta_z
3020
3021	DO j = j_start, j_end
3022	DO k = kts,ktf
3023	DO i = i_start, i_end
3024	titau1avg(i,k,j)=0.25*(titau1(i+1,k+1,j)+titau1(i,k+1,j)+ &
3025	titau1(i+1,k ,j)+titau1(i,k ,j))
3026	titau2avg(i,k,j)=0.25*(titau2(i,k+1,j+1)+titau2(i,k+1,j)+ &
3027	titau2(i,k ,j+1)+titau2(i,k ,j))
3028	! new
3029	tmpzx =0.25*( zx(i,k ,j)+zx(i+1,k ,j)+ &
3030	zx(i,k+1,j)+zx(i+1,k+1,j) )
3031	tmpzy =0.25*( zy(i,k ,j)+zy(i,k ,j+1)+ &
3032	zy(i,k+1,j)+zy(i,k+1,j+1) )
3033
3034	titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
3035	titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy
3036	! titau1avg(i,k,j)=titau1avg(i,k,j)tmpzxzeta_z(i,j)
3037	! titau2avg(i,k,j)=titau2avg(i,k,j)tmpzyzeta_z(i,j)
3038	ENDDO
3039	ENDDO
3040	ENDDO
3041
3042	DO j = j_start, j_end
3043	DO i = i_start, i_end
3044	titau1avg(i,ktf+1,j)=0.
3045	titau2avg(i,ktf+1,j)=0.
3046	ENDDO
3047	ENDDO
3048
3049	DO j = j_start, j_end
3050	DO k = kts+1,ktf
3051	DO i = i_start, i_end
3052
3053	mrdx=msftx(i,j)*rdx
3054	mrdy=msfty(i,j)*rdy
3055
3056	tendency(i,k,j)=tendency(i,k,j)- &
3057	(mrdx*(titau1(i+1,k,j)-titau1(i,k,j))+ &
3058	mrdy*(titau2(i,k,j+1)-titau2(i,k,j))- &
3059	msfty(i,j)rdz(i,k,j)(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
3060	titau2avg(i,k,j)-titau2avg(i,k-1,j) &
3061	) &
3062	)
3063	! msft(i,j)/dn(k)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
3064	! titau2avg(i,k,j)-titau2avg(i,k-1,j) &
3065	! ) &
3066	! )
3067	ENDDO
3068	ENDDO
3069	ENDDO
3070
3071	END SUBROUTINE horizontal_diffusion_w_2
3072
3073	!=======================================================================
3074	!=======================================================================
3075
3076	SUBROUTINE horizontal_diffusion_s (tendency, mu, config_flags, var, &
3077	msftx, msfty, msfux, msfuy, &
3078	msfvx, msfvy, xkhh, rdx, rdy, &
3079	fnm, fnp, cf1, cf2, cf3, &
3080	zx, zy, rdz, rdzw, dnw, dn, &
3081	doing_tke, &
3082	ids, ide, jds, jde, kds, kde, &
3083	ims, ime, jms, jme, kms, kme, &
3084	its, ite, jts, jte, kts, kte )
3085
3086	!-----------------------------------------------------------------------
3087	! Begin declarations.
3088
3089	IMPLICIT NONE
3090
3091	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3092
3093	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3094	ims, ime, jms, jme, kms, kme, &
3095	its, ite, jts, jte, kts, kte
3096
3097	LOGICAL, INTENT(IN ) :: doing_tke
3098
3099	REAL , INTENT(IN ) :: cf1, cf2, cf3
3100
3101	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3102	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3103	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
3104	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3105
3106	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux
3107	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfuy
3108	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvx
3109	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvy
3110	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msftx
3111	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfty
3112
3113	REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: mu
3114
3115	! REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), &
3116	! INTENT(IN ) :: xkhh
3117
3118	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
3119
3120	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: &
3121	xkhh, &
3122	rdz, &
3123	rdzw
3124
3125	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: var, &
3126	zx, &
3127	zy
3128
3129	REAL , INTENT(IN ) :: rdx, &
3130	rdy
3131
3132	! Local data
3133
3134	INTEGER :: i, j, k, ktf
3135
3136	INTEGER :: i_start, i_end, j_start, j_end
3137
3138	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: H1avg, &
3139	H2avg, &
3140	H1, &
3141	H2, &
3142	xkxavg
3143	! new
3144	! zxavg, &
3145	! zyavg
3146
3147	REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
3148
3149	REAL :: mrdx, mrdy, rcoup
3150	REAL :: tmpzx, tmpzy, tmpzeta_z, rdzu, rdzv
3151	INTEGER :: ktes1,ktes2
3152
3153	! End declarations.
3154	!-----------------------------------------------------------------------
3155
3156	ktf=MIN(kte,kde-1)
3157
3158	!-----------------------------------------------------------------------
3159	! scalars: t (.), u(\|), v(+), w(-)
3160	!
3161	! t u t u t v t v
3162	!
3163	! w - 3 - k+1 w - 3 - k+1
3164	!
3165	! t . 1 O 1 . k t . 2 O 2 . k
3166	!
3167	! w - 3 - k w - 3 - k
3168	!
3169	! t . \| . \| . k-1 t . + . + . k-1
3170	!
3171	! w - - - k-1 w - - - k-1
3172	!
3173	! t i-1 i i i+1 j-1 j j j+1
3174	!
3175
3176	ktes1=kte-1
3177	ktes2=kte-2
3178
3179	i_start = its
3180	i_end = MIN(ite,ide-1)
3181	j_start = jts
3182	j_end = MIN(jte,jde-1)
3183
3184	IF ( config_flags%open_xs .or. config_flags%specified .or. &
3185	config_flags%nested) i_start = MAX(ids+1,its)
3186	IF ( config_flags%open_xe .or. config_flags%specified .or. &
3187	config_flags%nested) i_end = MIN(ide-2,ite)
3188	IF ( config_flags%open_ys .or. config_flags%specified .or. &
3189	config_flags%nested) j_start = MAX(jds+1,jts)
3190	IF ( config_flags%open_ye .or. config_flags%specified .or. &
3191	config_flags%nested) j_end = MIN(jde-2,jte)
3192	IF ( config_flags%periodic_x ) i_start = its
3193	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
3194
3195	! diffusion of the TKE needs mutiple 2
3196
3197	IF ( doing_tke ) THEN
3198	DO j = j_start, j_end
3199	DO k = kts,ktf
3200	DO i = i_start, i_end
3201	tmptendf(i,k,j)=tendency(i,k,j)
3202	ENDDO
3203	ENDDO
3204	ENDDO
3205	ENDIF
3206
3207	! H1 = partial var over partial x
3208
3209	DO j = j_start, j_end
3210	DO k = kts, ktf
3211	DO i = i_start, i_end + 1
3212	! new
3213	! zxavg(i,k,j) =0.5*( zx(i-1,k,j)+ zx(i,k,j))
3214	xkxavg(i,k,j)=0.5*(xkhh(i-1,k,j)+xkhh(i,k,j))
3215	ENDDO
3216	ENDDO
3217	ENDDO
3218
3219	DO j = j_start, j_end
3220	DO k = kts+1, ktf
3221	DO i = i_start, i_end + 1
3222	H1avg(i,k,j)=0.5(fnm(k)(var(i-1,k ,j)+var(i,k ,j))+ &
3223	fnp(k)*(var(i-1,k-1,j)+var(i,k-1,j)))
3224	ENDDO
3225	ENDDO
3226	ENDDO
3227
3228	DO j = j_start, j_end
3229	DO i = i_start, i_end + 1
3230	H1avg(i,kts ,j)=0.5(cf1var(i ,1,j)+cf2*var(i ,2,j)+ &
3231	cf3var(i ,3,j)+cf1var(i-1,1,j)+ &
3232	cf2var(i-1,2,j)+cf3var(i-1,3,j))
3233	H1avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
3234	var(i,ktes2,j))0.5dnw(ktes1)/dn(ktes1)+ &
3235	var(i-1,ktes1,j)+(var(i-1,ktes1,j)- &
3236	var(i-1,ktes2,j))0.5dnw(ktes1)/dn(ktes1))
3237	ENDDO
3238	ENDDO
3239
3240	DO j = j_start, j_end
3241	DO k = kts, ktf
3242	DO i = i_start, i_end + 1
3243	! new
3244	tmpzx = 0.5*( zx(i,k,j)+ zx(i,k+1,j))
3245	rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
3246	H1(i,k,j)=-msfuy(i,j)xkxavg(i,k,j)( &
3247	rdx(var(i,k,j)-var(i-1,k,j)) - tmpzx &
3248	(H1avg(i,k+1,j)-H1avg(i,k,j))*rdzu )
3249
3250	! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
3251	! H1(i,k,j)=-msfuy(i,j)xkxavg(i,k,j)( &
3252	! rdx(var(i,k,j)-var(i-1,k,j)) - tmpzxtmpzeta_z* &
3253	! (H1avg(i,k+1,j)-H1avg(i,k,j))/dnw(k))
3254	ENDDO
3255	ENDDO
3256	ENDDO
3257
3258	! H2 = partial var over partial y
3259
3260	DO j = j_start, j_end + 1
3261	DO k = kts, ktf
3262	DO i = i_start, i_end
3263	! new
3264	! zyavg(i,k,j) =0.5*( zy(i,k,j-1)+ zy(i,k,j))
3265	xkxavg(i,k,j)=0.5*(xkhh(i,k,j-1)+xkhh(i,k,j))
3266	ENDDO
3267	ENDDO
3268	ENDDO
3269
3270	DO j = j_start, j_end + 1
3271	DO k = kts+1, ktf
3272	DO i = i_start, i_end
3273	! new
3274	H2avg(i,k,j)=0.5(fnm(k)(var(i,k ,j-1)+var(i,k ,j))+ &
3275	fnp(k)*(var(i,k-1,j-1)+var(i,k-1,j)))
3276	ENDDO
3277	ENDDO
3278	ENDDO
3279
3280	DO j = j_start, j_end + 1
3281	DO i = i_start, i_end
3282	H2avg(i,kts ,j)=0.5(cf1var(i,1,j )+cf2*var(i ,2,j)+ &
3283	cf3var(i,3,j )+cf1var(i,1,j-1)+ &
3284	cf2var(i,2,j-1)+cf3var(i,3,j-1))
3285	H2avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
3286	var(i,ktes2,j))0.5dnw(ktes1)/dn(ktes1)+ &
3287	var(i,ktes1,j-1)+(var(i,ktes1,j-1)- &
3288	var(i,ktes2,j-1))0.5dnw(ktes1)/dn(ktes1))
3289	ENDDO
3290	ENDDO
3291
3292	DO j = j_start, j_end + 1
3293	DO k = kts, ktf
3294	DO i = i_start, i_end
3295	! new
3296	tmpzy = 0.5*( zy(i,k,j)+ zy(i,k+1,j))
3297	rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
3298	H2(i,k,j)=-msfvy(i,j)xkxavg(i,k,j)( &
3299	rdy(var(i,k,j)-var(i,k,j-1)) - tmpzy &
3300	(H2avg(i ,k+1,j)-H2avg(i,k,j))*rdzv)
3301
3302	! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i,j-1))
3303	! H2(i,k,j)=-msfvy(i,j)xkxavg(i,k,j)( &
3304	! rdy(var(i,k,j)-var(i,k,j-1)) - tmpzytmpzeta_z* &
3305	! (H2avg(i ,k+1,j)-H2avg(i,k,j))/dnw(k))
3306	ENDDO
3307	ENDDO
3308	ENDDO
3309
3310	DO j = j_start, j_end
3311	DO k = kts+1, ktf
3312	DO i = i_start, i_end
3313	H1avg(i,k,j)=0.5(fnm(k)(H1(i+1,k ,j)+H1(i,k ,j))+ &
3314	fnp(k)*(H1(i+1,k-1,j)+H1(i,k-1,j)))
3315	H2avg(i,k,j)=0.5(fnm(k)(H2(i,k ,j+1)+H2(i,k ,j))+ &
3316	fnp(k)*(H2(i,k-1,j+1)+H2(i,k-1,j)))
3317	! new
3318	! zxavg(i,k,j)=fnm(k)zx(i,k,j)+fnp(k)zx(i,k-1,j)
3319	! zyavg(i,k,j)=fnm(k)zy(i,k,j)+fnp(k)zy(i,k-1,j)
3320
3321	! H1avg(i,k,j)=zxH1avgzeta_z
3322	! H2avg(i,k,j)=zyH2avgzeta_z
3323
3324	tmpzx = 0.5*( zx(i,k,j)+ zx(i+1,k,j ))
3325	tmpzy = 0.5*( zy(i,k,j)+ zy(i ,k,j+1))
3326
3327	H1avg(i,k,j)=H1avg(i,k,j)*tmpzx
3328	H2avg(i,k,j)=H2avg(i,k,j)*tmpzy
3329
3330	! H1avg(i,k,j)=H1avg(i,k,j)tmpzxzeta_z(i,j)
3331	! H2avg(i,k,j)=H2avg(i,k,j)tmpzyzeta_z(i,j)
3332	ENDDO
3333	ENDDO
3334	ENDDO
3335
3336	DO j = j_start, j_end
3337	DO i = i_start, i_end
3338	H1avg(i,kts ,j)=0.
3339	H1avg(i,ktf+1,j)=0.
3340	H2avg(i,kts ,j)=0.
3341	H2avg(i,ktf+1,j)=0.
3342	ENDDO
3343	ENDDO
3344
3345	DO j = j_start, j_end
3346	DO k = kts,ktf
3347	DO i = i_start, i_end
3348
3349	mrdx=msftx(i,j)*rdx
3350	mrdy=msfty(i,j)*rdy
3351
3352	tendency(i,k,j)=tendency(i,k,j)- &
3353	(mrdx0.5((mu(i+1,j)+mu(i,j))*H1(i+1,k,j)- &
3354	(mu(i-1,j)+mu(i,j))*H1(i ,k,j))+ &
3355	mrdy0.5((mu(i,j+1)+mu(i,j))*H2(i,k,j+1)- &
3356	(mu(i,j-1)+mu(i,j))*H2(i,k,j ))- &
3357	msfty(i,j)mu(i,j)(H1avg(i,k+1,j)-H1avg(i,k,j)+ &
3358	H2avg(i,k+1,j)-H2avg(i,k,j) &
3359	)*rdzw(i,k,j) &
3360	)
3361
3362	ENDDO
3363	ENDDO
3364	ENDDO
3365
3366	IF ( doing_tke ) THEN
3367	DO j = j_start, j_end
3368	DO k = kts,ktf
3369	DO i = i_start, i_end
3370	tendency(i,k,j)=tmptendf(i,k,j)+2.* &
3371	(tendency(i,k,j)-tmptendf(i,k,j))
3372	ENDDO
3373	ENDDO
3374	ENDDO
3375	ENDIF
3376
3377	END SUBROUTINE horizontal_diffusion_s
3378
3379	!=======================================================================
3380	!=======================================================================
3381
3382	SUBROUTINE vertical_diffusion_2 ( ru_tendf, rv_tendf, rw_tendf, rt_tendf, &
3383	tke_tendf, moist_tendf, n_moist, &
3384	chem_tendf, n_chem, &
3385	scalar_tendf, n_scalar, &
3386	tracer_tendf, n_tracer, &
3387	u_2, v_2, &
3388	thp,u_base,v_base,t_base,qv_base,mu,tke, &
3389	config_flags,defor13,defor23,defor33, &
3390	nba_mij, n_nba_mij, & !JDM
3391	div, &
3392	moist,chem,scalar,tracer,xkmv,xkhv,km_opt,&
3393	fnm, fnp, dn, dnw, rdz, rdzw, &
3394	hfx, qfx, ust, rho, &
3395	ids, ide, jds, jde, kds, kde, &
3396	ims, ime, jms, jme, kms, kme, &
3397	its, ite, jts, jte, kts, kte )
3398
3399	!-----------------------------------------------------------------------
3400	! Begin declarations.
3401
3402	IMPLICIT NONE
3403
3404	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3405
3406	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3407	ims, ime, jms, jme, kms, kme, &
3408	its, ite, jts, jte, kts, kte
3409
3410	INTEGER , INTENT(IN ) :: n_moist,n_chem,n_scalar,n_tracer,km_opt
3411
3412	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3413	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3414	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3415	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
3416	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
3417
3418	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: qv_base
3419	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base
3420	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: v_base
3421	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: t_base
3422
3423	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::ru_tendf,&
3424	rv_tendf,&
3425	rw_tendf,&
3426	tke_tendf,&
3427	rt_tendf
3428
3429	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
3430	INTENT(INOUT) :: moist_tendf
3431
3432	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
3433	INTENT(INOUT) :: chem_tendf
3434
3435	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
3436	INTENT(INOUT) :: scalar_tendf
3437	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_tracer) , &
3438	INTENT(INOUT) :: tracer_tendf
3439
3440	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
3441	INTENT(INOUT) :: moist
3442
3443	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), &
3444	INTENT(INOUT) :: chem
3445
3446	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
3447	INTENT(IN ) :: scalar
3448	REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_tracer) , &
3449	INTENT(IN ) :: tracer
3450
3451	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
3452	defor23, &
3453	defor33, &
3454	div, &
3455	xkmv, &
3456	xkhv, &
3457	tke, &
3458	rdz, &
3459	u_2, &
3460	v_2, &
3461	rdzw
3462
3463	INTEGER, INTENT( IN ) :: n_nba_mij !JDM
3464
3465	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
3466	:: nba_mij
3467
3468	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rho
3469	REAL , DIMENSION( ims:ime, jms:jme), INTENT(INOUT) :: hfx, &
3470	qfx
3471	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: ust
3472	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: thp
3473
3474	! LOCAL VAR
3475
3476	REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: var_mix
3477
3478	INTEGER :: im, i,j,k
3479	INTEGER :: i_start, i_end, j_start, j_end
3480
3481	! REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: xkhv
3482
3483	!***************************************************************************
3484	!***************************************************************************
3485	!MODIFICA VARIABILI PER I FLUSSI
3486	!
3487	REAL :: V0_u,V0_v,tao_xz,tao_yz,ustar,cd0
3488	REAL :: xsfc,psi1,vk2,zrough,lnz
3489	REAL :: heat_flux, moist_flux, heat_flux0
3490	REAL :: cpm
3491	!
3492	!FINE MODIFICA VARIABILI PER I FLUSSI
3493	!***************************************************************************
3494	!
3495
3496	! End declarations.
3497	!-----------------------------------------------------------------------
3498
3499	i_start = its
3500	i_end = MIN(ite,ide-1)
3501	j_start = jts
3502	j_end = MIN(jte,jde-1)
3503	!
3504	!-----------------------------------------------------------------------
3505
3506	CALL vertical_diffusion_u_2( ru_tendf, config_flags, mu, &
3507	defor13, xkmv, &
3508	nba_mij, n_nba_mij, & !JDM
3509	dnw, rdzw, fnm, fnp, &
3510	ids, ide, jds, jde, kds, kde, &
3511	ims, ime, jms, jme, kms, kme, &
3512	its, ite, jts, jte, kts, kte )
3513
3514
3515	CALL vertical_diffusion_v_2( rv_tendf, config_flags, mu, &
3516	defor23, xkmv, &
3517	nba_mij, n_nba_mij, & !JDM
3518	dnw, rdzw, fnm, fnp, &
3519	ids, ide, jds, jde, kds, kde, &
3520	ims, ime, jms, jme, kms, kme, &
3521	its, ite, jts, jte, kts, kte )
3522
3523	CALL vertical_diffusion_w_2( rw_tendf, config_flags, mu, &
3524	defor33, tke(ims,kms,jms), &
3525	nba_mij, n_nba_mij, & !JDM
3526	div, xkmv, &
3527	dn, rdz, &
3528	ids, ide, jds, jde, kds, kde, &
3529	ims, ime, jms, jme, kms, kme, &
3530	its, ite, jts, jte, kts, kte )
3531
3532	!*****************************************
3533	!*****************************************
3534	! MODIFICA al flusso di momento alla parete
3535	!
3536	vflux: SELECT CASE( config_flags%isfflx )
3537	CASE (0) ! Assume cd a constant, specified in namelist
3538	cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient
3539	! set in namelist.input
3540	!
3541	!calcolo del modulo della velocita
3542	DO j = j_start, j_end
3543	DO i = i_start, ite
3544	V0_u=0.
3545	tao_xz=0.
3546	V0_u= sqrt((u_2(i,kts,j)**2) + &
3547	(((v_2(i ,kts,j )+ &
3548	v_2(i ,kts,j+1)+ &
3549	v_2(i-1,kts,j )+ &
3550	v_2(i-1,kts,j+1))/4)**2))+epsilon
3551	tao_xz=cd0V0_uu_2(i,kts,j)
3552	ru_tendf(i,kts,j)=ru_tendf(i,kts,j) &
3553	-0.25(mu(i,j)+mu(i-1,j))tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
3554	ENDDO
3555	ENDDO
3556
3557	!
3558	DO j = j_start, jte
3559	DO i = i_start, i_end
3560	V0_v=0.
3561	tao_yz=0.
3562	V0_v= sqrt((v_2(i,kts,j)**2) + &
3563	(((u_2(i ,kts,j )+ &
3564	u_2(i ,kts,j-1)+ &
3565	u_2(i+1,kts,j )+ &
3566	u_2(i+1,kts,j-1))/4)**2))+epsilon
3567	tao_yz=cd0V0_vv_2(i,kts,j)
3568	rv_tendf(i,kts,j)=rv_tendf(i,kts,j) &
3569	-0.25(mu(i,j)+mu(i,j-1))tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
3570	ENDDO
3571	ENDDO
3572
3573	CASE (1,2) ! ustar computed from surface routine
3574	DO j = j_start, j_end
3575	DO i = i_start, ite
3576	V0_u=0.
3577	tao_xz=0.
3578	V0_u= sqrt((u_2(i,kts,j)**2) + &
3579	(((v_2(i ,kts,j )+ &
3580	v_2(i ,kts,j+1)+ &
3581	v_2(i-1,kts,j )+ &
3582	v_2(i-1,kts,j+1))/4)**2))+epsilon
3583	ustar=0.5*(ust(i,j)+ust(i-1,j))
3584	tao_xz=ustarustaru_2(i,kts,j)/V0_u
3585	ru_tendf(i,kts,j)=ru_tendf(i,kts,j) &
3586	-0.25(mu(i,j)+mu(i-1,j))tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
3587	ENDDO
3588	ENDDO
3589
3590	DO j = j_start, jte
3591	DO i = i_start, i_end
3592	V0_v=0.
3593	tao_yz=0.
3594	V0_v= sqrt((v_2(i,kts,j)**2) + &
3595	(((u_2(i ,kts,j )+ &
3596	u_2(i ,kts,j-1)+ &
3597	u_2(i+1,kts,j )+ &
3598	u_2(i+1,kts,j-1))/4)**2))+epsilon
3599	ustar=0.5*(ust(i,j)+ust(i,j-1))
3600	tao_yz=ustarustarv_2(i,kts,j)/V0_v
3601	rv_tendf(i,kts,j)=rv_tendf(i,kts,j) &
3602	-0.25(mu(i,j)+mu(i,j-1))tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
3603	ENDDO
3604	ENDDO
3605
3606	CASE DEFAULT
3607	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
3608	END SELECT vflux
3609
3610	!
3611	! FINE MODIFICA al flusso di momento alla parete
3612	!*****************************************
3613	!*****************************************
3614
3615	IF ( config_flags%mix_full_fields ) THEN
3616
3617	DO j=jts,min(jte,jde-1)
3618	DO k=kts,kte-1
3619	DO i=its,min(ite,ide-1)
3620	var_mix(i,k,j) = thp(i,k,j)
3621	ENDDO
3622	ENDDO
3623	ENDDO
3624
3625	ELSE
3626
3627	DO j=jts,min(jte,jde-1)
3628	DO k=kts,kte-1
3629	DO i=its,min(ite,ide-1)
3630	var_mix(i,k,j) = thp(i,k,j) - t_base(k)
3631	ENDDO
3632	ENDDO
3633	ENDDO
3634
3635	END IF
3636
3637	CALL vertical_diffusion_s( rt_tendf, config_flags, var_mix, mu, xkhv, &
3638	dn, dnw, rdz, rdzw, fnm, fnp, &
3639	.false., &
3640	ids, ide, jds, jde, kds, kde, &
3641	ims, ime, jms, jme, kms, kme, &
3642	its, ite, jts, jte, kts, kte )
3643
3644
3645	!*****************************************
3646	!*****************************************
3647	!MODIFICA al flusso di calore
3648	!
3649	!
3650	hflux: SELECT CASE( config_flags%isfflx )
3651	CASE (0,2) ! with fixed surface heat flux given in the namelist
3652	heat_flux = config_flags%tke_heat_flux ! constant heat flux value
3653	! set in namelist.input
3654	DO j = j_start, j_end
3655	DO i = i_start, i_end
3656	cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV))
3657	hfx(i,j)=heat_fluxcprho(i,1,j) ! provided for output only
3658	rt_tendf(i,kts,j)=rt_tendf(i,kts,j) &
3659	+mu(i,j)heat_fluxrdzw(i,kts,j)
3660	ENDDO
3661	ENDDO
3662
3663	CASE (1) ! use surface heat flux computed from surface routine
3664	DO j = j_start, j_end
3665	DO i = i_start, i_end
3666
3667	cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV))
3668	heat_flux = hfx(i,j)/cpm/rho(i,1,j)
3669	rt_tendf(i,kts,j)=rt_tendf(i,kts,j) &
3670	+mu(i,j)heat_fluxrdzw(i,kts,j)
3671
3672	ENDDO
3673	ENDDO
3674
3675	CASE DEFAULT
3676	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
3677	END SELECT hflux
3678
3679	!
3680	! FINE MODIFICA al flusso di calore
3681	!*****************************************
3682	!*****************************************
3683
3684	If (km_opt .eq. 2) then
3685	CALL vertical_diffusion_s( tke_tendf(ims,kms,jms), &
3686	config_flags, tke(ims,kms,jms), &
3687	mu, xkhv, &
3688	dn, dnw, rdz, rdzw, fnm, fnp, &
3689	.true., &
3690	ids, ide, jds, jde, kds, kde, &
3691	ims, ime, jms, jme, kms, kme, &
3692	its, ite, jts, jte, kts, kte )
3693	endif
3694
3695	IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN
3696
3697	moist_loop: do im = PARAM_FIRST_SCALAR, n_moist
3698
3699	IF ( (.not. config_flags%mix_full_fields) .and. (im == P_QV) ) THEN
3700
3701	DO j=jts,min(jte,jde-1)
3702	DO k=kts,kte-1
3703	DO i=its,min(ite,ide-1)
3704	var_mix(i,k,j) = moist(i,k,j,im) - qv_base(k)
3705	ENDDO
3706	ENDDO
3707	ENDDO
3708
3709	ELSE
3710
3711	DO j=jts,min(jte,jde-1)
3712	DO k=kts,kte-1
3713	DO i=its,min(ite,ide-1)
3714	var_mix(i,k,j) = moist(i,k,j,im)
3715	ENDDO
3716	ENDDO
3717	ENDDO
3718
3719	END IF
3720
3721
3722	CALL vertical_diffusion_s( moist_tendf(ims,kms,jms,im), &
3723	config_flags, var_mix, &
3724	mu, xkhv, &
3725	dn, dnw, rdz, rdzw, fnm, fnp, &
3726	.false., &
3727	ids, ide, jds, jde, kds, kde, &
3728	ims, ime, jms, jme, kms, kme, &
3729	its, ite, jts, jte, kts, kte )
3730
3731	!*****************************************
3732	!*****************************************
3733	!MODIFICATIONS for water vapor flux
3734	!
3735	!
3736	qflux: SELECT CASE( config_flags%isfflx )
3737	CASE (0)
3738	! do nothing
3739	CASE (1,2) ! with surface moisture flux
3740	IF ( im == P_QV ) THEN
3741	DO j = j_start, j_end
3742	DO i = i_start, i_end
3743	moist_flux = qfx(i,j)/rho(i,1,j)/(1.+moist(i,kts,j,P_QV))
3744	moist_tendf(i,kts,j,im)=moist_tendf(i,kts,j,im) &
3745	+mu(i,j)moist_fluxrdzw(i,kts,j)
3746	ENDDO
3747	ENDDO
3748	ENDIF
3749
3750	CASE DEFAULT
3751	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
3752	END SELECT qflux
3753	!
3754	! END MODIFICATIONS for water vapor flux
3755	!*****************************************
3756	!*****************************************
3757	ENDDO moist_loop
3758
3759	ENDIF
3760
3761	IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN
3762
3763	chem_loop: do im = PARAM_FIRST_SCALAR, n_chem
3764
3765	CALL vertical_diffusion_s( chem_tendf(ims,kms,jms,im), &
3766	config_flags, chem(ims,kms,jms,im), &
3767	mu, xkhv, &
3768	dn, dnw, rdz, rdzw, fnm, fnp, &
3769	.false., &
3770	ids, ide, jds, jde, kds, kde, &
3771	ims, ime, jms, jme, kms, kme, &
3772	its, ite, jts, jte, kts, kte )
3773	ENDDO chem_loop
3774
3775	ENDIF
3776
3777	IF (n_tracer .ge. PARAM_FIRST_SCALAR) THEN
3778
3779	tracer_loop: do im = PARAM_FIRST_SCALAR, n_tracer
3780
3781	CALL vertical_diffusion_s( tracer_tendf(ims,kms,jms,im), &
3782	config_flags, tracer(ims,kms,jms,im), &
3783	mu, xkhv, &
3784	dn, dnw, rdz, rdzw, fnm, fnp, &
3785	.false., &
3786	ids, ide, jds, jde, kds, kde, &
3787	ims, ime, jms, jme, kms, kme, &
3788	its, ite, jts, jte, kts, kte )
3789	ENDDO tracer_loop
3790
3791	ENDIF
3792
3793
3794	IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN
3795
3796	scalar_loop: do im = PARAM_FIRST_SCALAR, n_scalar
3797
3798	CALL vertical_diffusion_s( scalar_tendf(ims,kms,jms,im), &
3799	config_flags, scalar(ims,kms,jms,im), &
3800	mu, xkhv, &
3801	dn, dnw, rdz, rdzw, fnm, fnp, &
3802	.false., &
3803	ids, ide, jds, jde, kds, kde, &
3804	ims, ime, jms, jme, kms, kme, &
3805	its, ite, jts, jte, kts, kte )
3806	ENDDO scalar_loop
3807
3808	ENDIF
3809
3810	END SUBROUTINE vertical_diffusion_2
3811
3812	!=======================================================================
3813	!=======================================================================
3814
3815	SUBROUTINE vertical_diffusion_u_2( tendency, config_flags, mu, &
3816	defor13, xkmv, &
3817	nba_mij, n_nba_mij, & !JDM
3818	dnw, rdzw, fnm, fnp, &
3819	ids, ide, jds, jde, kds, kde, &
3820	ims, ime, jms, jme, kms, kme, &
3821	its, ite, jts, jte, kts, kte )
3822
3823	!-----------------------------------------------------------------------
3824	! Begin declarations.
3825
3826	IMPLICIT NONE
3827
3828	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3829
3830	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3831	ims, ime, jms, jme, kms, kme, &
3832	its, ite, jts, jte, kts, kte
3833
3834	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3835	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3836	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3837	! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
3838
3839	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3840
3841	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3842	INTENT(IN ) ::defor13, &
3843	xkmv, &
3844	rdzw
3845
3846	INTEGER, INTENT( IN ) :: n_nba_mij !JDM
3847
3848	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
3849	:: nba_mij
3850
3851	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
3852
3853	! LOCAL VARS
3854
3855	INTEGER :: i, j, k, ktf
3856
3857	INTEGER :: i_start, i_end, j_start, j_end
3858	INTEGER :: is_ext,ie_ext,js_ext,je_ext
3859
3860	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
3861
3862	REAL , DIMENSION( its:ite, jts:jte) :: zzavg
3863
3864	REAL :: rdzu
3865
3866	! End declarations.
3867	!-----------------------------------------------------------------------
3868
3869	ktf=MIN(kte,kde-1)
3870
3871	i_start = its
3872	i_end = ite
3873	j_start = jts
3874	j_end = MIN(jte,jde-1)
3875
3876	IF ( config_flags%open_xs .or. config_flags%specified .or. &
3877	config_flags%nested) i_start = MAX(ids+1,its)
3878	IF ( config_flags%open_xe .or. config_flags%specified .or. &
3879	config_flags%nested) i_end = MIN(ide-1,ite)
3880	IF ( config_flags%open_ys .or. config_flags%specified .or. &
3881	config_flags%nested) j_start = MAX(jds+1,jts)
3882	IF ( config_flags%open_ye .or. config_flags%specified .or. &
3883	config_flags%nested) j_end = MIN(jde-2,jte)
3884	IF ( config_flags%periodic_x ) i_start = its
3885	IF ( config_flags%periodic_x ) i_end = ite
3886
3887	! titau3 = titau13
3888	is_ext=0
3889	ie_ext=0
3890	js_ext=0
3891	je_ext=0
3892	CALL cal_titau_13_31( config_flags, titau3, defor13, &
3893	nba_mij(ims,kms,jms,P_m13), & !JDM
3894	mu, xkmv, fnm, fnp, &
3895	is_ext, ie_ext, js_ext, je_ext, &
3896	ids, ide, jds, jde, kds, kde, &
3897	ims, ime, jms, jme, kms, kme, &
3898	its, ite, jts, jte, kts, kte )
3899	!
3900	DO j = j_start, j_end
3901	DO k=kts+1,ktf
3902	DO i = i_start, i_end
3903
3904	rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
3905	tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j)-titau3(i,k,j))
3906
3907	ENDDO
3908	ENDDO
3909	ENDDO
3910
3911	! ******** MODIF...
3912	! we will pick up the surface drag (titau3(i,kts,j)) later
3913	!
3914	DO j = j_start, j_end
3915	k=kts
3916	DO i = i_start, i_end
3917
3918	rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
3919	tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j))
3920	ENDDO
3921	ENDDO
3922	! ******** MODIF...
3923
3924	END SUBROUTINE vertical_diffusion_u_2
3925
3926	!=======================================================================
3927	!=======================================================================
3928
3929	SUBROUTINE vertical_diffusion_v_2( tendency, config_flags, mu, &
3930	defor23, xkmv, &
3931	nba_mij, n_nba_mij, & !JDM
3932	dnw, rdzw, fnm, fnp, &
3933	ids, ide, jds, jde, kds, kde, &
3934	ims, ime, jms, jme, kms, kme, &
3935	its, ite, jts, jte, kts, kte )
3936	!-----------------------------------------------------------------------
3937	! Begin declarations.
3938
3939	IMPLICIT NONE
3940
3941	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3942
3943	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3944	ims, ime, jms, jme, kms, kme, &
3945	its, ite, jts, jte, kts, kte
3946	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
3947	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
3948	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
3949	! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
3950
3951	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
3952
3953	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
3954	INTENT(IN ) ::defor23, &
3955	xkmv, &
3956	rdzw
3957
3958	INTEGER, INTENT( IN ) :: n_nba_mij !JDM
3959
3960	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
3961	:: nba_mij
3962
3963	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
3964
3965	! LOCAL VARS
3966
3967	INTEGER :: i, j, k, ktf
3968
3969	INTEGER :: i_start, i_end, j_start, j_end
3970	INTEGER :: is_ext,ie_ext,js_ext,je_ext
3971
3972	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
3973
3974	REAL , DIMENSION( its:ite, jts:jte) :: zzavg
3975
3976	REAL :: rdzv
3977
3978	! End declarations.
3979	!-----------------------------------------------------------------------
3980
3981	ktf=MIN(kte,kde-1)
3982
3983	i_start = its
3984	i_end = MIN(ite,ide-1)
3985	j_start = jts
3986	j_end = jte
3987
3988	IF ( config_flags%open_xs .or. config_flags%specified .or. &
3989	config_flags%nested) i_start = MAX(ids+1,its)
3990	IF ( config_flags%open_xe .or. config_flags%specified .or. &
3991	config_flags%nested) i_end = MIN(ide-2,ite)
3992	IF ( config_flags%open_ys .or. config_flags%specified .or. &
3993	config_flags%nested) j_start = MAX(jds+1,jts)
3994	IF ( config_flags%open_ye .or. config_flags%specified .or. &
3995	config_flags%nested) j_end = MIN(jde-1,jte)
3996	IF ( config_flags%periodic_x ) i_start = its
3997	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
3998
3999	! titau3 = titau23
4000	is_ext=0
4001	ie_ext=0
4002	js_ext=0
4003	je_ext=0
4004	CALL cal_titau_23_32( config_flags, titau3, defor23, &
4005	nba_mij(ims,kms,jms,P_m23), & !JDM
4006	mu, xkmv, fnm, fnp, &
4007	is_ext, ie_ext, js_ext, je_ext, &
4008	ids, ide, jds, jde, kds, kde, &
4009	ims, ime, jms, jme, kms, kme, &
4010	its, ite, jts, jte, kts, kte )
4011
4012	DO j = j_start, j_end
4013	DO k = kts+1,ktf
4014	DO i = i_start, i_end
4015
4016	rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
4017	tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j)-titau3(i,k,j))
4018
4019	ENDDO
4020	ENDDO
4021	ENDDO
4022
4023	! ******** MODIF...
4024	! we will pick up the surface drag (titau3(i,kts,j)) later
4025	!
4026	DO j = j_start, j_end
4027	k=kts
4028	DO i = i_start, i_end
4029
4030	rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
4031	tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j))
4032
4033	ENDDO
4034	ENDDO
4035	! ******** MODIF...
4036
4037	END SUBROUTINE vertical_diffusion_v_2
4038
4039	!=======================================================================
4040	!=======================================================================
4041
4042	SUBROUTINE vertical_diffusion_w_2(tendency, config_flags, mu, &
4043	defor33, tke, &
4044	nba_mij, n_nba_mij, & !JDM
4045	div, xkmv, &
4046	dn, rdz, &
4047	ids, ide, jds, jde, kds, kde, &
4048	ims, ime, jms, jme, kms, kme, &
4049	its, ite, jts, jte, kts, kte )
4050
4051	!-----------------------------------------------------------------------
4052	! Begin declarations.
4053
4054	IMPLICIT NONE
4055
4056	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4057
4058	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4059	ims, ime, jms, jme, kms, kme, &
4060	its, ite, jts, jte, kts, kte
4061
4062	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
4063
4064	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
4065
4066	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
4067	INTENT(IN ) ::defor33, &
4068	tke, &
4069	div, &
4070	xkmv, &
4071	rdz
4072
4073	INTEGER, INTENT( IN ) :: n_nba_mij !JDM
4074
4075	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
4076	:: nba_mij
4077
4078	REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
4079
4080	! LOCAL VARS
4081
4082	INTEGER :: i, j, k, ktf
4083
4084	INTEGER :: i_start, i_end, j_start, j_end
4085	INTEGER :: is_ext,ie_ext,js_ext,je_ext
4086
4087	REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
4088
4089	! End declarations.
4090	!-----------------------------------------------------------------------
4091
4092	ktf=MIN(kte,kde-1)
4093
4094	i_start = its
4095	i_end = MIN(ite,ide-1)
4096	j_start = jts
4097	j_end = MIN(jte,jde-1)
4098
4099	IF ( config_flags%open_xs .or. config_flags%specified .or. &
4100	config_flags%nested) i_start = MAX(ids+1,its)
4101	IF ( config_flags%open_xe .or. config_flags%specified .or. &
4102	config_flags%nested) i_end = MIN(ide-2,ite)
4103	IF ( config_flags%open_ys .or. config_flags%specified .or. &
4104	config_flags%nested) j_start = MAX(jds+1,jts)
4105	IF ( config_flags%open_ye .or. config_flags%specified .or. &
4106	config_flags%nested) j_end = MIN(jde-2,jte)
4107	IF ( config_flags%periodic_x ) i_start = its
4108	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
4109
4110	! titau3 = titau33
4111	is_ext=0
4112	ie_ext=0
4113	js_ext=0
4114	je_ext=0
4115	CALL cal_titau_11_22_33( config_flags, titau3, &
4116	mu, tke, xkmv, defor33, &
4117	nba_mij(ims,kms,jms,P_m33), & !JDM
4118	is_ext, ie_ext, js_ext, je_ext, &
4119	ids, ide, jds, jde, kds, kde, &
4120	ims, ime, jms, jme, kms, kme, &
4121	its, ite, jts, jte, kts, kte )
4122
4123	! DO j = j_start, j_end
4124	! DO k = kts+1, ktf
4125	! DO i = i_start, i_end
4126	! titau3(i,k,j)=titau3(i,k,j)*zeta_z(i,j)
4127	! ENDDO
4128	! ENDDO
4129	! ENDDO
4130
4131	DO j = j_start, j_end
4132	DO k = kts+1, ktf
4133	DO i = i_start, i_end
4134	tendency(i,k,j)=tendency(i,k,j)-rdz(i,k,j)*(titau3(i,k,j)-titau3(i,k-1,j))
4135	ENDDO
4136	ENDDO
4137	ENDDO
4138
4139	END SUBROUTINE vertical_diffusion_w_2
4140
4141	!=======================================================================
4142	!=======================================================================
4143
4144	SUBROUTINE vertical_diffusion_s( tendency, config_flags, var, mu, xkhv, &
4145	dn, dnw, rdz, rdzw, fnm, fnp, &
4146	doing_tke, &
4147	ids, ide, jds, jde, kds, kde, &
4148	ims, ime, jms, jme, kms, kme, &
4149	its, ite, jts, jte, kts, kte )
4150
4151	!-----------------------------------------------------------------------
4152	! Begin declarations.
4153
4154	IMPLICIT NONE
4155
4156	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4157
4158	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4159	ims, ime, jms, jme, kms, kme, &
4160	its, ite, jts, jte, kts, kte
4161
4162	LOGICAL, INTENT(IN ) :: doing_tke
4163
4164	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
4165	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
4166	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
4167	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
4168
4169	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
4170
4171	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN) :: xkhv
4172
4173	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: mu
4174
4175	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
4176	INTENT(IN ) :: var, &
4177	rdz, &
4178	rdzw
4179	! LOCAL VARS
4180
4181	INTEGER :: i, j, k, ktf
4182
4183	INTEGER :: i_start, i_end, j_start, j_end
4184
4185	REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: H3, &
4186	xkxavg, &
4187	rravg
4188
4189	REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
4190
4191	! End declarations.
4192	!-----------------------------------------------------------------------
4193
4194	ktf=MIN(kte,kde-1)
4195
4196	i_start = its
4197	i_end = MIN(ite,ide-1)
4198	j_start = jts
4199	j_end = MIN(jte,jde-1)
4200
4201	IF ( config_flags%open_xs .or. config_flags%specified .or. &
4202	config_flags%nested) i_start = MAX(ids+1,its)
4203	IF ( config_flags%open_xe .or. config_flags%specified .or. &
4204	config_flags%nested) i_end = MIN(ide-2,ite)
4205	IF ( config_flags%open_ys .or. config_flags%specified .or. &
4206	config_flags%nested) j_start = MAX(jds+1,jts)
4207	IF ( config_flags%open_ye .or. config_flags%specified .or. &
4208	config_flags%nested) j_end = MIN(jde-2,jte)
4209	IF ( config_flags%periodic_x ) i_start = its
4210	IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
4211
4212	IF (doing_tke) THEN
4213	DO j = j_start, j_end
4214	DO k = kts,ktf
4215	DO i = i_start, i_end
4216	tmptendf(i,k,j)=tendency(i,k,j)
4217	ENDDO
4218	ENDDO
4219	ENDDO
4220	ENDIF
4221
4222	! H3
4223
4224	xkxavg = 0.
4225
4226	DO j = j_start, j_end
4227	DO k = kts+1,ktf
4228	DO i = i_start, i_end
4229	xkxavg(i,k,j)=fnm(k)xkhv(i,k,j)+fnp(k)xkhv(i,k-1,j)
4230	H3(i,k,j)=-xkxavg(i,k,j)(var(i,k,j)-var(i,k-1,j))rdz(i,k,j)
4231	! H3(i,k,j)=-xkxavg(i,k,j)zeta_z(i,j) &
4232	! (var(i,k,j)-var(i,k-1,j))/dn(k)
4233	ENDDO
4234	ENDDO
4235	ENDDO
4236
4237	DO j = j_start, j_end
4238	DO i = i_start, i_end
4239	H3(i,kts,j)=0.
4240	H3(i,ktf+1,j)=0.
4241	! H3(i,kts,j)=H3(i,kts+1,j)
4242	! H3(i,ktf+1,j)=H3(i,ktf,j)
4243	ENDDO
4244	ENDDO
4245
4246	DO j = j_start, j_end
4247	DO k = kts,ktf
4248	DO i = i_start, i_end
4249	tendency(i,k,j)=tendency(i,k,j) &
4250	-mu(i,j)(H3(i,k+1,j)-H3(i,k,j))rdzw(i,k,j)
4251	ENDDO
4252	ENDDO
4253	ENDDO
4254
4255	IF (doing_tke) THEN
4256	DO j = j_start, j_end
4257	DO k = kts,ktf
4258	DO i = i_start, i_end
4259	tendency(i,k,j)=tmptendf(i,k,j)+2.* &
4260	(tendency(i,k,j)-tmptendf(i,k,j))
4261	ENDDO
4262	ENDDO
4263	ENDDO
4264	ENDIF
4265
4266	END SUBROUTINE vertical_diffusion_s
4267
4268	!=======================================================================
4269	!=======================================================================
4270
4271	SUBROUTINE cal_titau_11_22_33( config_flags, titau, &
4272	mu, tke, xkx, defor, &
4273	mtau, & !JDM
4274	is_ext, ie_ext, js_ext, je_ext, &
4275	ids, ide, jds, jde, kds, kde, &
4276	ims, ime, jms, jme, kms, kme, &
4277	its, ite, jts, jte, kts, kte )
4278
4279	! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
4280	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4281	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4282
4283	! Purpose: This routine calculates stress terms (taus) for use in
4284	! the calculation of production of TKE by sheared wind
4285
4286	! References: Klemp and Wilhelmson (JAS 1978)
4287	! Deardorff (B-L Meteor 1980)
4288	! Chen and Dudhia (NCAR WRF physics report 2000)
4289
4290	! Key:
4291
4292	!-----------------------------------------------------------------------
4293	! Begin declarations.
4294
4295	IMPLICIT NONE
4296
4297	TYPE( grid_config_rec_type ), INTENT( IN ) &
4298	:: config_flags
4299
4300	INTEGER, INTENT( IN ) &
4301	:: ids, ide, jds, jde, kds, kde, &
4302	ims, ime, jms, jme, kms, kme, &
4303	its, ite, jts, jte, kts, kte
4304
4305	INTEGER, INTENT( IN ) &
4306	:: is_ext, ie_ext, js_ext, je_ext
4307
4308	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
4309	:: titau
4310
4311	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4312	:: defor, xkx, tke
4313
4314	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & !JDM
4315	:: mtau
4316
4317	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4318	:: mu
4319
4320	! Local variables.
4321
4322	INTEGER &
4323	:: i, j, k, ktf, i_start, i_end, j_start, j_end
4324
4325	! End declarations.
4326	!-----------------------------------------------------------------------
4327
4328	ktf = MIN( kte, kde-1 )
4329
4330	i_start = its
4331	i_end = ite
4332	j_start = jts
4333	j_end = jte
4334
4335	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4336	config_flags%nested) i_start = MAX( ids+1, its )
4337	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4338	config_flags%nested) i_end = MIN( ide-1, ite )
4339	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4340	config_flags%nested) j_start = MAX( jds+1, jts )
4341	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4342	config_flags%nested) j_end = MIN( jde-1, jte )
4343	IF ( config_flags%periodic_x ) i_start = its
4344	IF ( config_flags%periodic_x ) i_end = ite
4345
4346	i_start = i_start - is_ext
4347	i_end = i_end + ie_ext
4348	j_start = j_start - js_ext
4349	j_end = j_end + je_ext
4350
4351	IF ( config_flags%sfs_opt .GT. 0 ) THEN ! USE NBA MODEL SFS STRESSES
4352
4353	DO j = j_start, j_end
4354	DO k = kts, ktf
4355	DO i = i_start, i_end
4356
4357	titau(i,k,j) = mu(i,j) * mtau(i,k,j)
4358
4359	END DO
4360	END DO
4361	END DO
4362
4363	ELSE !NOT NBA
4364
4365	IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT
4366
4367	DO j = j_start, j_end
4368	DO k = kts, ktf
4369	DO i = i_start, i_end
4370
4371	titau(i,k,j) = - mu(i,j) * xkx(i,k,j) * defor(i,k,j)
4372	mtau(i,k,j) = - xkx(i,k,j) * defor(i,k,j)
4373
4374	END DO
4375	END DO
4376	END DO
4377
4378	ELSE !NO STRESS OUTPUT
4379
4380	DO j = j_start, j_end
4381	DO k = kts, ktf
4382	DO i = i_start, i_end
4383
4384	titau(i,k,j) = - mu(i,j) * xkx(i,k,j) * defor(i,k,j)
4385
4386	END DO
4387	END DO
4388	END DO
4389
4390	ENDIF
4391
4392	ENDIF
4393
4394	END SUBROUTINE cal_titau_11_22_33
4395
4396	!=======================================================================
4397	!=======================================================================
4398
4399	SUBROUTINE cal_titau_12_21( config_flags, titau, &
4400	mu, xkx, defor, &
4401	mtau, & !JDM
4402	is_ext, ie_ext, js_ext, je_ext, &
4403	ids, ide, jds, jde, kds, kde, &
4404	ims, ime, jms, jme, kms, kme, &
4405	its, ite, jts, jte, kts, kte )
4406
4407	! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR
4408	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4409	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4410
4411	! Pusrpose This routine calculates the stress terms (taus) for use in
4412	! the calculation of production of TKE by sheared wind
4413
4414	! References: Klemp and Wilhelmson (JAS 1978)
4415	! Deardorff (B-L Meteor 1980)
4416	! Chen and Dudhia (NCAR WRF physics report 2000)
4417
4418	! Key:
4419
4420	!-----------------------------------------------------------------------
4421	! Begin declarations.
4422
4423	IMPLICIT NONE
4424
4425	TYPE( grid_config_rec_type), INTENT( IN ) &
4426	:: config_flags
4427
4428	INTEGER, INTENT( IN ) &
4429	:: ids, ide, jds, jde, kds, kde, &
4430	ims, ime, jms, jme, kms, kme, &
4431	its, ite, jts, jte, kts, kte
4432
4433	INTEGER, INTENT( IN ) &
4434	:: is_ext, ie_ext, js_ext, je_ext
4435
4436	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
4437	:: titau
4438
4439	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4440	:: defor, xkx
4441
4442	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & !JDM
4443	:: mtau
4444
4445	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4446	:: mu
4447
4448	! Local variables.
4449
4450	INTEGER &
4451	:: i, j, k, ktf, i_start, i_end, j_start, j_end
4452
4453	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
4454	:: xkxavg
4455
4456	REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
4457	:: muavg
4458
4459	! End declarations.
4460	!-----------------------------------------------------------------------
4461
4462	ktf = MIN( kte, kde-1 )
4463
4464	! Needs one more point in the x and y directions.
4465
4466	i_start = its
4467	i_end = ite
4468	j_start = jts
4469	j_end = jte
4470
4471	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4472	config_flags%nested ) i_start = MAX( ids+1, its )
4473	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4474	config_flags%nested ) i_end = MIN( ide-1, ite )
4475	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4476	config_flags%nested ) j_start = MAX( jds+1, jts )
4477	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4478	config_flags%nested ) j_end = MIN( jde-1, jte )
4479	IF ( config_flags%periodic_x ) i_start = its
4480	IF ( config_flags%periodic_x ) i_end = ite
4481
4482	i_start = i_start - is_ext
4483	i_end = i_end + ie_ext
4484	j_start = j_start - js_ext
4485	j_end = j_end + je_ext
4486
4487	DO j = j_start, j_end
4488	DO k = kts, ktf
4489	DO i = i_start, i_end
4490	xkxavg(i,k,j) = 0.25 * ( xkx(i-1,k,j ) + xkx(i,k,j ) + &
4491	xkx(i-1,k,j-1) + xkx(i,k,j-1) )
4492	END DO
4493	END DO
4494	END DO
4495
4496	DO j = j_start, j_end
4497	DO i = i_start, i_end
4498	muavg(i,j) = 0.25 * ( mu(i-1,j ) + mu(i,j ) + &
4499	mu(i-1,j-1) + mu(i,j-1) )
4500	END DO
4501	END DO
4502
4503	! titau12 or titau21
4504
4505	IF ( config_flags%sfs_opt .GT. 0 ) THEN ! USE NBA MODEL SFS STRESSES
4506
4507	DO j = j_start, j_end
4508	DO k = kts, ktf
4509	DO i = i_start, i_end
4510
4511	titau(i,k,j) = muavg(i,j) * mtau(i,k,j)
4512
4513	END DO
4514	END DO
4515	END DO
4516
4517	ELSE ! NOT NBA
4518
4519	IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT
4520
4521	DO j = j_start, j_end
4522	DO k = kts, ktf
4523	DO i = i_start, i_end
4524
4525	titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4526	mtau(i,k,j) = - xkxavg(i,k,j) * defor(i,k,j)
4527
4528	END DO
4529	END DO
4530	END DO
4531
4532	ELSE ! NO STRESS OUTPUT
4533
4534	DO j = j_start, j_end
4535	DO k = kts, ktf
4536	DO i = i_start, i_end
4537
4538	titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4539
4540	END DO
4541	END DO
4542	END DO
4543
4544	ENDIF
4545
4546	ENDIF
4547
4548	END SUBROUTINE cal_titau_12_21
4549
4550	!=======================================================================
4551
4552	SUBROUTINE cal_titau_13_31( config_flags, titau, &
4553	defor, &
4554	mtau, & !JDM
4555	mu, xkx, fnm, fnp, &
4556	is_ext, ie_ext, js_ext, je_ext, &
4557	ids, ide, jds, jde, kds, kde, &
4558	ims, ime, jms, jme, kms, kme, &
4559	its, ite, jts, jte, kts, kte )
4560
4561	! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR
4562	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4563	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4564
4565	! Purpose: This routine calculates the stress terms (taus) for use in
4566	! the calculation of production of TKE by sheared wind
4567
4568	! References: Klemp and Wilhelmson (JAS 1978)
4569	! Deardorff (B-L Meteor 1980)
4570	! Chen and Dudhia (NCAR WRF physics report 2000)
4571
4572	! Key:
4573
4574	!-----------------------------------------------------------------------
4575	! Begin declarations.
4576
4577	IMPLICIT NONE
4578
4579	TYPE( grid_config_rec_type), INTENT( IN ) &
4580	:: config_flags
4581
4582	INTEGER, INTENT( IN ) &
4583	:: ids, ide, jds, jde, kds, kde, &
4584	ims, ime, jms, jme, kms, kme, &
4585	its, ite, jts, jte, kts, kte
4586
4587	INTEGER, INTENT( IN ) &
4588	:: is_ext, ie_ext, js_ext, je_ext
4589
4590	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
4591	:: fnm, fnp
4592
4593	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
4594	:: titau
4595
4596	REAL, DIMENSION( ims:ime, kms:kme, jms:jme), INTENT( IN ) &
4597	:: defor, xkx
4598
4599	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & !JDM
4600	:: mtau
4601
4602	REAL, DIMENSION( ims:ime, jms:jme), INTENT( IN ) &
4603	:: mu
4604
4605	! Local variables.
4606
4607	INTEGER &
4608	:: i, j, k, ktf, i_start, i_end, j_start, j_end
4609
4610	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
4611	:: xkxavg
4612
4613	REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
4614	:: muavg
4615
4616	! End declarations.
4617	!-----------------------------------------------------------------------
4618
4619	ktf = MIN( kte, kde-1 )
4620
4621	! Find ide-1 and jde-1 for averaging to p point.
4622
4623	i_start = its
4624	i_end = ite
4625	j_start = jts
4626	j_end = MIN( jte, jde-1 )
4627
4628	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4629	config_flags%nested) i_start = MAX( ids+1, its )
4630	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4631	config_flags%nested) i_end = MIN( ide-1, ite )
4632	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4633	config_flags%nested) j_start = MAX( jds+1, jts )
4634	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4635	config_flags%nested) j_end = MIN( jde-2, jte )
4636	IF ( config_flags%periodic_x ) i_start = its
4637	IF ( config_flags%periodic_x ) i_end = ite
4638
4639	i_start = i_start - is_ext
4640	i_end = i_end + ie_ext
4641	j_start = j_start - js_ext
4642	j_end = j_end + je_ext
4643
4644	DO j = j_start, j_end
4645	DO k = kts+1, ktf
4646	DO i = i_start, i_end
4647	xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i-1,k ,j) ) + &
4648	fnp(k) * ( xkx(i,k-1,j) + xkx(i-1,k-1,j) ) )
4649	END DO
4650	END DO
4651	END DO
4652
4653	DO j = j_start, j_end
4654	DO i = i_start, i_end
4655	muavg(i,j) = 0.5 * ( mu(i,j) + mu(i-1,j) )
4656	END DO
4657	END DO
4658
4659	IF ( config_flags%sfs_opt .GT. 0 ) THEN ! USE NBA MODEL SFS STRESSES
4660
4661	DO j = j_start, j_end
4662	DO k = kts+1, ktf
4663	DO i = i_start, i_end
4664	titau(i,k,j) = muavg(i,j) * mtau(i,k,j)
4665	ENDDO
4666	ENDDO
4667	ENDDO
4668
4669	ELSE ! NOT NBA
4670
4671	IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT
4672
4673	DO j = j_start, j_end
4674	DO k = kts+1, ktf
4675	DO i = i_start, i_end
4676
4677	titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4678	mtau(i,k,j) = - xkxavg(i,k,j) * defor(i,k,j)
4679
4680	ENDDO
4681	ENDDO
4682	ENDDO
4683
4684	ELSE ! NO STRESS OUTPUT
4685
4686	DO j = j_start, j_end
4687	DO k = kts+1, ktf
4688	DO i = i_start, i_end
4689
4690	titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4691
4692	ENDDO
4693	ENDDO
4694	ENDDO
4695
4696	ENDIF
4697
4698	ENDIF
4699
4700	DO j = j_start, j_end
4701	DO i = i_start, i_end
4702	titau(i,kts ,j) = 0.0
4703	titau(i,ktf+1,j) = 0.0
4704	ENDDO
4705	ENDDO
4706
4707	END SUBROUTINE cal_titau_13_31
4708
4709	!=======================================================================
4710	!=======================================================================
4711
4712	SUBROUTINE cal_titau_23_32( config_flags, titau, defor, &
4713	mtau, & !JDM
4714	mu, xkx, fnm, fnp, &
4715	is_ext, ie_ext, js_ext, je_ext, &
4716	ids, ide, jds, jde, kds, kde, &
4717	ims, ime, jms, jme, kms, kme, &
4718	its, ite, jts, jte, kts, kte )
4719
4720	! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
4721	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
4722	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
4723
4724	! Purpose: This routine calculates stress terms (taus) for use in
4725	! the calculation of production of TKE by sheared wind
4726
4727	! References: Klemp and Wilhelmson (JAS 1978)
4728	! Deardorff (B-L Meteor 1980)
4729	! Chen and Dudhia (NCAR WRF physics report 2000)
4730
4731	! Key:
4732
4733	!-----------------------------------------------------------------------
4734	! Begin declarations.
4735
4736	IMPLICIT NONE
4737
4738	TYPE( grid_config_rec_type ), INTENT( IN ) &
4739	:: config_flags
4740
4741	INTEGER, INTENT( IN ) &
4742	:: ids, ide, jds, jde, kds, kde, &
4743	ims, ime, jms, jme, kms, kme, &
4744	its, ite, jts, jte, kts, kte
4745
4746	INTEGER, INTENT( IN ) &
4747	:: is_ext,ie_ext,js_ext,je_ext
4748
4749	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
4750	:: fnm, fnp
4751
4752	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) &
4753	:: titau
4754
4755	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
4756	:: defor, xkx
4757
4758	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & !JDM
4759	:: mtau
4760
4761	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
4762	:: mu
4763
4764	! Local variables.
4765
4766	INTEGER &
4767	:: i, j, k, ktf, i_start, i_end, j_start, j_end
4768
4769	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
4770	:: xkxavg
4771
4772	REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) &
4773	:: muavg
4774
4775	! End declarations.
4776	!-----------------------------------------------------------------------
4777
4778	ktf = MIN( kte, kde-1 )
4779
4780	! Find ide-1 and jde-1 for averaging to p point.
4781
4782	i_start = its
4783	i_end = MIN( ite, ide-1 )
4784	j_start = jts
4785	j_end = jte
4786
4787	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
4788	config_flags%nested) i_start = MAX( ids+1, its )
4789	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
4790	config_flags%nested) i_end = MIN( ide-2, ite )
4791	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
4792	config_flags%nested) j_start = MAX( jds+1, jts )
4793	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
4794	config_flags%nested) j_end = MIN( jde-1, jte )
4795	IF ( config_flags%periodic_x ) i_start = its
4796	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
4797
4798	i_start = i_start - is_ext
4799	i_end = i_end + ie_ext
4800	j_start = j_start - js_ext
4801	j_end = j_end + je_ext
4802
4803	DO j = j_start, j_end
4804	DO k = kts+1, ktf
4805	DO i = i_start, i_end
4806	xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i,k ,j-1) ) + &
4807	fnp(k) * ( xkx(i,k-1,j) + xkx(i,k-1,j-1) ) )
4808	END DO
4809	END DO
4810	END DO
4811
4812	DO j = j_start, j_end
4813	DO i = i_start, i_end
4814	muavg(i,j) = 0.5 * ( mu(i,j) + mu(i,j-1) )
4815	END DO
4816	END DO
4817
4818	IF ( config_flags%sfs_opt .EQ. 1 ) THEN ! USE NBA MODEL SFS STRESSES
4819
4820	DO j = j_start, j_end
4821	DO k = kts+1, ktf
4822	DO i = i_start, i_end
4823
4824	titau(i,k,j) = muavg(i,j) * mtau(i,k,j)
4825
4826	END DO
4827	END DO
4828	END DO
4829
4830	ELSE ! NOT NBA
4831
4832	IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT
4833
4834	DO j = j_start, j_end
4835	DO k = kts+1, ktf
4836	DO i = i_start, i_end
4837
4838	titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4839	mtau(i,k,j) = - xkxavg(i,k,j) * defor(i,k,j)
4840
4841	END DO
4842	END DO
4843	END DO
4844
4845	ELSE ! NO STRESS OUTPUT
4846
4847	DO j = j_start, j_end
4848	DO k = kts+1, ktf
4849	DO i = i_start, i_end
4850
4851	titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j)
4852
4853	END DO
4854	END DO
4855	END DO
4856
4857	ENDIF
4858
4859	ENDIF
4860
4861	DO j = j_start, j_end
4862	DO i = i_start, i_end
4863	titau(i,kts ,j) = 0.0
4864	titau(i,ktf+1,j) = 0.0
4865	END DO
4866	END DO
4867
4868	END SUBROUTINE cal_titau_23_32
4869
4870	!=======================================================================
4871	!=======================================================================
4872
4873	SUBROUTINE phy_bc ( config_flags,div,defor11,defor22,defor33, &
4874	defor12,defor13,defor23,xkmh,xkmv,xkhh,xkhv,tke, &
4875	RUBLTEN, RVBLTEN, &
4876	RUCUTEN, RVCUTEN, &
4877	RUSHTEN, RVSHTEN, &
4878	ids, ide, jds, jde, kds, kde, &
4879	ims, ime, jms, jme, kms, kme, &
4880	ips, ipe, jps, jpe, kps, kpe, &
4881	its, ite, jts, jte, kts, kte )
4882
4883	!------------------------------------------------------------------------------
4884	! Begin declarations.
4885
4886	IMPLICIT NONE
4887
4888	TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
4889
4890	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4891	ims, ime, jms, jme, kms, kme, &
4892	ips, ipe, jps, jpe, kps, kpe, &
4893	its, ite, jts, jte, kts, kte
4894
4895	REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::RUBLTEN, &
4896	RVBLTEN, &
4897	RUCUTEN, &
4898	RVCUTEN, &
4899	RUSHTEN, &
4900	RVSHTEN, &
4901	defor11, &
4902	defor22, &
4903	defor33, &
4904	defor12, &
4905	defor13, &
4906	defor23, &
4907	xkmh, &
4908	xkmv, &
4909	xkhh, &
4910	xkhv, &
4911	tke, &
4912	div
4913
4914	! End declarations.
4915	!-----------------------------------------------------------------------
4916
4917	IF(config_flags%bl_pbl_physics .GT. 0) THEN
4918
4919	CALL set_physical_bc3d( RUBLTEN , 't', config_flags, &
4920	ids, ide, jds, jde, kds, kde, &
4921	ims, ime, jms, jme, kms, kme, &
4922	ips, ipe, jps, jpe, kps, kpe, &
4923	its, ite, jts, jte, kts, kte )
4924
4925	CALL set_physical_bc3d( RVBLTEN , 't', config_flags, &
4926	ids, ide, jds, jde, kds, kde, &
4927	ims, ime, jms, jme, kms, kme, &
4928	ips, ipe, jps, jpe, kps, kpe, &
4929	its, ite, jts, jte, kts, kte )
4930
4931	ENDIF
4932
4933	IF(config_flags%cu_physics .GT. 0) THEN
4934
4935	CALL set_physical_bc3d( RUCUTEN , 't', config_flags, &
4936	ids, ide, jds, jde, kds, kde, &
4937	ims, ime, jms, jme, kms, kme, &
4938	ips, ipe, jps, jpe, kps, kpe, &
4939	its, ite, jts, jte, kts, kte )
4940
4941	CALL set_physical_bc3d( RVCUTEN , 't', config_flags, &
4942	ids, ide, jds, jde, kds, kde, &
4943	ims, ime, jms, jme, kms, kme, &
4944	ips, ipe, jps, jpe, kps, kpe, &
4945	its, ite, jts, jte, kts, kte )
4946
4947	ENDIF
4948
4949	IF(config_flags%shcu_physics .GT. 0) THEN
4950
4951	CALL set_physical_bc3d( RUSHTEN , 't', config_flags, &
4952	ids, ide, jds, jde, kds, kde, &
4953	ims, ime, jms, jme, kms, kme, &
4954	ips, ipe, jps, jpe, kps, kpe, &
4955	its, ite, jts, jte, kts, kte )
4956
4957	CALL set_physical_bc3d( RVSHTEN , 't', config_flags, &
4958	ids, ide, jds, jde, kds, kde, &
4959	ims, ime, jms, jme, kms, kme, &
4960	ips, ipe, jps, jpe, kps, kpe, &
4961	its, ite, jts, jte, kts, kte )
4962
4963	ENDIF
4964
4965	! move out of the conditional, below; horiz coeffs needed for
4966	! all diff_opt cases. JM
4967
4968	CALL set_physical_bc3d( xkmh , 't', config_flags, &
4969	ids, ide, jds, jde, kds, kde, &
4970	ims, ime, jms, jme, kms, kme, &
4971	ips, ipe, jps, jpe, kps, kpe, &
4972	its, ite, jts, jte, kts, kte )
4973
4974	CALL set_physical_bc3d( xkhh , 't', config_flags, &
4975	ids, ide, jds, jde, kds, kde, &
4976	ims, ime, jms, jme, kms, kme, &
4977	ips, ipe, jps, jpe, kps, kpe, &
4978	its, ite, jts, jte, kts, kte )
4979
4980	IF(config_flags%diff_opt .eq. 2) THEN
4981
4982	CALL set_physical_bc3d( xkmv , 't', config_flags, &
4983	ids, ide, jds, jde, kds, kde, &
4984	ims, ime, jms, jme, kms, kme, &
4985	ips, ipe, jps, jpe, kps, kpe, &
4986	its, ite, jts, jte, kts, kte )
4987
4988	CALL set_physical_bc3d( xkhv , 't', config_flags, &
4989	ids, ide, jds, jde, kds, kde, &
4990	ims, ime, jms, jme, kms, kme, &
4991	ips, ipe, jps, jpe, kps, kpe, &
4992	its, ite, jts, jte, kts, kte )
4993
4994	CALL set_physical_bc3d( div , 't', config_flags, &
4995	ids, ide, jds, jde, kds, kde, &
4996	ims, ime, jms, jme, kms, kme, &
4997	ips, ipe, jps, jpe, kps, kpe, &
4998	its, ite, jts, jte, kts, kte )
4999
5000	CALL set_physical_bc3d( defor11 , 't', config_flags, &
5001	ids, ide, jds, jde, kds, kde, &
5002	ims, ime, jms, jme, kms, kme, &
5003	ips, ipe, jps, jpe, kps, kpe, &
5004	its, ite, jts, jte, kts, kte )
5005
5006	CALL set_physical_bc3d( defor22 , 't', config_flags, &
5007	ids, ide, jds, jde, kds, kde, &
5008	ims, ime, jms, jme, kms, kme, &
5009	ips, ipe, jps, jpe, kps, kpe, &
5010	its, ite, jts, jte, kts, kte )
5011
5012	CALL set_physical_bc3d( defor33 , 't', config_flags, &
5013	ids, ide, jds, jde, kds, kde, &
5014	ims, ime, jms, jme, kms, kme, &
5015	ips, ipe, jps, jpe, kps, kpe, &
5016	its, ite, jts, jte, kts, kte )
5017
5018	CALL set_physical_bc3d( defor12 , 'd', config_flags, &
5019	ids, ide, jds, jde, kds, kde, &
5020	ims, ime, jms, jme, kms, kme, &
5021	ips, ipe, jps, jpe, kps, kpe, &
5022	its, ite, jts, jte, kts, kte )
5023
5024	CALL set_physical_bc3d( defor13 , 'e', config_flags, &
5025	ids, ide, jds, jde, kds, kde, &
5026	ims, ime, jms, jme, kms, kme, &
5027	ips, ipe, jps, jpe, kps, kpe, &
5028	its, ite, jts, jte, kts, kte )
5029
5030	CALL set_physical_bc3d( defor23 , 'f', config_flags, &
5031	ids, ide, jds, jde, kds, kde, &
5032	ims, ime, jms, jme, kms, kme, &
5033	ips, ipe, jps, jpe, kps, kpe, &
5034	its, ite, jts, jte, kts, kte )
5035
5036	ENDIF
5037
5038	END SUBROUTINE phy_bc
5039
5040	!=======================================================================
5041	!=======================================================================
5042
5043	SUBROUTINE tke_rhs( tendency, BN2, config_flags, &
5044	defor11, defor22, defor33, &
5045	defor12, defor13, defor23, &
5046	u, v, w, div, tke, mu, &
5047	theta, p, p8w, t8w, z, fnm, fnp, &
5048	cf1, cf2, cf3, msftx, msfty, &
5049	xkmh, xkmv, xkhv, &
5050	rdx, rdy, dx, dy, dt, zx, zy, &
5051	rdz, rdzw, dn, dnw, isotropic, &
5052	hfx, qfx, qv, ust, rho, &
5053	ids, ide, jds, jde, kds, kde, &
5054	ims, ime, jms, jme, kms, kme, &
5055	its, ite, jts, jte, kts, kte )
5056
5057	!-----------------------------------------------------------------------
5058	! Begin declarations.
5059
5060	IMPLICIT NONE
5061
5062	TYPE( grid_config_rec_type ), INTENT( IN ) &
5063	:: config_flags
5064
5065	INTEGER, INTENT( IN ) &
5066	:: ids, ide, jds, jde, kds, kde, &
5067	ims, ime, jms, jme, kms, kme, &
5068	its, ite, jts, jte, kts, kte
5069
5070	INTEGER, INTENT( IN ) :: isotropic
5071	REAL, INTENT( IN ) &
5072	:: cf1, cf2, cf3, dt, rdx, rdy, dx, dy
5073
5074	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
5075	:: fnm, fnp, dnw, dn
5076
5077	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5078	:: msftx, msfty
5079
5080	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5081	:: tendency
5082
5083	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5084	:: defor11, defor22, defor33, defor12, defor13, defor23, &
5085	div, BN2, tke, xkmh, xkmv, xkhv, zx, zy, u, v, w, theta, &
5086	p, p8w, t8w, z, rdz, rdzw
5087
5088	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5089	:: mu
5090
5091	REAL, DIMENSION ( ims:ime, jms:jme ), INTENT( IN ) &
5092	:: hfx, ust, qfx
5093	REAL, DIMENSION ( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) &
5094	:: qv, rho
5095
5096	! Local variables.
5097
5098	INTEGER &
5099	:: i, j, k, ktf, i_start, i_end, j_start, j_end
5100
5101	! End declarations.
5102	!-----------------------------------------------------------------------
5103
5104	CALL tke_shear( tendency, config_flags, &
5105	defor11, defor22, defor33, &
5106	defor12, defor13, defor23, &
5107	u, v, w, tke, ust, mu, fnm, fnp, &
5108	cf1, cf2, cf3, msftx, msfty, &
5109	xkmh, xkmv, &
5110	rdx, rdy, zx, zy, rdz, rdzw, dnw, dn, &
5111	ids, ide, jds, jde, kds, kde, &
5112	ims, ime, jms, jme, kms, kme, &
5113	its, ite, jts, jte, kts, kte )
5114
5115	CALL tke_buoyancy( tendency, config_flags, mu, &
5116	tke, xkhv, BN2, theta, dt, &
5117	hfx, qfx, qv, rho, &
5118	ids, ide, jds, jde, kds, kde, &
5119	ims, ime, jms, jme, kms, kme, &
5120	its, ite, jts, jte, kts, kte )
5121
5122	CALL tke_dissip( tendency, config_flags, &
5123	mu, tke, bn2, theta, p8w, t8w, z, &
5124	dx, dy,rdz, rdzw, isotropic, &
5125	msftx, msfty, &
5126	ids, ide, jds, jde, kds, kde, &
5127	ims, ime, jms, jme, kms, kme, &
5128	its, ite, jts, jte, kts, kte )
5129
5130	! Set a lower limit on TKE.
5131
5132	ktf = MIN( kte, kde-1 )
5133	i_start = its
5134	i_end = MIN( ite, ide-1 )
5135	j_start = jts
5136	j_end = MIN( jte, jde-1 )
5137
5138	IF ( config_flags%open_xs .or. config_flags%specified .or. &
5139	config_flags%nested) i_start = MAX(ids+1,its)
5140	IF ( config_flags%open_xe .or. config_flags%specified .or. &
5141	config_flags%nested) i_end = MIN(ide-2,ite)
5142	IF ( config_flags%open_ys .or. config_flags%specified .or. &
5143	config_flags%nested) j_start = MAX(jds+1,jts)
5144	IF ( config_flags%open_ye .or. config_flags%specified .or. &
5145	config_flags%nested) j_end = MIN(jde-2,jte)
5146	IF ( config_flags%periodic_x ) i_start = its
5147	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
5148
5149	DO j = j_start, j_end
5150	DO k = kts, ktf
5151	DO i = i_start, i_end
5152	tendency(i,k,j) = max( tendency(i,k,j), -mu(i,j) * max( 0.0 , tke(i,k,j) ) / dt )
5153	END DO
5154	END DO
5155	END DO
5156
5157	END SUBROUTINE tke_rhs
5158
5159	!=======================================================================
5160	!=======================================================================
5161
5162	SUBROUTINE tke_buoyancy( tendency, config_flags, mu, &
5163	tke, xkhv, BN2, theta, dt, &
5164	hfx, qfx, qv, rho, &
5165	ids, ide, jds, jde, kds, kde, &
5166	ims, ime, jms, jme, kms, kme, &
5167	its, ite, jts, jte, kts, kte )
5168
5169	!-----------------------------------------------------------------------
5170	! Begin declarations.
5171
5172	IMPLICIT NONE
5173
5174	TYPE( grid_config_rec_type ), INTENT( IN ) &
5175	:: config_flags
5176
5177	INTEGER, INTENT( IN ) &
5178	:: ids, ide, jds, jde, kds, kde, &
5179	ims, ime, jms, jme, kms, kme, &
5180	its, ite, jts, jte, kts, kte
5181
5182	REAL, INTENT( IN ) &
5183	:: dt
5184
5185	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5186	:: tendency
5187
5188	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5189	:: xkhv, tke, BN2, theta
5190
5191	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5192	:: mu
5193
5194	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) &
5195	:: qv, rho
5196
5197	REAL, DIMENSION(ims:ime, jms:jme ), INTENT ( IN ) :: hfx, qfx
5198
5199	! Local variables.
5200
5201	INTEGER &
5202	:: i, j, k, ktf
5203
5204	INTEGER &
5205	:: i_start, i_end, j_start, j_end
5206
5207	REAL :: heat_flux, heat_flux0
5208
5209	REAL :: cpm
5210
5211	! End declarations.
5212	!-----------------------------------------------------------------------
5213
5214	!-----------------------------------------------------------------------
5215	! Add to the TKE tendency the term that accounts for production of TKE
5216	! due to buoyant motions.
5217
5218	ktf = MIN( kte, kde-1 )
5219	i_start = its
5220	i_end = MIN( ite, ide-1 )
5221	j_start = jts
5222	j_end = MIN( jte, jde-1 )
5223
5224	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
5225	config_flags%nested ) i_start = MAX( ids+1, its )
5226	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
5227	config_flags%nested ) i_end = MIN( ide-2, ite )
5228	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
5229	config_flags%nested ) j_start = MAX( jds+1, jts )
5230	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
5231	config_flags%nested ) j_end = MIN( jde-2, jte )
5232	IF ( config_flags%periodic_x ) i_start = its
5233	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
5234
5235	DO j = j_start, j_end
5236	DO k = kts+1, ktf
5237	DO i = i_start, i_end
5238	tendency(i,k,j) = tendency(i,k,j) - mu(i,j) * xkhv(i,k,j) * BN2(i,k,j)
5239	END DO
5240	END DO
5241	END DO
5242
5243	! MARTA: change in the computation of the tke's tendency at the surface.
5244	! the buoyancy flux is the average of the surface heat flux (0.06) and the
5245	! flux at the first w level
5246	!
5247	! WCS 040331
5248
5249	hflux: SELECT CASE( config_flags%isfflx )
5250	CASE (0,2) ! with fixed surface heat flux given in the namelist
5251	heat_flux0 = config_flags%tke_heat_flux ! constant heat flux value
5252	! set in namelist.input
5253	! LES mods
5254	K=KTS
5255	DO j = j_start, j_end
5256	DO i = i_start, i_end
5257	heat_flux = heat_flux0
5258	tendency(i,k,j)= tendency(i,k,j) - &
5259	mu(i,j)((xkhv(i,k,j)BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2.
5260
5261	ENDDO
5262	ENDDO
5263
5264	CASE (1) ! use surface heat flux computed from surface routine
5265	K=KTS
5266	DO j = j_start, j_end
5267	DO i = i_start, i_end
5268	cpm = cp * (1. + 0.8*qv(i,k,j))
5269	heat_flux = (hfx(i,j)/cpm)/rho(i,k,j)
5270
5271	tendency(i,k,j)= tendency(i,k,j) - &
5272	mu(i,j)((xkhv(i,k,j)BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2.
5273
5274	ENDDO
5275	ENDDO
5276
5277	CASE DEFAULT
5278	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
5279	END SELECT hflux
5280
5281	! end of MARTA/WCS change
5282
5283	! The tendency array now includes production of TKE from buoyant
5284	! motions.
5285	!-----------------------------------------------------------------------
5286
5287	END SUBROUTINE tke_buoyancy
5288
5289	!=======================================================================
5290	!=======================================================================
5291
5292	SUBROUTINE tke_dissip( tendency, config_flags, &
5293	mu, tke, bn2, theta, p8w, t8w, z, &
5294	dx, dy, rdz, rdzw, isotropic, &
5295	msftx, msfty, &
5296	ids, ide, jds, jde, kds, kde, &
5297	ims, ime, jms, jme, kms, kme, &
5298	its, ite, jts, jte, kts, kte )
5299
5300	! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR
5301	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
5302	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
5303
5304	! Purpose: This routine calculates dissipation of turbulent kinetic
5305	! energy.
5306
5307	! References: Klemp and Wilhelmson (JAS 1978)
5308	! Deardorff (B-L Meteor 1980)
5309	! Chen and Dudhia (NCAR WRF physics report 2000)
5310
5311	!-----------------------------------------------------------------------
5312	! Begin declarations.
5313
5314	IMPLICIT NONE
5315
5316	TYPE( grid_config_rec_type ), INTENT( IN ) &
5317	:: config_flags
5318
5319	INTEGER, INTENT( IN ) &
5320	:: ids, ide, jds, jde, kds, kde, &
5321	ims, ime, jms, jme, kms, kme, &
5322	its, ite, jts, jte, kts, kte
5323
5324	INTEGER, INTENT( IN ) :: isotropic
5325	REAL, INTENT( IN ) &
5326	:: dx, dy
5327
5328	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5329	:: tendency
5330
5331	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5332	:: tke, bn2, theta, p8w, t8w, z, rdz, rdzw
5333
5334	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5335	:: mu
5336
5337	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5338	:: msftx, msfty
5339	! Local variables.
5340
5341	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
5342	:: dthrdn
5343
5344	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
5345	:: l_scale
5346
5347	REAL, DIMENSION( its:ite ) &
5348	:: sumtke, sumtkez
5349
5350	INTEGER &
5351	:: i, j, k, ktf, i_start, i_end, j_start, j_end
5352
5353	REAL &
5354	:: disp_len, deltas, coefc, tmpdz, len_s, thetasfc, &
5355	thetatop, len_0, tketmp, tmp, ce1, ce2, c_k
5356
5357	! End declarations.
5358	!-----------------------------------------------------------------------
5359	c_k = config_flags%c_k
5360
5361	ce1 = ( c_k / 0.10 ) * 0.19
5362	ce2 = max( 0.0 , 0.93 - ce1 )
5363
5364	ktf = MIN( kte, kde-1 )
5365	i_start = its
5366	i_end = MIN(ite,ide-1)
5367	j_start = jts
5368	j_end = MIN(jte,jde-1)
5369
5370	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
5371	config_flags%nested) i_start = MAX( ids+1, its )
5372	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
5373	config_flags%nested) i_end = MIN( ide-2, ite )
5374	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
5375	config_flags%nested) j_start = MAX( jds+1, jts )
5376	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
5377	config_flags%nested) j_end = MIN( jde-2, jte )
5378	IF ( config_flags%periodic_x ) i_start = its
5379	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
5380
5381	CALL calc_l_scale( config_flags, tke, BN2, l_scale, &
5382	i_start, i_end, ktf, j_start, j_end, &
5383	dx, dy, rdzw, msftx, msfty, &
5384	ids, ide, jds, jde, kds, kde, &
5385	ims, ime, jms, jme, kms, kme, &
5386	its, ite, jts, jte, kts, kte )
5387	DO j = j_start, j_end
5388	DO k = kts, ktf
5389	DO i = i_start, i_end
5390	deltas = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333
5391	tketmp = MAX( tke(i,k,j), 1.0e-6 )
5392
5393	! Apply Deardorff's (1980) "wall effect" at the bottom of the domain.
5394	! For LES with fine grid, no need for this wall effect!
5395
5396	IF ( k .eq. kts .or. k .eq. ktf ) then
5397	coefc = 3.9
5398	ELSE
5399	coefc = ce1 + ce2 * l_scale(i,k,j) / deltas
5400	END IF
5401
5402	tendency(i,k,j) = tendency(i,k,j) - &
5403	mu(i,j) * coefc * tketmp**1.5 / l_scale(i,k,j)
5404	END DO
5405	END DO
5406	END DO
5407
5408	END SUBROUTINE tke_dissip
5409
5410	!=======================================================================
5411	!=======================================================================
5412
5413	SUBROUTINE tke_shear( tendency, config_flags, &
5414	defor11, defor22, defor33, &
5415	defor12, defor13, defor23, &
5416	u, v, w, tke, ust, mu, fnm, fnp, &
5417	cf1, cf2, cf3, msftx, msfty, &
5418	xkmh, xkmv, &
5419	rdx, rdy, zx, zy, rdz, rdzw, dn, dnw, &
5420	ids, ide, jds, jde, kds, kde, &
5421	ims, ime, jms, jme, kms, kme, &
5422	its, ite, jts, jte, kts, kte )
5423
5424	! History: Sep 2003 Rewritten by George Bryan and Jason Knievel,
5425	! NCAR
5426	! Oct 2001 Converted to mass core by Bill Skamarock, NCAR
5427	! Aug 2000 Original code by Shu-Hua Chen, UC-Davis
5428
5429	! Purpose: This routine calculates the production of turbulent
5430	! kinetic energy by stresses due to sheared wind.
5431
5432	! References: Klemp and Wilhelmson (JAS 1978)
5433	! Deardorff (B-L Meteor 1980)
5434	! Chen and Dudhia (NCAR WRF physics report 2000)
5435
5436	! Key:
5437
5438	! avg temporary working array
5439	! cf1
5440	! cf2
5441	! cf3
5442	! defor11 deformation term ( du/dx + du/dx )
5443	! defor12 deformation term ( dv/dx + du/dy ); same as defor21
5444	! defor13 deformation term ( dw/dx + du/dz ); same as defor31
5445	! defor22 deformation term ( dv/dy + dv/dy )
5446	! defor23 deformation term ( dw/dy + dv/dz ); same as defor32
5447	! defor33 deformation term ( dw/dz + dw/dz )
5448	! div 3-d divergence
5449	! dn
5450	! dnw
5451	! fnm
5452	! fnp
5453	! msftx
5454	! msfty
5455	! rdx
5456	! rdy
5457	! tendency
5458	! titau tau (stress tensor) with a tilde, indicating division by
5459	! a map-scale factor and the fraction of the total modeled
5460	! atmosphere beneath a given altitude (titau = tau/m/zeta)
5461	! tke turbulent kinetic energy
5462
5463	!-----------------------------------------------------------------------
5464	! Begin declarations.
5465
5466	IMPLICIT NONE
5467
5468	TYPE( grid_config_rec_type ), INTENT( IN ) &
5469	:: config_flags
5470
5471	INTEGER, INTENT( IN ) &
5472	:: ids, ide, jds, jde, kds, kde, &
5473	ims, ime, jms, jme, kms, kme, &
5474	its, ite, jts, jte, kts, kte
5475
5476	REAL, INTENT( IN ) &
5477	:: cf1, cf2, cf3, rdx, rdy
5478
5479	REAL, DIMENSION( kms:kme ), INTENT( IN ) &
5480	:: fnm, fnp, dn, dnw
5481
5482	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5483	:: msftx, msfty
5484
5485	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5486	:: tendency
5487
5488	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5489	:: defor11, defor22, defor33, defor12, defor13, defor23, &
5490	tke, xkmh, xkmv, zx, zy, u, v, w, rdz, rdzw
5491
5492	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5493	:: mu
5494
5495	REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5496	:: ust
5497
5498	! Local variables.
5499
5500	INTEGER &
5501	:: i, j, k, ktf, ktes1, ktes2, &
5502	i_start, i_end, j_start, j_end, &
5503	is_ext, ie_ext, js_ext, je_ext
5504
5505	REAL &
5506	:: mtau
5507
5508	REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) &
5509	:: avg, titau, tmp2
5510
5511	REAL, DIMENSION( its:ite, kts:kte, jts:jte ) &
5512	:: titau12, tmp1, zxavg, zyavg
5513
5514	REAL :: absU, cd0, Cd
5515
5516	! End declarations.
5517	!-----------------------------------------------------------------------
5518
5519	ktf = MIN( kte, kde-1 )
5520	ktes1 = kte-1
5521	ktes2 = kte-2
5522
5523	i_start = its
5524	i_end = MIN( ite, ide-1 )
5525	j_start = jts
5526	j_end = MIN( jte, jde-1 )
5527
5528	IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
5529	config_flags%nested ) i_start = MAX( ids+1, its )
5530	IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
5531	config_flags%nested ) i_end = MIN( ide-2, ite )
5532	IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
5533	config_flags%nested ) j_start = MAX( jds+1, jts )
5534	IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
5535	config_flags%nested ) j_end = MIN( jde-2, jte )
5536	IF ( config_flags%periodic_x ) i_start = its
5537	IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
5538
5539	DO j = j_start, j_end
5540	DO k = kts, ktf
5541	DO i = i_start, i_end
5542	zxavg(i,k,j) = 0.25 * ( zx(i,k ,j) + zx(i+1,k ,j) + &
5543	zx(i,k+1,j) + zx(i+1,k+1,j) )
5544	zyavg(i,k,j) = 0.25 * ( zy(i,k ,j) + zy(i,k ,j+1) + &
5545	zy(i,k+1,j) + zy(i,k+1,j+1) )
5546	END DO
5547	END DO
5548	END DO
5549
5550	! Begin calculating production of turbulence due to shear. The approach
5551	! is to add together contributions from six terms, each of which is the
5552	! square of a deformation that is then multiplied by an exchange
5553	! coefficiant. The same exchange coefficient is assumed for horizontal
5554	! and vertical coefficients for some of the terms (the vertical value is
5555	! the one used).
5556
5557	! For defor11.
5558
5559	DO j = j_start, j_end
5560	DO k = kts, ktf
5561	DO i = i_start, i_end
5562	tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
5563	mu(i,j) * xkmh(i,k,j) * ( ( defor11(i,k,j) )**2 )
5564	END DO
5565	END DO
5566	END DO
5567
5568	! For defor22.
5569
5570	DO j = j_start, j_end
5571	DO k = kts, ktf
5572	DO i = i_start, i_end
5573	tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
5574	mu(i,j) * xkmh(i,k,j) * ( ( defor22(i,k,j) )**2 )
5575	END DO
5576	END DO
5577	END DO
5578
5579	! For defor33.
5580
5581	DO j = j_start, j_end
5582	DO k = kts, ktf
5583	DO i = i_start, i_end
5584	tendency(i,k,j) = tendency(i,k,j) + 0.5 * &
5585	mu(i,j) * xkmv(i,k,j) * ( ( defor33(i,k,j) )**2 )
5586	END DO
5587	END DO
5588	END DO
5589
5590	! For defor12.
5591
5592	DO j = j_start, j_end
5593	DO k = kts, ktf
5594	DO i = i_start, i_end
5595	avg(i,k,j) = 0.25 * &
5596	( ( defor12(i ,k,j)2 ) + ( defor12(i ,k,j+1)2 ) + &
5597	( defor12(i+1,k,j)2 ) + ( defor12(i+1,k,j+1)2 ) )
5598	END DO
5599	END DO
5600	END DO
5601
5602	DO j = j_start, j_end
5603	DO k = kts, ktf
5604	DO i = i_start, i_end
5605	tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmh(i,k,j) * avg(i,k,j)
5606	END DO
5607	END DO
5608	END DO
5609
5610	! For defor13.
5611
5612	DO j = j_start, j_end
5613	DO k = kts+1, ktf
5614	DO i = i_start, i_end+1
5615	tmp2(i,k,j) = defor13(i,k,j)
5616	END DO
5617	END DO
5618	END DO
5619
5620	DO j = j_start, j_end
5621	DO i = i_start, i_end+1
5622	tmp2(i,kts ,j) = 0.0
5623	tmp2(i,ktf+1,j) = 0.0
5624	END DO
5625	END DO
5626
5627	DO j = j_start, j_end
5628	DO k = kts, ktf
5629	DO i = i_start, i_end
5630	avg(i,k,j) = 0.25 * &
5631	( ( tmp2(i ,k+1,j)2 ) + ( tmp2(i ,k,j)2 ) + &
5632	( tmp2(i+1,k+1,j)2 ) + ( tmp2(i+1,k,j)2 ) )
5633	END DO
5634	END DO
5635	END DO
5636
5637	DO j = j_start, j_end
5638	DO k = kts, ktf
5639	DO i = i_start, i_end
5640	tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
5641	END DO
5642	END DO
5643	END DO
5644
5645	!MARTA: add the drag at the surface; WCS 040331
5646	K=KTS
5647
5648	uflux: SELECT CASE( config_flags%isfflx )
5649	CASE (0) ! Assume cd a constant, specified in namelist
5650
5651	cd0 = config_flags%tke_drag_coefficient ! drag coefficient set
5652	! in namelist.input
5653	DO j = j_start, j_end
5654	DO i = i_start, i_end
5655
5656	absU=0.5sqrt((u(i,k,j)+u(i+1,k,j))2+(v(i,k,j)+v(i,k,j+1))*2)
5657	Cd = cd0
5658	tendency(i,k,j) = tendency(i,k,j) + &
5659	mu(i,j)( (u(i,k,j)+u(i+1,k,j))0.5* &
5660	CdabsU(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 )
5661
5662	END DO
5663	END DO
5664
5665	CASE (1,2) ! ustar computed from surface routine
5666
5667	DO j = j_start, j_end
5668	DO i = i_start, i_end
5669
5670	absU=0.5sqrt((u(i,k,j)+u(i+1,k,j))2+(v(i,k,j)+v(i,k,j+1))*2)+epsilon
5671	Cd = (ust(i,j)2)/(absU2)
5672	tendency(i,k,j) = tendency(i,k,j) + &
5673	mu(i,j)( (u(i,k,j)+u(i+1,k,j))0.5* &
5674	CdabsU(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 )
5675
5676	END DO
5677	END DO
5678
5679	CASE DEFAULT
5680	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
5681	END SELECT uflux
5682	! end of MARTA/WCS change
5683
5684	! For defor23.
5685
5686	DO j = j_start, j_end+1
5687	DO k = kts+1, ktf
5688	DO i = i_start, i_end
5689	tmp2(i,k,j) = defor23(i,k,j)
5690	END DO
5691	END DO
5692	END DO
5693
5694	DO j = j_start, j_end+1
5695	DO i = i_start, i_end
5696	tmp2(i,kts, j) = 0.0
5697	tmp2(i,ktf+1,j) = 0.0
5698	END DO
5699	END DO
5700
5701	DO j = j_start, j_end
5702	DO k = kts, ktf
5703	DO i = i_start, i_end
5704	avg(i,k,j) = 0.25 * &
5705	( ( tmp2(i,k+1,j )2 ) + ( tmp2(i,k,j )2) + &
5706	( tmp2(i,k+1,j+1)2 ) + ( tmp2(i,k,j+1)2) )
5707	END DO
5708	END DO
5709	END DO
5710
5711	DO j = j_start, j_end
5712	DO k = kts, ktf
5713	DO i = i_start, i_end
5714	tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j)
5715	END DO
5716	END DO
5717	END DO
5718
5719	!MARTA: add the drag at the surface; WCS 040331
5720	K=KTS
5721
5722	vflux: SELECT CASE( config_flags%isfflx )
5723	CASE (0) ! Assume cd a constant, specified in namelist
5724
5725	cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient
5726	! set in namelist.input
5727	DO j = j_start, j_end
5728	DO i = i_start, i_end
5729
5730	absU=0.5sqrt((u(i,k,j)+u(i+1,k,j))2+(v(i,k,j)+v(i,k,j+1))*2)
5731	Cd = cd0
5732	tendency(i,k,j) = tendency(i,k,j) + &
5733	mu(i,j)( (v(i,k,j)+v(i,k,j+1))0.5* &
5734	CdabsU(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 )
5735
5736	END DO
5737	END DO
5738
5739	CASE (1,2) ! ustar computed from surface routine
5740
5741	DO j = j_start, j_end
5742	DO i = i_start, i_end
5743
5744	absU=0.5sqrt((u(i,k,j)+u(i+1,k,j))2+(v(i,k,j)+v(i,k,j+1))*2)+epsilon
5745	Cd = (ust(i,j)2)/(absU2)
5746	tendency(i,k,j) = tendency(i,k,j) + &
5747	mu(i,j)( (v(i,k,j)+v(i,k,j+1))0.5* &
5748	CdabsU(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 )
5749
5750	END DO
5751	END DO
5752
5753	CASE DEFAULT
5754	CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' )
5755	END SELECT vflux
5756	! end of MARTA/WCS change
5757
5758	END SUBROUTINE tke_shear
5759
5760	!=======================================================================
5761	!=======================================================================
5762
5763	SUBROUTINE compute_diff_metrics( config_flags, ph, phb, z, rdz, rdzw, &
5764	zx, zy, rdx, rdy, &
5765	ids, ide, jds, jde, kds, kde, &
5766	ims, ime, jms, jme, kms, kme, &
5767	its, ite, jts, jte, kts, kte )
5768
5769	!-----------------------------------------------------------------------
5770	! Begin declarations.
5771
5772	IMPLICIT NONE
5773
5774	TYPE( grid_config_rec_type ), INTENT( IN ) &
5775	:: config_flags
5776
5777	INTEGER, INTENT( IN ) &
5778	:: ids, ide, jds, jde, kds, kde, &
5779	ims, ime, jms, jme, kms, kme, &
5780	its, ite, jts, jte, kts, kte
5781
5782	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5783	:: ph, phb
5784
5785	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) &
5786	:: rdz, rdzw, zx, zy, z
5787
5788	REAL, INTENT( IN ) &
5789	:: rdx, rdy
5790
5791	! Local variables.
5792
5793	INTEGER &
5794	:: i, j, k, i_start, i_end, j_start, j_end, ktf
5795
5796	! End declarations.
5797	!-----------------------------------------------------------------------
5798
5799	ktf = MIN( kte, kde-1 )
5800
5801	! Bug fix, WCS, 22 april 2002.
5802
5803	! We need rdzw in halo for average to u and v points.
5804
5805	j_start = jts-1
5806	j_end = jte
5807
5808	! Begin with dz computations.
5809
5810	DO j = j_start, j_end
5811
5812	IF ( ( j_start >= jts ) .AND. ( j_end <= MIN( jte, jde-1 ) ) ) THEN
5813	i_start = its-1
5814	i_end = ite
5815	ELSE
5816	i_start = its
5817	i_end = MIN( ite, ide-1 )
5818	END IF
5819
5820	! Compute z at w points for rdz and rdzw computations. We'll switch z
5821	! to z at p points before returning
5822
5823	DO k = 1, kte
5824
5825	! Bug fix, WCS, 22 april 2002
5826
5827	DO i = i_start, i_end
5828	z(i,k,j) = ( ph(i,k,j) + phb(i,k,j) ) / g
5829	END DO
5830	END DO
5831
5832	DO k = 1, ktf
5833	DO i = i_start, i_end
5834	rdzw(i,k,j) = 1.0 / ( z(i,k+1,j) - z(i,k,j) )
5835	END DO
5836	END DO
5837
5838	DO k = 2, ktf
5839	DO i = i_start, i_end
5840	rdz(i,k,j) = 2.0 / ( z(i,k+1,j) - z(i,k-1,j) )
5841	END DO
5842	END DO
5843
5844	! Bug fix, WCS, 22 april 2002; added the following code
5845
5846	DO i = i_start, i_end
5847	rdz(i,1,j) = 2./(z(i,2,j)-z(i,1,j))
5848	END DO
5849
5850	END DO
5851
5852	! End bug fix.
5853
5854	! Now compute zx and zy; we'll assume that the halo for ph and phb is
5855	! properly filled.
5856
5857	i_start = its
5858	i_end = MIN( ite, ide-1 )
5859	j_start = jts
5860	j_end = MIN( jte, jde-1 )
5861
5862	DO j = j_start, j_end
5863	DO k = 1, kte
5864	DO i = MAX( ids+1, its ), i_end
5865	zx(i,k,j) = rdx * ( phb(i,k,j) - phb(i-1,k,j) ) / g
5866	END DO
5867	END DO
5868	END DO
5869
5870	DO j = j_start, j_end
5871	DO k = 1, kte
5872	DO i = MAX( ids+1, its ), i_end
5873	zx(i,k,j) = zx(i,k,j) + rdx * ( ph(i,k,j) - ph(i-1,k,j) ) / g
5874	END DO
5875	END DO
5876	END DO
5877
5878	DO j = MAX( jds+1, jts ), j_end
5879	DO k = 1, kte
5880	DO i = i_start, i_end
5881	zy(i,k,j) = rdy * ( phb(i,k,j) - phb(i,k,j-1) ) / g
5882	END DO
5883	END DO
5884	END DO
5885
5886	DO j = MAX( jds+1, jts ), j_end
5887	DO k = 1, kte
5888	DO i = i_start, i_end
5889	zy(i,k,j) = zy(i,k,j) + rdy * ( ph(i,k,j) - ph(i,k,j-1) ) / g
5890	END DO
5891	END DO
5892	END DO
5893
5894	! Some b.c. on zx and zy.
5895
5896	IF ( .NOT. config_flags%periodic_x ) THEN
5897
5898	IF ( ite == ide ) THEN
5899	DO j = j_start, j_end
5900	DO k = 1, ktf
5901	zx(ide,k,j) = 0.0
5902	END DO
5903	END DO
5904	END IF
5905
5906	IF ( its == ids ) THEN
5907	DO j = j_start, j_end
5908	DO k = 1, ktf
5909	zx(ids,k,j) = 0.0
5910	END DO
5911	END DO
5912	END IF
5913
5914	ELSE
5915
5916	IF ( ite == ide ) THEN
5917	DO j=j_start,j_end
5918	DO k=1,ktf
5919	zx(ide,k,j) = rdx * ( phb(ide,k,j) - phb(ide-1,k,j) ) / g
5920	END DO
5921	END DO
5922
5923	DO j = j_start, j_end
5924	DO k = 1, ktf
5925	zx(ide,k,j) = zx(ide,k,j) + rdx * ( ph(ide,k,j) - ph(ide-1,k,j) ) / g
5926	END DO
5927	END DO
5928	END IF
5929
5930	IF ( its == ids ) THEN
5931	DO j = j_start, j_end
5932	DO k = 1, ktf
5933	zx(ids,k,j) = rdx * ( phb(ids,k,j) - phb(ids-1,k,j) ) / g
5934	END DO
5935	END DO
5936
5937	DO j =j_start,j_end
5938	DO k =1,ktf
5939	zx(ids,k,j) = zx(ids,k,j) + rdx * ( ph(ids,k,j) - ph(ids-1,k,j) ) / g
5940	END DO
5941	END DO
5942	END IF
5943
5944	END IF
5945
5946	IF ( .NOT. config_flags%periodic_y ) THEN
5947
5948	IF ( jte == jde ) THEN
5949	DO k =1, ktf
5950	DO i =i_start, i_end
5951	zy(i,k,jde) = 0.0
5952	END DO
5953	END DO
5954	END IF
5955
5956	IF ( jts == jds ) THEN
5957	DO k =1, ktf
5958	DO i =i_start, i_end
5959	zy(i,k,jds) = 0.0
5960	END DO
5961	END DO
5962	END IF
5963
5964	ELSE
5965
5966	IF ( jte == jde ) THEN
5967	DO j=j_start, j_end
5968	DO k=1, ktf
5969	zy(i,k,jde) = rdy * ( phb(i,k,jde) - phb(i,k,jde-1) ) / g
5970	END DO
5971	END DO
5972
5973	DO j = j_start, j_end
5974	DO k = 1, ktf
5975	zy(i,k,jde) = zy(i,k,jde) + rdy * ( ph(i,k,jde) - ph(i,k,jde-1) ) / g
5976	END DO
5977	END DO
5978	END IF
5979
5980	IF ( jts == jds ) THEN
5981	DO j = j_start, j_end
5982	DO k = 1, ktf
5983	zy(i,k,jds) = rdy * ( phb(i,k,jds) - phb(i,k,jds-1) ) / g
5984	END DO
5985	END DO
5986
5987	DO j = j_start, j_end
5988	DO k = 1, ktf
5989	zy(i,k,jds) = zy(i,k,jds) + rdy * ( ph(i,k,jds) - ph(i,k,jds-1) ) / g
5990	END DO
5991	END DO
5992	END IF
5993
5994	END IF
5995
5996	! Calculate z at p points.
5997
5998	DO j = j_start, j_end
5999	DO k = 1, ktf
6000	DO i = i_start, i_end
6001	z(i,k,j) = 0.5 * &
6002	( ph(i,k,j) + phb(i,k,j) + ph(i,k+1,j) + phb(i,k+1,j) ) / g
6003	END DO
6004	END DO
6005	END DO
6006
6007	END SUBROUTINE compute_diff_metrics
6008
6009	!=======================================================================
6010	!=======================================================================
6011
6012	END MODULE module_diffusion_em
6013
6014	!=======================================================================
6015	!=======================================================================
6016
6017

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

Download in other formats:

Original Format