Context Navigation

module_diffusion_em.F @ 3085

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

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

Download in other formats:

Original Format