dyn_em/module_ieva_em.F

   1 MODULE module_ieva_em
   2
   3   USE module_bc
   4   USE module_model_constants
   5   USE module_wrf_error
   6
   7   REAL(KIND=8), PARAMETER :: eps        = 1.0d-08
   8   REAL,         PARAMETER :: alpha_max  = 1.1
   9   REAL,         PARAMETER :: alpha_min  = 0.8
  10
  11 CONTAINS
  12
  13 LOGICAL FUNCTION CHK_IEVA( config_flags, rk_step )
  14
  15    IMPLICIT NONE
  16
  17    TYPE(grid_config_rec_type), INTENT(IN) :: config_flags
  18    INTEGER,                    INTENT(IN) :: rk_step
  19
  20    INTEGER :: zadvect_implicit
  21    INTEGER :: rk_order
  22
  23    rk_order         = config_flags%rk_ord
  24    zadvect_implicit = config_flags%zadvect_implicit
  25
  26    CHK_IEVA = .FALSE.
  27
  28 ! These option to run on the sub-steps are off - they dont make much of a difference
  29
  30 !  IF( zadvect_implicit .eq. 3 ) THEN   ! DO IEVA on all sub-steps of RK integrator
  31 !      CHK_IEVA = .TRUE.
  32 !  ENDIF
  33 !
  34 !  IF( zadvect_implicit .eq. 2 .and. rk_step .ge. rk_order-1 ) THEN ! DO IEVA on last two sub-steps of RK integrator
  35 !      CHK_IEVA = .TRUE.
  36 !  ENDIF
  37
  38    IF( zadvect_implicit .gt. 0 .and. rk_step .eq. rk_order ) THEN  ! DO IEVA on last sub-step of RK integrator
  39        CHK_IEVA = .TRUE.
  40    ENDIF
  41
  42 RETURN
  43 END FUNCTION CHK_IEVA
  44
  45 !-------------------------------------------------------------------------------
  46 ! IEVA code for splitting vertical velocity into pieces for explicit and implicit
  47 ! advection by the routines based on Shchepetkin (2015)
  48
  49 SUBROUTINE WW_SPLIT(wwE, wwI,                         &
  50                     u, v, ww,                         &
  51                     mut, rdnw, msfty,                 &
  52                     c1f, c2f,                         &
  53                     rdx, rdy, msfux, msfuy,           &
  54                     msfvx, msfvy, dt,                 &
  55                     config_flags, rk_step,            &
  56                     ids, ide, jds, jde, kds, kde,     &
  57                     ims, ime, jms, jme, kms, kme,     &
  58                     its, ite, jts, jte, kts, kte )
  59
  60    TYPE( grid_config_rec_type ) ,   INTENT(IN   ) :: config_flags
  61
  62    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
  63                                        ims, ime, jms, jme, kms, kme, &
  64                                        its, ite, jts, jte, kts, kte, &
  65                                        rk_step
  66
  67    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::   u, v, ww
  68
  69    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) ::  wwE, wwI
  70
  71    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: mut
  72
  73    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw
  74    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: c1f, c2f
  75
  76    REAL, INTENT(IN)    :: dt
  77    REAL, INTENT(IN)    :: rdx, rdy
  78    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux, msfuy
  79    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfvx, msfvy, msfty
  80
  81 ! Local variables
  82
  83    REAL(KIND=4)    :: cr, cx, cy, c0, c2d, cw_max2, cw_min, cw, cff, rdz
  84    INTEGER         :: i, j, k, ii, jj, i_start, i_end, j_start, j_end, ktf
  85    LOGICAL         :: print_flag = .true.
  86    SAVE            :: print_flag
  87    LOGICAL         :: degrade_xs, degrade_ys, degrade_xe, degrade_ye
  88
  89 ! Implicit advection parameters based on Shchepetkin (2015)
  90 ! Moved these to the top of the module.
  91
  92 !  REAL, PARAMETER :: alpha_max  = 1.1
  93 !  REAL, PARAMETER :: alpha_min  = 0.9
  94    REAL, PARAMETER :: cmnx_ratio = alpha_min/alpha_max
  95    REAL, PARAMETER :: cutoff     = 2.0 - cmnx_ratio
  96    REAL, PARAMETER :: r4cmx      = 1.0/(4.0 - 4.0*cmnx_ratio )
  97    REAL, PARAMETER :: Ceps       = 0.9
  98
  99    LOGICAL         :: ieva
 100
 101 ! Debug declarations
 102
 103 !  REAL,    PARAMETER :: binfrac = 0.1
 104 !  INTEGER, PARAMETER :: binsize = 11
 105 !  INTEGER            :: count(binsize)
 106 !  real(KIND=4)       :: binavg(binsize), crs(binsize)
 107
 108 ! Check to see if IEVA is on for this RK3 substep
 109
 110    ieva = CHK_IEVA(config_flags, rk_step)
 111
 112 ! =========================================================================
 113 ! Loop limits:  need to compute wwE & wwI one halo ring out for staggered variables.
 114
 115    ktf     = kte
 116    i_start = its-1
 117    i_end   = MIN(ite,ide-1)+1
 118    j_start = jts-1
 119    j_end   = MIN(jte,jde-1)+1
 120
 121
 122 ! =========================================================================
 123 ! PRINT OUT IEVA INFO
 124
 125    IF( print_flag .and. ieva ) THEN
 126      write(wrf_err_message,*) '-- BEGIN WW_SPLIT INFO ----------------------------------------'
 127      CALL wrf_debug( 0, wrf_err_message )
 128      write(wrf_err_message,*) 'RK_STEP ##', rk_step
 129      CALL wrf_debug( 0, wrf_err_message )
 130      write(wrf_err_message,*) 'IMS:IME     /  JME:JME     ', ims, ime, jms, jme
 131      CALL wrf_debug( 0, wrf_err_message )
 132      write(wrf_err_message,*) 'ISTART:IEND / JSTART:JEND  ', i_start,i_end,j_start,j_end
 133      CALL wrf_debug( 0, wrf_err_message )
 134      write(wrf_err_message,*) 'KMS:KME     /   KDS:KDE    KTS:KTE   ', kms,kme,kds,kde,kts,kte
 135      CALL wrf_debug( 0, wrf_err_message )
 136      write(wrf_err_message,*) '----------------------------------------'
 137      CALL wrf_debug( 0, wrf_err_message )
 138      write(wrf_err_message,*) 'WW_SPLIT:  ZADVECT_IMPLICT = ', config_flags%zadvect_implicit
 139      CALL wrf_debug( 0, wrf_err_message )
 140      write(wrf_err_message,*) 'WW_SPLIT:  IEVA = ', ieva
 141      CALL wrf_debug( 0, wrf_err_message )
 142      IF( config_flags%zadvect_implicit > 1 ) THEN
 143        write(wrf_err_message,*) 'WW_SPLIT:  WARNING: IEVA IS CONFIG TO ONLY RUN ON LAST RK-substep'
 144        CALL wrf_debug( 0, wrf_err_message )
 145      ENDIF
 146      write(wrf_err_message,*) 'WW_SPLIT:  dt = ', dt
 147      CALL wrf_debug( 0, wrf_err_message )
 148      write(wrf_err_message,*) 'WW_SPLIT:  alpha_max/min = ', alpha_max, alpha_min
 149      CALL wrf_debug( 0, wrf_err_message )
 150      write(wrf_err_message,*) '-- END WW_SPLIT INFO ----------------------------------------'
 151      CALL wrf_debug( 0, wrf_err_message )
 152      print_flag = .false.
 153    ENDIF
 154
 155 ! =========================================================================
 156 ! Loop limits
 157
 158    IF( ieva ) THEN   ! Implicit adaptive advection
 159
 160 !=========================================================================
 161 ! DEBUG CODE BLOCK
 162 !
 163 !    binavg(:) = 0.0
 164 !    crs(:)    = 0.0
 165 !    count(:)  = 0
 166 !
 167 !=========================================================================
 168
 169      DO j = j_start, j_end
 170        DO k = kts+1, ktf-1
 171          DO i = i_start, i_end
 172            wwE(i,k,j) = 0.
 173            wwI(i,k,j) = 0.
 174          ENDDO
 175        ENDDO
 176      ENDDO
 177
 178 ! Set upper and lower boundary values to be explicit advection
 179
 180      DO j = j_start, j_end
 181        DO i = i_start,i_end
 182
 183          wwE(i,1,j) = ww(i,    1,j)
 184          wwI(i,1,j) = 0.0
 185
 186          wwE(i,ktf,j) = ww(i,ktf,j)
 187          wwI(i,ktf,j) = 0.0
 188
 189        ENDDO
 190      ENDDO
 191
 192      DO j = j_start, j_end
 193
 194        DO k = 2, ktf-1
 195
 196          DO i = i_start,i_end
 197
 198 !===================================================================================================
 199 !
 200 ! Translating Shchepetkin (2015, hereafter S2015), all the volumes and surface areas cancel out (not obvious
 201 ! unless you read the definitions of U, V, and W right after equation 3.7; they are u*dy*dz, v*dx*dz, w*dx*dy)...
 202 ! You can write the algorithm using the adjusted maximum vertical courant number (alpha_max - eps*hor_div)
 203 ! versus the actual vertical courant number.  This code is based on the appendix C of Shchepetkin (2015).
 204 ! The CX and CY measure the maximum possible divergence for a zone, e.g., the Lipschlitz number.  If
 205 ! this number ~O(1), then the maximum allowed vertical courant number for the explicit advection is close
 206 ! to zeros, the scheme switches to fully implicit.  Ceps is a fudge factor for when the vertical and
 207 ! horizontal discretizations are of different order.  It is set to 0.9, which allows the some headroom for
 208 ! the maximum horizontal Courant number since we are using an RK3 integrator which has a larger stability
 209 ! region than the integrators used in Shchepetkin (2015).
 210 !
 211 !===================================================================================================
 212
 213 ! S2015 Lipshlitz courant number.  Sometimes it can be touchy so we use simpler estimate below
 214 !
 215 !          IF( ww(i,k,j) < 0.0 ) THEN    !! Omega has opposite sign of W
 216 !            cx = rdx*(max(u(i+1,k,j  ), 0.0) - min(u(i,k,j), 0.0))
 217 !            cy = rdy*(max(v(i,  k,j+1), 0.0) - min(v(i,k,j), 0.0))
 218 !          ELSE
 219 !            cx = rdx*(max(u(i+1,k-1,j  ), 0.0) - min(u(i,k-1,j), 0.0))
 220 !            cy = rdy*(max(v(i,  k-1,j+1), 0.0) - min(v(i,k-1,j), 0.0))
 221 !          ENDIF
 222
 223 !          cw_max = max(alpha_max - dt*Ceps*(cx + cy),0.0)                 ! S2015 adjusted max vertical courant
 224
 225 ! Compute a mean flow courant number to adjust the max vertical courant number
 226
 227            cx = 0.25*rdx*(u(i+1,k,j  ) + u(i,k,j) + u(i+1,k-1,j  ) + u(i,k-1,j))
 228            cy = 0.25*rdy*(v(i,  k,j+1) + v(i,k,j) + v(i  ,k-1,j+1) + v(i,k-1,j))
 229
 230            cw_max = max(alpha_max - dt*Ceps*sqrt(cx**2 + cy**2),0.0)       ! simpler adjusted max vertical courant
 231
 232            cr     = ww(i,k,j) * dt * rdnw(k) / (c1f(k)*mut(i,j)+c2f(k))    ! the actual vertical courant number
 233
 234            IF( cw_max > 0.0 ) THEN
 235
 236              wfrac   = 1.0   ! initialize wfrac
 237
 238              cw_max2 = cw_max**2
 239              cw_min  = cw_max*cmnx_ratio
 240              cw      = abs(cr)
 241
 242              if ( cw < cw_min ) then
 243                cff = cw_max2
 244              elseif ( cw < cutoff*cw_min ) then
 245                cff = cw_max2 + r4cmx*(cw-cw_min)**2
 246              else
 247                cff = cw_max*cw
 248              endif
 249
 250              wfrac = cw_max2 / cff
 251              wfrac = amax1(amin1(wfrac, 1.0), 0.0)
 252
 253              wwE(i,k,j) = ww(i,k,j) * wfrac
 254              wwI(i,k,j) = ww(i,k,j) * (1.0 - wfrac)
 255
 256 !=========================================================================
 257 ! DEBUG CODE BLOCK: compute distribution of wfrac
 258 !
 259 !            DO mm = 1, binsize
 260 !              IF( (mm-1)*binfrac / wfrac <= 1.0 .and. &
 261 !                      mm*binfrac / wfrac  > 1.0 ) THEN
 262 !                count(mm)  = count(mm) + 1
 263 !                binavg(mm) = binavg(mm) + wfrac
 264 !                crs(mm)    = crs(mm) + cw
 265 !              ENDIF
 266 !            ENDDO
 267 !
 268 !=========================================================================
 269
 270            ELSE   ! If cw_max <= 0 - then the horizontal flow is very fast
 271                   ! so do all vertical transport implicitly
 272
 273 !=========================================================================
 274 ! DEBUG CODE BLOCK:  write out information when fully implicit
 275 !
 276 !            write(wrf_err_message,*) '-- BEGIN WW_SPLIT ----------------------------------------'
 277 !            CALL wrf_debug(100, wrf_err_message )
 278 !            write(wrf_err_message,*) ' MAX VERTICAL CFL ~ 0: i,j,k, cw_max, dt*(cx+cy) ',i,j,k,cw_max,Ceps*dt*(cx+cy)
 279 !            CALL wrf_debug(100, wrf_err_message )
 280 !            write(wrf_err_message,*) '-- END WW_SPLIT ----------------------------------------'
 281 !            CALL wrf_debug( 100, wrf_err_message )
 282 !
 283 !=========================================================================
 284
 285              wwE(i,k,j) = 0.0
 286              wwI(i,k,j) = ww(i,k,j)
 287
 288            ENDIF
 289
 290          ENDDO
 291        ENDDO
 292      ENDDO
 293
 294    ELSE  ! NOT doing IEVA, pure explicit advection
 295
 296      DO j = j_start, j_end
 297        DO k = kts, ktf
 298          DO i = i_start,i_end
 299            wwE(i,k,j) = ww(i,k,j)
 300            wwI(i,k,j) = 0.0
 301          ENDDO
 302        ENDDO
 303      ENDDO
 304
 305    ENDIF
 306
 307 !=========================================================================
 308 ! DEBUG - print distribution of wfrac
 309 !
 310 !  write(wrf_err_message,*) '-- BEGIN WW_SPLIT ----------------------------------------'
 311 !  CALL wrf_debug( 100, wrf_err_message )
 312 !  DO mm = 1,11
 313 !    write(wrf_err_message,FMT='(a, 2(2x,f4.2),2x,i9,2(2x,f6.4))' ) &
 314 !     'WW_SPLIT: BIN #, count', (mm-1)*binfrac, mm*binfrac, count(mm), binavg(mm)/(count(mm)+1), crs(mm)/(count(mm)+1)
 315 !    CALL wrf_debug( 100, wrf_err_message )
 316 !  ENDDO
 317 !  write(wrf_err_message,*) '-- END WW_SPLIT ----------------------------------------'
 318 !  CALL wrf_debug( 100, wrf_err_message )
 319 !
 320 !=========================================================================
 321
 322 RETURN
 323 END SUBROUTINE WW_SPLIT
 324
 325
 326  SUBROUTINE TRIDIAG2D(a, b, c, r, bb, is, ie, istart, iend, ks, ke, kstart, kend)
 327 !--------------------------------------------------------------------------------------------------
 328 !     Solves for a vector u of length n the tridiagonal linear set      ! from numerical recipes
 329 !     m u = r, where a, b and c are the three main diagonals of matrix  !
 330 !     m(n,n), the other terms are 0. r is the right side vector.        !
 331 !
 332 !     a(1) and c(n) are not used in the calculation
 333 !--------------------------------------------------------------------------------------------------
 334
 335     integer,                              intent(in)    :: is, ie, ks, ke, istart, iend, kstart, kend
 336     real(kind=8), dimension(is:ie,ks:ke), intent(in)    :: a, b, c, r
 337     real(kind=8), dimension(is:ie,ks:ke), intent(inout) :: bb
 338 !
 339 ! Local Variables
 340 !
 341 !   real(kind=8), dimension(ks:ke) :: gam
 342 !   real(kind=8), dimension(is:ie) :: bet
 343     real(kind=8), dimension(ks:ke) :: gam
 344     real(kind=8), dimension(is:ie) :: bet
 345     integer                        :: i, k
 346 !
 347 ! Copy into temp storage...
 348 !
 349     DO i = istart, iend
 350
 351       bet(i) = b(i,kstart)
 352
 353     ENDDO
 354
 355 !   IF ( ANY( bet < 1.0e-15 ) ) THEN
 356 !     write(wrf_err_message,*) '----------------------------------------'
 357 !     CALL wrf_debug( 0, wrf_err_message )
 358 !     write(wrf_err_message,*) "TriDiag2D Error!!! Min value of diagonal is ", minval(bet)
 359 !     CALL wrf_debug( 0, wrf_err_message )
 360 !     write(wrf_err_message,*) '----------------------------------------'
 361 !     CALL wrf_debug( 0, wrf_err_message )
 362 !     STOP
 363 !   ENDIF
 364
 365     DO i = istart, iend
 366
 367        bb(i,kstart) = r(i,kstart)/bet(i)
 368
 369        DO k = kstart+1, kend
 370
 371          gam(k) = c(i,k-1) / bet(i)
 372          bet(i) = b(i,k) - a(i,k) * gam(k)
 373
 374          bb(i,k) = ( r(i,k) - a(i,k) * bb(i,k-1) ) / bet(i)
 375
 376        ENDDO
 377
 378        DO k = kend-1, kstart, -1
 379
 380          bb(i,k) = bb(i,k) - gam(k+1) * bb(i,k+1)
 381
 382        ENDDO
 383
 384      ENDDO
 385
 386  END SUBROUTINE TRIDIAG2D
 387 !--------------------------------------------------------------------------------------------------
 388  SUBROUTINE TRIDIAG(a, b, c, r, bb, ks, ke, kstart, kend)
 389 !--------------------------------------------------------------------------------------------------
 390 !     Solves for a vector u of length n the tridiagonal linear set      ! from numerical recipes
 391 !     m u = r, where a, b and c are the three main diagonals of matrix  !
 392 !     m(n,n), the other terms are 0. r is the right side vector.        !
 393 !
 394 !     a(1) and c(n) are not used in the calculation
 395 !--------------------------------------------------------------------------------------------------
 396     IMPLICIT NONE
 397
 398     integer,                        intent(in)  :: ks, ke, kstart, kend
 399     real(kind=4), dimension(ks:ke), intent(in ) :: a, b, c, r
 400     real(kind=4), dimension(ks:ke), intent(out) :: bb
 401
 402     !
 403     ! Local Variables
 404     !
 405
 406     real(kind=4), dimension(size(b,1)) :: gam
 407     real(kind=4)                       :: bet
 408     integer                            :: k
 409
 410     bet = b(kstart)
 411
 412     bb(kstart) = r(kstart)/bet
 413
 414     DO k = kstart+1, kend
 415
 416       gam(k) = c(k-1) / bet
 417       bet    = b(k) - a(k) * gam(k)
 418
 419       bb(k) = ( r(k) - a(k) * bb(k-1) ) / bet
 420
 421     ENDDO
 422
 423     DO k = kend-1, kstart, -1
 424       bb(k) = bb(k) - gam(k+1) * bb(k+1)
 425     ENDDO
 426
 427 END SUBROUTINE TRIDIAG
 428 !-------------------------------------------------------------------------------
 429 !
 430
 431 SUBROUTINE CALC_MUT_NEW( u, v, c1h, c2h,                     &
 432                          mut_old, muu, muv, mut_new,         &
 433                          dt, rdx, rdy, msftx, msfty,         &
 434                          msfux, msfuy, msfvx, msfvx_inv,     &
 435                          msfvy, rdnw,                        &
 436                          ids, ide, jds, jde, kds, kde,       &
 437                          ims, ime, jms, jme, kms, kme,       &
 438                          its, ite, jts, jte, kts, kte    )
 439
 440    IMPLICIT NONE
 441
 442    ! Input data
 443
 444
 445    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 446                                  ims, ime, jms, jme, kms, kme, &
 447                                  its, ite, jts, jte, kts, kte
 448
 449    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: u, v
 450    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mut_old, muu, muv, &
 451                                                             msftx, msfty, &
 452                                                             msfux, msfuy, &
 453                                                             msfvx, msfvy, &
 454                                                             msfvx_inv
 455    REAL , DIMENSION( kms:kme ) , INTENT(IN   ) :: rdnw
 456    REAL , DIMENSION( kms:kme ) , INTENT(IN   ) :: c1h, c2h
 457
 458    REAL , DIMENSION( ims:ime, jms:jme ) , INTENT(OUT  ) :: mut_new
 459    REAL , INTENT(IN)  :: dt, rdx, rdy
 460
 461    ! Local data
 462
 463    INTEGER :: i, j, k, ktf, i_start, i_end, j_start, j_end
 464    REAL , DIMENSION( its-1:ite ) :: dmdt
 465    REAL , DIMENSION( its-1:ite, kts:kte ) :: divv
 466    REAL :: bigerr
 467
 468 !<DESCRIPTION>
 469 !
 470 !  The algorithm integrates the continuity equation through the column
 471 !  to find a new column mass needed for IEVA calculations
 472 !
 473 !</DESCRIPTION>
 474
 475
 476     ktf=MIN(kte,kde-1)
 477
 478     i_start = its-1
 479     i_end   = MIN(ite,ide-1)
 480     j_start = jts-1
 481     j_end   = MIN(jte,jde-1)
 482
 483
 484     DO j = j_start, j_end
 485
 486       DO i = i_start, i_end
 487         dmdt(i) = 0.
 488       ENDDO
 489
 490 !       Comments on the modifications for map scale factors
 491 !       ADT eqn 47 / my (putting rho -> 'mu') is:
 492 !       (1/my) partial d mu/dt = -mx partial d/dx(mu u/my)
 493 !                                -mx partial d/dy(mu v/mx)
 494 !                                -partial d/dz(mu w/my)
 495 !
 496 !       Using nu instead of z the last term becomes:
 497 !                                -partial d/dnu((c1(k)*mu(dnu/dt))/my)
 498 !
 499 !       Integrating with respect to nu over ALL levels, with dnu/dt=0 at top
 500 !       and bottom, the last term becomes = 0
 501 !
 502 !       Integral|bot->top[(1/my) partial d mu/dt]dnu =
 503 !       Integral|bot->top[-mx partial d/dx(mu u/my)
 504 !                         -mx partial d/dy(mu v/mx)]dnu
 505 !
 506 !       muu='mu'[on u]/my, muv='mu'[on v]/mx
 507 !       (1/my) partial d mu/dt is independent of nu
 508 !         => LHS = Integral|bot->top[con]dnu = conservation*(-1) = -dmdt
 509 !
 510 !         => dmdt = mx*Integral|bot->top[partial d/dx(mu u/my) +
 511 !                                        partial d/dy(mu v/mx)]dnu
 512 !         => dmdt = sum_bot->top[divv]
 513 !       where
 514 !         divv=mx*[partial d/dx(mu u/my) + partial d/dy(mu v/mx)]*delta nu
 515
 516       DO k=kts,ktf
 517         DO i = i_start, i_end
 518
 519         divv(i,k) = (msftx(i,j) / rdnw(k)) *                                                                        &
 520             (rdx*((c1h(k)*muu(i+1,j)+c2h(k))*u(i+1,k,j)/msfuy(i+1,j)-(c1h(k)*muu(i,j)+c2h(k))*u(i,k,j)/msfuy(i,j))  &
 521             +rdy*((c1h(k)*muv(i,j+1)+c2h(k))*v(i,k,j+1)*msfvx_inv(i,j+1)-(c1h(k)*muv(i,j)+c2h(k))*v(i,k,j)*msfvx_inv(i,j)) )
 522
 523         dmdt(i) = dmdt(i) + divv(i,k)
 524
 525        ENDDO
 526       ENDDO
 527
 528 !       Further map scale factor notes:
 529 !       Now integrate from bottom to top, level by level:
 530 !       mu dnu/dt/my [k+1] = mu dnu/dt/my [k] + [-(1/my) partial d mu/dt
 531 !                           -mx partial d/dx(mu u/my)
 532 !                           -mx partial d/dy(mu v/mx)]*dnu[k->k+1]
 533 !       ww [k+1] = ww [k] -(1/my) partial d mu/dt * dnu[k->k+1] - divv[k]
 534 !                = ww [k] -dmdt * dnw[k] - divv[k]
 535
 536
 537       DO i = i_start, i_end
 538
 539          mut_new(i,j) = mut_old(i,j) - dt*dmdt(i)
 540
 541 ! DEBUG code.
 542 !          bigerr = (mut_new(i,j) / mut_old(i,j))
 543 !
 544 !          IF( abs(bigerr) > 2.0 .or. mut_old(i,j) < 0.0 ) THEN
 545 !             write(wrf_err_message,*) '------------- CALC_MU_TMP ---------------------------'
 546 !             CALL wrf_debug( 0, wrf_err_message )
 547 !             write(wrf_err_message,*) 'MU_ERR: ', bigerr, mut_new(i,j), mut_old(i,j), dt, dt*dmdt(i)
 548 !             CALL wrf_debug( 0, wrf_err_message )
 549 !         ENDIF
 550
 551       ENDDO   ! END I
 552
 553     ENDDO     ! END J
 554
 555
 556 RETURN
 557 END SUBROUTINE CALC_MUT_NEW
 558 !-------------------------------------------------------------------------------
 559 !
 560 SUBROUTINE advect_ph_implicit( ph, pho, tendency, phb,        &
 561                                ru, rv, wwE, wwI, w,           &
 562                                c1, c2,                        &
 563                                mut, config_flags,             &
 564                                msfux, msfuy, msfvx, msfvy,    &
 565                                msftx, msfty,                  &
 566                                fzm, fzp,                      &
 567                                dt_rk,                         &
 568                                rdx, rdy, rdzw,                &
 569                                ids, ide, jds, jde, kds, kde,  &
 570                                ims, ime, jms, jme, kms, kme,  &
 571                                its, ite, jts, jte, kts, kte  )
 572
 573    IMPLICIT NONE
 574
 575    ! Input data
 576
 577    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
 578
 579    INTEGER ,                 INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 580                                               ims, ime, jms, jme, kms, kme, &
 581                                               its, ite, jts, jte, kts, kte
 582
 583    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: ph,    &
 584                                                                       pho,   &
 585                                                                       phb,   &
 586                                                                       ru,    &
 587                                                                       rv,    &
 588                                                                       wwI,   &
 589                                                                       wwE,   &
 590                                                                       w
 591
 592    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mut
 593    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
 594
 595    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,  &
 596                                                                     msfuy,  &
 597                                                                     msfvx,  &
 598                                                                     msfvy,  &
 599                                                                     msftx,  &
 600                                                                     msfty
 601
 602    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm,  &
 603                                                                   fzp,  &
 604                                                                   rdzw, &
 605                                                                   c1,   &
 606                                                                   c2
 607
 608    REAL ,                                        INTENT(IN   ) :: rdx,  &
 609                                                                   rdy
 610    REAL ,                                        INTENT(IN   ) :: dt_rk
 611
 612    ! Local data
 613
 614    INTEGER :: i, j, k, itf, jtf, ktf
 615    INTEGER :: i_start, i_end, j_start, j_end
 616    INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
 617    INTEGER :: jmin, jmax, jp, jm, imin, imax, im, ip
 618    REAL    :: wiL, wiR, wiC, weL, weR, dz
 619
 620    REAL(KIND=8), DIMENSION(its:ite,kts:kte) :: at, bt, ct, rt, bb
 621    REAL(KIND=8), DIMENSION(its:ite)         :: btmp(ims:ime)
 622
 623    LOGICAL :: specified
 624    INTEGER :: iup, jup, kup, idn, jdn, kdn, kp1, km1, valid_ik
 625
 626    ktf     = kte-1
 627    i_start = its
 628    i_end   = MIN(ite,ide-1)
 629    j_start = jts
 630    j_end   = MIN(jte,jde-1)
 631
 632 ! Vertical loop to do implicit advection.....
 633 ! Different from the other schemes - the advection here is not in flux form...
 634
 635    DO j = j_start, j_end
 636
 637      DO k = kts+1, ktf
 638        DO i = i_start, i_end
 639
 640 ! We are choosing to do the vertical advection here the same as in RHS_PH,
 641 ! e.g., 2nd order box averaging instead of the usual upwinding...
 642 ! Divide my mu(i,j) to make sure the diagnonal term is ~O(1)
 643
 644 ! Upwind
 645          wiC     = 0.5*wwI(i,k,j)*(rdzw(k-1)+rdzw(k)) * msfty(i,j) / (c1(k)*mut(i,j)+c2(k))
 646          at(i,k) = - dt_rk*max(wiC,0.0)
 647          ct(i,k) =   dt_rk*min(wiC,0.0)
 648          btmp(i) =   - at(i,k) - ct(i,k)
 649
 650 ! 2nd order centered for both implicit piece of advection
 651 !          wiL     = 0.5*(wwI(i,k-1,j)+wwI(i,k,j)) * rdzw(k-1) * msfty(i,j) / mut(i,j)
 652 !          wiR     = 0.5*(wwI(i,k+1,j)+wwI(i,k,j)) * rdzw(k)   * msfty(i,j) / mut(i,j)
 653 !          at(i,k) = - dt_rk*wiL*fzp(k)
 654 !          ct(i,k) =   dt_rk*wiR*fzm(k)
 655 !          btmp(i) =   - at(i,k) - ct(i,k)
 656
 657 ! Diagonal term
 658
 659          bt(i,k) =   1.0 + btmp(i)
 660
 661 ! 2nd order centered for explict piece of advection  (TURNED OFF AND COMPUTED IN RHS_PH)
 662 !
 663 !          weL   = -dt_rk*fzp(k)*0.5*(wwE(i,k-1,j)+wwE(i,k,j)) * rdzw(k-1) * msfty(i,j) / mut(i,j)
 664 !          weR   =  dt_rk*fzm(k)*0.5*(wwE(i,k+1,j)+wwE(i,k,j)) * rdzw(k)   * msfty(i,j) / mut(i,j)
 665
 666 ! Eq here is more complicated because the solution to the tridiagonal is in phi', not total
 667 ! phi.  So we retain the vertical advection of phb on RHS.
 668
 669          rt(i,k) = tendency(i,k,j) * dt_rk * msfty(i,j) / (c1(k)*mut(i,j)+c2(k))        &
 670 !                  - weL*ph (i,k-1,j)     + (weL + weR)*ph (i,k,j) - weR*ph (i,k+1,j)     &
 671 !                  - weL*phb(i,k-1,j)     + (weL + weR)*phb(i,k,j) - weR*phb(i,k+1,j)     &
 672                  - at(i,k)*pho(i,k-1,j) -     btmp(i)*pho(i,k,j) - ct(i,k)*pho(i,k+1,j) &
 673                  - at(i,k)*phb(i,k-1,j) -     btmp(i)*phb(i,k,j) - ct(i,k)*phb(i,k+1,j)
 674
 675 ! BC's need because tendencies of phi, at kts and kte are not zero, so add them to RHS of Tridiagonal eqs
 676 ! NOT USED!
 677 !          IF( k == kts+1 ) THEN
 678 !             rt(i,k) = rt(i,k) - at(i,k)*(ph(i,k-1,j) - pho(i,k-1,j))       ! lower boundary condition
 679 !          ENDIF
 680 !
 681 !          IF( k == ktf ) THEN    ! upper boundary condition
 682 !           rt(i,k)   = rt(i,k) - ct(i,k) * (ph(i,k+1,j) - pho(i,k+1,j))
 683 !         ENDIF
 684
 685        ENDDO
 686      ENDDO
 687
 688      CALL tridiag2D(at, bt, ct, rt, bb, its, ite, i_start, i_end, kts, kte, kts+1, ktf)
 689
 690      DO k = kts+1, ktf
 691        DO i = i_start, i_end
 692
 693          tendency(i,k,j) = sngl(bb(i,k)) * (c1(k)*mut(i,j)+c2(k)) / dt_rk / msfty(i,j)
 694
 695        ENDDO
 696      ENDDO
 697
 698    ENDDO
 699
 700
 701 RETURN
 702 END SUBROUTINE advect_ph_implicit
 703 !-------------------------------------------------------------------------------
 704 !
 705 SUBROUTINE advect_s_implicit( s, s_old, tendency,            &
 706                               ru, rv, rom,                   &
 707                               c1, c2,                        &
 708                               mut_old, mut, mut_new,         &
 709                               config_flags,                  &
 710                               msfux, msfuy, msfvx, msfvy,    &
 711                               msftx, msfty,                  &
 712                               fzm, fzp,                      &
 713                               dt_rk,                         &
 714                               rdx, rdy, rdzw,                &
 715                               ids, ide, jds, jde, kds, kde,  &
 716                               ims, ime, jms, jme, kms, kme,  &
 717                               its, ite, jts, jte, kts, kte  )
 718
 719    IMPLICIT NONE
 720
 721    ! Input data
 722
 723    TYPE(grid_config_rec_type), INTENT(IN) :: config_flags
 724
 725    INTEGER ,                   INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
 726                                              ims, ime, jms, jme, kms, kme, &
 727                                              its, ite, jts, jte, kts, kte
 728
 729    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: s,     &
 730                                                                    s_old, &
 731                                                                    ru,    &
 732                                                                    rv,    &
 733                                                                    rom
 734 !-------------------------------------------------------------------------------
 735 ! LWJ: definitions of various column masses
 736 ! mut     ==> current column mass from sub-step
 737 ! mut_old ==> "n" time level column mass
 738 ! mut_new ==> "n+1*" estimated column needed for dynamical variables where we dont
 739 !            have time-avg column mass.  For scalars (not theta) mut_new == mut
 740
 741    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: mut, mut_old, mut_new
 742
 743 !-------------------------------------------------------------------------------
 744
 745    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
 746
 747    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,  &
 748                                                                     msfuy,  &
 749                                                                     msfvx,  &
 750                                                                     msfvy,  &
 751                                                                     msftx,  &
 752                                                                     msfty
 753
 754    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm,  &
 755                                                                   fzp,  &
 756                                                                   rdzw, &
 757                                                                   c1,   &
 758                                                                   c2
 759
 760    REAL ,                                        INTENT(IN   ) :: rdx,  &
 761                                                                   rdy
 762    REAL ,                                        INTENT(IN   ) :: dt_rk
 763
 764 ! Local data
 765
 766    INTEGER :: i, j, k, itf, jtf, ktf
 767    INTEGER :: i_start, i_end, j_start, j_end
 768    INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
 769    INTEGER :: jmin, jmax, jp, jm, imin, imax, im, ip
 770    INTEGER :: jp1, jp0, jtmp
 771
 772    INTEGER :: horz_order, vert_order
 773
 774    REAL    :: mrdx, mrdy, ub, vb, uw, vw, dvm, dvp, wiL, wiR, dz
 775
 776    REAL(KIND=8), DIMENSION(its:ite,kts:kte) :: at, bt, ct, rt, bb
 777    REAL(KIND=8), DIMENSION(its:ite)         :: btmp
 778
 779    INTEGER :: iup, jup, kup, idn, jdn, kdn, kp1, km1
 780    INTEGER :: valid_ik, skip
 781
 782 !
 783
 784    ktf     = MIN(kte,kde-1)
 785
 786 !-- horizontal loop bounds copied from module_advect_em.F/advect_scalar
 787
 788    i_start = its
 789    i_end   = MIN(ite,ide-1)
 790    j_start = jts
 791    j_end   = MIN(jte,jde-1)
 792
 793 ! J outer loop to do implicit advection.....
 794
 795    DO j = j_start, j_end
 796
 797      DO k = kts, ktf
 798
 799        km1 = k - 1
 800        kp1 = k + 1
 801        IF( k .eq. ktf ) kp1 = ktf
 802        IF( k .eq. kts ) km1 = kts
 803
 804        DO i = i_start, i_end
 805
 806          ! Redefine mass fluxes with new temporary column mass
 807
 808          wiL   = rom(i,k,  j) * rdzw(k) / (c1(k)*mut_new(i,j)+c2(k))
 809          wiR   = rom(i,k+1,j) * rdzw(k) / (c1(k)*mut_new(i,j)+c2(k))
 810
 811          at(i,k) = - dt_rk*max(wiL,0.0)
 812          ct(i,k) =   dt_rk*min(wiR,0.0)
 813          btmp(i) =   dt_rk*(max(wiR,0.0) - min(wiL,0.0))
 814          bt(i,k) = 1.0 + btmp(i)
 815          rt(i,k) = dt_rk*tendency(i,k,j)  &
 816                    - (c1(k)*mut_old(i,j)+c2(k))*(at(i,k)*s_old(i,km1,j) + btmp(i)*s_old(i,k,j) + ct(i,k)*s_old(i,kp1,j))
 817        ENDDO
 818      ENDDO
 819
 820      CALL tridiag2D(at, bt, ct, rt, bb, its, ite, i_start, i_end, kts, kte, kts, ktf)
 821
 822      DO k = kts, ktf
 823        DO i = i_start, i_end
 824
 825          ! Rescale tendency with old column mass for consistency in update step.
 826
 827          tendency(i,k,j) = sngl(bb(i,k)) / dt_rk
 828
 829        ENDDO
 830      ENDDO
 831
 832    ENDDO  ! J-LOOP
 833
 834 RETURN
 835 END SUBROUTINE advect_s_implicit
 836 !-------------------------------------------------------------------------------
 837 !
 838 SUBROUTINE advect_u_implicit( u, u_old, tendency,            &
 839                               ru, rv, rom,                   &
 840                               c1, c2,                        &
 841                               muu_old, muu, muu_new,         &
 842                               config_flags,                  &
 843                               msfux, msfuy, msfvx, msfvy,    &
 844                               msftx, msfty,                  &
 845                               fzm, fzp,                      &
 846                               dt_rk,                         &
 847                               rdx, rdy, rdzw,                &
 848                               ids, ide, jds, jde, kds, kde,  &
 849                               ims, ime, jms, jme, kms, kme,  &
 850                               its, ite, jts, jte, kts, kte  )
 851
 852    IMPLICIT NONE
 853
 854 ! Input data
 855
 856    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
 857
 858    INTEGER ,                 INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 859                                               ims, ime, jms, jme, kms, kme, &
 860                                               its, ite, jts, jte, kts, kte
 861
 862    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: u,     &
 863                                                                       u_old, &
 864                                                                       ru,    &
 865                                                                       rv,    &
 866                                                                       rom
 867
 868 !-------------------------------------------------------------------------------
 869 ! LWJ: definitions of various column masses
 870 ! mut    ==> current column mass from sub-step
 871 ! mu_old ==> "n" time level column mass
 872 ! mu_new ==> "n+1*" estimated column needed for dynamical variables where we dont
 873 !            have time-avg column mass.  For scalars (not theta) mut_new == mut
 874
 875    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: muu, muu_old, muu_new
 876
 877 !-------------------------------------------------------------------------------
 878
 879    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
 880
 881    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,  &
 882                                                                     msfuy,  &
 883                                                                     msfvx,  &
 884                                                                     msfvy,  &
 885                                                                     msftx,  &
 886                                                                     msfty
 887
 888    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm,  &
 889                                                                   fzp,  &
 890                                                                   rdzw, &
 891                                                                   c1,   &
 892                                                                   c2
 893
 894    REAL ,                                        INTENT(IN   ) :: rdx,  &
 895                                                                   rdy
 896    REAL ,                                        INTENT(IN   ) :: dt_rk
 897
 898    ! Local data
 899
 900    INTEGER :: i, j, k, itf, jtf, ktf
 901    INTEGER :: i_start, i_end, j_start, j_end
 902    INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
 903    INTEGER :: jmin, jmax, jp, jm, imin, imax, im, ip
 904    INTEGER :: jp1, jp0, jtmp
 905
 906    INTEGER :: horz_order, vert_order
 907
 908    REAL    :: wiL, wiR, dz
 909
 910    REAL(KIND=8), DIMENSION(its:ite,kts:kte) :: at, bt, ct, rt, bb
 911    REAL(KIND=8), DIMENSION(its:ite)         :: btmp
 912
 913    LOGICAL :: specified
 914    INTEGER :: iup, jup, kup, idn, jdn, kdn, kp1, km1
 915    INTEGER :: valid_ik
 916
 917    ktf     = MIN(kte,kde-1)
 918
 919 !-- loop bounds copied from module_advect_em.F/advect_u
 920
 921    specified = .false.
 922    if(config_flags%specified .or. config_flags%nested) specified = .true.
 923
 924    i_start = its
 925    i_end   = ite
 926    j_start = jts
 927    j_end   = min(jte,jde-1)
 928
 929 !  IF ( config_flags%open_xs ) i_start = MAX(ids+1,its)
 930 !  IF ( config_flags%open_xe ) i_end   = MIN(ide-1,ite)
 931
 932    IF ( config_flags%open_ys .or. specified ) i_start = MAX(ids+1,its)
 933    IF ( config_flags%open_ye .or. specified ) i_end   = MIN(ide-1,ite)
 934    IF ( config_flags%periodic_x ) i_start = its
 935    IF ( config_flags%periodic_x ) i_end = ite
 936
 937 ! Vertical loop to do implicit advection.....
 938
 939    DO j = j_start, j_end
 940
 941      DO k = kts, ktf
 942
 943        km1 = k - 1
 944        kp1 = k + 1
 945        IF( k .eq. ktf ) kp1 = ktf
 946        IF( k .eq. kts ) km1 = kts
 947
 948        DO i = i_start, i_end
 949
 950 ! SO IMPORTANT to * 1/DZ HERE CAUSE IT CHANGES Wi* SIGN TO BE CORRECT FOR UPWINDING!!!!
 951          wiL   = 0.5*(rom(i-1,k,  j)+rom(i,k,  j)) * rdzw(k) * msfuy(i,j) / (c1(k)*muu_new(i,j)+c2(k))
 952          wiR   = 0.5*(rom(i-1,k+1,j)+rom(i,k+1,j)) * rdzw(k) * msfuy(i,j) / (c1(k)*muu_new(i,j)+c2(k))
 953
 954          at(i,k) = - dt_rk*max(wiL,0.0)
 955          ct(i,k) =   dt_rk*min(wiR,0.0)
 956          btmp(i) =   dt_rk*(max(wiR,0.0) - min(wiL,0.0))
 957          bt(i,k) = 1.0 + btmp(i)
 958          rt(i,k) = dt_rk*tendency(i,k,j)*msfuy(i,j)  &
 959                  - (c1(k)*muu_old(i,j)+c2(k))*(at(i,k)*u_old(i,km1,j) + btmp(i)*u_old(i,k,j) + ct(i,k)*u_old(i,kp1,j))
 960
 961        ENDDO
 962      ENDDO
 963
 964      CALL tridiag2D(at, bt, ct, rt, bb, its, ite, i_start, i_end, kts, kte, kts, ktf)
 965
 966      DO k = kts, ktf
 967        DO i = i_start, i_end
 968
 969          tendency(i,k,j) = sngl(bb(i,k)) / dt_rk / msfuy(i,j)
 970
 971        ENDDO
 972      ENDDO
 973
 974    ENDDO ! J-LOOP
 975
 976 RETURN
 977 END SUBROUTINE advect_u_implicit
 978 !-------------------------------------------------------------------------------
 979 !
 980 SUBROUTINE advect_v_implicit( v, v_old, tendency,            &
 981                               ru, rv, rom,                   &
 982                               c1, c2,                        &
 983                               muv_old, muv, muv_new,         &
 984                               config_flags,                  &
 985                               msfux, msfuy, msfvx, msfvy,    &
 986                               msftx, msfty,                  &
 987                               fzm, fzp,                      &
 988                               dt_rk,                         &
 989                               rdx, rdy, rdzw,                &
 990                               ids, ide, jds, jde, kds, kde,  &
 991                               ims, ime, jms, jme, kms, kme,  &
 992                               its, ite, jts, jte, kts, kte  )
 993
 994    IMPLICIT NONE
 995
 996    ! Input data
 997
 998    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
 999
1000    INTEGER ,                 INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1001                                               ims, ime, jms, jme, kms, kme, &
1002                                               its, ite, jts, jte, kts, kte
1003
1004    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: v,     &
1005                                                                       v_old, &
1006                                                                       ru,    &
1007                                                                       rv,    &
1008                                                                       rom
1009
1010 !-------------------------------------------------------------------------------
1011 ! LWJ: definitions of various column masses
1012 ! mut    ==> current column mass from sub-step
1013 ! mu_old ==> "n" time level column mass
1014 ! mu_new ==> "n+1*" estimated column needed for dynamical variables where we dont
1015 !            have time-avg column mass.  For scalars (not theta) mut_new == mut
1016
1017    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: muv, muv_old, muv_new
1018
1019 !-------------------------------------------------------------------------------
1020
1021    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
1022
1023    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,  &
1024                                                                     msfuy,  &
1025                                                                     msfvx,  &
1026                                                                     msfvy,  &
1027                                                                     msftx,  &
1028                                                                     msfty
1029
1030    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm,  &
1031                                                                   fzp,  &
1032                                                                   rdzw, &
1033                                                                   c1,   &
1034                                                                   c2
1035
1036    REAL ,                                        INTENT(IN   ) :: rdx,  &
1037                                                                   rdy
1038    REAL ,                                        INTENT(IN   ) :: dt_rk
1039
1040 ! Local data
1041
1042    INTEGER :: i, j, k, itf, jtf, ktf
1043    INTEGER :: i_start, i_end, j_start, j_end
1044    INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
1045    INTEGER :: jmin, jmax, jp, jm, imin, imax, im, ip
1046    INTEGER :: jp1, jp0, jtmp
1047
1048    INTEGER :: horz_order, vert_order
1049
1050    REAL    :: mrdx, mrdy, ub, vb, uw, vw, dvm, dvp, wiL, wiR, dz
1051
1052    REAL(KIND=8), DIMENSION(its:ite,kts:kte) :: at, bt, ct, rt, bb
1053    REAL(KIND=8), DIMENSION(its:ite)         :: btmp
1054
1055    LOGICAL :: specified
1056    INTEGER :: iup, jup, kup, idn, jdn, kdn, kp1, km1
1057    INTEGER :: valid_ik
1058
1059    ktf     = MIN(kte,kde-1)
1060
1061 !-- loop bounds copied from module_advect_em.F/advect_v
1062
1063    specified = .false.
1064    if(config_flags%specified .or. config_flags%nested) specified = .true.
1065
1066    i_start = its
1067    i_end   = MIN(ite,ide-1)
1068    j_start = jts
1069    j_end   = jte
1070
1071 ! Polar boundary conditions are like open or specified
1072 ! We don't want to calculate vertical v tendencies at the N or S pole
1073    IF ( config_flags%open_ys .or. specified .or. config_flags%polar ) j_start = MAX(jds+1,jts)
1074    IF ( config_flags%open_ye .or. specified .or. config_flags%polar ) j_end   = MIN(jde-1,jte)
1075
1076 ! Main Loop
1077
1078    DO j = j_start, j_end
1079
1080      DO k = kts, ktf
1081
1082        km1 = k - 1
1083        kp1 = k + 1
1084        IF( k .eq. ktf ) kp1 = ktf
1085        IF( k .eq. kts ) km1 = kts
1086
1087        DO i = i_start, i_end
1088
1089
1090          wiL   = 0.5*(rom(i,k,  j-1)+rom(i,k,  j)) * rdzw(k) * msfvy(i,j) / (c1(k)*muv_new(i,j)+c2(k))
1091          wiR   = 0.5*(rom(i,k+1,j-1)+rom(i,k+1,j)) * rdzw(k) * msfvy(i,j) / (c1(k)*muv_new(i,j)+c2(k))
1092
1093          at(i,k) = - dt_rk*max(wiL,0.0)
1094          ct(i,k) =   dt_rk*min(wiR,0.0)
1095          btmp(i) =   dt_rk*(max(wiR,0.0) - min(wiL,0.0))
1096          bt(i,k) = 1.0 + btmp(i)
1097          rt(i,k) = dt_rk*tendency(i,k,j) * msfvx(i,j) &
1098                  - (c1(k)*muv_old(i,j)+c2(k))*(at(i,k)*v_old(i,km1,j) + btmp(i)*v_old(i,k,j) + ct(i,k)*v_old(i,kp1,j))
1099
1100        ENDDO
1101      ENDDO
1102
1103      CALL tridiag2D(at, bt, ct, rt, bb, its, ite, i_start, i_end, kts, kte, kts, ktf)
1104
1105      DO k = kts, ktf
1106        DO i = i_start, i_end
1107
1108        tendency(i,k,j) = sngl(bb(i,k)) / dt_rk / msfvx(i,j)
1109
1110        ENDDO
1111      ENDDO
1112
1113    ENDDO
1114
1115 RETURN
1116 END SUBROUTINE advect_v_implicit
1117
1118 !-------------------------------------------------------------------------------
1119 !
1120 SUBROUTINE advect_w_implicit( w, w_old, tendency,            &
1121                               utend, vtend, ht, rom,         &
1122                               ph_new, ph_old, ph_tend,       &
1123                               c1, c2,                        &
1124                               cf1, cf2, cf3,                 &
1125                               mut_old, mut, mut_new,         &
1126                               config_flags,                  &
1127                               msfux, msfuy, msfvx, msfvy,    &
1128                               msftx, msfty,                  &
1129                               fzm, fzp,                      &
1130                               dt_rk,                         &
1131                               rdx, rdy, rdzu,                &
1132                               ids, ide, jds, jde, kds, kde,  &
1133                               ims, ime, jms, jme, kms, kme,  &
1134                               its, ite, jts, jte, kts, kte  )
1135
1136    IMPLICIT NONE
1137
1138 ! Input data
1139
1140    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
1141
1142    INTEGER ,                 INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1143                                               ims, ime, jms, jme, kms, kme, &
1144                                               its, ite, jts, jte, kts, kte
1145
1146    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: w,     &
1147                                                                       w_old, &
1148                                                                       utend, &
1149                                                                       vtend, &
1150                                                                      ph_old, &
1151                                                                      ph_new, &
1152                                                                     ph_tend, &
1153                                                                       rom
1154
1155 !-------------------------------------------------------------------------------
1156 ! LWJ: definitions of various column masses
1157 ! mut    ==> current column mass from sub-step
1158 ! mu_old ==> "n" time level column mass
1159 ! mu_new ==> "n+1*" estimated column needed for dynamical variables where we dont
1160 !            have time-avg column mass.  For scalars (not theta) mut_new == mut
1161
1162    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: mut, mut_old, mut_new
1163
1164 !-------------------------------------------------------------------------------
1165
1166    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
1167
1168    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,  &
1169                                                                     msfuy,  &
1170                                                                     msfvx,  &
1171                                                                     msfvy,  &
1172                                                                     msftx,  &
1173                                                                     msfty,  &
1174                                                                     ht
1175
1176    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm,  &
1177                                                                   fzp,  &
1178                                                                   rdzu, &
1179                                                                   c1,   &
1180                                                                   c2
1181
1182    REAL ,                                        INTENT(IN   ) :: rdx,  &
1183                                                                   rdy,  &
1184                                                                   cf1,  &
1185                                                                   cf2,  &
1186                                                                   cf3
1187
1188    REAL ,                                        INTENT(IN   ) :: dt_rk
1189
1190 ! Local data
1191
1192    INTEGER :: i, j, k, itf, jtf, ktf
1193    INTEGER :: i_start, i_end, j_start, j_end
1194    INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f
1195    REAL    :: wiL, wiR, dz, dw
1196
1197    REAL(KIND=8), DIMENSION(its:ite,kts:kte) :: at, bt, ct, rt, bb
1198    REAL(KIND=8), DIMENSION(its:ite)         :: btmp
1199
1200    LOGICAL :: specified
1201    INTEGER :: iup, jup, kup, idn, jdn, kdn, kp1, km1
1202    INTEGER :: valid_ik
1203
1204    ktf     = MIN(kte,kde-1)
1205    i_start = its
1206    i_end   = MIN(ite,ide-1)
1207    j_start = jts
1208    j_end   = MIN(jte,jde-1)
1209
1210 ! Main loop
1211
1212    DO j = j_start, j_end
1213
1214      DO k = kts+1, ktf
1215        DO i = i_start, i_end
1216
1217 ! SO IMPORTANT to * 1/DZ HERE CAUSE IT CHANGES Wi* SIGN TO BE CORRECT FOR UPWINDING!!!!
1218          wiL   = 0.5*(rom(i,k-1,j)+rom(i,k,j)) * rdzu(k) * msfty(i,j) / (c1(k)*mut_new(i,j)+c2(k))
1219          wiR   = 0.5*(rom(i,k+1,j)+rom(i,k,j)) * rdzu(k) * msfty(i,j) / (c1(k)*mut_new(i,j)+c2(k))
1220
1221          at(i,k) = - dt_rk*max(wiL,0.0)
1222          ct(i,k) =   dt_rk*min(wiR,0.0)
1223          btmp(i) =   dt_rk*(max(wiR,0.0) - min(wiL,0.0))
1224          bt(i,k) =   1.0 + btmp(i)
1225          rt(i,k) = dt_rk*tendency(i,k,j) * msfty(i,j)  &
1226                  - (c1(k)*mut_old(i,j)+c2(k))*(at(i,k)*w_old(i,k-1,j) + btmp(i)*w_old(i,k,j) + ct(i,k)*w_old(i,k+1,j))
1227
1228 ! BC at lower boundary needed because (W(i,kts,j)  .ne. 0) and value need for RHS of tridiagonal at kts+1
1229 ! Since we already included w_old, just need to compute dw increment from UTEND, VTEND at kts
1230
1231         IF( k == kts+1 ) THEN
1232          dw =  msfty(i,j)*.5*rdy*(                                       &
1233                            (ht(i,j+1)-ht(i,j  ))                         &
1234           *(cf1*vtend(i,1,j+1)+cf2*vtend(i,2,j+1)+cf3*vtend(i,3,j+1))    &
1235                           +(ht(i,j  )-ht(i,j-1))                         &
1236           *(cf1*vtend(i,1,j  )+cf2*vtend(i,2,j  )+cf3*vtend(i,3,j  ))  ) &
1237                  +msftx(i,j)*.5*rdx*(                                    &
1238                            (ht(i+1,j)-ht(i,j  ))                         &
1239           *(cf1*utend(i+1,1,j)+cf2*utend(i+1,2,j)+cf3*utend(i+1,3,j))    &
1240                           +(ht(i,j  )-ht(i-1,j))                         &
1241           *(cf1*utend(i  ,1,j)+cf2*utend(i  ,2,j)+cf3*utend(i  ,3,j))  )
1242
1243           rt(i,k) = rt(i,k) - (c1(k)*mut(i,j)+c2(k))*at(i,k)*dt_rk*dw
1244         ENDIF
1245
1246 ! W-increment at upper boundary is computed from phi
1247
1248         IF( k == ktf ) THEN    ! upper boundary condition
1249           dw = msfty(i,j)*(  (ph_new(i,k+1,j)-ph_old(i,k+1,j))/dt_rk     &
1250                             - ph_tend(i,k+1,j)/(c1(k)*mut(i,j)+c2(k))/g)
1251
1252           rt(i,k)   = rt(i,k) - (c1(k)*mut(i,j)+c2(k))*ct(i,k) * (dw - w_old(i,k+1,j))
1253         ENDIF
1254
1255        ENDDO
1256      ENDDO
1257
1258      CALL tridiag2D(at, bt, ct, rt, bb, its, ite, i_start, i_end, kts, kte, kts+1, ktf)
1259
1260      DO k = kts+1, ktf
1261        DO i = i_start, i_end
1262
1263          tendency(i,k,j) = sngl(bb(i,k)) / dt_rk / msfty(i,j)
1264
1265        ENDDO
1266      ENDDO
1267
1268    ENDDO
1269
1270 RETURN
1271 END SUBROUTINE advect_w_implicit
1272
1273 END MODULE module_ieva_em
1274
1275 !-----------------------------------