Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_diffusion_em.F @ 2757

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

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

Download in other formats:

Original Format