source: trunk/WRF.COMMON/WRFV2/dyn_em/module_diffusion_em.F @ 2756

Last change on this file since 2756 was 1585, checked in by aslmd, 8 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
1018SUBROUTINE 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(degrad90*tmp)*kmmax*dampcoef
1088         dampkv(i,k,j)=cos(degrad90*tmp)*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(degrad90*tmp)*kmmax*dampcoef
1106         dampkv(i,k,j)=cos(degrad90*tmp)*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
1126END 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
1364SUBROUTINE 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
1420END SUBROUTINE isotropic_km
1421
1422!=======================================================================
1423!=======================================================================
1424
1425SUBROUTINE 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 = (dx*dy*dz)**(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_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_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_s*c_s*mlen_v*mlen_v*tmp, 1.0E-6*mlen_v*mlen_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=(dx*dy/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_s*c_s*deltas*deltas*tmp, 1.0E-6*deltas*deltas )
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
1575END SUBROUTINE smag_km
1576
1577!=======================================================================
1578!=======================================================================
1579
1580SUBROUTINE 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_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_h )
1661         xkmh(i,k,j)=c_s*c_s*mlen_h*mlen_h*tmp
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
1671END 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 = (dx*dy*dz)**(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_bound*pr_inv, xkmh(i,k,j) * pr_inv )
1840        xkhv(i,k,j)  = MIN( kv_tke_upper_bound*pr_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
1919SUBROUTINE 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
2125SUBROUTINE 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
2298END SUBROUTINE horizontal_diffusion_u_2
2299
2300!=======================================================================
2301!=======================================================================
2302
2303SUBROUTINE 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
2473END SUBROUTINE horizontal_diffusion_v_2
2474
2475!=======================================================================
2476!=======================================================================
2477
2478SUBROUTINE 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)*tmpzx*zeta_z(i,j)
2619!      titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy*zeta_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
2655END SUBROUTINE horizontal_diffusion_w_2
2656
2657!=======================================================================
2658!=======================================================================
2659
2660SUBROUTINE 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*(cf1*var(i  ,1,j)+cf2*var(i  ,2,j)+ &
2811                            cf3*var(i  ,3,j)+cf1*var(i-1,1,j)+  &
2812                            cf2*var(i-1,2,j)+cf3*var(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.5*dnw(ktes1)/dn(ktes1)+ &
2815                            var(i-1,ktes1,j)+(var(i-1,ktes1,j)- &
2816                            var(i-1,ktes2,j))*0.5*dnw(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)) - tmpzx*tmpzeta_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*(cf1*var(i,1,j  )+cf2*var(i  ,2,j)+ &
2862                            cf3*var(i,3,j  )+cf1*var(i,1,j-1)+  &
2863                            cf2*var(i,2,j-1)+cf3*var(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.5*dnw(ktes1)/dn(ktes1)+ &
2866                            var(i,ktes1,j-1)+(var(i,ktes1,j-1)- &
2867                            var(i,ktes2,j-1))*0.5*dnw(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)) - tmpzy*tmpzeta_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)=zx*H1avg*zeta_z
2901! H2avg(i,k,j)=zy*H2avg*zeta_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)*tmpzx*zeta_z(i,j)
2910!      H2avg(i,k,j)=H2avg(i,k,j)*tmpzy*zeta_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           (mrdx*0.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            mrdy*0.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
2956END SUBROUTINE horizontal_diffusion_s
2957
2958!=======================================================================
2959!=======================================================================
2960
2961SUBROUTINE 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_u*u_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_v*v_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_flux*rdzw(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
3273END SUBROUTINE vertical_diffusion_2
3274
3275!=======================================================================
3276!=======================================================================
3277
3278SUBROUTINE 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
3379END SUBROUTINE vertical_diffusion_u_2
3380
3381!=======================================================================
3382!=======================================================================
3383
3384SUBROUTINE 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
3485END SUBROUTINE vertical_diffusion_v_2
3486
3487!=======================================================================
3488!=======================================================================
3489
3490SUBROUTINE 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
3579END SUBROUTINE vertical_diffusion_w_2
3580
3581!=======================================================================
3582!=======================================================================
3583
3584SUBROUTINE 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
3706END 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
4168SUBROUTINE 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
4308END 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 = (dx*dy*dz)**(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.5*sqrt((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                     cd0*absU*defor13(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.5*sqrt((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                     cd0*absU*defor23(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.