wrfv2_fire/dyn_em/module_diffusion_em.F

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