source: trunk/WRF.COMMON/WRFV3/dyn_em/module_diffusion_em.F @ 2761

Last change on this file since 2761 was 2761, checked in by aslmd, 2 years ago

Applied planetary adaptation changes to WRFV3. job done previously by LMD_LES_MARS_install. Moved Registry.EM.newphys to put it simply in Registry.EM

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