dyn_em/module_small_step_em.F

   1
   2 !WRF:MODEL_LAYER:DYNAMICS
   3 !
   4
   5 !  SMALL_STEP code for the geometric height coordinate model
   6 !
   7 !---------------------------------------------------------------------------
   8
   9 MODULE module_small_step_em
  10
  11    USE module_configure
  12    USE module_model_constants
  13
  14 CONTAINS
  15
  16 SUBROUTINE small_step_prep( u_1, u_2, v_1, v_2, w_1, w_2, &
  17                             t_1, t_2, ph_1, ph_2,         &
  18                             mub, mu_1, mu_2,              &
  19                             muu, muus, muv, muvs,         &
  20                             mut, muts, mudf,              &
  21                             c1h, c2h, c1f, c2f,           &
  22                             c3h, c4h, c3f, c4f,           &
  23                             u_save, v_save, w_save,       &
  24                             t_save, ph_save, mu_save,     &
  25                             ww, ww_save,                  &
  26                             c2a, pb, p, alt,              &
  27                             msfux, msfuy, msfvx,          &
  28                             msfvx_inv,                    &
  29                             msfvy, msftx, msfty,          &
  30                             rdx, rdy,                     &
  31                             rk_step,                      &
  32                             ids,ide, jds,jde, kds,kde,    &
  33                             ims,ime, jms,jme, kms,kme,    &
  34                             its,ite, jts,jte, kts,kte    )
  35
  36   IMPLICIT NONE  ! religion first
  37
  38 ! declarations for the stuff coming in
  39
  40   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
  41   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
  42   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
  43
  44   INTEGER,      INTENT(IN   )    :: rk_step
  45
  46   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: u_1,   &
  47                                                               v_1,   &
  48                                                               w_1,   &
  49                                                               t_1,   &
  50                                                               ph_1
  51
  52   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(  OUT) :: u_save,   &
  53                                                               v_save,   &
  54                                                               w_save,   &
  55                                                               t_save,   &
  56                                                               ph_save
  57
  58   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: u_2,   &
  59                                                               v_2,   &
  60                                                               w_2,   &
  61                                                               t_2,   &
  62                                                               ph_2
  63
  64   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(  OUT) :: c2a, &
  65                                                                ww_save
  66
  67   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) ::  pb,  &
  68                                                                 p,   &
  69                                                                 alt, &
  70                                                                 ww
  71
  72   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(INOUT) :: mu_1,mu_2
  73
  74   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(IN   ) :: mub,  &
  75                                                                muu,  &
  76                                                                muv,  &
  77                                                                mut,  &
  78                                                                msfux,&
  79                                                                msfuy,&
  80                                                                msfvx,&
  81                                                                msfvx_inv,&
  82                                                                msfvy,&
  83                                                                msftx,&
  84                                                                msfty
  85
  86   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(  OUT) :: muus, &
  87                                                                muvs, &
  88                                                                muts, &
  89                                                                mudf
  90
  91   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(  OUT) :: mu_save
  92
  93   REAL, INTENT(IN) :: rdx,rdy
  94
  95   REAL, DIMENSION( kms:kme ),INTENT(IN   ) :: c1h, c2h, c1f, c2f, &
  96                                               c3h, c4h, c3f, c4f
  97
  98 ! local variables
  99
 100   INTEGER :: i, j, k
 101   INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
 102   INTEGER :: i_endu, j_endv
 103
 104
 105 !<DESCRIPTION>
 106 !
 107 !  small_step_prep prepares the prognostic variables for the small timestep.
 108 !  This includes switching time-levels in the arrays and computing coupled
 109 !  perturbation variables for the small timestep
 110 !  (i.e. mu*u" = mu(t)*u(t)-mu(*)*u(*); mu*u" is advanced during the small
 111 !  timesteps
 112 !
 113 !</DESCRIPTION>
 114
 115     i_start = its
 116     i_end   = min(ite,ide-1)
 117     j_start = jts
 118     j_end   = min(jte,jde-1)
 119     k_start = kts
 120     k_end = min(kte,kde-1)
 121
 122     i_endu = ite
 123     j_endv = jte
 124
 125     !  if this is the first RK step, reset *_1 to *_2
 126     !  (we are replacing the t-dt fields with the time t fields)
 127
 128     IF ((rk_step == 1) ) THEN
 129
 130       DO j=j_start, j_end
 131       DO i=i_start, i_end
 132         MU_1(i,j)=MU_2(i,j)
 133         ww_save(i,kde,j) = 0.
 134         ww_save(i,1,j) = 0.
 135         MUDF(i,j) = 0.  !  initialize external mode div damp to zero
 136       ENDDO
 137       ENDDO
 138
 139       DO j=j_start, j_end
 140       DO k=k_start, k_end
 141       DO i=i_start, i_endu
 142         u_1(i,k,j) = u_2(i,k,j)
 143       ENDDO
 144       ENDDO
 145       ENDDO
 146
 147       DO j=j_start, j_endv
 148       DO k=k_start, k_end
 149       DO i=i_start, i_end
 150         v_1(i,k,j) = v_2(i,k,j)
 151       ENDDO
 152       ENDDO
 153       ENDDO
 154
 155       DO j=j_start, j_end
 156       DO k=k_start, k_end
 157       DO i=i_start, i_end
 158         t_1(i,k,j) = t_2(i,k,j)
 159       ENDDO
 160       ENDDO
 161       ENDDO
 162
 163       DO j=j_start, j_end
 164       DO k=k_start, min(kde,kte)
 165       DO i=i_start, i_end
 166         w_1(i,k,j)  = w_2(i,k,j)
 167         ph_1(i,k,j) = ph_2(i,k,j)
 168       ENDDO
 169       ENDDO
 170       ENDDO
 171
 172     DO j=j_start, j_end
 173       DO i=i_start, i_end
 174         MUTS(i,j)=MUB(i,j)+MU_2(i,j)
 175       ENDDO
 176       DO i=i_start, i_endu
 177         MUUS(i,j) = MUU(i,j)
 178       ENDDO
 179     ENDDO
 180
 181     DO j=j_start, j_endv
 182     DO i=i_start, i_end
 183         MUVS(i,j) = MUV(i,j)
 184     ENDDO
 185     ENDDO
 186
 187     DO j=j_start, j_end
 188       DO i=i_start, i_end
 189         MU_SAVE(i,j)=MU_2(i,j)
 190         MU_2(i,j)=0.
 191       ENDDO
 192     ENDDO
 193
 194     ELSE
 195
 196     DO j=j_start, j_end
 197       DO i=i_start, i_end
 198         MUTS(i,j)=MUB(i,j)+MU_1(i,j)
 199       ENDDO
 200       DO i=i_start, i_endu
 201         MUUS(i,j)=0.5*(MUB(i,j)+MU_1(i,j)+MUB(i-1,j)+MU_1(i-1,j))
 202       ENDDO
 203     ENDDO
 204
 205     DO j=j_start, j_endv
 206     DO i=i_start, i_end
 207       MUVS(i,j)=0.5*(MUB(i,j)+MU_1(i,j)+MUB(i,j-1)+MU_1(i,j-1))
 208     ENDDO
 209     ENDDO
 210
 211     DO j=j_start, j_end
 212       DO i=i_start, i_end
 213         MU_SAVE(i,j)=MU_2(i,j)
 214         MU_2(i,j)=MU_1(i,j)-MU_2(i,j)
 215       ENDDO
 216     ENDDO
 217
 218
 219     END IF
 220
 221     ! set up the small timestep variables
 222
 223       DO j=j_start, j_end
 224       DO i=i_start, i_end
 225         ww_save(i,kde,j) = 0.
 226         ww_save(i,1,j) = 0.
 227       ENDDO
 228       ENDDO
 229
 230     DO j=j_start, j_end
 231     DO k=k_start, k_end
 232     DO i=i_start, i_end
 233       c2a(i,k,j) = cpovcv*(pb(i,k,j)+p(i,k,j))/alt(i,k,j)
 234     ENDDO
 235     ENDDO
 236     ENDDO
 237
 238     DO j=j_start, j_end
 239     DO k=k_start, k_end
 240     DO i=i_start, i_endu
 241       u_save(i,k,j) = u_2(i,k,j)
 242       ! u coupled with my
 243       u_2(i,k,j) = ((c1h(k)*muus(i,j)+c2h(k))*u_1(i,k,j)-(c1h(k)*muu(i,j)+c2h(k))*u_2(i,k,j))/msfuy(i,j)
 244     ENDDO
 245     ENDDO
 246     ENDDO
 247
 248     DO j=j_start, j_endv
 249     DO k=k_start, k_end
 250     DO i=i_start, i_end
 251       v_save(i,k,j) = v_2(i,k,j)
 252       ! v coupled with mx
 253 !      v_2(i,k,j) = ((c1h(k)*muvs(i,j)+c2h(k))*v_1(i,k,j)-(c1h(k)*muv(i,j)+c2h(k))*v_2(i,k,j))/msfvx(i,j)
 254       v_2(i,k,j) = ((c1h(k)*muvs(i,j)+c2h(k))*v_1(i,k,j)-(c1h(k)*muv(i,j)+c2h(k))*v_2(i,k,j))*msfvx_inv(i,j)
 255     ENDDO
 256     ENDDO
 257     ENDDO
 258
 259     DO j=j_start, j_end
 260     DO k=k_start, k_end
 261     DO i=i_start, i_end
 262       t_save(i,k,j) = t_2(i,k,j)
 263       t_2(i,k,j) = (c1h(k)*muts(i,j)+c2h(k))*t_1(i,k,j)-(c1h(k)*mut(i,j)+c2h(k))*t_2(i,k,j)
 264     ENDDO
 265     ENDDO
 266     ENDDO
 267
 268     DO j=j_start, j_end
 269 !    DO k=k_start, min(kde,kte)
 270     DO k=k_start, kde
 271     DO i=i_start, i_end
 272       w_save(i,k,j) = w_2(i,k,j)
 273       ! w coupled with my
 274       w_2(i,k,j)  = ((c1f(k)*muts(i,j)+c2f(k))* w_1(i,k,j)-(c1f(k)*mut(i,j)+c2f(k))* w_2(i,k,j))/msfty(i,j)
 275       ph_save(i,k,j) = ph_2(i,k,j)
 276       ph_2(i,k,j) = ph_1(i,k,j)-ph_2(i,k,j)
 277     ENDDO
 278     ENDDO
 279     ENDDO
 280
 281       DO j=j_start, j_end
 282 !      DO k=k_start, min(kde,kte)
 283       DO k=k_start, kde
 284       DO i=i_start, i_end
 285         ww_save(i,k,j) = ww(i,k,j)
 286       ENDDO
 287       ENDDO
 288       ENDDO
 289
 290 END SUBROUTINE small_step_prep
 291
 292 !-------------------------------------------------------------------------
 293
 294
 295 SUBROUTINE small_step_finish( u_2, u_1, v_2, v_1, w_2, w_1,    &
 296                               t_2, t_1, ph_2, ph_1, ww, ww1,   &
 297                               mu_2, mu_1,                      &
 298                               mut, muts, muu, muus, muv, muvs, &
 299                               c1h, c2h, c1f, c2f,              &
 300                               c3h, c4h, c3f, c4f,              &
 301                               u_save, v_save, w_save,          &
 302                               t_save, ph_save, mu_save,        &
 303                               msfux, msfuy, msfvx, msfvy,      &
 304                               msftx, msfty,                    &
 305                               h_diabatic,                      &
 306                               number_of_small_timesteps,dts,   &
 307                               rk_step, rk_order,               &
 308                               ids,ide, jds,jde, kds,kde,       &
 309                               ims,ime, jms,jme, kms,kme,       &
 310                               its,ite, jts,jte, kts,kte       )
 311
 312
 313   IMPLICIT NONE  ! religion first
 314
 315 !  stuff passed in
 316
 317   INTEGER,                  INTENT(IN   ) :: ids,ide, jds,jde, kds,kde
 318   INTEGER,                  INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
 319   INTEGER,                  INTENT(IN   ) :: its,ite, jts,jte, kts,kte
 320   INTEGER,                  INTENT(IN   ) :: number_of_small_timesteps
 321   INTEGER,                  INTENT(IN   ) :: rk_step, rk_order
 322   REAL,                     INTENT(IN   ) :: dts
 323
 324   REAL,   DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) :: u_1, &
 325                                                                  v_1, &
 326                                                                  w_1, &
 327                                                                  t_1, &
 328                                                                  ww1, &
 329                                                                  ph_1
 330
 331   REAL,   DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: u_2, &
 332                                                                  v_2, &
 333                                                                  w_2, &
 334                                                                  t_2, &
 335                                                                  ww,  &
 336                                                                  ph_2
 337
 338   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN   ) :: u_save,   &
 339                                                               v_save,   &
 340                                                               w_save,   &
 341                                                               t_save,   &
 342                                                               ph_save,  &
 343                                                              h_diabatic
 344
 345   REAL,   DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: muus, muvs
 346   REAL,   DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2, mu_1
 347   REAL,   DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mut, muts, &
 348                                                         muu, muv, mu_save
 349   REAL,   DIMENSION(ims:ime, jms:jme), INTENT(IN   ) :: msfux, msfuy, &
 350                                                         msfvx, msfvy, &
 351                                                         msftx, msfty
 352
 353   REAL, DIMENSION( kms:kme ),INTENT(IN   ) :: c1h, c2h, c1f, c2f, &
 354                                               c3h, c4h, c3f, c4f
 355
 356
 357 ! local stuff
 358
 359   INTEGER         :: i,j,k
 360   INTEGER :: i_start, i_end, j_start, j_end, i_endu, j_endv
 361
 362 !<DESCRIPTION>
 363 !
 364 !  small_step_finish reconstructs the full uncoupled prognostic variables
 365 !  from the coupled perturbation variables used in the small timesteps.
 366 !
 367 !</DESCRIPTION>
 368
 369   i_start = its
 370   i_end   = min(ite,ide-1)
 371   j_start = jts
 372   j_end   = min(jte,jde-1)
 373
 374   i_endu = ite
 375   j_endv = jte
 376
 377 !    addition of time level t back into variables
 378
 379   DO j = j_start, j_endv
 380   DO k = kds, kde-1
 381   DO i = i_start, i_end
 382     ! v coupled with mx
 383     v_2(i,k,j) = (msfvx(i,j)*v_2(i,k,j) + v_save(i,k,j)*(c1h(k)*muv(i,j)+c2h(k)))/(c1h(k)*muvs(i,j)+c2h(k))
 384   ENDDO
 385   ENDDO
 386   ENDDO
 387
 388   DO j = j_start, j_end
 389   DO k = kds, kde-1
 390   DO i = i_start, i_endu
 391     ! u coupled with my
 392     u_2(i,k,j) = (msfuy(i,j)*u_2(i,k,j) + u_save(i,k,j)*(c1h(k)*muu(i,j)+c2h(k)))/(c1h(k)*muus(i,j)+c2h(k))
 393   ENDDO
 394   ENDDO
 395   ENDDO
 396
 397   DO j = j_start, j_end
 398   DO k = kds, kde
 399   DO i = i_start, i_end
 400     ! w coupled with my
 401     w_2(i,k,j) = (msfty(i,j)*w_2(i,k,j) + w_save(i,k,j)*(c1f(k)*mut(i,j)+c2f(k)))/(c1f(k)*muts(i,j)+c2f(k))
 402     ph_2(i,k,j) = ph_2(i,k,j) + ph_save(i,k,j)
 403     ww(i,k,j) = ww(i,k,j) + ww1(i,k,j)
 404   ENDDO
 405   ENDDO
 406   ENDDO
 407
 408   IF ( rk_step < rk_order ) THEN
 409     DO j = j_start, j_end
 410     DO k = kds, kde-1
 411     DO i = i_start, i_end
 412       t_2(i,k,j) = (t_2(i,k,j) + t_save(i,k,j)*(c1h(k)*mut(i,j)+c2h(k)))/(c1h(k)*muts(i,j)+c2h(k))
 413     ENDDO
 414     ENDDO
 415     ENDDO
 416   ELSE
 417
 418     DO j = j_start, j_end
 419     DO k = kds, kde-1
 420     DO i = i_start, i_end
 421       t_2(i,k,j) = (t_2(i,k,j) - dts*number_of_small_timesteps*(c1h(k)*mut(i,j)+c2h(k))*h_diabatic(i,k,j) &
 422                     + t_save(i,k,j)*(c1h(k)*mut(i,j)+c2h(k)))/(c1h(k)*muts(i,j)+c2h(k))
 423     ENDDO
 424     ENDDO
 425     ENDDO
 426   ENDIF
 427
 428   DO j = j_start, j_end
 429   DO i = i_start, i_end
 430     mu_2(i,j) = mu_2(i,j) + mu_save(i,j)
 431   ENDDO
 432   ENDDO
 433
 434 END SUBROUTINE small_step_finish
 435
 436 !-----------------------------------------------------------------------
 437
 438 SUBROUTINE calc_p_rho( al, p, ph,                    &
 439                        alt, t_2, t_1, c2a, pm1,      &
 440                        mu, mut,                      &
 441                        c1h, c2h, c1f, c2f,           &
 442                        c3h, c4h, c3f, c4f,           &
 443                        znu, t0,                      &
 444                        rdnw, dnw, smdiv,             &
 445                        non_hydrostatic, step,        &
 446                        ids, ide, jds, jde, kds, kde, &
 447                        ims, ime, jms, jme, kms, kme, &
 448                        its,ite, jts,jte, kts,kte    )
 449
 450   IMPLICIT NONE  ! religion first
 451
 452 ! declarations for the stuff coming in
 453
 454   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
 455   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
 456   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
 457
 458   INTEGER,      INTENT(IN   )    :: step
 459
 460   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(  OUT) :: al,   &
 461                                                                p
 462
 463   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN   ) :: alt,   &
 464                                                               t_2,   &
 465                                                               t_1,   &
 466                                                               c2a
 467
 468   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph, pm1
 469
 470   REAL, DIMENSION(ims:ime, jms:jme)         , INTENT(IN   ) :: mu,   &
 471                                                                mut
 472
 473   REAL, DIMENSION(kms:kme)         , INTENT(IN   ) :: dnw,  &
 474                                                       rdnw, &
 475                                                       znu
 476
 477   REAL,                                       INTENT(IN   ) :: t0, smdiv
 478
 479   REAL, DIMENSION( kms:kme ),INTENT(IN   ) :: c1h, c2h, c1f, c2f, &
 480                                               c3h, c4h, c3f, c4f
 481
 482   LOGICAL, INTENT(IN   )  :: non_hydrostatic
 483
 484 ! local variables
 485
 486   INTEGER :: i, j, k
 487   INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
 488   REAL    :: ptmp
 489
 490 !<DESCRIPTION>
 491 !
 492 !  For the nonhydrostatic option,
 493 !  calc_p_rho computes the perturbation inverse density and
 494 !  perturbation pressure from the hydrostatic relation and the
 495 !  linearized equation of state, respectively.
 496 !
 497 !  For the hydrostatic option,
 498 !  calc_p_rho computes the perturbation pressure, perturbation density,
 499 !  and perturbation geopotential
 500 !  from the vertical coordinate definition, linearized equation of state
 501 !  and the hydrostatic relation, respectively.
 502 !
 503 !  forward weighting of the pressure (divergence damping) is also
 504 !  computed here.
 505 !
 506 !</DESCRIPTION>
 507
 508    i_start = its
 509    i_end   = min(ite,ide-1)
 510    j_start = jts
 511    j_end   = min(jte,jde-1)
 512    k_start = kts
 513    k_end = min(kte,kde-1)
 514
 515    IF (non_hydrostatic) THEN
 516      DO j=j_start, j_end
 517      DO k=k_start, k_end
 518      DO i=i_start, i_end
 519
 520 !  al computation is all dry, so ok with moisture
 521
 522       al(i,k,j)=-1./(c1h(k)*Mut(i,j)+c2h(k))*(alt(i,k,j)*(c1h(k)*mu(i,j))  &
 523              +rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
 524
 525 !  this is temporally linearized p, no moisture correction needed
 526
 527       p(i,k,j)=c2a(i,k,j)*(alt(i,k,j)*(t_2(i,k,j)-(c1h(k)*mu(i,j))*t_1(i,k,j))  &
 528                        /((c1h(k)*Mut(i,j)+c2h(k))*(t0+t_1(i,k,j)))-al (i,k,j))
 529
 530      ENDDO
 531      ENDDO
 532      ENDDO
 533
 534    ELSE  ! hydrostatic calculation
 535
 536        DO j=j_start, j_end
 537        DO k=k_start, k_end
 538        DO i=i_start, i_end
 539          p(i,k,j)=MU(i,j)*c3h(k)
 540          al(i,k,j)=alt(i,k,j)*(t_2(i,k,j)-(c1h(k)*mu(i,j))*t_1(i,k,j))           &
 541                       /((c1h(k)*Mut(i,j)+c2h(k))*(t0+t_1(i,k,j)))-p(i,k,j)/c2a(i,k,j)
 542          ph(i,k+1,j)=ph(i,k,j)-dnw(k)*((c1h(k)*Mut(i,j)+c2h(k))*al (i,k,j)              &
 543                           +(c1h(k)*mu(i,j))*alt(i,k,j))
 544        ENDDO
 545        ENDDO
 546        ENDDO
 547    END IF
 548 !  divergence damping setup
 549      IF (step == 0) then   ! we're initializing small timesteps
 550        DO j=j_start, j_end
 551        DO k=k_start, k_end
 552        DO i=i_start, i_end
 553          pm1(i,k,j)=p(i,k,j)
 554        ENDDO
 555        ENDDO
 556        ENDDO
 557      ELSE                     ! we're in the small timesteps
 558        DO j=j_start, j_end    ! and adding div damping component
 559        DO k=k_start, k_end
 560        DO i=i_start, i_end
 561          ptmp = p(i,k,j)
 562          p(i,k,j) = p(i,k,j) + smdiv*(p(i,k,j)-pm1(i,k,j))
 563          pm1(i,k,j) = ptmp
 564        ENDDO
 565        ENDDO
 566        ENDDO
 567      END IF
 568 END SUBROUTINE calc_p_rho
 569 !----------------------------------------------------------------------
 570 SUBROUTINE calc_coef_w( a,alpha,gamma,              &
 571                         mut,                        &
 572                         c1h, c2h, c1f, c2f,         &
 573                         c3h, c4h, c3f, c4f,         &
 574                         cqw,                        &
 575                         rdn, rdnw, c2a,             &
 576                         dts, g, epssm, top_lid,     &
 577                         ids,ide, jds,jde, kds,kde,  & ! domain dims
 578                         ims,ime, jms,jme, kms,kme,  & ! memory dims
 579                         its,ite, jts,jte, kts,kte  )  ! tile   dims
 580
 581       IMPLICIT NONE  ! religion first
 582 !  passed in through the call
 583   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
 584   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
 585   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
 586   LOGICAL,      INTENT(IN   )    :: top_lid
 587   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN   ) :: c2a,  &
 588                                                                cqw
 589   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: alpha, &
 590                                                                gamma, &
 591                                                                a
 592   REAL, DIMENSION(ims:ime, jms:jme),          INTENT(IN   ) :: mut
 593   REAL, DIMENSION(kms:kme),                   INTENT(IN   ) :: rdn,   &
 594                                                                rdnw
 595   REAL,                                       INTENT(IN   ) :: epssm, &
 596                                                                dts,   &
 597                                                                g
 598   REAL, DIMENSION( kms:kme ),INTENT(IN   ) :: c1h, c2h, c1f, c2f, &
 599                                               c3h, c4h, c3f, c4f
 600 !  Local stack data.
 601   REAL, DIMENSION(ims:ime)                         :: cof
 602   REAL  :: b, c
 603   REAL  :: muthmutf_kk, muthmutf_km1k, muthmutf_kkp1
 604   INTEGER :: i, j, k, kk, i_start, i_end, j_start, j_end, k_start, k_end
 605   INTEGER :: ij, ijp, ijm, lid_flag
 606 !<DESCRIPTION>
 607 !
 608 !  calc_coef_w calculates the coefficients needed for the
 609 !  implicit solution of the vertical momentum and geopotential equations.
 610 !  This requires solution of a tri-diagonal equation.
 611 !
 612 !</DESCRIPTION>
 613       i_start = its
 614       i_end   = min(ite,ide-1)
 615       j_start = jts
 616       j_end   = min(jte,jde-1)
 617       k_start = kts
 618       k_end   = kte-1
 619       lid_flag=1
 620       IF(top_lid)lid_flag=0
 621      outer_j_loop:  DO j = j_start, j_end
 622         k = kde-1
 623         DO i = i_start, i_end
 624           cof(i)  = (.5*dts*g*(1.+epssm))**2
 625           a(i, 2 ,j) = 0.
 626           a(i,kde,j) =-2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)*lid_flag/((c1h(k)*MUT(i,j)+c2h(k))*(c1f(k)*MUT(i,j)+c2f(k)))
 627           gamma(i,1 ,j) = 0.
 628         ENDDO
 629         DO kk=3,kde-1
 630         k=kk-1
 631         DO i=i_start, i_end
 632           a(i,kk,j) =   -cqw(i,kk,j)*cof(i)*rdn(kk)* rdnw(kk-1)*c2a(i,kk-1,j)/((c1h(k)*MUT(i,j)+c2h(k))*(c1f(k)*MUT(i,j)+c2f(k)))
 633         ENDDO
 634         ENDDO
 635         DO k=2,kde-1
 636         DO i=i_start, i_end
 637           b = 1.+cqw(i,k,j)*cof(i)*rdn(k)*(rdnw(k  )*c2a(i,k,  j)/((c1h(k)*MUT(i,j)+c2h(k))*(c1f(k)*MUT(i,j)+c2f(k)))  &
 638                                           +rdnw(k-1)*c2a(i,k-1,j)/((c1h(k-1)*MUT(i,j)+c2h(k-1))*(c1f(k)*MUT(i,j)+c2f(k))) )
 639           c =   -cqw(i,k,j)*cof(i)*rdn(k)*rdnw(k  )*c2a(i,k,j  )/((c1h(k)*MUT(i,j)+c2h(k))*(c1f(k+1)*MUT(i,j)+c2f(k+1)))
 640           alpha(i,k,j) = 1./(b-a(i,k,j)*gamma(i,k-1,j))
 641           gamma(i,k,j) = c*alpha(i,k,j)
 642         ENDDO
 643         ENDDO
 644         k=kde
 645         DO i=i_start, i_end
 646           b = 1.+2.*cof(i)*rdnw(kde-1)**2*c2a(i,kde-1,j)/((c1h(k-1)*MUT(i,j)+c2h(k-1))*(c1f(k)*MUT(i,j)+c2f(k)))
 647           c = 0.
 648           alpha(i,kde,j) = 1./(b-a(i,kde,j)*gamma(i,kde-1,j))
 649           gamma(i,kde,j) = c*alpha(i,kde,j)
 650         ENDDO
 651       ENDDO outer_j_loop
 652 END SUBROUTINE calc_coef_w
 653 !----------------------------------------------------------------------
 654 SUBROUTINE advance_uv ( u, ru_tend, v, rv_tend,        &
 655                         p, pb,                         &
 656                         ph, php, alt, al, mu,          &
 657                         muu, cqu, muv, cqv, mudf,      &
 658                         c1h, c2h, c1f, c2f,            &
 659                         c3h, c4h, c3f, c4f,            &
 660                         msfux, msfuy, msfvx,           &
 661                         msfvx_inv, msfvy,              &
 662                         rdx, rdy, dts,                 &
 663                         cf1, cf2, cf3, fnm, fnp,       &
 664                         emdiv,                         &
 665                         rdnw, config_flags, spec_zone, &
 666                         non_hydrostatic, top_lid,      &
 667                         ids, ide, jds, jde, kds, kde,  &
 668                         ims, ime, jms, jme, kms, kme,  &
 669                         its, ite, jts, jte, kts, kte  )
 670       IMPLICIT NONE  ! religion first
 671 ! stuff coming in
 672       TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
 673       LOGICAL, INTENT(IN   ) :: non_hydrostatic, top_lid
 674       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
 675       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
 676       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
 677       INTEGER,      INTENT(IN   )    :: spec_zone
 678       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),  &
 679             INTENT(INOUT) ::                          &
 680                                                   u,  &
 681                                                   v
 682       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
 683             INTENT(IN   ) ::                          &
 684                                              ru_tend, &
 685                                              rv_tend, &
 686                                              ph,      &
 687                                              php,     &
 688                                              p,       &
 689                                              pb,      &
 690                                              alt,     &
 691                                              al,      &
 692                                              cqu,     &
 693                                              cqv
 694       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: muu,  &
 695                                                                 muv,  &
 696                                                                 mu,   &
 697                                                                 mudf
 698       REAL, DIMENSION( kms:kme ),              INTENT(IN   ) :: fnm,    &
 699                                                                 fnp ,   &
 700                                                                 rdnw
 701       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msfux,  &
 702                                                                 msfuy,  &
 703                                                                 msfvx,  &
 704                                                                 msfvy,  &
 705                                                                 msfvx_inv
 706       REAL,                                    INTENT(IN   ) :: rdx,    &
 707                                                                 rdy,    &
 708                                                                 dts,    &
 709                                                                 cf1,    &
 710                                                                 cf2,    &
 711                                                                 cf3,    &
 712                                                                 emdiv
 713       REAL, DIMENSION( kms:kme ),INTENT(IN   ) :: c1h, c2h, c1f, c2f, &
 714                                                   c3h, c4h, c3f, c4f
 715
 716 !  Local 3d array from the stack (note tile size)
 717       REAL, DIMENSION (its:ite, kts:kte) :: dpn, dpxy
 718       REAL, DIMENSION (its:ite) :: mudf_xy
 719       REAL                      :: dx, dy
 720       INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
 721       INTEGER :: i_endu, j_endv, k_endw
 722       INTEGER :: i_start_up, i_end_up, j_start_up, j_end_up
 723       INTEGER :: i_start_vp, i_end_vp, j_start_vp, j_end_vp
 724       INTEGER :: i_start_u_tend, i_end_u_tend, j_start_v_tend, j_end_v_tend
 725 !<DESCRIPTION>
 726 !
 727 !  advance_uv advances the explicit perturbation horizontal momentum
 728 !  equations (u,v) by adding in the large-timestep tendency along with
 729 !  the small timestep pressure gradient tendency.
 730 !
 731 !</DESCRIPTION>
 732 !  now, the real work.
 733 !  set the loop bounds taking into account boundary conditions.
 734     IF( config_flags%nested .or. config_flags%specified ) THEN
 735       i_start = max( its,ids+spec_zone )
 736       i_end   = min( ite,ide-spec_zone-1 )
 737       j_start = max( jts,jds+spec_zone )
 738       j_end   = min( jte,jde-spec_zone-1 )
 739       k_start = kts
 740       k_end   = min( kte, kde-1 )
 741       i_endu = min( ite,ide-spec_zone )
 742       j_endv = min( jte,jde-spec_zone )
 743       k_endw = k_end
 744       IF( config_flags%periodic_x) THEN
 745         i_start = its
 746         i_end   = min(ite,ide-1)
 747         i_endu = ite
 748       ENDIF
 749     ELSE
 750       i_start = its
 751       i_end   = min(ite,ide-1)
 752       j_start = jts
 753       j_end   = min(jte,jde-1)
 754       k_start = kts
 755       k_end   = kte-1
 756       i_endu = ite
 757       j_endv = jte
 758       k_endw = k_end
 759     ENDIF
 760       i_start_up = i_start
 761       i_end_up   = i_endu
 762       j_start_up = j_start
 763       j_end_up   = j_end
 764       i_start_vp = i_start
 765       i_end_vp   = i_end
 766       j_start_vp = j_start
 767       j_end_vp   = j_endv
 768       IF ( (config_flags%open_xs   .or.     &
 769             config_flags%symmetric_xs   )   &
 770             .and. (its == ids) )            &
 771                  i_start_up = i_start_up + 1
 772       IF ( (config_flags%open_xe    .or.  &
 773             config_flags%symmetric_xe   ) &
 774              .and. (ite == ide) )         &
 775                  i_end_up   = i_end_up - 1
 776       IF ( (config_flags%open_ys    .or.   &
 777             config_flags%symmetric_ys  .or.   &
 778             config_flags%polar   )  &
 779                      .and. (jts == jds) )  &
 780                  j_start_vp = j_start_vp + 1
 781       IF ( (config_flags%open_ye     .or. &
 782             config_flags%symmetric_ye  .or.   &
 783             config_flags%polar   )  &
 784             .and. (jte == jde) )          &
 785                  j_end_vp   = j_end_vp - 1
 786       i_start_u_tend = i_start
 787       i_end_u_tend   = i_endu
 788       j_start_v_tend = j_start
 789       j_end_v_tend   = j_endv
 790       IF ( config_flags%symmetric_xs .and. (its == ids) ) &
 791                      i_start_u_tend = i_start_u_tend+1
 792       IF ( config_flags%symmetric_xe .and. (ite == ide) ) &
 793                      i_end_u_tend = i_end_u_tend-1
 794       IF ( config_flags%symmetric_ys .and. (jts == jds) ) &
 795                      j_start_v_tend = j_start_v_tend+1
 796       IF ( config_flags%symmetric_ye .and. (jte == jde) ) &
 797                      j_end_v_tend = j_end_v_tend-1
 798    dx = 1./rdx
 799    dy = 1./rdy
 800 !  start real calculations.
 801 !  first, u
 802   u_outer_j_loop: DO j = j_start, j_end
 803    DO k = k_start, k_end
 804    DO i = i_start_u_tend, i_end_u_tend
 805      u(i,k,j) = u(i,k,j) + dts*ru_tend(i,k,j)
 806    ENDDO
 807    ENDDO
 808    DO i = i_start_up, i_end_up
 809      MUDF_XY(i)= -emdiv*dx*(MUDF(i,j)-MUDF(i-1,j))/msfuy(i,j)
 810    ENDDO
 811    DO k = k_start, k_end
 812    DO i = i_start_up, i_end_up
 813 !     Comments on map scale factors:
 814 !     x pressure gradient: ADT eqn 44, penultimate term on RHS
 815 !        = -(mx/my)*(mu/rho)*partial dp/dx
 816 !     [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
 817 !     Klemp et al. splits into 2 terms:
 818 !        mu alpha partial dp/dx + partial dp/dnu * partial dphi/dx
 819 !     then into 4 terms:
 820 !        mu alpha partial dp'/dx + nu mu alpha' partial dmubar/dx +
 821 !      + mu partial dphi/dx + partial dphi'/dx * (partial dp'/dnu - mu')
 822 !
 823 !     first 3 terms:
 824 !     ph, alt, p, al, pb not coupled
 825 !     since we want tendency to fit ADT eqn 44 (coupled) we need to
 826 !     multiply by (mx/my):
 827 !
 828      dpxy(i,k)= (msfux(i,j)/msfuy(i,j))*.5*rdx*(c1h(k)*muu(i,j)+c2h(k))*(         &
 829        ((ph (i,k+1,j)-ph (i-1,k+1,j))+(ph (i,k,j)-ph (i-1,k,j)))  &
 830       +(alt(i,k  ,j)+alt(i-1,k  ,j))*(p  (i,k,j)-p  (i-1,k,j))  &
 831       +(al (i,k  ,j)+al (i-1,k  ,j))*(pb (i,k,j)-pb (i-1,k,j)) )
 832    ENDDO
 833    ENDDO
 834    IF (non_hydrostatic) THEN
 835      DO i = i_start_up, i_end_up
 836        dpn(i,1) = .5*( cf1*(p(i,1,j)+p(i-1,1,j))  &
 837                       +cf2*(p(i,2,j)+p(i-1,2,j))  &
 838                       +cf3*(p(i,3,j)+p(i-1,3,j)) )
 839        dpn(i,kde) = 0.
 840      ENDDO
 841      IF (top_lid) THEN
 842        DO i = i_start_up, i_end_up
 843          dpn(i,kde) =.5*( cf1*(p(i-1,kde-1,j)+p(i,kde-1,j))   &
 844                          +cf2*(p(i-1,kde-2,j)+p(i,kde-2,j))   &
 845                          +cf3*(p(i-1,kde-3,j)+p(i,kde-3,j))  )
 846        ENDDO
 847      ENDIF
 848      DO k = k_start+1, k_end
 849        DO i = i_start_up, i_end_up
 850          dpn(i,k) = .5*( fnm(k)*(p(i,k  ,j)+p(i-1,k  ,j))   &
 851                         +fnp(k)*(p(i,k-1,j)+p(i-1,k-1,j)) )
 852        ENDDO
 853      ENDDO
 854 !    Comments on map scale factors:
 855 !    4th term:
 856 !    php, dpn, mu not coupled
 857 !    since we want tendency to fit ADT eqn 44 (coupled) we need to
 858 !    multiply by (mx/my):
 859      DO k = k_start, k_end
 860        DO i = i_start_up, i_end_up
 861          dpxy(i,k)=dpxy(i,k) + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))*        &
 862            (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*((c1h(k)*mu(i-1,j))+(c1h(k)*mu(i,j))))
 863        ENDDO
 864      ENDDO
 865    END IF
 866    DO k = k_start, k_end
 867      DO i = i_start_up, i_end_up
 868        u(i,k,j)=u(i,k,j)-dts*cqu(i,k,j)*dpxy(i,k)+(c1h(k)*mudf_xy(i))
 869      ENDDO
 870    ENDDO
 871    ENDDO u_outer_j_loop
 872 ! now v
 873   v_outer_j_loop: DO j = j_start_v_tend, j_end_v_tend
 874    DO k = k_start, k_end
 875    DO i = i_start, i_end
 876      v(i,k,j) = v(i,k,j) + dts*rv_tend(i,k,j)
 877    ENDDO
 878    ENDDO
 879    DO i = i_start, i_end
 880      MUDF_XY(i)= -emdiv*dy*(MUDF(i,j)-MUDF(i,j-1))*msfvx_inv(i,j)
 881    ENDDO
 882      IF (     ( j >= j_start_vp)  &
 883        .and.( j <= j_end_vp  ) )  THEN
 884          DO k = k_start, k_end
 885   DO i = i_start, i_end
 886 !      Comments on map scale factors:
 887 !      y pressure gradient: ADT eqn 45, penultimate term on RHS
 888 !         = -(my/mx)*(mu/rho)*partial dp/dy
 889 !      [i.e., first rho->mu; 2nd still rho; alpha=1/rho]
 890 !      Klemp et al. splits into 2 terms:
 891 !         mu alpha partial dp/dy + partial dp/dnu * partial dphi/dy
 892 !      then into 4 terms:
 893 !         mu alpha partial dp'/dy + nu mu alpha' partial dmubar/dy +
 894 !       + mu partial dphi/dy + partial dphi'/dy * (partial dp'/dnu - mu')
 895 !
 896 !      first 3 terms:
 897 !      ph, alt, p, al, pb not coupled
 898 !      since we want tendency to fit ADT eqn 45 (coupled) we need to
 899 !      multiply by (my/mx):
 900 !      mudf_xy is NOT a map scale factor coupling
 901 !      it is some sort of divergence damping
 902        dpxy(i,k)= (msfvy(i,j)/msfvx(i,j))*.5*rdy*(c1h(k)*muv(i,j)+c2h(k))*(       &
 903         ((ph(i,k+1,j)-ph(i,k+1,j-1))+(ph (i,k,j)-ph (i,k,j-1)))   &
 904         +(alt(i,k  ,j)+alt(i,k  ,j-1))*(p  (i,k,j)-p  (i,k,j-1))  &
 905         +(al (i,k  ,j)+al (i,k  ,j-1))*(pb (i,k,j)-pb (i,k,j-1)) )
 906      ENDDO
 907      ENDDO
 908      IF (non_hydrostatic) THEN
 909        DO i = i_start, i_end
 910          dpn(i,1) = .5*( cf1*(p(i,1,j)+p(i,1,j-1))  &
 911                         +cf2*(p(i,2,j)+p(i,2,j-1))  &
 912                         +cf3*(p(i,3,j)+p(i,3,j-1)) )
 913          dpn(i,kde) = 0.
 914        ENDDO
 915      IF (top_lid) THEN
 916        DO i = i_start, i_end
 917          dpn(i,kde) =.5*( cf1*(p(i,kde-1,j-1)+p(i,kde-1,j))   &
 918                          +cf2*(p(i,kde-2,j-1)+p(i,kde-2,j))   &
 919                          +cf3*(p(i,kde-3,j-1)+p(i,kde-3,j))  )
 920        ENDDO
 921      ENDIF
 922        DO k = k_start+1, k_end
 923          DO i = i_start, i_end
 924            dpn(i,k) = .5*( fnm(k)*(p(i,k  ,j)+p(i,k  ,j-1))  &
 925                           +fnp(k)*(p(i,k-1,j)+p(i,k-1,j-1)) )
 926          ENDDO
 927        ENDDO
 928 !      Comments on map scale factors:
 929 !      4th term:
 930 !      php, dpn, mu not coupled
 931 !      since we want tendency to fit ADT eqn 45 (coupled) we need to
 932 !      multiply by (my/mx):
 933        DO k = k_start, k_end
 934          DO i = i_start, i_end
 935            dpxy(i,k)=dpxy(i,k) + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))*    &
 936              (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*((c1h(k)*mu(i,j-1))+(c1h(k)*mu(i,j))))
 937          ENDDO
 938        ENDDO
 939      END IF
 940      DO k = k_start, k_end
 941        DO i = i_start, i_end
 942          v(i,k,j)=v(i,k,j)-dts*cqv(i,k,j)*dpxy(i,k)+(c1h(k)*mudf_xy(i))
 943        ENDDO
 944      ENDDO
 945    END IF
 946   ENDDO  v_outer_j_loop
 947   ! The check for j_start_vp and j_end_vp makes sure that the edges in v
 948   ! are not updated.  Let's add a polar boundary condition check here for
 949   ! safety to ensure that v at the poles never gets to be non-zero in the
 950   ! small time steps.
 951   IF (config_flags%polar) THEN
 952      IF (jts == jds) THEN
 953         DO k = k_start, k_end
 954         DO i = i_start, i_end
 955            v(i,k,jds) = 0.
 956         ENDDO
 957         ENDDO
 958      END IF
 959      IF (jte == jde) THEN
 960         DO k = k_start, k_end
 961         DO i = i_start, i_end
 962            v(i,k,jde) = 0.
 963         ENDDO
 964         ENDDO
 965      END IF
 966   END IF
 967 END SUBROUTINE advance_uv
 968 !---------------------------------------------------------------------
 969 SUBROUTINE advance_mu_t( ww, ww_1, u, u_1, v, v_1,            &
 970                          mu, mut, muave, muts, muu, muv, mudf,&
 971                          c1h, c2h, c1f, c2f,                  &
 972                          c3h, c4h, c3f, c4f,                  &
 973                          uam, vam, wwam, t, t_1,              &
 974                          t_ave, ft, mu_tend,                  &
 975                          rdx, rdy, dts, epssm,                &
 976                          dnw, fnm, fnp, rdnw,                 &
 977                          msfux, msfuy, msfvx, msfvx_inv,      &
 978                          msfvy, msftx, msfty,                 &
 979                          step, config_flags,                  &
 980                          ids, ide, jds, jde, kds, kde,        &
 981                          ims, ime, jms, jme, kms, kme,        &
 982                          its, ite, jts, jte, kts, kte        )
 983   IMPLICIT NONE  ! religion first
 984 ! stuff coming in
 985   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
 986   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
 987   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
 988   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
 989   INTEGER,      INTENT(IN   )    :: step
 990   REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),   &
 991             INTENT(IN   ) ::                       &
 992                                               u,   &
 993                                               v,   &
 994                                               u_1, &
 995                                               v_1, &
 996                                               t_1, &
 997                                               ft
 998   REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),      &
 999             INTENT(INOUT) ::                          &
1000                                               ww,     &
1001                                               ww_1,   &
1002                                               t,      &
1003                                               t_ave,  &
1004                                               uam,    &
1005                                               vam,    &
1006                                               wwam
1007
1008   REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: muu,  &
1009                                                             muv,  &
1010                                                             mut,  &
1011                                                             msfux,&
1012                                                             msfuy,&
1013                                                             msfvx,&
1014                                                             msfvx_inv,&
1015                                                             msfvy,&
1016                                                             msftx,&
1017                                                             msfty,&
1018                                                             mu_tend
1019   REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(  OUT) :: muave, &
1020                                                             muts,  &
1021                                                             mudf
1022   REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: mu
1023   REAL, DIMENSION( kms:kme ),              INTENT(IN   ) :: fnm,    &
1024                                                             fnp,    &
1025                                                             dnw,    &
1026                                                             rdnw
1027   REAL,                                    INTENT(IN   ) :: rdx,    &
1028                                                             rdy,    &
1029                                                             dts,    &
1030                                                             epssm
1031   REAL, DIMENSION( kms:kme ),INTENT(IN   ) :: c1h, c2h, c1f, c2f, &
1032                                               c3h, c4h, c3f, c4f
1033 !  Local arrays from the stack (note tile size)
1034   REAL, DIMENSION (its:ite, kts:kte) :: wdtn, dvdxi
1035   REAL, DIMENSION (its:ite) :: dmdt
1036   INTEGER :: i, j, k, kk, i_start, i_end, j_start, j_end, k_start, k_end
1037   INTEGER :: i_endu, j_endv
1038   REAL    :: acc
1039 !<DESCRIPTION>
1040 !
1041 !  advance_mu_t advances the explicit perturbation theta equation and the mass
1042 !  conservation equation.  In addition, the small timestep omega is updated,
1043 !  and some quantities needed in other places are squirrelled away.
1044 !
1045 !</DESCRIPTION>
1046 !  now, the real work.
1047 !  set the loop bounds taking into account boundary conditions.
1048   i_start = its
1049   i_end   = min(ite,ide-1)
1050   j_start = jts
1051   j_end   = min(jte,jde-1)
1052   k_start = kts
1053   k_end   = kte-1
1054   IF ( .NOT. config_flags%periodic_x )THEN
1055      IF ( config_flags%specified .or. config_flags%nested ) then
1056        i_start = max(its,ids+1)
1057        i_end   = min(ite,ide-2)
1058      ENDIF
1059   ENDIF
1060   IF ( config_flags%specified .or. config_flags%nested ) then
1061      j_start = max(jts,jds+1)
1062      j_end   = min(jte,jde-2)
1063   ENDIF
1064   i_endu = ite
1065   j_endv = jte
1066 !        CALCULATION OF WW (dETA/dt)
1067    DO j = j_start, j_end
1068      DO i=i_start, i_end
1069             DMDT(i) = 0.
1070      ENDDO
1071 !  NOTE:  mu is not coupled with the map scale factor.
1072 !         ww (omega) IS coupled with the map scale factor.
1073 !         Being coupled with the map scale factor means
1074 !           multiplication by (1/msft) in this case.
1075 !  Comments on map scale factors
1076 !  ADT eqn 47:
1077 !  partial drho/dt = -mx*my[partial d/dx(rho u/my) + partial d/dy(rho v/mx)]
1078 !                    -partial d/dz(rho w)
1079 !  with rho -> mu, dividing by my, and with partial d/dnu(rho nu/my [=ww])
1080 !  as the final term (because we're looking for d_nu_/dt)
1081 !
1082 !  begin by integrating with respect to nu from bottom to top
1083 !  BCs are ww=0 at both
1084 !  final term gives 0
1085 !  first term gives Integral([1/my]partial d mu/dt) over total column = dm/dt
1086 !  RHS remaining is Integral(-mx[partial d/dx(mu u/my) +
1087 !                                partial d/dy(mu v/mx)]) over column
1088 !  lines below find RHS terms at each level then set dmdt = sum over all levels
1089 !
1090 !  [don't divide the below by msfty until find ww, since dmdt is used in
1091 !   the meantime]
1092      DO k=k_start, k_end
1093      DO i=i_start, i_end
1094          dvdxi(i,k) = msftx(i,j)*msfty(i,j)*(                                      &
1095                      rdy*( (v(i,k,j+1)+(c1h(k)*muv(i,j+1)+c2h(k))*v_1(i,k,j+1)*msfvx_inv(i,j+1))   &
1096                           -(v(i,k,j  )+(c1h(k)*muv(i,j  )+c2h(k))*v_1(i,k,j  )*msfvx_inv(i,j  )) ) &
1097                     +rdx*( (u(i+1,k,j)+(c1h(k)*muu(i+1,j)+c2h(k))*u_1(i+1,k,j)/msfuy(i+1,j))       &
1098                           -(u(i,k,j  )+(c1h(k)*muu(i  ,j)+c2h(k))*u_1(i,k,j  )/msfuy(i  ,j)) ))
1099         DMDT(i)    = DMDT(i) + dnw(k)*dvdxi(i,k)
1100      ENDDO
1101      ENDDO
1102      DO i=i_start, i_end
1103        MUAVE(i,j) = MU(i,j)
1104        MU(i,j) = MU(i,j)+dts*(DMDT(i)+MU_TEND(i,j))
1105        MUDF(i,j) = (DMDT(i)+MU_TEND(i,j)) ! save tendency for div damp filter
1106        MUTS(i,j) = MUT(i,j)+MU(i,j)
1107        MUAVE(i,j) =.5*((1.+epssm)*MU(i,j)+(1.-epssm)*MUAVE(i,j))
1108      ENDDO
1109      DO kk=2,k_end
1110      k=kk-1
1111      DO i=i_start, i_end
1112        ww(i,kk,j)=ww(i,kk-1,j)-dnw(kk-1)*((c1h(k)*dmdt(i))+dvdxi(i,kk-1)+(c1h(k)*mu_tend(i,j)))/msfty(i,j)
1113      ENDDO
1114      END DO
1115 !  NOTE:  ww_1 (large timestep ww) is already coupled with the
1116 !         map scale factor
1117      DO k=1,k_end
1118      DO i=i_start, i_end
1119        ww(i,k,j)=ww(i,k,j)-ww_1(i,k,j)
1120      END DO
1121      END DO
1122    ENDDO
1123 ! CALCULATION OF THETA
1124 ! NOTE: theta'' is not coupled with the map-scale factor,
1125 !       while the theta'' tendency is coupled (i.e., mult by 1/msft)
1126 ! Comments on map scale factors
1127 ! BUT NOTE THAT both are mass coupled
1128 ! in flux form equations (Klemp et al.) Theta = mu*theta
1129 !
1130 ! scalar eqn: partial d/dt(rho q/my) = -mx[partial d/dx(q rho u/my) +
1131 !                                          partial d/dy(q rho v/mx)]
1132 !                                      - partial d/dz(q rho w/my)
1133 ! with rho -> mu, and with partial d/dnu(q rho nu/my) as the final term
1134 !
1135 ! adding previous tendency contribution which was map scale factor coupled
1136 ! (had been divided by msfty)
1137 ! need to uncouple before updating uncoupled Theta (by adding)
1138    DO j=j_start, j_end
1139      DO k=1,k_end
1140      DO i=i_start, i_end
1141        t_ave(i,k,j) = t(i,k,j)
1142        t   (i,k,j) = t(i,k,j) + msfty(i,j)*dts*ft(i,k,j)
1143      END DO
1144      END DO
1145    ENDDO
1146    DO j=j_start, j_end
1147      DO i=i_start, i_end
1148        wdtn(i,1  )=0.
1149        wdtn(i,kde)=0.
1150      ENDDO
1151      DO k=2,k_end
1152      DO i=i_start, i_end
1153         ! for scalar eqn RHS term 3
1154         wdtn(i,k)= ww(i,k,j)*(fnm(k)*t_1(i,k  ,j)+fnp(k)*t_1(i,k-1,j))
1155      ENDDO
1156      ENDDO
1157 ! scalar eqn, RHS terms 1, 2 and 3
1158 ! multiply by msfty to uncouple result for Theta from map scale factor
1159      DO k=1,k_end
1160      DO i=i_start, i_end
1161        ! multiplication by msfty uncouples result for Theta
1162        t(i,k,j) = t(i,k,j) - dts*msfty(i,j)*(              &
1163                           ! multiplication by mx needed for RHS terms 1 & 2
1164                           msftx(i,j)*(                     &
1165                .5*rdy*                                     &
1166               ( v(i,k,j+1)*(t_1(i,k,j+1)+t_1(i,k, j ))     &
1167                -v(i,k,j  )*(t_1(i,k, j )+t_1(i,k,j-1)) )   &
1168              + .5*rdx*                                     &
1169               ( u(i+1,k,j)*(t_1(i+1,k,j)+t_1(i  ,k,j))     &
1170                -u(i  ,k,j)*(t_1(i  ,k,j)+t_1(i-1,k,j)) ) ) &
1171              + rdnw(k)*( wdtn(i,k+1)-wdtn(i,k) ) )
1172      ENDDO
1173      ENDDO
1174    ENDDO
1175 END SUBROUTINE advance_mu_t
1176
1177 !------------------------------------------------------------
1178 SUBROUTINE advance_w( w, rw_tend, ww, w_save, u, v,  &
1179                       mu1, mut, muave, muts,      &
1180                       c1h, c2h, c1f, c2f,         &
1181                       c3h, c4h, c3f, c4f,         &
1182                       t_2ave, t_2, t_1,           &
1183                       ph, ph_1, phb, ph_tend,     &
1184                       ht, c2a, cqw, alt, alb,     &
1185                       a, alpha, gamma,            &
1186                       rdx, rdy, dts, t0, epssm,   &
1187                       dnw, fnm, fnp, rdnw, rdn,   &
1188                       cf1, cf2, cf3, msftx, msfty,&
1189                       config_flags, top_lid,      &
1190                       ids,ide, jds,jde, kds,kde,  & ! domain dims
1191                       ims,ime, jms,jme, kms,kme,  & ! memory dims
1192                       its,ite, jts,jte, kts,kte  )  ! tile   dims
1193 ! We have used msfty for msft_inv but have not thought through w equation
1194 ! pieces properly yet, so we will have to hope that it is okay
1195 ! We think we have found a slight error in surface w calculation
1196   IMPLICIT NONE ! religion first
1197
1198 ! stuff coming in
1199   TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
1200   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1201   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1202   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1203   LOGICAL,      INTENT(IN   )    :: top_lid
1204       REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
1205             INTENT(INOUT) ::                          &
1206                                              t_2ave,  &
1207                                              w,       &
1208                                              ph
1209       REAL, DIMENSION(  ims:ime , kms:kme, jms:jme ), &
1210             INTENT(IN   ) ::                          &
1211                                              rw_tend, &
1212                                              ww,      &
1213                                              w_save,  &
1214                                              u,       &
1215                                              v,       &
1216                                              t_2,     &
1217                                              t_1,     &
1218                                              ph_1,    &
1219                                              phb,     &
1220                                              ph_tend, &
1221                                              alpha,   &
1222                                              gamma,   &
1223                                              a,       &
1224                                              c2a,     &
1225                                              cqw,     &
1226                                              alb,     &
1227                                              alt
1228       REAL, DIMENSION( ims:ime , jms:jme ), &
1229             INTENT(IN   )  ::               &
1230                                    mu1,     &
1231                                    mut,     &
1232                                    muave,   &
1233                                    muts,    &
1234                                    ht,      &
1235                                    msftx,   &
1236                                    msfty
1237       REAL, DIMENSION( kms:kme ),  INTENT(IN   )  :: fnp,     &
1238                                                      fnm,     &
1239                                                      rdnw,    &
1240                                                      rdn,     &
1241                                                      dnw
1242       REAL,   INTENT(IN   )  :: rdx,     &
1243                                 rdy,     &
1244                                 dts,     &
1245                                 cf1,     &
1246                                 cf2,     &
1247                                 cf3,     &
1248                                 t0,      &
1249                                 epssm
1250       REAL, DIMENSION( kms:kme ),INTENT(IN   ) :: c1h, c2h, c1f, c2f, &
1251                                                   c3h, c4h, c3f, c4f
1252 !  Stack based 3d data, tile size.
1253       REAL, DIMENSION( its:ite ) :: mut_inv, msft_inv
1254       REAL, DIMENSION( its:ite, kts:kte ) :: rhs, wdwn
1255       INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
1256       REAL, DIMENSION (kts:kte) :: dampwt
1257       real :: htop,hbot,hdepth,hk
1258       real    :: pi,dampmag
1259       REAL :: muthk, muthkm1
1260 !<DESCRIPTION>
1261 !
1262 !  advance_w advances the implicit w and geopotential equations.
1263 !
1264 !</DESCRIPTION>
1265 !  set loop limits.
1266 !  Currently set for periodic boundary conditions
1267       i_start = its
1268       i_end   = min(ite,ide-1)
1269       j_start = jts
1270       j_end   = min(jte,jde-1)
1271       k_start = kts
1272       k_end   = kte-1
1273       IF ( .NOT. config_flags%periodic_x )THEN
1274          IF ( config_flags%specified .or. config_flags%nested ) then
1275            i_start = max(its,ids+1)
1276            i_end   = min(ite,ide-2)
1277          ENDIF
1278       ENDIF
1279       IF ( config_flags%specified .or. config_flags%nested ) then
1280          j_start = max(jts,jds+1)
1281          j_end   = min(jte,jde-2)
1282       ENDIF
1283    pi = 4.*atan(1.)
1284    dampmag = dts*config_flags%dampcoef
1285    hdepth=config_flags%zdamp
1286 ! calculation of phi and w equations
1287 !   Comments on map scale factors:
1288 !   phi equation is:
1289 !    partial d/dt(rho phi/my) = -mx partial d/dx(phi rho u/my)
1290 !                               -mx partial d/dy(phi rho v/mx)
1291 !                               - partial d/dz(phi rho w/my) + rho g w/my
1292 !   as with scalar equation, use uncoupled value (here phi) to find the
1293 !   coupled tendency (rho phi/my)
1294 !   here as usual rho -> ~'mu'
1295 !
1296 !   w eqn [divided by my] is:
1297 !     partial d/dt(rho w/my) = -mx partial d/dx(w rho u/my)
1298 !                              -mx partial d/dy(v rho v/mx)
1299 !                              - partial d/dz(w rho w/my)
1300 !                              +rho[(u*u+v*v)/a + 2 u omega cos(lat) -
1301 !                                   (1/rho) partial dp/dz - g + Fz]/my
1302 !   here as usual rho -> ~'mu'
1303 !
1304 !  'u,v,w' sent here must be coupled variables (= rho w/my etc.) by their usage
1305     DO i=i_start, i_end
1306       rhs(i,1) = 0.
1307     ENDDO
1308   j_loop_w:  DO j = j_start, j_end
1309     DO i=i_start, i_end
1310        msft_inv(i) = 1./msfty(i,j)
1311     ENDDO
1312     DO k=1, k_end
1313     DO i=i_start, i_end
1314       t_2ave(i,k,j)=.5*((1.+epssm)*t_2(i,k,j)       &
1315                     +(1.-epssm)*t_2ave(i,k,j))
1316       t_2ave(i,k,j)=(t_2ave(i,k,j) + (c1h(k)*Muave(i,j))*t0) &
1317                     /((c1h(k)*Muts(i,j)+c2h(k))*(t0+t_1(i,k,j)))
1318       rhs(i,k+1) = dts*(ph_tend(i,k+1,j) + .5*g*(1.-epssm)*w(i,k+1,j))
1319     ENDDO
1320     ENDDO
1321 !   building up RHS of phi equation
1322 !   ph_tend contains terms 1 and 2; now adding 3 and 4 in stages:
1323 !   here rhs = delta t [ph_tend + ~g*w/2 - ~ww * partial d phi/dz]
1324     IF (config_flags%phi_adv_z == 2) THEN
1325       !  First get staggered partial d/dnu(phi) and then multiply with omega
1326       !  to avoid double staggering of omega
1327
1328       DO k=1, k_end
1329       DO i=i_start, i_end
1330          wdwn(i,k+1)=rdnw(k)*(ph_1(i,k+1,j)-ph_1(i,k,j)+phb(i,k+1,j)-phb(i,k,j))
1331       ENDDO
1332       ENDDO
1333
1334       DO k=2,k_end
1335       DO i=i_start, i_end
1336          rhs(i,k) = rhs(i,k)-dts*ww(i,k,j)*( fnm(k)*wdwn(i,k+1)  &
1337                                             +fnp(k)*wdwn(i,k  ) )
1338       ENDDO
1339       ENDDO
1340     ELSE
1341       !  First destagger omega and multiply with partial d/dnu(phi), then stagger the product
1342
1343       DO k=1, k_end
1344       DO i=i_start, i_end
1345          wdwn(i,k+1)=.5*(ww(i,k+1,j)+ww(i,k,j))*rdnw(k)    &
1346               *(ph_1(i,k+1,j)-ph_1(i,k,j)+phb(i,k+1,j)-phb(i,k,j))
1347       ENDDO
1348       ENDDO
1349
1350       DO k=2,k_end
1351       DO i=i_start, i_end
1352          rhs(i,k) = rhs(i,k)-dts*( fnm(k)*wdwn(i,k+1)  &
1353                                   +fnp(k)*wdwn(i,k  ) )
1354       ENDDO
1355       ENDDO
1356     ENDIF
1357 !  NOTE:  phi'' is not coupled with the map-scale factor  (1/m),
1358 !         but it's tendency is, so must multiply by msft here
1359 !   Comments on map scale factors:
1360 !   building up RHS of phi equation
1361 !   ph_tend contains terms 1 and 2; now adding 3 and 4 in stages:
1362 !   here rhs = ph_previous + (msft/mu)*[rhs_previous - ~ww * delta t *
1363 !                                                      partial d phi/dz]
1364 !            = ph_previous + (msft*delta t/mu)*[ph_tend + ~g*w/2 - ~ww *
1365 !                                                         partial d phi/dz]
1366     DO k=2,k_end+1
1367     DO i=i_start, i_end
1368       rhs(i,k) = ph(i,k,j) + msfty(i,j)*rhs(i,k)/(c1f(k)*mut(i,j)+c2f(k))
1369       if(top_lid .and. k.eq.k_end+1)rhs(i,k)=0.
1370     ENDDO
1371     ENDDO
1372 !  lower boundary condition on w
1373 !   Comments on map scale factors:
1374 !   Chain rule: if Z=Z(X,Y) [true at the surface] then
1375 !      dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1376 !   using capitals to denote actual values
1377 !   In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1378 !      w = dz/dt = mx u dz/dx + my v dz/dy
1379 !   [where dz/dx is just the surface height change between x
1380 !    gridpoints, and dz/dy is the change between y gridpoints]
1381 !   [cf1, cf2 and cf3 do vertical weighting of u or v values nearest
1382 !    the surface]
1383     DO i=i_start, i_end
1384        w(i,1,j)=                                           &
1385                 msfty(i,j)*.5*rdy*(                        &
1386                          (ht(i,j+1)-ht(i,j  ))             &
1387         *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1))    &
1388                         +(ht(i,j  )-ht(i,j-1))             &
1389         *(cf1*v(i,1,j  )+cf2*v(i,2,j  )+cf3*v(i,3,j  ))  ) &
1390                +msftx(i,j)*.5*rdx*(                        &
1391                          (ht(i+1,j)-ht(i,j  ))             &
1392         *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j))    &
1393                         +(ht(i,j  )-ht(i-1,j))             &
1394         *(cf1*u(i  ,1,j)+cf2*u(i  ,2,j)+cf3*u(i  ,3,j))  )
1395      ENDDO
1396 !
1397 ! Jammed 3 doubly nested loops over k/i into 1 for slight improvement
1398 ! in efficiency.  No change in results (bit-for-bit). JM 20040514
1399 ! (left a blank line where the other two k/i-loops were)
1400 !
1401 !   above surface, begin by adding delta t * previous (coupled) w tendency
1402     DO k=2,k_end
1403       DO i=i_start, i_end
1404
1405         w(i,k,j)=w(i,k,j)+dts*rw_tend(i,k,j)                       &
1406                  + msft_inv(i)*cqw(i,k,j)*(                        &
1407             +.5*dts*g*rdn(k)*                                      &
1408              (c2a(i,k  ,j)*rdnw(k  )/(c1h(k)*MUT(i,j)+c2h(k))                         &
1409         *((1.+epssm)*(rhs(i,k+1  )-rhs(i,k    ))                   &
1410          +(1.-epssm)*(ph(i,k+1,j)-ph(i,k  ,j)))                    &
1411              -c2a(i,k-1,j)*rdnw(k-1)/(c1h(k-1)*MUT(i,j)+c2h(k-1))                       &
1412         *((1.+epssm)*(rhs(i,k    )-rhs(i,k-1  ))                   &
1413          +(1.-epssm)*(ph(i,k  ,j)-ph(i,k-1,j)))))                  &
1414                 +dts*g*msft_inv(i)*(rdn(k)*                        &
1415              (c2a(i,k  ,j)*alt(i,k  ,j)*t_2ave(i,k  ,j)            &
1416              -c2a(i,k-1,j)*alt(i,k-1,j)*t_2ave(i,k-1,j))           &
1417               - (c1f(k)*muave(i,j)))
1418       ENDDO
1419     ENDDO
1420     K=k_end+1
1421     DO i=i_start, i_end
1422        w(i,k,j)=w(i,k,j)+dts*rw_tend(i,k,j)                         &
1423            +msft_inv(i)*(                                        &
1424          -.5*dts*g/(c1h(k-1)*MUT(i,j)+c2h(k-1))*rdnw(k-1)**2*2.*c2a(i,k-1,j)            &
1425              *((1.+epssm)*(rhs(i,k  )-rhs(i,k-1  ))                 &
1426               +(1.-epssm)*(ph(i,k,j)-ph(i,k-1,j)))                  &
1427          -dts*g*(2.*rdnw(k-1)*                                      &
1428                    c2a(i,k-1,j)*alt(i,k-1,j)*t_2ave(i,k-1,j)        &
1429               + (c1f(k)*muave(i,j))) )
1430        if(top_lid)w(i,k,j) = 0.
1431     ENDDO
1432
1433     DO k=2,k_end+1
1434     DO i=i_start, i_end
1435        w(i,k,j)=(w(i,k,j)-a(i,k,j)*w(i,k-1,j))*alpha(i,k,j)
1436     ENDDO
1437     ENDDO
1438
1439     DO k=k_end,2,-1
1440     DO i=i_start, i_end
1441        w (i,k,j)=w (i,k,j)-gamma(i,k,j)*w(i,k+1,j)
1442     ENDDO
1443     ENDDO
1444
1445     IF (config_flags%damp_opt .eq. 3) THEN
1446       DO k=k_end+1,2,-1
1447       DO i=i_start, i_end
1448           htop=(ph_1(i,k_end+1,j)+phb(i,k_end+1,j))/g
1449           hk=(ph_1(i,k,j)+phb(i,k,j))/g
1450           hbot=htop-hdepth
1451           dampwt(k) = 0.
1452           if(hk .ge. hbot)then
1453             dampwt(k) = dampmag*sin(0.5*pi*(hk-hbot)/hdepth)*sin(0.5*pi*(hk-hbot)/hdepth)
1454           endif
1455           w(i,k,j) = (w(i,k,j) - dampwt(k)*(c1f(k)*mut(i,j)+c2f(k))*w_save(i,k,j))/(1.+dampwt(k))
1456       ENDDO
1457       ENDDO
1458     ENDIF
1459
1460     DO k=k_end+1,2,-1
1461     DO i=i_start, i_end
1462        ph(i,k,j) = rhs(i,k)+msfty(i,j)*.5*dts*g*(1.+epssm) &
1463                       *w(i,k,j)/(c1f(k)*muts(i,j)+c2f(k))
1464     ENDDO
1465     ENDDO
1466
1467   ENDDO j_loop_w
1468
1469 END SUBROUTINE advance_w
1470
1471 !---------------------------------------------------------------------
1472
1473 SUBROUTINE sumflux ( ru, rv, ww,                             &
1474                      u_lin, v_lin, ww_lin,                   &
1475                      muu, muv,                               &
1476                      c1h, c2h, c1f, c2f,                     &
1477                      c3h, c4h, c3f, c4f,                     &
1478                      ru_m, rv_m, ww_m, epssm,                &
1479                      msfux, msfuy, msfvx, msfvx_inv, msfvy,  &
1480                      iteration , number_of_small_timesteps,  &
1481                      ids,ide, jds,jde, kds,kde,              &
1482                      ims,ime, jms,jme, kms,kme,              &
1483                      its,ite, jts,jte, kts,kte              )
1484
1485
1486   IMPLICIT NONE  ! religion first
1487
1488 ! declarations for the stuff coming in
1489
1490   INTEGER,      INTENT(IN   )    :: number_of_small_timesteps
1491   INTEGER,      INTENT(IN   )    :: iteration
1492   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1493   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1494   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1495
1496   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),  INTENT(IN   ) :: ru, &
1497                                                                 rv, &
1498                                                                 ww, &
1499                                                                 u_lin,  &
1500                                                                 v_lin,  &
1501                                                                 ww_lin
1502
1503
1504   REAL, DIMENSION(ims:ime, kms:kme, jms:jme) , INTENT(INOUT) :: ru_m, &
1505                                                                 rv_m, &
1506                                                                 ww_m
1507   REAL, DIMENSION(ims:ime, jms:jme) , INTENT(IN   ) :: muu, muv,      &
1508                                                        msfux, msfuy,  &
1509                                                        msfvx, msfvy, msfvx_inv
1510
1511   REAL, DIMENSION( kms:kme ),INTENT(IN   ) :: c1h, c2h, c1f, c2f, &
1512                                               c3h, c4h, c3f, c4f
1513
1514   INTEGER :: mini, minj, mink
1515
1516
1517   REAL, INTENT(IN   )  ::  epssm
1518   INTEGER   :: i,j,k
1519
1520
1521 !<DESCRIPTION>
1522 !
1523 !  update the small-timestep time-averaged mass fluxes;  these
1524 !  are needed for consistent mass-conserving scalar advection.
1525 !
1526 !</DESCRIPTION>
1527
1528     IF (iteration == 1 )THEN
1529       DO  j = jts, jte
1530       DO  k = kts, kte
1531       DO  i = its, ite
1532         ru_m(i,k,j)  = 0.
1533         rv_m(i,k,j)  = 0.
1534         ww_m(i,k,j)  = 0.
1535       ENDDO
1536       ENDDO
1537       ENDDO
1538     ENDIF
1539
1540   mini = min(ide-1,ite)
1541   minj = min(jde-1,jte)
1542   mink = min(kde-1,kte)
1543
1544
1545     DO  j = jts, minj
1546     DO  k = kts, mink
1547     DO  i = its, mini
1548       ru_m(i,k,j)  = ru_m(i,k,j) + ru(i,k,j)
1549       rv_m(i,k,j)  = rv_m(i,k,j) + rv(i,k,j)
1550       ww_m(i,k,j)  = ww_m(i,k,j) + ww(i,k,j)
1551     ENDDO
1552     ENDDO
1553     ENDDO
1554
1555     IF (ite .GT. mini) THEN
1556       DO  j = jts, minj
1557       DO  k = kts, mink
1558       DO  i = mini+1, ite
1559         ru_m(i,k,j)  = ru_m(i,k,j) + ru(i,k,j)
1560       ENDDO
1561       ENDDO
1562       ENDDO
1563     END IF
1564     IF (jte .GT. minj) THEN
1565       DO  j = minj+1, jte
1566       DO  k = kts, mink
1567       DO  i = its, mini
1568  rv_m(i,k,j)  = rv_m(i,k,j) + rv(i,k,j)
1569       ENDDO
1570       ENDDO
1571       ENDDO
1572     END IF
1573     IF ( kte .GT. mink) THEN
1574       DO  j = jts, minj
1575       DO  k = mink+1, kte
1576       DO  i = its, mini
1577         ww_m(i,k,j)  = ww_m(i,k,j) + ww(i,k,j)
1578       ENDDO
1579       ENDDO
1580       ENDDO
1581     END IF
1582
1583   IF (iteration == number_of_small_timesteps) THEN
1584
1585     DO  j = jts, minj
1586     DO  k = kts, mink
1587     DO  i = its, mini
1588       ru_m(i,k,j)  = ru_m(i,k,j) / number_of_small_timesteps   &
1589                      + (c1h(k)*muu(i,j)+c2h(k))*u_lin(i,k,j)/msfuy(i,j)
1590       rv_m(i,k,j)  = rv_m(i,k,j) / number_of_small_timesteps   &
1591                      + (c1h(k)*muv(i,j)+c2h(k))*v_lin(i,k,j)*msfvx_inv(i,j)
1592       ww_m(i,k,j)  = ww_m(i,k,j) / number_of_small_timesteps   &
1593                      + ww_lin(i,k,j)
1594     ENDDO
1595     ENDDO
1596     ENDDO
1597
1598
1599     IF (ite .GT. mini) THEN
1600       DO  j = jts, minj
1601       DO  k = kts, mink
1602       DO  i = mini+1, ite
1603         ru_m(i,k,j)  = ru_m(i,k,j) / number_of_small_timesteps   &
1604                      + (c1h(k)*muu(i,j)+c2h(k))*u_lin(i,k,j)/msfuy(i,j)
1605       ENDDO
1606       ENDDO
1607       ENDDO
1608     END IF
1609     IF (jte .GT. minj) THEN
1610       DO  j = minj+1, jte
1611       DO  k = kts, mink
1612       DO  i = its, mini
1613  rv_m(i,k,j)  = rv_m(i,k,j) / number_of_small_timesteps   &
1614                      + (c1h(k)*muv(i,j)+c2h(k))*v_lin(i,k,j)*msfvx_inv(i,j)
1615       ENDDO
1616       ENDDO
1617       ENDDO
1618     END IF
1619     IF ( kte .GT. mink) THEN
1620       DO  j = jts, minj
1621       DO  k = mink+1, kte
1622       DO  i = its, mini
1623         ww_m(i,k,j)  = ww_m(i,k,j) / number_of_small_timesteps   &
1624                      + ww_lin(i,k,j)
1625       ENDDO
1626       ENDDO
1627       ENDDO
1628     END IF
1629
1630   ENDIF
1631
1632
1633 END SUBROUTINE sumflux
1634
1635 !---------------------------------------------------------------------
1636
1637   SUBROUTINE init_module_small_step
1638   END SUBROUTINE init_module_small_step
1639
1640 #if ( WRFPLUS == 1 )
1641 !----------------------------------------------------------------------
1642 SUBROUTINE advance_all (u, ru_tend,  v, rv_tend,       &
1643                         w, rw_tend,  t,  t_tend,       &
1644                         mu, mu_tend, ph, ph_tend,      &
1645                         muu, muv, mut,                 &
1646                         msfuy, msfvx, msfty,           &
1647                         dts,                           &
1648                         config_flags, spec_zone,       &
1649                         ids, ide, jds, jde, kds, kde,  &
1650                         ims, ime, jms, jme, kms, kme,  &
1651                         ips, ipe, jps, jpe, kps, kpe,  &
1652                         its, ite, jts, jte, kts, kte  )
1653
1654       IMPLICIT NONE  ! religion first
1655
1656       TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
1657
1658       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1659       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1660       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1661       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1662       INTEGER,      INTENT(IN   )    :: spec_zone
1663
1664       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),  &
1665             INTENT(INOUT) ::                          &
1666                                                   u,  &
1667                                                   v,  &
1668                                                   w,  &
1669                                                   t,  &
1670                                                   ph
1671
1672       REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
1673             INTENT(IN   ) ::                          &
1674                                              ru_tend, &
1675                                              rv_tend, &
1676                                              rw_tend, &
1677                                               t_tend, &
1678                                              ph_tend
1679
1680       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(INOUT) :: mu
1681       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: mu_tend
1682
1683       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: muu,  &
1684                                                                 muv,  &
1685                                                                 mut
1686
1687       REAL, DIMENSION( ims:ime , jms:jme ),    INTENT(IN   ) :: msfuy,  &
1688                                                                 msfvx,  &
1689                                                                 msfty
1690
1691       REAL,                                    INTENT(IN   ) :: dts
1692
1693
1694       INTEGER :: i,j,k, i_start, i_end, j_start, j_end, k_start, k_end
1695       INTEGER :: i_endu, j_endv, k_endw
1696
1697       INTEGER :: i_start_u_tend, i_end_u_tend, j_start_v_tend, j_end_v_tend
1698
1699 !<DESCRIPTION>
1700 !
1701 !  advance_all advances the explicit perturbation horizontal momentum
1702 !  equations (u,v) by adding in the large-timestep tendency along with
1703 !  the small timestep pressure gradient tendency.
1704 !
1705 !</DESCRIPTION>
1706
1707 !  now, the real work.
1708 !  set the loop bounds taking into account boundary conditions.
1709
1710    IF( config_flags%nested .or. config_flags%specified ) THEN
1711       i_start = max( its,ids+spec_zone )
1712       i_end   = min( ite,ide-spec_zone-1 )
1713       j_start = max( jts,jds+spec_zone )
1714       j_end   = min( jte,jde-spec_zone-1 )
1715       k_start = kts
1716       k_end   = min( kte, kde-1 )
1717
1718       i_endu = min( ite,ide-spec_zone )
1719       j_endv = min( jte,jde-spec_zone )
1720       k_endw = k_end
1721
1722       IF( config_flags%periodic_x) THEN
1723         i_start = its
1724         i_end   = min(ite,ide-1)
1725         i_endu = ite
1726       ENDIF
1727    ELSE
1728       i_start = its
1729       i_end   = min(ite,ide-1)
1730       j_start = jts
1731       j_end   = min(jte,jde-1)
1732       k_start = kts
1733       k_end   = kte-1
1734
1735       i_endu = ite
1736       j_endv = jte
1737       k_endw = k_end
1738    ENDIF
1739
1740    i_start_u_tend = i_start
1741    i_end_u_tend   = i_endu
1742    j_start_v_tend = j_start
1743    j_end_v_tend   = j_endv
1744
1745    IF ( config_flags%symmetric_xs .and. (its == ids) ) &
1746                   i_start_u_tend = i_start_u_tend+1
1747    IF ( config_flags%symmetric_xe .and. (ite == ide) ) &
1748                   i_end_u_tend = i_end_u_tend-1
1749    IF ( config_flags%symmetric_ys .and. (jts == jds) ) &
1750                   j_start_v_tend = j_start_v_tend+1
1751    IF ( config_flags%symmetric_ye .and. (jte == jde) ) &
1752                   j_end_v_tend = j_end_v_tend-1
1753
1754    u_outer_j_loop: DO j = j_start, j_end
1755       DO k = k_start, k_end
1756       DO i = i_start_u_tend, i_end_u_tend
1757         u(i,k,j) = u(i,k,j) + dts*ru_tend(i,k,j)*msfuy(i,j)/muu(i,j)
1758       ENDDO
1759       ENDDO
1760    ENDDO u_outer_j_loop
1761
1762    v_outer_j_loop: DO j = j_start_v_tend, j_end_v_tend
1763       DO k = k_start, k_end
1764       DO i = i_start, i_end
1765         v(i,k,j) = v(i,k,j) + dts*rv_tend(i,k,j)*msfvx(i,j)/muv(i,j)
1766       ENDDO
1767       ENDDO
1768    ENDDO v_outer_j_loop
1769
1770    i_start = its
1771    i_end   = min(ite,ide-1)
1772    j_start = jts
1773    j_end   = min(jte,jde-1)
1774    k_start = kts
1775    k_end   = kte-1
1776    IF ( .NOT. config_flags%periodic_x )THEN
1777      IF ( config_flags%specified .or. config_flags%nested ) then
1778        i_start = max(its,ids+1)
1779        i_end   = min(ite,ide-2)
1780      ENDIF
1781    ENDIF
1782    IF ( config_flags%specified .or. config_flags%nested ) then
1783      j_start = max(jts,jds+1)
1784      j_end   = min(jte,jde-2)
1785    ENDIF
1786
1787    DO j = j_start, j_end
1788    DO i=i_start, i_end
1789      mu(i,j) = mu(i,j)+dts*mu_tend(i,j)
1790    ENDDO
1791    ENDDO
1792
1793    DO j=j_start, j_end
1794      DO k=1,k_end
1795      DO i=i_start, i_end
1796        t(i,k,j) = t(i,k,j) + msfty(i,j)*dts*t_tend(i,k,j)/mut(i,j)
1797      END DO
1798      END DO
1799    ENDDO
1800
1801    i_start = its
1802    i_end   = min(ite,ide-1)
1803    j_start = jts
1804    j_end   = min(jte,jde-1)
1805    k_start = kts
1806    k_end   = kte-1
1807    IF ( .NOT. config_flags%periodic_x )THEN
1808       IF ( config_flags%specified .or. config_flags%nested ) then
1809         i_start = max(its,ids+1)
1810         i_end   = min(ite,ide-2)
1811       ENDIF
1812    ENDIF
1813    IF ( config_flags%specified .or. config_flags%nested ) then
1814       j_start = max(jts,jds+1)
1815       j_end   = min(jte,jde-2)
1816    ENDIF
1817
1818 IF ( config_flags%non_hydrostatic ) THEN
1819   j_loop_w:  DO j = j_start, j_end
1820     DO k=2,k_end+1
1821       DO i=i_start, i_end
1822         w(i,k,j)=w(i,k,j)+dts*rw_tend(i,k,j)*msfty(i,j)/mut(i,j)
1823         ph(i,k,j) = ph(i,k,j)+dts*ph_tend(i,k,j)*msfty(i,j)/mut(i,j)
1824       ENDDO
1825     ENDDO
1826   ENDDO j_loop_w
1827 ENDIF
1828
1829 END SUBROUTINE advance_all
1830
1831 !----------------------------------------------------------------------
1832
1833 SUBROUTINE save_ph_mu ( ph_1, ph_2, ph_save, &
1834                         mu_1, mu_2, mu_save, &
1835                         rk_step,                      &
1836                         ids,ide, jds,jde, kds,kde,    &
1837                         ims,ime, jms,jme, kms,kme,    &
1838                         its,ite, jts,jte, kts,kte )
1839
1840   IMPLICIT NONE  ! religion first
1841
1842 ! declarations for the stuff coming in
1843
1844   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1845   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1846   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1847
1848   INTEGER,      INTENT(IN   )    :: rk_step
1849
1850   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph_1
1851   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(  OUT) :: ph_save
1852   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph_2
1853   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_1,mu_2
1854   REAL, DIMENSION(ims:ime, jms:jme), INTENT(  OUT) :: mu_save
1855
1856 ! local variables
1857
1858   INTEGER :: i, j, k
1859   INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
1860
1861
1862     i_start = its
1863     i_end   = min(ite,ide-1)
1864     j_start = jts
1865     j_end   = min(jte,jde-1)
1866     k_start = kts
1867     k_end = min(kte,kde-1)
1868
1869     !  if this is the first RK step, reset *_1 to *_2
1870     !  (we are replacing the t-dt fields with the time t fields)
1871
1872     IF ((rk_step == 1) ) THEN
1873
1874       DO j=j_start, j_end
1875       DO k=k_start, min(kde,kte)
1876       DO i=i_start, i_end
1877         ph_1(i,k,j) = ph_2(i,k,j)
1878       ENDDO
1879       ENDDO
1880       ENDDO
1881
1882       DO j=j_start, j_end
1883       DO i=i_start, i_end
1884         mu_1(i,j)=mu_2(i,j)
1885       ENDDO
1886       ENDDO
1887
1888      END IF
1889
1890     ! set up the small timestep variables
1891
1892     DO j=j_start, j_end
1893 !    DO k=k_start, min(kde,kte)
1894     DO k=k_start, kde
1895     DO i=i_start, i_end
1896       ph_save(i,k,j) = ph_2(i,k,j)
1897       ph_2(i,k,j) = ph_1(i,k,j)-ph_2(i,k,j)
1898     ENDDO
1899     ENDDO
1900     ENDDO
1901
1902     DO j=j_start, j_end
1903     DO i=i_start, i_end
1904       mu_save(i,j)=mu_2(i,j)
1905       mu_2(i,j)=mu_1(i,j)-mu_2(i,j)
1906     ENDDO
1907     ENDDO
1908
1909 END SUBROUTINE save_ph_mu
1910
1911 !----------------------------------------------------------------------
1912
1913 SUBROUTINE restore_ph_mu ( ph_2, ph_save, &
1914                            mu_2, mu_save, &
1915                            ids,ide, jds,jde, kds,kde,    &
1916                            ims,ime, jms,jme, kms,kme,    &
1917                            its,ite, jts,jte, kts,kte )
1918
1919   IMPLICIT NONE  ! religion first
1920
1921 ! declarations for the stuff coming in
1922
1923   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1924   INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1925   INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1926
1927   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(IN   ) :: ph_save
1928   REAL, DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(INOUT) :: ph_2
1929   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: mu_2
1930   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN   ) :: mu_save
1931
1932 ! local variables
1933
1934   INTEGER :: i, j, k
1935   INTEGER :: i_start, i_end, j_start, j_end
1936
1937
1938   i_start = its
1939   i_end   = min(ite,ide-1)
1940   j_start = jts
1941   j_end   = min(jte,jde-1)
1942
1943   DO j = j_start, j_end
1944   DO k = kds, kde
1945   DO i = i_start, i_end
1946     ph_2(i,k,j) = ph_2(i,k,j) + ph_save(i,k,j)
1947   ENDDO
1948   ENDDO
1949   ENDDO
1950
1951   DO j = j_start, j_end
1952   DO i = i_start, i_end
1953     mu_2(i,j) = mu_2(i,j) + mu_save(i,j)
1954   ENDDO
1955   ENDDO
1956
1957 END SUBROUTINE restore_ph_mu
1958 !----------------------------------------------------------------------
1959 #endif
1960
1961 END MODULE module_small_step_em