dyn_em/module_sfs_nba.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2
   3 !==============================================================================
   4 !
   5 ! Copyright 2009. Lawrence Livermore National Security, LLC. All rights reserved.
   6 ! This work was produced at the Lawrence Livermore National Laboratory (LLNL) under
   7 ! contract no. DE-AC52-07NA27344 (Contract 44) between the U.S. Department of Energy (DOE)
   8 ! and Lawrence Livermore National Security, LLC (LLNS) for the operation of LLNL. Copyright
   9 ! is reserved to Lawrence Livermore National Security, LLC for purposes of controlled
  10 ! dissemination, commercialization through formal licensing, or other disposition under
  11 ! terms of Contract 44; DOE policies, regulations and orders; and U.S. statutes. The rights
  12 ! of the Federal Government are reserved under Contract 44.
  13 !
  14 ! DISCLAIMER
  15 ! This work was prepared as an account of work sponsored by an agency of the United States
  16 ! Government. Neither the United States Government nor Lawrence Livermore National
  17 ! Security, LLC nor any of their employees, makes any warranty, express or implied, or
  18 ! assumes any liability or responsibility for the accuracy, completeness, or usefulness of
  19 ! any information, apparatus, product, or process disclosed, or represents that its use
  20 ! would not infringe privately-owned rights. Reference herein to any specific commercial
  21 ! products, process, or service by trade name, trademark, manufacturer or otherwise does
  22 ! not necessarily constitute or imply its endorsement, recommendation, or favoring by the
  23 ! United States Government or Lawrence Livermore National Security, LLC. The views and
  24 ! opinions of authors expressed herein do not necessarily state or reflect those of the
  25 ! United States Government or Lawrence Livermore National Security, LLC, and shall not be
  26 ! used for advertising or product endorsement purposes.
  27 !
  28 ! LICENSING REQUIREMENTS
  29 ! Any use, reproduction, modification, or distribution of this software or documentation
  30 ! for commercial purposes requires a license from Lawrence Livermore National Security,
  31 ! LLC. Contact: Lawrence Livermore National Laboratory, Industrial Partnerships Office,
  32 ! P.O. Box 808, L-795, Livermore, CA 94551
  33 !
  34 !=============================================================================
  35 !
  36 ! Modification History:
  37 !
  38 ! Implemented 12/2009 by Jeff Mirocha, jmirocha@llnl.gov
  39 !
  40 !=============================================================================
  41
  42 MODULE module_sfs_nba
  43
  44   USE module_configure, ONLY : grid_config_rec_type
  45
  46   IMPLICIT NONE
  47
  48   REAL :: c1, c2, c3, ce, cb, cs ! global model parameters
  49
  50 CONTAINS
  51
  52 !=============================================================================
  53
  54 SUBROUTINE calc_mij_constants( )
  55
  56 !-----------------------------------------------------------------------------
  57 !
  58 ! PURPOSE: Compute constants for Mij calculations
  59 !
  60 !-----------------------------------------------------------------------------
  61
  62   IMPLICIT NONE
  63
  64   REAL :: sk, pi ! local model parameters
  65
  66 !-----------------------------------------------------------------------------
  67
  68   sk = 0.5
  69   pi = 3.1415927
  70   cb = 0.36
  71
  72   cs = ( ( 8.0*( 1.0+cb ) )/( 27.0*pi**2 ) )**0.5
  73   c1 = ( ( 960.0**0.5 )*cb )/( 7.0*( 1.0+cb )*sk )
  74   c2 = -c1
  75   ce = ( ( 8.0*pi/27.0 )**( 1.0/3.0 ) )*cs**( 4.0/3.0 )
  76   c3 = ( ( 27.0/( 8.0*pi ) )**( 1.0/3.0 ) )*cs**( 2.0/3.0 )
  77
  78   RETURN
  79
  80 END SUBROUTINE calc_mij_constants
  81
  82 !=============================================================================
  83
  84 SUBROUTINE calc_smnsmn( smnsmn,                       &
  85                         s11, s22, s33,                &
  86                         s12, s13, s23,                &
  87                         config_flags,                 &
  88                         ids, ide, jds, jde, kds, kde, &
  89                         ims, ime, jms, jme, kms, kme, &
  90                         ips, ipe, jps, jpe, kps, kpe, &
  91                         its, ite, jts, jte, kts, kte  )
  92
  93 !-----------------------------------------------------------------------------
  94 !
  95 ! PURPOSE: Compute Smn*Smn = S11^2 + S22^2 + S33^2 + 2*(S12^2 + S13^2 +S23^2)
  96 !
  97 !-----------------------------------------------------------------------------
  98
  99   IMPLICIT NONE
 100
 101   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT( OUT ) &
 102   :: smnsmn        ! Smn*Smn                             (s-2)
 103
 104   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT( IN  ) &
 105   :: s11         & ! 2*deformation element 11            (s-1)
 106    , s22         & ! 2*deformation element 22            (s-1)
 107    , s33         & ! 2*deformation element 33            (s-1)
 108    , s12         & ! 2*deformation element 12            (s-1)
 109    , s13         & ! 2*deformation element 13            (s-1)
 110    , s23           ! 2*deformation element 23            (s-1)
 111
 112   TYPE (grid_config_rec_type),                           INTENT( IN  ) &
 113   :: config_flags
 114
 115   INTEGER,                                               INTENT( IN  ) &
 116   :: ids, ide, jds, jde, kds, kde, &
 117      ims, ime, jms, jme, kms, kme, &
 118      ips, ipe, jps, jpe, kps, kpe, &
 119      its, ite, jts, jte, kts, kte
 120
 121 ! LOCAL VARIABLES ------------------------------------------------------------
 122
 123   REAL :: tmp
 124
 125   INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf
 126
 127 !-----------------------------------------------------------------------------
 128 !
 129 ! Set loop limits,
 130 ! taken from /dyn_em/module_diffusion_em.F SUBROUTINE smag_km
 131 !
 132 !-----------------------------------------------------------------------------
 133
 134   ktf = min(kte,kde-1)
 135
 136   i_start = its
 137   i_end   = MIN(ite,ide-1)
 138   j_start = jts
 139   j_end   = MIN(jte,jde-1)
 140
 141   IF ( config_flags%open_xs .or. config_flags%specified .or. &
 142        config_flags%nested) i_start = MAX(ids+1,its)
 143   IF ( config_flags%open_xe .or. config_flags%specified .or. &
 144        config_flags%nested) i_end   = MIN(ide-2,ite)
 145   IF ( config_flags%open_ys .or. config_flags%specified .or. &
 146        config_flags%nested) j_start = MAX(jds+1,jts)
 147   IF ( config_flags%open_ye .or. config_flags%specified .or. &
 148        config_flags%nested) j_end   = MIN(jde-2,jte)
 149   IF ( config_flags%periodic_x ) i_start = its
 150   IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
 151
 152 !-----------------------------------------------------------------------------
 153
 154 ! Below the 0.25 factor divides the incoming WRF deformations,
 155 ! which are multiplied by a factor of 2, by 2
 156
 157   DO j=j_start,j_end
 158   DO k=kts,ktf
 159   DO i=i_start,i_end
 160
 161     smnsmn(i,k,j) = 0.25*( s11(i,k,j)*s11(i,k,j) + &
 162                            s22(i,k,j)*s22(i,k,j) + &
 163                            s33(i,k,j)*s33(i,k,j) )
 164
 165   END DO
 166   END DO
 167   END DO
 168
 169 ! Below the 0.125 factor accounts for the four-point averaging (0.25)
 170 ! and divides the incoming WRF deformation elements by 2 (0.5)
 171
 172   DO j=j_start,j_end
 173   DO k=kts,ktf
 174   DO i=i_start,i_end
 175
 176     tmp = 0.125*( s12(i  ,k,j) + s12(i  ,k,j+1) + &
 177                   s12(i+1,k,j) + s12(i+1,k,j+1) )
 178     smnsmn(i,k,j) = smnsmn(i,k,j) + 2.0*tmp*tmp
 179
 180   END DO
 181   END DO
 182   END DO
 183
 184   DO j=j_start,j_end
 185   DO k=kts,ktf
 186   DO i=i_start,i_end
 187
 188     tmp = 0.125*( s13(i  ,k+1,j) + s13(i  ,k,j) + &
 189                   s13(i+1,k+1,j) + s13(i+1,k,j) )
 190     smnsmn(i,k,j) = smnsmn(i,k,j) + 2.0*tmp*tmp
 191
 192   END DO
 193   END DO
 194   END DO
 195
 196   DO j=j_start,j_end
 197   DO k=kts,ktf
 198   DO i=i_start,i_end
 199
 200     tmp = 0.125*( s23(i,k+1,j  ) + s23(i,k,j  ) + &
 201                   s23(i,k+1,j+1) + s23(i,k,j+1) )
 202     smnsmn(i,k,j) = smnsmn(i,k,j) + 2.0*tmp*tmp
 203
 204   END DO
 205   END DO
 206   END DO
 207
 208   RETURN
 209
 210 END SUBROUTINE calc_smnsmn
 211
 212 !=============================================================================
 213
 214 SUBROUTINE calc_mii( m11, m22, m33,                &
 215                      s11, s22, s33,                &
 216                      s12, s13, s23,                &
 217                      r12, r13, r23, smnsmn,        &
 218                      tke, rdzw, dx, dy,            &
 219                      config_flags,                 &
 220                      ids, ide, jds, jde, kds, kde, &
 221                      ims, ime, jms, jme, kms, kme, &
 222                      ips, ipe, jps, jpe, kps, kpe, &
 223                      its, ite, jts, jte, kts, kte  )
 224
 225 !-----------------------------------------------------------------------------
 226 !
 227 ! PURPOSE: Compute Mij for i = j
 228 !
 229 !-----------------------------------------------------------------------------
 230
 231   IMPLICIT NONE
 232
 233   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT( OUT ) &
 234   :: m11         & ! NBA stress element 11               (m2 s-2)
 235    , m22         & ! NBA stress element 22               (m2 s-2)
 236    , m33           ! NBA stress element 33               (m2 s-2)
 237
 238   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT( IN  ) &
 239   :: s11         & ! 2*deformation element 11            (s-1)
 240    , s22         & ! 2*deformation element 22            (s-1)
 241    , s33         & ! 2*deformation element 33            (s-1)
 242    , s12         & ! 2*deformation element 12            (s-1)
 243    , s13         & ! 2*deformation element 13            (s-1)
 244    , s23         & ! 2*deformation element 23            (s-1)
 245    , r12         & ! 2*rotation element 12               (s-1)
 246    , r13         & ! 2*rotation element 13               (s-1)
 247    , r23         & ! 2*rotation element 23               (s-1)
 248    , smnsmn      & ! Smn*Smn                             (s-2)
 249    , tke         & ! tke                                 (m2 s-2)
 250    , rdzw          ! 1/dz at w-levels                    (m-1)
 251
 252   REAL,                                                  INTENT( IN  ) &
 253   :: dx          & ! grid spacing in x                   (m)
 254    , dy            ! grid spacing in y                   (m)
 255
 256   TYPE (grid_config_rec_type),                           INTENT( IN  ) &
 257   :: config_flags
 258
 259   INTEGER,                                               INTENT( IN  ) &
 260   :: ids, ide, jds, jde, kds, kde, &
 261      ims, ime, jms, jme, kms, kme, &
 262      ips, ipe, jps, jpe, kps, kpe, &
 263      its, ite, jts, jte, kts, kte
 264
 265 ! LOCAL VARIABLES ------------------------------------------------------------
 266
 267   REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! sij/2, rij/2
 268   :: ss11        &
 269    , ss22        &
 270    , ss33        &
 271    , ss12        &
 272    , ss13        &
 273    , ss23        &
 274    , rr12        &
 275    , rr13        &
 276    , rr23
 277
 278   REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! projected to c
 279   :: ss12c       &
 280    , rr12c       &
 281    , ss13c       &
 282    , rr13c       &
 283    , ss23c       &
 284    , rr23c
 285
 286   REAL :: delta, a, b
 287
 288   INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf, is_ext, js_ext
 289
 290 !-----------------------------------------------------------------------------
 291 !
 292 ! Set loop limits,
 293 ! taken from /dyn_em/module_diffusion_em.F SUBROUTINE cal_titau_11_22_33
 294 !
 295 !-----------------------------------------------------------------------------
 296
 297   ktf = MIN( kte, kde-1 )
 298
 299   i_start = its
 300   i_end   = ite
 301   j_start = jts
 302   j_end   = jte
 303
 304     IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
 305          config_flags%nested) i_start = MAX( ids+1, its )
 306     IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
 307          config_flags%nested) i_end   = MIN( ide-1, ite )
 308     IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
 309          config_flags%nested) j_start = MAX( jds+1, jts )
 310     IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
 311          config_flags%nested) j_end   = MIN( jde-1, jte )
 312       IF ( config_flags%periodic_x ) i_start = its
 313       IF ( config_flags%periodic_x ) i_end = ite
 314
 315   is_ext = 1
 316   js_ext = 1
 317
 318   i_start = i_start - is_ext
 319   j_start = j_start - js_ext
 320
 321 !-----------------------------------------------------------------------------
 322 !
 323 ! Divide WRF deformations, which are multiplied by 2, by 2
 324 !
 325 !-----------------------------------------------------------------------------
 326
 327   DO j=j_start,j_end+1
 328   DO k=kts,ktf
 329   DO i=i_start,i_end+1
 330
 331     ss11(i,k,j)=s11(i,k,j)/2.0
 332     ss22(i,k,j)=s22(i,k,j)/2.0
 333     ss33(i,k,j)=s33(i,k,j)/2.0
 334     ss12(i,k,j)=s12(i,k,j)/2.0
 335     ss13(i,k,j)=s13(i,k,j)/2.0
 336     ss23(i,k,j)=s23(i,k,j)/2.0
 337     rr12(i,k,j)=r12(i,k,j)/2.0
 338     rr13(i,k,j)=r13(i,k,j)/2.0
 339     rr23(i,k,j)=r23(i,k,j)/2.0
 340
 341   END DO
 342   END DO
 343   END DO
 344
 345   DO j=j_start,j_end+1
 346   DO i=i_start,i_end+1
 347
 348     ss13(i,kde,j) = 0.0
 349     ss23(i,kde,j) = 0.0
 350     rr13(i,kde,j) = 0.0
 351     rr23(i,kde,j) = 0.0
 352
 353   END DO
 354   END DO
 355
 356 !-----------------------------------------------------------------------------
 357 !
 358 ! Project variables to c
 359 !
 360 !-----------------------------------------------------------------------------
 361
 362   DO j = j_start, j_end
 363   DO k = kts, ktf
 364   DO i = i_start, i_end
 365
 366     ss12c(i,k,j) = 0.25*( ss12(i  ,k  ,j  ) + ss12(i  ,k  ,j+1) + &
 367                           ss12(i+1,k  ,j  ) + ss12(i+1,k  ,j+1) )
 368
 369     rr12c(i,k,j) = 0.25*( rr12(i  ,k  ,j  ) + rr12(i  ,k  ,j+1) + &
 370                           rr12(i+1,k  ,j  ) + rr12(i+1,k  ,j+1) )
 371
 372     ss13c(i,k,j) = 0.25*( ss13(i  ,k+1,j  ) + ss13(i  ,k  ,j  ) + &
 373                           ss13(i+1,k+1,j  ) + ss13(i+1,k  ,j  ) )
 374
 375     rr13c(i,k,j) = 0.25*( rr13(i  ,k+1,j  ) + rr13(i  ,k  ,j  ) + &
 376                           rr13(i+1,k+1,j  ) + rr13(i+1,k  ,j  ) )
 377
 378     ss23c(i,k,j) = 0.25*( ss23(i  ,k+1,j  ) + ss23(i  ,k  ,j  ) + &
 379                           ss23(i  ,k+1,j+1) + ss23(i  ,k  ,j+1) )
 380
 381     rr23c(i,k,j) = 0.25*( rr23(i  ,k+1,j  ) + rr23(i  ,k  ,j  ) + &
 382                           rr23(i  ,k+1,j+1) + rr23(i  ,k  ,j+1) )
 383
 384   ENDDO
 385   ENDDO
 386   ENDDO
 387
 388 !-----------------------------------------------------------------------------
 389 !
 390 ! Calculate M11, M22 and M33
 391 !
 392 !-----------------------------------------------------------------------------
 393
 394   IF ( config_flags%sfs_opt .EQ. 1 ) THEN !Do not use TKE
 395
 396     DO j=j_start,j_end
 397     DO k=kts,ktf
 398     DO i=i_start,i_end
 399
 400       delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
 401       a = -1.0*( cs*delta )**2
 402
 403       m11(i,k,j) = a*(   2.0*sqrt( 2.0*smnsmn(i,k,j) )*ss11(i,k,j) &
 404                        + c1*(   ss11(i,k,j) *ss11(i,k,j)           &
 405                               + ss12c(i,k,j)*ss12c(i,k,j)          &
 406                               + ss13c(i,k,j)*ss13c(i,k,j)          &
 407                               - smnsmn(i,k,j)/3.0                  &
 408                             )                                      &
 409                        + c2*( -2.0*(   ss12c(i,k,j)*rr12c(i,k,j)   &
 410                                      + ss13c(i,k,j)*rr13c(i,k,j)   &
 411                                    )                               &
 412                             )                                      &
 413                      )
 414
 415       m22(i,k,j) = a*(   2.0*sqrt( 2.0*smnsmn(i,k,j) )*ss22(i,k,j) &
 416                        + c1*(   ss22(i,k,j) *ss22(i,k,j)           &
 417                               + ss12c(i,k,j)*ss12c(i,k,j)          &
 418                               + ss23c(i,k,j)*ss23c(i,k,j)          &
 419                               - smnsmn(i,k,j)/3.0                  &
 420                             )                                      &
 421                        + c2*(  2.0*(   ss12c(i,k,j)*rr12c(i,k,j)   &
 422                                      - ss23c(i,k,j)*rr23c(i,k,j)   &
 423                                    )                               &
 424                             )                                      &
 425                      )
 426
 427       m33(i,k,j) = a*(   2.0*sqrt( 2.0*smnsmn(i,k,j) )*ss33(i,k,j) &
 428                        + c1*(   ss33(i,k,j) *ss33(i,k,j)            &
 429                               + ss13c(i,k,j)*ss13c(i,k,j)          &
 430                               + ss23c(i,k,j)*ss23c(i,k,j)          &
 431                               - smnsmn(i,k,j)/3.0                  &
 432                             )                                      &
 433                        + c2*(  2.0*(   ss13c(i,k,j)*rr13c(i,k,j)   &
 434                                      + ss23c(i,k,j)*rr23c(i,k,j)   &
 435                                    )                               &
 436                             )                                      &
 437                      )
 438
 439     ENDDO
 440     ENDDO
 441     ENDDO
 442
 443   ELSE !(config_flags%sfs_opt .EQ. 2) Use TKE
 444
 445     DO j=j_start,j_end
 446     DO k=kts,ktf
 447     DO i=i_start,i_end
 448
 449       delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
 450       a = -1.0*ce*delta
 451       b = c3*delta
 452
 453       m11(i,k,j) = a*(   2.0*sqrt( tke(i,k,j) )*ss11(i,k,j)            &
 454                        + b*(                                           &
 455                                c1*(   ss11(i,k,j) *ss11(i,k,j)         &
 456                                     + ss12c(i,k,j)*ss12c(i,k,j)        &
 457                                     + ss13c(i,k,j)*ss13c(i,k,j)        &
 458                                     - smnsmn(i,k,j)/3.0                &
 459                                           )                            &
 460                              + c2*( -2.0*(   ss12c(i,k,j)*rr12c(i,k,j) &
 461                                            + ss13c(i,k,j)*rr13c(i,k,j) &
 462                                          )                             &
 463                                   )                                    &
 464                            )                                           &
 465                      )
 466
 467       m22(i,k,j) = a*(   2.0*sqrt( tke(i,k,j) )*ss22(i,k,j)            &
 468                        + b*(                                           &
 469                                c1*(   ss22(i,k,j) *ss22(i,k,j)         &
 470                                     + ss12c(i,k,j)*ss12c(i,k,j)        &
 471                                     + ss23c(i,k,j)*ss23c(i,k,j)        &
 472                                     - smnsmn(i,k,j)/3.0                &
 473                                             )                          &
 474                              + c2*(  2.0*(   ss12c(i,k,j)*rr12c(i,k,j) &
 475                                            - ss23c(i,k,j)*rr23c(i,k,j) &
 476                                          )                             &
 477                                   )                                    &
 478                            )                                           &
 479                      )
 480
 481       m33(i,k,j) = a*(   2.0*sqrt( tke(i,k,j) )*ss33(i,k,j)            &
 482                        + b*(                                           &
 483                                c1*(   ss33(i,k,j) *ss33(i,k,j)         &
 484                                     + ss13c(i,k,j)*ss13c(i,k,j)        &
 485                                     + ss23c(i,k,j)*ss23c(i,k,j)        &
 486                                     - smnsmn(i,k,j)/3.0                &
 487                                   )                                    &
 488                              + c2*(  2.0*(   ss13c(i,k,j)*rr13c(i,k,j) &
 489                                            + ss23c(i,k,j)*rr23c(i,k,j) &
 490                                          )                             &
 491                                   )                                    &
 492                            )                                           &
 493                      )
 494
 495
 496     ENDDO
 497     ENDDO
 498     ENDDO
 499
 500   ENDIF
 501
 502   RETURN
 503
 504 END SUBROUTINE calc_mii
 505
 506 !=============================================================================
 507
 508 SUBROUTINE calc_m12( m12,                          &
 509                      s11, s22,                     &
 510                      s12, s13, s23,                &
 511                      r12, r13, r23, smnsmn,        &
 512                      tke, rdzw, dx, dy,            &
 513                      config_flags,                 &
 514                      ids, ide, jds, jde, kds, kde, &
 515                      ims, ime, jms, jme, kms, kme, &
 516                      ips, ipe, jps, jpe, kps, kpe, &
 517                      its, ite, jts, jte, kts, kte  )
 518
 519 !-----------------------------------------------------------------------------
 520 !
 521 ! PURPOSE: Compute M12
 522 !
 523 !-----------------------------------------------------------------------------
 524
 525   IMPLICIT NONE
 526
 527   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT( OUT ) &
 528   :: m12           ! NBA stress element 12               (m2 s-2)
 529
 530   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT( IN  ) &
 531   :: s11         & ! 2*deformation element 11            (s-1)
 532    , s22         & ! 2*deformation element 22            (s-1)
 533    , s12         & ! 2*deformation element 12            (s-1)
 534    , s13         & ! 2*deformation element 13            (s-1)
 535    , s23         & ! 2*deformation element 23            (s-1)
 536    , r12         & ! 2*rotation element 12               (s-1)
 537    , r13         & ! 2*rotation element 13               (s-1)
 538    , r23         & ! 2*rotation element 23               (s-1)
 539    , smnsmn      & ! Smn*Smn                             (s-2)
 540    , tke         & ! tke                                 (m2 s-2)
 541    , rdzw          ! 1/dz at w-levels                    (m-1)
 542
 543   REAL,                                                  INTENT( IN  ) &
 544   :: dx          & ! grid spacing in x                   (m)
 545    , dy            ! grid spacing in y                   (m)
 546
 547   TYPE (grid_config_rec_type),                           INTENT( IN  ) &
 548   :: config_flags
 549
 550   INTEGER,                                               INTENT( IN  ) &
 551   :: ids, ide, jds, jde, kds, kde, &
 552      ims, ime, jms, jme, kms, kme, &
 553      ips, ipe, jps, jpe, kps, kpe, &
 554      its, ite, jts, jte, kts, kte
 555
 556 ! LOCAL VARIABLES ------------------------------------------------------------
 557
 558   REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! sij/2, rij/2
 559   :: ss11        &
 560    , ss22        &
 561    , ss12        &
 562    , ss13        &
 563    , ss23        &
 564    , rr12        &
 565    , rr13        &
 566    , rr23
 567
 568
 569   REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! projected to d
 570   :: tked        &
 571    , ss11d       &
 572    , ss22d       &
 573    , ss13d       &
 574    , ss23d       &
 575    , rr13d       &
 576    , rr23d       &
 577    , smnsmnd
 578
 579   REAL :: delta, a, b
 580
 581   INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf, je_ext, ie_ext
 582
 583 !-----------------------------------------------------------------------------
 584 !
 585 ! Set loop limits,
 586 ! taken from /dyn_em/module_diffusion_em.F SUBROUTINE cal_titau_12_21
 587 !
 588 !-----------------------------------------------------------------------------
 589
 590   ktf = MIN( kte, kde-1 )
 591
 592 ! Needs one more point in the x and y directions.
 593
 594   i_start = its
 595   i_end   = ite
 596   j_start = jts
 597   j_end   = jte
 598
 599     IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
 600          config_flags%nested ) i_start = MAX( ids+1, its )
 601     IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
 602          config_flags%nested ) i_end   = MIN( ide-1, ite )
 603     IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
 604          config_flags%nested ) j_start = MAX( jds+1, jts )
 605     IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
 606          config_flags%nested ) j_end   = MIN( jde-1, jte )
 607       IF ( config_flags%periodic_x ) i_start = its
 608       IF ( config_flags%periodic_x ) i_end = ite
 609
 610   je_ext = 1
 611   ie_ext = 1
 612
 613   i_end = i_end + ie_ext
 614   j_end = j_end + je_ext
 615
 616 !-----------------------------------------------------------------------------
 617 !
 618 ! Divide WRF deformations, which are multiplied by 2, by 2
 619 !
 620 !-----------------------------------------------------------------------------
 621
 622   DO j=j_start-1,j_end
 623   DO k=kts,ktf
 624   DO i=i_start-1,i_end
 625
 626     ss11(i,k,j)=s11(i,k,j)/2.0
 627     ss22(i,k,j)=s22(i,k,j)/2.0
 628     ss12(i,k,j)=s12(i,k,j)/2.0
 629     ss13(i,k,j)=s13(i,k,j)/2.0
 630     ss23(i,k,j)=s23(i,k,j)/2.0
 631     rr12(i,k,j)=r12(i,k,j)/2.0
 632     rr13(i,k,j)=r13(i,k,j)/2.0
 633     rr23(i,k,j)=r23(i,k,j)/2.0
 634
 635   END DO
 636   END DO
 637   END DO
 638
 639   DO j=j_start-1,j_end
 640   DO i=i_start-1,i_end
 641
 642     ss13(i,kde,j) = 0.0
 643     ss23(i,kde,j) = 0.0
 644     rr13(i,kde,j) = 0.0
 645     rr23(i,kde,j) = 0.0
 646
 647   END DO
 648   END DO
 649
 650 !-----------------------------------------------------------------------------
 651 !
 652 ! Project variables to d
 653 !
 654 !-----------------------------------------------------------------------------
 655
 656   DO j = j_start, j_end
 657   DO k = kts, ktf
 658   DO i = i_start, i_end
 659
 660     tked(i,k,j) = 0.25*( tke(i-1,k  ,j  ) + tke(i  ,k  ,j  ) + &
 661                          tke(i-1,k  ,j-1) + tke(i  ,k  ,j-1) )
 662
 663     smnsmnd(i,k,j) = 0.25*( smnsmn(i-1,k  ,j  ) + smnsmn(i  ,k  ,j  ) + &
 664                             smnsmn(i-1,k  ,j-1) + smnsmn(i  ,k  ,j-1) )
 665
 666     ss11d(i,k,j) = 0.25*( ss11(i-1,k  ,j  ) + ss11(i  ,k  ,j  ) + &
 667                           ss11(i-1,k  ,j-1) + ss11(i  ,k  ,j-1) )
 668
 669     ss22d(i,k,j) = 0.25*( ss22(i-1,k  ,j  ) + ss22(i  ,k  ,j  ) + &
 670                           ss22(i-1,k  ,j-1) + ss22(i  ,k  ,j-1) )
 671
 672     ss13d(i,k,j) = 0.25*( ss13(i  ,k+1,j  ) + ss13(i  ,k+1,j-1) + &
 673                           ss13(i  ,k  ,j  ) + ss13(i  ,k  ,j-1) )
 674
 675     rr13d(i,k,j) = 0.25*( rr13(i  ,k+1,j  ) + rr13(i  ,k+1,j-1) + &
 676                           rr13(i  ,k  ,j  ) + rr13(i  ,k  ,j-1) )
 677
 678     ss23d(i,k,j) = 0.25*( ss23(i  ,k+1,j  ) + ss23(i-1,k+1,j  ) + &
 679                           ss23(i  ,k  ,j  ) + ss23(i-1,k  ,j  ) )
 680
 681     rr23d(i,k,j) = 0.25*( rr23(i  ,k+1,j  ) + rr23(i-1,k+1,j  ) + &
 682                           rr23(i  ,k  ,j  ) + rr23(i-1,k  ,j  ) )
 683
 684   END DO
 685   END DO
 686   END DO
 687
 688 !-----------------------------------------------------------------------------
 689 !
 690 ! Calculate M12
 691 !
 692 !-----------------------------------------------------------------------------
 693
 694   IF ( config_flags%sfs_opt .EQ. 1 ) THEN !Do not use TKE
 695
 696     DO j=j_start,j_end
 697     DO k=kts,ktf
 698     DO i=i_start,i_end
 699
 700       delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
 701       a = -1.0*( cs*delta )**2
 702
 703       m12(i,k,j) = a*(   2.0*sqrt( 2.0*smnsmnd(i,k,j) )*ss12(i,k,j) &
 704                        + c1*(   ss11d(i,k,j)*ss12(i,k,j)            &
 705                               + ss22d(i,k,j)*ss12(i,k,j)            &
 706                               + ss13d(i,k,j)*ss23d(i,k,j)           &
 707                             )                                       &
 708                        + c2*(   ss11d(i,k,j)*rr12(i,k,j)            &
 709                               - ss13d(i,k,j)*rr23d(i,k,j)           &
 710                               - ss22d(i,k,j)*rr12(i,k,j)            &
 711                               - ss23d(i,k,j)*rr13d(i,k,j)           &
 712                             )                                       &
 713                       )
 714
 715     ENDDO
 716     ENDDO
 717     ENDDO
 718
 719   ELSE !(config_flags%sfs_opt .EQ. 2) Use TKE
 720
 721     DO j=j_start,j_end
 722     DO k=kts,ktf
 723     DO i=i_start,i_end
 724
 725       delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
 726       a = -1.0*ce*delta
 727       b = c3*delta
 728
 729       m12(i,k,j) = a*(   2.0*sqrt( tked(i,k,j) )*s12(i,k,j)     &
 730                        + b*(                                    &
 731                                c1*(   ss11d(i,k,j)*ss12(i,k,j)  &
 732                                     + ss22d(i,k,j)*ss12(i,k,j)  &
 733                                     + ss13d(i,k,j)*ss23d(i,k,j) &
 734                                   )                             &
 735                              + c2*(   ss11d(i,k,j)*rr12(i,k,j)  &
 736                                     - ss13d(i,k,j)*rr23d(i,k,j) &
 737                                     - ss22d(i,k,j)*rr12(i,k,j)  &
 738                                     - ss23d(i,k,j)*rr13d(i,k,j) &
 739                                   )                             &
 740                            )                                    &
 741                      )
 742     ENDDO
 743     ENDDO
 744     ENDDO
 745
 746   ENDIF
 747
 748   RETURN
 749
 750 END SUBROUTINE calc_m12
 751
 752 !=============================================================================
 753
 754 SUBROUTINE calc_m13( m13,                          &
 755                      s11, s33,                     &
 756                      s12, s13, s23,                &
 757                      r12, r13, r23, smnsmn,        &
 758                      tke, rdzw, dx, dy,            &
 759                      fnm, fnp,                     &
 760                      config_flags,                 &
 761                      ids, ide, jds, jde, kds, kde, &
 762                      ims, ime, jms, jme, kms, kme, &
 763                      ips, ipe, jps, jpe, kps, kpe, &
 764                      its, ite, jts, jte, kts, kte  )
 765
 766 !-----------------------------------------------------------------------------
 767 !
 768 ! PURPOSE: Compute M13
 769 !
 770 !-----------------------------------------------------------------------------
 771
 772   IMPLICIT NONE
 773
 774   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT( OUT ) &
 775   :: m13           !                                     (m2 s-2)
 776
 777   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT( IN  ) &
 778   :: s11         & ! 2*deformation element 11            (s-1)
 779    , s33         & ! 2*deformation element 33            (s-1)
 780    , s12         & ! 2*deformation element 12            (s-1)
 781    , s13         & ! 2*deformation element 13            (s-1)
 782    , s23         & ! 2*deformation element 23            (s-1)
 783    , r12         & ! 2*rotation element 12               (s-1)
 784    , r13         & ! 2*rotation element 13               (s-1)
 785    , r23         & ! 2*rotation element 23               (s-1)
 786    , smnsmn      & ! Smn*Smn                             (s-2)
 787    , tke         & ! tke                                 (m2 s-2)
 788    , rdzw          ! 1/dz at w-levels                    (m-1)
 789
 790   REAL,                                                  INTENT( IN  ) &
 791   :: dx          & ! grid spacing in x                   (m)
 792    , dy            ! grid spacing in y                   (m)
 793
 794   REAL, DIMENSION( kms:kme ),                            INTENT( IN  ) &
 795   :: fnm         & ! vertical interpolation coefficients
 796    , fnp           !
 797
 798   TYPE (grid_config_rec_type),                           INTENT( IN  ) &
 799   :: config_flags
 800
 801   INTEGER,                                               INTENT( IN  ) &
 802   :: ids, ide, jds, jde, kds, kde, &
 803      ims, ime, jms, jme, kms, kme, &
 804      ips, ipe, jps, jpe, kps, kpe, &
 805      its, ite, jts, jte, kts, kte
 806
 807 ! LOCAL VARIABLES ------------------------------------------------------------
 808
 809   REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! sij/2, rij/2
 810   :: ss11        &
 811    , ss33        &
 812    , ss12        &
 813    , ss13        &
 814    , ss23        &
 815    , rr12        &
 816    , rr13        &
 817    , rr23
 818
 819   REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! projected to e
 820   :: tkee        &
 821    , ss11e       &
 822    , ss33e       &
 823    , ss12e       &
 824    , ss23e       &
 825    , rr12e       &
 826    , rr23e       &
 827    , smnsmne
 828
 829   REAL :: delta, a, b
 830
 831   INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf, ie_ext
 832
 833 !-----------------------------------------------------------------------------
 834 !
 835 ! Set loop limits,
 836 ! taken from /dyn_em/module_diffusion_em.F SUBROUTINE cal_titau_13_31
 837 !
 838 !-----------------------------------------------------------------------------
 839
 840   ktf = MIN( kte, kde-1 )
 841
 842 ! Find ide-1 and jde-1 for averaging to p point.
 843
 844   i_start = its
 845   i_end   = ite
 846   j_start = jts
 847   j_end   = MIN( jte, jde-1 )
 848
 849     IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
 850          config_flags%nested) i_start = MAX( ids+1, its )
 851     IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
 852          config_flags%nested) i_end   = MIN( ide-1, ite )
 853     IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
 854          config_flags%nested) j_start = MAX( jds+1, jts )
 855     IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
 856          config_flags%nested) j_end   = MIN( jde-2, jte )
 857       IF ( config_flags%periodic_x ) i_start = its
 858       IF ( config_flags%periodic_x ) i_end = ite
 859
 860   ie_ext = 1
 861   i_end = i_end + ie_ext
 862
 863 !-----------------------------------------------------------------------------
 864 !
 865 ! Divide WRF deformations, which are multiplied by 2, by 2
 866 !
 867 !-----------------------------------------------------------------------------
 868
 869   DO j=j_start,j_end+1
 870   DO k=kts,ktf
 871   DO i=i_start-1,i_end
 872
 873     ss11(i,k,j)=s11(i,k,j)/2.0
 874     ss33(i,k,j)=s33(i,k,j)/2.0
 875     ss12(i,k,j)=s12(i,k,j)/2.0
 876     ss13(i,k,j)=s13(i,k,j)/2.0
 877     ss23(i,k,j)=s23(i,k,j)/2.0
 878     rr12(i,k,j)=r12(i,k,j)/2.0
 879     rr13(i,k,j)=r13(i,k,j)/2.0
 880     rr23(i,k,j)=r23(i,k,j)/2.0
 881
 882   END DO
 883   END DO
 884   END DO
 885
 886 !-----------------------------------------------------------------------------
 887 !
 888 ! Project variables to e
 889 !
 890 !-----------------------------------------------------------------------------
 891
 892   DO j = j_start, j_end
 893   DO k = kts+1, ktf
 894   DO i = i_start, i_end
 895
 896     tkee(i,k,j) = 0.5*( fnm(k)*( tke(i,k  ,j) + tke(i-1,k  ,j) ) + &
 897                         fnp(k)*( tke(i,k-1,j) + tke(i-1,k-1,j) ) )
 898
 899     smnsmne(i,k,j) = 0.5*( fnm(k)*( smnsmn(i,k  ,j) + smnsmn(i-1,k  ,j) ) + &
 900                            fnp(k)*( smnsmn(i,k-1,j) + smnsmn(i-1,k-1,j) ) )
 901
 902     ss11e(i,k,j) = 0.5*( fnm(k)*( ss11(i  ,k  ,j  ) + ss11(i-1,k  ,j  ) ) + &
 903                          fnp(k)*( ss11(i  ,k-1,j  ) + ss11(i-1,k-1,j  ) ) )
 904
 905     ss33e(i,k,j) = 0.5*( fnm(k)*( ss33(i  ,k  ,j  ) + ss33(i-1,k  ,j  ) ) + &
 906                          fnp(k)*( ss33(i  ,k-1,j  ) + ss33(i-1,k-1,j  ) ) )
 907
 908     ss12e(i,k,j) = 0.5*( fnm(k)*( ss12(i  ,k  ,j  ) + ss12(i  ,k  ,j+1) ) + &
 909                          fnp(k)*( ss12(i  ,k-1,j  ) + ss12(i  ,k-1,j+1) ) )
 910
 911     rr12e(i,k,j) = 0.5*( fnm(k)*( rr12(i  ,k  ,j  ) + rr12(i  ,k  ,j+1) ) + &
 912                          fnp(k)*( rr12(i  ,k-1,j  ) + rr12(i  ,k-1,j+1) ) )
 913
 914     ss23e(i,k,j) = 0.25*( ss23(i  ,k  ,j) + ss23(i  ,k  ,j+1) + &
 915                           ss23(i-1,k  ,j) + ss23(i-1,k  ,j+1) )
 916
 917     rr23e(i,k,j) = 0.25*( rr23(i  ,k  ,j) + rr23(i  ,k  ,j+1) + &
 918                           rr23(i-1,k  ,j) + rr23(i-1,k  ,j+1) )
 919
 920   END DO
 921   END DO
 922   END DO
 923
 924 !-----------------------------------------------------------------------------
 925 !
 926 ! Calculate M_13
 927 !
 928 !-----------------------------------------------------------------------------
 929
 930
 931   IF ( config_flags%sfs_opt .EQ. 1 ) THEN !Do not use TKE
 932
 933     DO j=j_start,j_end
 934     DO k=kts+1,ktf
 935     DO i=i_start,i_end
 936
 937       delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
 938       a = -1.0*( cs*delta )**2
 939
 940       m13(i,k,j) = a*(   2.0*sqrt( 2.0*smnsmne(i,k,j) )*ss13(i,k,j) &
 941                        + c1*(   ss11e(i,k,j)*ss13(i,k,j)            &
 942                               + ss12e(i,k,j)*ss23e(i,k,j)           &
 943                               + ss13(i,k,j)*ss33e(i,k,j)            &
 944                             )                                       &
 945                        + c2*(   ss11e(i,k,j)*rr13(i,k,j)            &
 946                               + ss12e(i,k,j)*rr23e(i,k,j)           &
 947                               - ss23e(i,k,j)*rr12e(i,k,j)           &
 948                               - ss33e(i,k,j)*rr13(i,k,j)            &
 949                             )                                       &
 950                      )
 951
 952     ENDDO
 953     ENDDO
 954     ENDDO
 955
 956   ELSE !(config_flags%sfs_opt .EQ. 2) Use TKE
 957
 958     DO j=j_start,j_end
 959     DO k=kts+1,ktf
 960     DO i=i_start,i_end
 961
 962       delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
 963       a = -1.0*ce*delta
 964       b = c3*delta
 965
 966       m13(i,k,j) = a*(   2.0*sqrt( tkee(i,k,j) )*ss13(i,k,j)    &
 967                        + b*(                                    &
 968                                c1*(   ss11e(i,k,j)*ss13(i,k,j)  &
 969                                     + ss12e(i,k,j)*ss23e(i,k,j) &
 970                                     + ss13(i,k,j)*ss33e(i,k,j)  &
 971                                   )                             &
 972                              + c2*(   ss11e(i,k,j)*rr13(i,k,j)  &
 973                                     + ss12e(i,k,j)*rr23e(i,k,j) &
 974                                     - ss23e(i,k,j)*rr12e(i,k,j) &
 975                                     - ss33e(i,k,j)*rr13(i,k,j)  &
 976                                   )                             &
 977                            )                                    &
 978                      )
 979
 980     ENDDO
 981     ENDDO
 982     ENDDO
 983
 984   ENDIF
 985
 986   RETURN
 987
 988 END SUBROUTINE calc_m13
 989
 990 !=============================================================================
 991
 992 SUBROUTINE calc_m23( m23,                          &
 993                      s22, s33,                     &
 994                      s12, s13, s23,                &
 995                      r12, r13, r23, smnsmn,        &
 996                      tke, rdzw, dx, dy,            &
 997                      fnm, fnp,                     &
 998                      config_flags,                 &
 999                      ids, ide, jds, jde, kds, kde, &
1000                      ims, ime, jms, jme, kms, kme, &
1001                      ips, ipe, jps, jpe, kps, kpe, &
1002                      its, ite, jts, jte, kts, kte  )
1003
1004 !-----------------------------------------------------------------------------
1005 !
1006 ! PURPOSE: Compute M23
1007 !
1008 !-----------------------------------------------------------------------------
1009
1010   IMPLICIT NONE
1011
1012   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT( OUT ) &
1013   :: m23           !                                     (m2 s-2)
1014
1015   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),          INTENT(  IN ) &
1016   :: s22         & ! 2*deformation element 22            (s-1)
1017    , s33         & ! 2*deformation element 33            (s-1)
1018    , s12         & ! 2*deformation element 12            (s-1)
1019    , s13         & ! 2*deformation element 13            (s-1)
1020    , s23         & ! 2*deformation element 23            (s-1)
1021    , r12         & ! 2*rotation element 12               (s-1)
1022    , r13         & ! 2*rotation element 13               (s-1)
1023    , r23         & ! 2*rotation element 23               (s-1)
1024    , smnsmn      & ! Smn*Smn                             (s-2)
1025    , tke         & ! tke                                 (m2 s-2)
1026    , rdzw          ! 1/dz at w-levels                    (m-1)
1027
1028   REAL,                                                  INTENT( IN  ) &
1029   :: dx          & ! grid spacing in x                   (m)
1030    , dy            ! grid spacing in y                   (m)
1031
1032   REAL, DIMENSION( kms:kme ),                            INTENT( IN  ) &
1033   :: fnm         & ! vertical interpolation coefficients
1034    , fnp           !
1035
1036   TYPE (grid_config_rec_type),                           INTENT( IN  ) &
1037   :: config_flags
1038
1039   INTEGER,                                               INTENT( IN  ) &
1040   :: ids, ide, jds, jde, kds, kde, &
1041      ims, ime, jms, jme, kms, kme, &
1042      ips, ipe, jps, jpe, kps, kpe, &
1043      its, ite, jts, jte, kts, kte
1044
1045 ! LOCAL VARIABLES ------------------------------------------------------------
1046
1047   REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! sij/2, rij/2
1048   :: ss22        &
1049    , ss33        &
1050    , ss12        &
1051    , ss13        &
1052    , ss23        &
1053    , rr12        &
1054    , rr13        &
1055    , rr23
1056
1057   REAL, DIMENSION( its-1:ite+1, kms:kme, jts-1:jte+1 ) & ! projected to f
1058   :: tkef        &
1059    , ss22f       &
1060    , ss33f       &
1061    , ss12f       &
1062    , ss13f       &
1063    , rr12f       &
1064    , rr13f       &
1065    , smnsmnf
1066
1067   REAL :: delta, a, b
1068
1069   INTEGER :: i, j, k, i_start, i_end, j_start, j_end, ktf, je_ext
1070
1071 !-----------------------------------------------------------------------------
1072 !
1073 ! Set loop limits,
1074 ! taken from /dyn_em/module_diffusion_em.F SUBROUTINE cal_titau_23_32
1075 !
1076 !-----------------------------------------------------------------------------
1077
1078   ktf = MIN( kte, kde-1 )
1079
1080 ! Find ide-1 and jde-1 for averaging to p point.
1081
1082   i_start = its
1083   i_end   = MIN( ite, ide-1 )
1084   j_start = jts
1085   j_end   = jte
1086
1087     IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
1088          config_flags%nested) i_start = MAX( ids+1, its )
1089     IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
1090          config_flags%nested) i_end   = MIN( ide-2, ite )
1091     IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
1092          config_flags%nested) j_start = MAX( jds+1, jts )
1093     IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
1094          config_flags%nested) j_end   = MIN( jde-1, jte )
1095       IF ( config_flags%periodic_x ) i_start = its
1096       IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 )
1097
1098   je_ext = 1
1099   j_end = j_end + je_ext
1100
1101 !-----------------------------------------------------------------------------
1102 !
1103 ! Divide WRF deformations, which are multiplied by 2, by 2
1104 !
1105 !-----------------------------------------------------------------------------
1106
1107   DO j=j_start-1,j_end
1108   DO k=kts,ktf
1109   DO i=i_start,i_end+1
1110
1111     ss22(i,k,j)=s22(i,k,j)/2.0
1112     ss33(i,k,j)=s33(i,k,j)/2.0
1113     ss12(i,k,j)=s12(i,k,j)/2.0
1114     ss13(i,k,j)=s13(i,k,j)/2.0
1115     ss23(i,k,j)=s23(i,k,j)/2.0
1116     rr12(i,k,j)=r12(i,k,j)/2.0
1117     rr13(i,k,j)=r13(i,k,j)/2.0
1118     rr23(i,k,j)=r23(i,k,j)/2.0
1119
1120   END DO
1121   END DO
1122   END DO
1123
1124 !-----------------------------------------------------------------------------
1125 !
1126 ! Project variables to f
1127 !
1128 !-----------------------------------------------------------------------------
1129
1130   DO j = j_start, j_end
1131   DO k = kts+1, ktf
1132   DO i = i_start, i_end
1133
1134     tkef(i,k,j) = 0.5*( fnm(k)*( tke(i  ,k  ,j  ) + tke(i  ,k  ,j-1) ) + &
1135                         fnp(k)*( tke(i  ,k-1,j  ) + tke(i  ,k-1,j-1) ) )
1136
1137     smnsmnf(i,k,j) = 0.5*( fnm(k)*( smnsmn(i  ,k  ,j  ) + smnsmn(i  ,k  ,j-1) ) + &
1138                            fnp(k)*( smnsmn(i  ,k-1,j  ) + smnsmn(i  ,k-1,j-1) ) )
1139
1140     ss22f(i,k,j) = 0.5*( fnm(k)*( ss22(i  ,k  ,j  ) + ss22(i  ,k  ,j-1) ) + &
1141                          fnp(k)*( ss22(i  ,k-1,j  ) + ss22(i  ,k-1,j-1) ) )
1142
1143     ss33f(i,k,j) = 0.5*( fnm(k)*( ss33(i  ,k  ,j  ) + ss33(i  ,k  ,j-1) ) + &
1144                          fnp(k)*( ss33(i  ,k-1,j  ) + ss33(i  ,k-1,j-1) ) )
1145
1146     ss12f(i,k,j) = 0.5*( fnm(k)*( ss12(i  ,k  ,j  ) + ss12(i+1,k  ,j  ) ) + &
1147                          fnp(k)*( ss12(i  ,k-1,j  ) + ss12(i+1,k-1,j  ) ) )
1148
1149     rr12f(i,k,j) = 0.5*( fnm(k)*( rr12(i  ,k  ,j  ) + rr12(i+1,k  ,j  ) ) + &
1150                          fnp(k)*( rr12(i  ,k-1,j  ) + rr12(i+1,k-1,j  ) ) )
1151
1152     ss13f(i,k,j) = 0.25*( ss13(i  ,k  ,j  ) + ss13(i  ,k  ,j-1) + &
1153                           ss13(i+1,k  ,j-1) + ss13(i+1,k  ,j  ) )
1154
1155     rr13f(i,k,j) = 0.25*( rr13(i  ,k  ,j  ) + rr13(i  ,k  ,j-1) + &
1156                           rr13(i+1,k  ,j-1) + rr13(i+1,k  ,j  ) )
1157
1158   END DO
1159   END DO
1160   END DO
1161
1162 !-----------------------------------------------------------------------------
1163 !
1164 ! Calculate M23
1165 !
1166 !-----------------------------------------------------------------------------
1167
1168   IF ( config_flags%sfs_opt .EQ. 1 ) THEN !Do not use TKE
1169
1170     DO j=j_start,j_end
1171     DO k=kts+1,ktf
1172     DO i=i_start,i_end
1173
1174       delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
1175       a = -1.0*( cs*delta )**2
1176
1177       m23(i,k,j) = a*(   2.0*sqrt( 2.0*smnsmnf(i,k,j) )*ss23(i,k,j) &
1178                        + c1*(   ss12f(i,k,j)*ss13f(i,k,j)           &
1179                               + ss22f(i,k,j)*ss23(i,k,j)            &
1180                               + ss23(i,k,j) *ss33f(i,k,j)           &
1181                              )                                      &
1182                        + c2*(   ss12f(i,k,j)*rr13f(i,k,j)           &
1183                               + ss22f(i,k,j)*rr23(i,k,j)            &
1184                               + ss13f(i,k,j)*rr12f(i,k,j)           &
1185                               - ss33f(i,k,j)*rr23(i,k,j)            &
1186                             )                                       &
1187                      )
1188
1189     ENDDO
1190     ENDDO
1191     ENDDO
1192
1193   ELSE !(config_flags%sfs_opt .EQ. 2) Use TKE
1194
1195     DO j=j_start,j_end
1196     DO k=kts+1,ktf
1197     DO i=i_start,i_end
1198
1199       delta = ( dx * dy / rdzw(i,k,j) )**0.33333333
1200       a = -1.0*ce*delta
1201       b = c3*delta
1202
1203       m23(i,k,j) = a*(   2.0*sqrt( tkef(i,k,j) )*ss23(i,k,j)    &
1204                        + b*(                                    &
1205                                c1*(   ss12f(i,k,j)*ss13f(i,k,j) &
1206                                     + ss22f(i,k,j)*ss23(i,k,j)  &
1207                                     + ss23(i,k,j) *ss33f(i,k,j) &
1208                                   )                             &
1209                              + c2*(   ss12f(i,k,j)*rr13f(i,k,j) &
1210                                     + ss22f(i,k,j)*rr23(i,k,j)  &
1211                                     + ss13f(i,k,j)*rr12f(i,k,j) &
1212                                     - ss33f(i,k,j)*rr23(i,k,j)  &
1213                                   )                             &
1214                            )                                    &
1215                      )
1216
1217     ENDDO
1218     ENDDO
1219     ENDDO
1220
1221   ENDIF
1222
1223   RETURN
1224
1225 END SUBROUTINE calc_m23
1226
1227 !=============================================================================
1228
1229 END MODULE module_sfs_nba