wrfv2_fire/dyn_em/module_em.F

   1 !WRF:MODEL_LAYER:DYNAMICS
   2 !
   3
   4 MODULE module_em
   5
   6    USE module_model_constants
   7    USE module_advect_em
   8    USE module_big_step_utilities_em
   9    USE module_state_description
  10    USE module_damping_em
  11
  12 CONTAINS
  13
  14 !------------------------------------------------------------------------
  15
  16 SUBROUTINE rk_step_prep  ( config_flags, rk_step,           &
  17                            u, v, w, t, ph, mu,              &
  18                            moist,                           &
  19                            ru, rv, rw, ww, php, alt,        &
  20                            muu, muv,                        &
  21                            mub, mut, phb, pb, p, al, alb,   &
  22                            cqu, cqv, cqw,                   &
  23                            msfux, msfuy,                    &
  24                            msfvx, msfvx_inv, msfvy,         &
  25                            msftx, msfty,                    &
  26                            fnm, fnp, dnw, rdx, rdy,         &
  27                            n_moist,                         &
  28                            ids, ide, jds, jde, kds, kde,    &
  29                            ims, ime, jms, jme, kms, kme,    &
  30                            its, ite, jts, jte, kts, kte    )
  31
  32    IMPLICIT NONE
  33
  34
  35    !  Input data.
  36
  37    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
  38
  39    INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
  40                                     ims, ime, jms, jme, kms, kme, &
  41                                     its, ite, jts, jte, kts, kte
  42
  43    INTEGER ,       INTENT(IN   ) :: n_moist, rk_step
  44
  45    REAL ,          INTENT(IN   ) :: rdx, rdy
  46
  47    REAL , DIMENSION(  ims:ime , kms:kme, jms:jme ) ,                      &
  48                                                INTENT(IN   ) ::  u,       &
  49                                                                  v,       &
  50                                                                  w,       &
  51                                                                  t,       &
  52                                                                  ph,      &
  53                                                                  phb,     &
  54                                                                  pb,      &
  55                                                                  al,      &
  56                                                                  alb
  57
  58    REAL , DIMENSION( ims:ime , kms:kme , jms:jme  ) ,                     &
  59                                                INTENT(  OUT) ::  ru,      &
  60                                                                  rv,      &
  61                                                                  rw,      &
  62                                                                  ww,      &
  63                                                                  php,     &
  64                                                                  cqu,     &
  65                                                                  cqv,     &
  66                                                                  cqw,     &
  67                                                                  alt
  68
  69    REAL , DIMENSION(  ims:ime , kms:kme, jms:jme ) ,                      &
  70                                                INTENT(IN   ) ::  p
  71
  72
  73
  74
  75    REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist ), INTENT(   IN) :: &
  76                                                            moist
  77
  78    REAL , DIMENSION( ims:ime , jms:jme ) ,    INTENT(IN   ) :: msftx,     &
  79                                                                msfty,     &
  80                                                                msfux,     &
  81                                                                msfuy,     &
  82                                                                msfvx,     &
  83                                                                msfvx_inv, &
  84                                                                msfvy,     &
  85                                                                mu,        &
  86                                                                mub
  87
  88    REAL , DIMENSION( ims:ime , jms:jme ) ,    INTENT(  OUT) :: muu,    &
  89                                                                muv,    &
  90                                                                mut
  91
  92    REAL , DIMENSION( kms:kme ) ,    INTENT(IN   ) :: fnm, fnp, dnw
  93
  94    integer :: k
  95
  96
  97 !<DESCRIPTION>
  98 !
  99 !  rk_step_prep prepares a number of diagnostic quantities
 100 !  in preperation for a Runge-Kutta timestep.  subroutines called
 101 !  by rk_step_prep calculate
 102 !
 103 !  (1) total column dry air mass (mut, call to calculate_full)
 104 !
 105 !  (2) total column dry air mass at u and v points
 106 !      (muu, muv, call to calculate_mu_uv)
 107 !
 108 !  (3) mass-coupled velocities for advection
 109 !      (ru, rv, and rw, call to couple_momentum)
 110 !
 111 !  (4) omega (call to calc_ww_cp)
 112 !
 113 !  (5) moisture coefficients (cqu, cqv, cqw, call to calc_cq)
 114 !
 115 !  (6) inverse density (alt, call to calc_alt)
 116 !
 117 !  (7) geopotential at pressure points (php, call to calc_php)
 118 !
 119 !</DESCRIPTION>
 120
 121    CALL calculate_full( mut, mub, mu,             &
 122                         ids, ide, jds, jde, 1, 2, &
 123                         ims, ime, jms, jme, 1, 1, &
 124                         its, ite, jts, jte, 1, 1 )
 125
 126    CALL calc_mu_uv ( config_flags,                  &
 127                      mu, mub, muu, muv,             &
 128                      ids, ide, jds, jde, kds, kde,  &
 129                      ims, ime, jms, jme, kms, kme,  &
 130                      its, ite, jts, jte, kts, kte  )
 131
 132    CALL couple_momentum( muu, ru, u, msfuy,             &
 133                          muv, rv, v, msfvx, msfvx_inv,  &
 134                          mut, rw, w, msfty,             &
 135                          ids, ide, jds, jde, kds, kde,  &
 136                          ims, ime, jms, jme, kms, kme,  &
 137                          its, ite, jts, jte, kts, kte  )
 138
 139 !  new call, couples V with mu, also has correct map factors.  WCS, 3 june 2001
 140    CALL calc_ww_cp ( u, v, mu, mub, ww,               &
 141                      rdx, rdy, msftx, msfty,          &
 142                      msfux, msfuy, msfvx, msfvx_inv,  &
 143                      msfvy, dnw,                      &
 144                      ids, ide, jds, jde, kds, kde,    &
 145                      ims, ime, jms, jme, kms, kme,    &
 146                      its, ite, jts, jte, kts, kte    )
 147
 148    CALL calc_cq ( moist, cqu, cqv, cqw, n_moist, &
 149                   ids, ide, jds, jde, kds, kde,  &
 150                   ims, ime, jms, jme, kms, kme,  &
 151                   its, ite, jts, jte, kts, kte  )
 152
 153    CALL calc_alt ( alt, al, alb,                 &
 154                    ids, ide, jds, jde, kds, kde, &
 155                    ims, ime, jms, jme, kms, kme, &
 156                    its, ite, jts, jte, kts, kte )
 157
 158    CALL calc_php ( php, ph, phb,                 &
 159                    ids, ide, jds, jde, kds, kde, &
 160                    ims, ime, jms, jme, kms, kme, &
 161                    its, ite, jts, jte, kts, kte )
 162
 163 END SUBROUTINE rk_step_prep
 164
 165 !-------------------------------------------------------------------------------
 166
 167 SUBROUTINE rk_tendency ( config_flags, rk_step,                           &
 168                          ru_tend, rv_tend, rw_tend, ph_tend, t_tend,      &
 169                          ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf, &
 170                          mu_tend, u_save, v_save, w_save, ph_save,        &
 171                          t_save, mu_save, RTHFTEN,                        &
 172                          ru, rv, rw, ww,                                  &
 173                          u, v, w, t, ph,                                  &
 174                          u_old, v_old, w_old, t_old, ph_old,              &
 175                          h_diabatic, phb,t_init,                          &
 176                          mu, mut, muu, muv, mub,                          &
 177                          al, alt, p, pb, php, cqu, cqv, cqw,              &
 178                          u_base, v_base, t_base, qv_base, z_base,         &
 179                          msfux, msfuy, msfvx, msfvx_inv,                  &
 180                          msfvy, msftx, msfty,                             &
 181                          xlat, f, e, sina, cosa,                          &
 182                          fnm, fnp, rdn, rdnw,                             &
 183                          dt, rdx, rdy, khdif, kvdif, xkmhd, xkhh,         &
 184                          diff_6th_opt, diff_6th_factor,                   &
 185                          dampcoef,zdamp,damp_opt,                         &
 186                          cf1, cf2, cf3, cfn, cfn1, n_moist,               &
 187                          non_hydrostatic, top_lid,                        &
 188                          u_frame, v_frame,                                &
 189                          ids, ide, jds, jde, kds, kde,                    &
 190                          ims, ime, jms, jme, kms, kme,                    &
 191                          its, ite, jts, jte, kts, kte,                    &
 192                          max_vert_cfl, max_horiz_cfl)
 193
 194    IMPLICIT NONE
 195
 196    !  Input data.
 197
 198    TYPE(grid_config_rec_type)    ,           INTENT(IN   ) :: config_flags
 199
 200    INTEGER ,               INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 201                                             ims, ime, jms, jme, kms, kme, &
 202                                             its, ite, jts, jte, kts, kte
 203
 204    LOGICAL ,               INTENT(IN   ) :: non_hydrostatic, top_lid
 205
 206    INTEGER ,               INTENT(IN   ) :: n_moist, rk_step
 207
 208    REAL , DIMENSION( ims:ime , kms:kme, jms:jme  ) ,              &
 209                                         INTENT(IN   ) :: ru,      &
 210                                                          rv,      &
 211                                                          rw,      &
 212                                                          ww,      &
 213                                                          u,       &
 214                                                          v,       &
 215                                                          w,       &
 216                                                          t,       &
 217                                                          ph,      &
 218                                                          u_old,   &
 219                                                          v_old,   &
 220                                                          w_old,   &
 221                                                          t_old,   &
 222                                                          ph_old,  &
 223                                                          phb,     &
 224                                                          al,      &
 225                                                          alt,     &
 226                                                          p,       &
 227                                                          pb,      &
 228                                                          php,     &
 229                                                          cqu,     &
 230                                                          cqv,     &
 231                                                          t_init,  &
 232                                                          xkmhd,   &
 233                                                          xkhh,    &
 234                                                          h_diabatic
 235
 236    REAL , DIMENSION( ims:ime , kms:kme, jms:jme  ) ,              &
 237                                         INTENT(OUT  ) :: ru_tend, &
 238                                                          rv_tend, &
 239                                                          rw_tend, &
 240                                                          t_tend,  &
 241                                                          ph_tend, &
 242                                                          RTHFTEN, &
 243                                                           u_save, &
 244                                                           v_save, &
 245                                                           w_save, &
 246                                                          ph_save, &
 247                                                           t_save
 248
 249    REAL , DIMENSION( ims:ime , kms:kme, jms:jme  ) ,               &
 250                                         INTENT(INOUT) :: ru_tendf, &
 251                                                          rv_tendf, &
 252                                                          rw_tendf, &
 253                                                          t_tendf,  &
 254                                                          ph_tendf, &
 255                                                          cqw
 256
 257    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(  OUT) :: mu_tend, &
 258                                                                     mu_save
 259
 260    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,   &
 261                                                                     msfuy,   &
 262                                                                     msfvx,   &
 263                                                                     msfvx_inv,   &
 264                                                                     msfvy,   &
 265                                                                     msftx,   &
 266                                                                     msfty,   &
 267                                                                     xlat,    &
 268                                                                     f,       &
 269                                                                     e,       &
 270                                                                     sina,    &
 271                                                                     cosa,    &
 272                                                                     mu,      &
 273                                                                     mut,     &
 274                                                                     mub,     &
 275                                                                     muu,     &
 276                                                                     muv
 277
 278    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fnm,     &
 279                                                                   fnp,     &
 280                                                                   rdn,     &
 281                                                                   rdnw,    &
 282                                                                   u_base,  &
 283                                                                   v_base,  &
 284                                                                   t_base,  &
 285                                                                   qv_base, &
 286                                                                   z_base
 287
 288    REAL ,                                      INTENT(IN   ) :: rdx,     &
 289                                                                 rdy,     &
 290                                                                 dt,      &
 291                                                                 u_frame, &
 292                                                                 v_frame, &
 293                                                                 khdif,   &
 294                                                                 kvdif
 295    INTEGER, INTENT( IN ) :: diff_6th_opt
 296    REAL,    INTENT( IN ) :: diff_6th_factor
 297
 298    INTEGER, INTENT( IN ) :: damp_opt
 299
 300    REAL, INTENT( IN ) :: zdamp, dampcoef
 301
 302    REAL, INTENT( OUT ) :: max_horiz_cfl
 303    REAL, INTENT( OUT ) :: max_vert_cfl
 304
 305    REAL    :: kdift, khdq, kvdq, cfn, cfn1, cf1, cf2, cf3
 306    INTEGER :: i,j,k
 307    INTEGER :: time_step
 308
 309 !<DESCRIPTION>
 310 !
 311 !  rk_tendency computes the large-timestep tendency terms in the
 312 !  momentum, thermodynamic (theta), and geopotential equations.
 313 !  These terms include:
 314 !
 315 !  (1) advection (for u, v, w, theta - calls to advect_u, advect_v,
 316 !                 advect_w, and advact_scalar).
 317 !
 318 !  (2) geopotential equation terms (advection and "gw" - call to rhs_ph).
 319 !
 320 !  (3) buoyancy term in vertical momentum equation (call to pg_buoy_w).
 321 !
 322 !  (4) Coriolis and curvature terms in u,v,w momentum equations
 323 !      (calls to subroutines coriolis, curvature)
 324 !
 325 !  (5) 3D diffusion on coordinate surfaces.
 326 !
 327 !</DESCRIPTION>
 328
 329    CALL zero_tend ( ru_tend,                      &
 330                     ids, ide, jds, jde, kds, kde, &
 331                     ims, ime, jms, jme, kms, kme, &
 332                     its, ite, jts, jte, kts, kte )
 333
 334    CALL zero_tend ( rv_tend,                      &
 335                     ids, ide, jds, jde, kds, kde, &
 336                     ims, ime, jms, jme, kms, kme, &
 337                     its, ite, jts, jte, kts, kte )
 338
 339    CALL zero_tend ( rw_tend,                      &
 340                     ids, ide, jds, jde, kds, kde, &
 341                     ims, ime, jms, jme, kms, kme, &
 342                     its, ite, jts, jte, kts, kte )
 343
 344    CALL zero_tend ( t_tend,                       &
 345                     ids, ide, jds, jde, kds, kde, &
 346                     ims, ime, jms, jme, kms, kme, &
 347                     its, ite, jts, jte, kts, kte )
 348
 349    CALL zero_tend ( ph_tend,                      &
 350                     ids, ide, jds, jde, kds, kde, &
 351                     ims, ime, jms, jme, kms, kme, &
 352                     its, ite, jts, jte, kts, kte )
 353
 354    CALL zero_tend ( u_save,                       &
 355                     ids, ide, jds, jde, kds, kde, &
 356                     ims, ime, jms, jme, kms, kme, &
 357                     its, ite, jts, jte, kts, kte )
 358
 359    CALL zero_tend ( v_save,                       &
 360                     ids, ide, jds, jde, kds, kde, &
 361                     ims, ime, jms, jme, kms, kme, &
 362                     its, ite, jts, jte, kts, kte )
 363
 364    CALL zero_tend ( w_save,                       &
 365                     ids, ide, jds, jde, kds, kde, &
 366                     ims, ime, jms, jme, kms, kme, &
 367                     its, ite, jts, jte, kts, kte )
 368
 369    CALL zero_tend ( ph_save,                       &
 370                     ids, ide, jds, jde, kds, kde, &
 371                     ims, ime, jms, jme, kms, kme, &
 372                     its, ite, jts, jte, kts, kte )
 373
 374    CALL zero_tend ( t_save,                       &
 375                     ids, ide, jds, jde, kds, kde, &
 376                     ims, ime, jms, jme, kms, kme, &
 377                     its, ite, jts, jte, kts, kte )
 378
 379    CALL zero_tend ( mu_tend,                  &
 380                     ids, ide, jds, jde, 1, 1, &
 381                     ims, ime, jms, jme, 1, 1, &
 382                     its, ite, jts, jte, 1, 1 )
 383
 384    CALL zero_tend ( mu_save,                  &
 385                     ids, ide, jds, jde, 1, 1, &
 386                     ims, ime, jms, jme, 1, 1, &
 387                     its, ite, jts, jte, 1, 1 )
 388
 389      !  advection tendencies
 390      CALL nl_get_time_step ( 1, time_step )
 391
 392      CALL advect_u ( u, u , ru_tend, ru, rv, ww,   &
 393                      mut, time_step, config_flags, &
 394                      msfux, msfuy, msfvx, msfvy,   &
 395                      msftx, msfty,                 &
 396                      fnm, fnp, rdx, rdy, rdnw,     &
 397                      ids, ide, jds, jde, kds, kde, &
 398                      ims, ime, jms, jme, kms, kme, &
 399                      its, ite, jts, jte, kts, kte )
 400
 401      CALL advect_v ( v, v , rv_tend, ru, rv, ww,   &
 402                      mut, time_step, config_flags, &
 403                      msfux, msfuy, msfvx, msfvy,   &
 404                      msftx, msfty,                 &
 405                      fnm, fnp, rdx, rdy, rdnw,     &
 406                      ids, ide, jds, jde, kds, kde, &
 407                      ims, ime, jms, jme, kms, kme, &
 408                      its, ite, jts, jte, kts, kte )
 409
 410      IF (non_hydrostatic)                          &
 411      CALL advect_w ( w, w, rw_tend, ru, rv, ww,    &
 412                      mut, time_step, config_flags, &
 413                      msfux, msfuy, msfvx, msfvy,   &
 414                      msftx, msfty,                 &
 415                      fnm, fnp, rdx, rdy, rdn,      &
 416                      ids, ide, jds, jde, kds, kde, &
 417                      ims, ime, jms, jme, kms, kme, &
 418                      its, ite, jts, jte, kts, kte )
 419
 420 !  theta flux divergence
 421
 422      CALL advect_scalar ( t, t, t_tend, ru, rv, ww,     &
 423                           mut, time_step, config_flags, &
 424                           msfux, msfuy, msfvx, msfvy,   &
 425                           msftx, msfty, fnm, fnp,       &
 426                           rdx, rdy, rdnw,               &
 427                           ids, ide, jds, jde, kds, kde, &
 428                           ims, ime, jms, jme, kms, kme, &
 429                           its, ite, jts, jte, kts, kte )
 430
 431      IF ( config_flags%cu_physics == GDSCHEME  .OR.     &
 432           config_flags%cu_physics == G3SCHEME ) THEN
 433
 434      ! theta advection only:
 435
 436          CALL set_tend( RTHFTEN, t_tend, msfty,          &
 437                         ids, ide, jds, jde, kds, kde,    &
 438                         ims, ime, jms, jme, kms, kme,    &
 439                         its, ite, jts, jte, kts, kte     )
 440
 441      END IF
 442
 443      CALL rhs_ph( ph_tend, u, v, ww, ph, ph, phb, w, &
 444                   mut, muu, muv,                     &
 445                   fnm, fnp,                          &
 446                   rdnw, cfn, cfn1, rdx, rdy,         &
 447                   msfux, msfuy, msfvx,               &
 448                   msfvx_inv, msfvy,                  &
 449                   msftx, msfty,                      &
 450                   non_hydrostatic,                   &
 451                   config_flags,                      &
 452                   ids, ide, jds, jde, kds, kde,      &
 453                   ims, ime, jms, jme, kms, kme,      &
 454                   its, ite, jts, jte, kts, kte      )
 455
 456   CALL horizontal_pressure_gradient( ru_tend,rv_tend,                &
 457                                      ph,alt,p,pb,al,php,cqu,cqv,     &
 458                                      muu,muv,mu,fnm,fnp,rdnw,        &
 459                                      cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
 460                                      msfvx,msfvy,msftx,msfty,        &
 461                                      config_flags, non_hydrostatic,  &
 462                                      top_lid,                        &
 463                                      ids, ide, jds, jde, kds, kde,   &
 464                                      ims, ime, jms, jme, kms, kme,   &
 465                                      its, ite, jts, jte, kts, kte   )
 466
 467   IF (non_hydrostatic)                            &
 468   CALL pg_buoy_w( rw_tend, p, cqw, mu, mub,       &
 469                   rdnw, rdn, g, msftx, msfty,     &
 470                   ids, ide, jds, jde, kds, kde,   &
 471                   ims, ime, jms, jme, kms, kme,   &
 472                   its, ite, jts, jte, kts, kte   )
 473
 474   CALL w_damp   ( rw_tend, max_vert_cfl,            &
 475                     max_horiz_cfl,                  &
 476                     u, v, ww, w, mut, rdnw,         &
 477                     rdx, rdy, msfux, msfuy, msfvx,  &
 478                     msfvy, dt, config_flags,        &
 479                     ids, ide, jds, jde, kds, kde,   &
 480                     ims, ime, jms, jme, kms, kme,   &
 481                     its, ite, jts, jte, kts, kte   )
 482
 483   IF(config_flags%pert_coriolis) THEN
 484
 485     CALL perturbation_coriolis ( ru, rv, rw,                   &
 486                                  ru_tend,  rv_tend,  rw_tend,  &
 487                                  config_flags,                 &
 488                                  u_base, v_base, z_base,       &
 489                                  muu, muv, phb, ph,            &
 490                                  msftx, msfty, msfux, msfuy,   &
 491                                  msfvx, msfvy,                 &
 492                                  f, e, sina, cosa, fnm, fnp,   &
 493                                  ids, ide, jds, jde, kds, kde, &
 494                                  ims, ime, jms, jme, kms, kme, &
 495                                  its, ite, jts, jte, kts, kte )
 496   ELSE
 497     CALL coriolis ( ru, rv, rw,                   &
 498                     ru_tend,  rv_tend,  rw_tend,  &
 499                     config_flags,                 &
 500                     msftx, msfty, msfux, msfuy,   &
 501                     msfvx, msfvy,                 &
 502                     f, e, sina, cosa, fnm, fnp,   &
 503                     ids, ide, jds, jde, kds, kde, &
 504                     ims, ime, jms, jme, kms, kme, &
 505                     its, ite, jts, jte, kts, kte )
 506
 507   END IF
 508
 509   CALL curvature ( ru, rv, rw, u, v, w,            &
 510                    ru_tend,  rv_tend,  rw_tend,    &
 511                    config_flags,                   &
 512                    msfux, msfuy, msfvx, msfvy,     &
 513                    msftx, msfty,                   &
 514                    xlat, fnm, fnp, rdx, rdy,       &
 515                    ids, ide, jds, jde, kds, kde,   &
 516                    ims, ime, jms, jme, kms, kme,   &
 517                    its, ite, jts, jte, kts, kte   )
 518
 519 ! Damping option added for Held-Suarez test (also uses lw option HELDSUAREZ)
 520   IF (config_flags%ra_lw_physics == HELDSUAREZ) THEN
 521      CALL held_suarez_damp ( ru_tend, rv_tend,               &
 522                              ru,rv,p,pb,                     &
 523                              ids, ide, jds, jde, kds, kde,   &
 524                              ims, ime, jms, jme, kms, kme,   &
 525                              its, ite, jts, jte, kts, kte   )
 526   END IF
 527
 528 !**************************************************************
 529 !
 530 !  Next, the terms that we integrate only with forward-in-time
 531 !  (evaluate with time t variables).
 532 !
 533 !**************************************************************
 534
 535   forward_step: IF( rk_step == 1 ) THEN
 536
 537     diff_opt1 : IF (config_flags%diff_opt .eq. 1) THEN
 538
 539         CALL horizontal_diffusion ('u', u, ru_tendf, mut, config_flags, &
 540                                         msfux, msfuy, msfvx, msfvx_inv, &
 541                                         msfvy,msftx, msfty,             &
 542                                         khdif, xkmhd, rdx, rdy,         &
 543                                         ids, ide, jds, jde, kds, kde,   &
 544                                         ims, ime, jms, jme, kms, kme,   &
 545                                         its, ite, jts, jte, kts, kte   )
 546
 547         CALL horizontal_diffusion ('v', v, rv_tendf, mut, config_flags, &
 548                                         msfux, msfuy, msfvx, msfvx_inv, &
 549                                         msfvy,msftx, msfty,             &
 550                                         khdif, xkmhd, rdx, rdy,         &
 551                                         ids, ide, jds, jde, kds, kde,   &
 552                                         ims, ime, jms, jme, kms, kme,   &
 553                                         its, ite, jts, jte, kts, kte   )
 554
 555         CALL horizontal_diffusion ('w', w, rw_tendf, mut, config_flags, &
 556                                         msfux, msfuy, msfvx, msfvx_inv, &
 557                                         msfvy,msftx, msfty,             &
 558                                         khdif, xkmhd, rdx, rdy,         &
 559                                         ids, ide, jds, jde, kds, kde,   &
 560                                         ims, ime, jms, jme, kms, kme,   &
 561                                         its, ite, jts, jte, kts, kte   )
 562
 563         khdq = 3.*khdif
 564         CALL horizontal_diffusion_3dmp ( 'm', t, t_tendf, mut,            &
 565                                          config_flags, t_init,            &
 566                                          msfux, msfuy, msfvx, msfvx_inv,  &
 567                                          msfvy, msftx, msfty,             &
 568                                          khdq , xkhh, rdx, rdy,           &
 569                                          ids, ide, jds, jde, kds, kde,    &
 570                                          ims, ime, jms, jme, kms, kme,    &
 571                                          its, ite, jts, jte, kts, kte    )
 572
 573         pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
 574
 575           CALL vertical_diffusion_u ( u, ru_tendf, config_flags,      &
 576                                       u_base,                         &
 577                                       alt, muu, rdn, rdnw, kvdif,     &
 578                                       ids, ide, jds, jde, kds, kde,   &
 579                                       ims, ime, jms, jme, kms, kme,   &
 580                                       its, ite, jts, jte, kts, kte   )
 581
 582           CALL vertical_diffusion_v ( v, rv_tendf, config_flags,      &
 583                                       v_base,                         &
 584                                       alt, muv, rdn, rdnw, kvdif,     &
 585                                       ids, ide, jds, jde, kds, kde,   &
 586                                       ims, ime, jms, jme, kms, kme,   &
 587                                       its, ite, jts, jte, kts, kte   )
 588
 589           IF (non_hydrostatic)                                           &
 590           CALL vertical_diffusion ( 'w', w, rw_tendf, config_flags,      &
 591                                     alt, mut, rdn, rdnw, kvdif,          &
 592                                     ids, ide, jds, jde, kds, kde,        &
 593                                     ims, ime, jms, jme, kms, kme,        &
 594                                     its, ite, jts, jte, kts, kte        )
 595
 596           kvdq = 3.*kvdif
 597           CALL vertical_diffusion_3dmp ( t, t_tendf, config_flags, t_init,     &
 598                                          alt, mut, rdn, rdnw, kvdq ,           &
 599                                          ids, ide, jds, jde, kds, kde,         &
 600                                          ims, ime, jms, jme, kms, kme,         &
 601                                          its, ite, jts, jte, kts, kte         )
 602
 603         ENDIF pbl_test
 604
 605    !  Theta tendency computations.
 606
 607     END IF diff_opt1
 608
 609     IF ( diff_6th_opt .NE. 0 ) THEN
 610
 611       CALL sixth_order_diffusion( 'u', u, ru_tendf, mut, dt,          &
 612                                        config_flags,                  &
 613                                        diff_6th_opt, diff_6th_factor, &
 614                                        ids, ide, jds, jde, kds, kde,  &
 615                                        ims, ime, jms, jme, kms, kme,  &
 616                                        its, ite, jts, jte, kts, kte )
 617
 618       CALL sixth_order_diffusion( 'v', v, rv_tendf, mut, dt,          &
 619                                        config_flags,                  &
 620                                        diff_6th_opt, diff_6th_factor, &
 621                                        ids, ide, jds, jde, kds, kde,  &
 622                                        ims, ime, jms, jme, kms, kme,  &
 623                                        its, ite, jts, jte, kts, kte )
 624
 625       IF (non_hydrostatic)                                            &
 626       CALL sixth_order_diffusion( 'w', w, rw_tendf, mut, dt,          &
 627                                        config_flags,                  &
 628                                        diff_6th_opt, diff_6th_factor, &
 629                                        ids, ide, jds, jde, kds, kde,  &
 630                                        ims, ime, jms, jme, kms, kme,  &
 631                                        its, ite, jts, jte, kts, kte )
 632
 633       CALL sixth_order_diffusion( 'm', t,  t_tendf, mut, dt,          &
 634                                        config_flags,                  &
 635                                        diff_6th_opt, diff_6th_factor, &
 636                                        ids, ide, jds, jde, kds, kde,  &
 637                                        ims, ime, jms, jme, kms, kme,  &
 638                                        its, ite, jts, jte, kts, kte )
 639
 640     ENDIF
 641
 642     IF( damp_opt .eq. 2 )                                      &
 643        CALL rk_rayleigh_damp( ru_tendf, rv_tendf,              &
 644                               rw_tendf, t_tendf,               &
 645                               u, v, w, t, t_init,              &
 646                               mut, muu, muv, ph, phb,          &
 647                               u_base, v_base, t_base, z_base,  &
 648                               dampcoef, zdamp,                 &
 649                               ids, ide, jds, jde, kds, kde,    &
 650                               ims, ime, jms, jme, kms, kme,    &
 651                               its, ite, jts, jte, kts, kte   )
 652
 653   END IF forward_step
 654
 655 END SUBROUTINE rk_tendency
 656
 657 !-------------------------------------------------------------------------------
 658
 659 SUBROUTINE rk_addtend_dry ( ru_tend, rv_tend, rw_tend, ph_tend, t_tend,      &
 660                             ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf, &
 661                             u_save, v_save, w_save, ph_save, t_save,         &
 662                             mu_tend, mu_tendf, rk_step,                      &
 663                             h_diabatic, mut, msftx, msfty, msfux, msfuy,     &
 664                             msfvx, msfvx_inv, msfvy,                         &
 665                             ids,ide, jds,jde, kds,kde,                       &
 666                             ims,ime, jms,jme, kms,kme,                       &
 667                             ips,ipe, jps,jpe, kps,kpe,                       &
 668                             its,ite, jts,jte, kts,kte                       )
 669
 670    IMPLICIT NONE
 671
 672    !  Input data.
 673
 674    INTEGER ,               INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 675                                             ims, ime, jms, jme, kms, kme, &
 676                                             ips, ipe, jps, jpe, kps, kpe, &
 677                                             its, ite, jts, jte, kts, kte
 678    INTEGER ,               INTENT(IN   ) :: rk_step
 679
 680    REAL , DIMENSION( ims:ime , kms:kme, jms:jme  ) , INTENT(INOUT) :: ru_tend, &
 681                                                                       rv_tend, &
 682                                                                       rw_tend, &
 683                                                                       ph_tend, &
 684                                                                       t_tend,  &
 685                                                                       ru_tendf, &
 686                                                                       rv_tendf, &
 687                                                                       rw_tendf, &
 688                                                                       ph_tendf, &
 689                                                                       t_tendf
 690
 691    REAL , DIMENSION( ims:ime , jms:jme  ) , INTENT(INOUT) :: mu_tend, &
 692                                                              mu_tendf
 693
 694    REAL , DIMENSION( ims:ime , kms:kme, jms:jme  ) , INTENT(IN   ) ::  u_save,  &
 695                                                                        v_save,  &
 696                                                                        w_save,  &
 697                                                                       ph_save,  &
 698                                                                        t_save,  &
 699                                                                       h_diabatic
 700
 701    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: mut,       &
 702                                                                     msftx,     &
 703                                                                     msfty,     &
 704                                                                     msfux,     &
 705                                                                     msfuy,     &
 706                                                                     msfvx,     &
 707                                                                     msfvx_inv, &
 708                                                                     msfvy
 709
 710
 711 ! Local
 712    INTEGER :: i, j, k
 713
 714 !<DESCRIPTION>
 715 !
 716 ! rk_addtend_dry constructs the full large-timestep tendency terms for
 717 ! momentum (u,v,w), theta and geopotential equations.   This is accomplished
 718 ! by combining the physics tendencies (in *tendf; these are computed
 719 ! the first RK substep, held fixed thereafter) with the RK tendencies
 720 ! (in *tend, these include advection, pressure gradient, etc;
 721 ! these change each rk substep).  Output is in *tend.
 722 !
 723 !</DESCRIPTION>
 724
 725 !  Finally, add the forward-step tendency to the rk_tendency
 726
 727 ! u/v/w/save contain bc tendency that needs to be multiplied by msf
 728 ! (u by msfuy, v by msfvx)
 729 !  before adding it to physics tendency (*tendf)
 730 ! For momentum we need the final tendency to include an inverse msf
 731 ! physics/bc tendency needs to be divided, advection tendency already has it
 732
 733 ! For scalars we need the final tendency to include an inverse msf (msfty)
 734 ! advection tendency is OK, physics/bc tendency needs to be divided by msf
 735
 736    DO j = jts,MIN(jte,jde-1)
 737    DO k = kts,kte-1
 738    DO i = its,ite
 739      ! multiply by my to uncouple u
 740      IF(rk_step == 1)ru_tendf(i,k,j) = ru_tendf(i,k,j) +  u_save(i,k,j)*msfuy(i,j)
 741      ! divide by my to couple u
 742      ru_tend(i,k,j) = ru_tend(i,k,j) + ru_tendf(i,k,j)/msfuy(i,j)
 743    ENDDO
 744    ENDDO
 745    ENDDO
 746
 747    DO j = jts,jte
 748    DO k = kts,kte-1
 749    DO i = its,MIN(ite,ide-1)
 750      ! multiply by mx to uncouple v
 751      IF(rk_step == 1)rv_tendf(i,k,j) = rv_tendf(i,k,j) +  v_save(i,k,j)*msfvx(i,j)
 752      ! divide by mx to couple v
 753      rv_tend(i,k,j) = rv_tend(i,k,j) + rv_tendf(i,k,j)*msfvx_inv(i,j)
 754    ENDDO
 755    ENDDO
 756    ENDDO
 757
 758    DO j = jts,MIN(jte,jde-1)
 759    DO k = kts,kte
 760    DO i = its,MIN(ite,ide-1)
 761      ! multiply by my to uncouple w
 762      IF(rk_step == 1)rw_tendf(i,k,j) = rw_tendf(i,k,j) +  w_save(i,k,j)*msfty(i,j)
 763      ! divide by my to couple w
 764      rw_tend(i,k,j) = rw_tend(i,k,j) + rw_tendf(i,k,j)/msfty(i,j)
 765      IF(rk_step == 1)ph_tendf(i,k,j) = ph_tendf(i,k,j) +  ph_save(i,k,j)
 766      ! divide by my to couple scalar
 767      ph_tend(i,k,j) = ph_tend(i,k,j) + ph_tendf(i,k,j)/msfty(i,j)
 768    ENDDO
 769    ENDDO
 770    ENDDO
 771
 772    DO j = jts,MIN(jte,jde-1)
 773    DO k = kts,kte-1
 774    DO i = its,MIN(ite,ide-1)
 775      IF(rk_step == 1)t_tendf(i,k,j) = t_tendf(i,k,j) +  t_save(i,k,j)
 776      ! divide by my to couple theta
 777       t_tend(i,k,j) =  t_tend(i,k,j) +  t_tendf(i,k,j)/msfty(i,j)  &
 778                                      +  mut(i,j)*h_diabatic(i,k,j)/msfty(i,j)
 779      ! divide by my to couple heating
 780    ENDDO
 781    ENDDO
 782    ENDDO
 783
 784    DO j = jts,MIN(jte,jde-1)
 785    DO i = its,MIN(ite,ide-1)
 786 ! mu tendencies not coupled with 1/msf
 787       mu_tend(i,j) =  mu_tend(i,j) +  mu_tendf(i,j)
 788    ENDDO
 789    ENDDO
 790
 791 END SUBROUTINE rk_addtend_dry
 792
 793 !-------------------------------------------------------------------------------
 794
 795 SUBROUTINE rk_scalar_tend ( scs, sce, config_flags,          &
 796                             rk_step, dt,                     &
 797                             ru, rv, ww, mut, mub, mu_old,    &
 798                             alt,                             &
 799                             scalar_old, scalar,              &
 800                             scalar_tends, advect_tend,       &
 801                             RQVFTEN,                         &
 802                             base, moist_step, fnm, fnp,      &
 803                             msfux, msfuy, msfvx, msfvx_inv,  &
 804                             msfvy, msftx, msfty,             &
 805                             rdx, rdy, rdn, rdnw,             &
 806                             khdif, kvdif, xkmhd,             &
 807                             diff_6th_opt, diff_6th_factor,   &
 808                             pd_advection,                    &
 809                             ids, ide, jds, jde, kds, kde,    &
 810                             ims, ime, jms, jme, kms, kme,    &
 811                             its, ite, jts, jte, kts, kte    )
 812
 813    IMPLICIT NONE
 814
 815    !  Input data.
 816
 817    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
 818
 819    INTEGER ,                INTENT(IN   ) :: rk_step, scs, sce
 820    INTEGER ,                INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 821                                              ims, ime, jms, jme, kms, kme, &
 822                                              its, ite, jts, jte, kts, kte
 823
 824    LOGICAL , INTENT(IN   ) :: moist_step
 825
 826    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ),                &
 827                                          INTENT(IN   )  :: scalar, scalar_old
 828
 829    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ),                      &
 830                                          INTENT(INOUT)  :: scalar_tends
 831
 832    REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(INOUT) :: advect_tend
 833
 834    REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(OUT  ) :: RQVFTEN
 835
 836    REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(IN   ) ::     ru,  &
 837                                                                       rv,  &
 838                                                                       ww,  &
 839                                                                       xkmhd,  &
 840                                                                       alt
 841
 842
 843    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fnm,  &
 844                                                                   fnp,  &
 845                                                                   rdn,  &
 846                                                                   rdnw, &
 847                                                                   base
 848
 849    REAL , DIMENSION( ims:ime , jms:jme ) ,       INTENT(IN   ) :: msfux,    &
 850                                                                   msfuy,    &
 851                                                                   msfvx,    &
 852                                                                   msfvx_inv,    &
 853                                                                   msfvy,    &
 854                                                                   msftx,    &
 855                                                                   msfty,    &
 856                                                                   mub,     &
 857                                                                   mut,     &
 858                                                                   mu_old
 859
 860    REAL ,                                        INTENT(IN   ) :: rdx,     &
 861                                                                   rdy,     &
 862                                                                   khdif,   &
 863                                                                   kvdif
 864
 865    INTEGER, INTENT( IN ) :: diff_6th_opt
 866    REAL,    INTENT( IN ) :: diff_6th_factor
 867
 868    REAL ,                                        INTENT(IN   ) :: dt
 869
 870    LOGICAL, INTENT(IN   ) :: pd_advection
 871
 872    ! Local data
 873
 874    INTEGER :: im, i,j,k
 875    INTEGER :: time_step
 876
 877    REAL    :: khdq, kvdq, tendency
 878
 879 !<DESCRIPTION>
 880 !
 881 ! rk_scalar_tend calls routines that computes scalar tendency from advection
 882 ! and 3D mixing (TKE or fixed eddy viscosities).
 883 !
 884 !</DESCRIPTION>
 885
 886
 887    khdq = khdif/prandtl
 888    kvdq = kvdif/prandtl
 889
 890    scalar_loop : DO im = scs, sce
 891
 892      CALL zero_tend ( advect_tend(ims,kms,jms),     &
 893                       ids, ide, jds, jde, kds, kde, &
 894                       ims, ime, jms, jme, kms, kme, &
 895                       its, ite, jts, jte, kts, kte )
 896
 897      CALL nl_get_time_step ( 1, time_step )
 898
 899       IF( (rk_step == 3) .and. pd_advection ) THEN
 900
 901        CALL advect_scalar_pd       ( scalar(ims,kms,jms,im),             &
 902                                      scalar_old(ims,kms,jms,im),         &
 903                                      advect_tend(ims,kms,jms),           &
 904                                      ru, rv, ww, mut, mub, mu_old,       &
 905                                      config_flags,                       &
 906                                      msfux, msfuy, msfvx, msfvy,         &
 907                                      msftx, msfty, fnm, fnp,             &
 908                                      rdx, rdy, rdnw,dt,                  &
 909                                      ids, ide, jds, jde, kds, kde,       &
 910                                      ims, ime, jms, jme, kms, kme,       &
 911                                      its, ite, jts, jte, kts, kte     )
 912
 913       ELSE
 914
 915        CALL advect_scalar     ( scalar(ims,kms,jms,im),        &
 916                                 scalar(ims,kms,jms,im),        &
 917                                 advect_tend(ims,kms,jms),      &
 918                                 ru, rv, ww, mut, time_step,    &
 919                                 config_flags,                  &
 920                                 msfux, msfuy, msfvx, msfvy,    &
 921                                 msftx, msfty, fnm, fnp,        &
 922                                 rdx, rdy, rdnw,                &
 923                                 ids, ide, jds, jde, kds, kde,  &
 924                                 ims, ime, jms, jme, kms, kme,  &
 925                                 its, ite, jts, jte, kts, kte  )
 926       END IF
 927
 928      IF((config_flags%cu_physics == GDSCHEME .OR. config_flags%cu_physics == G3SCHEME) &
 929                      .and. moist_step .and. ( im == P_QV) ) THEN
 930
 931         CALL set_tend( RQVFTEN, advect_tend, msfty,    &
 932                        ids, ide, jds, jde, kds, kde,   &
 933                        ims, ime, jms, jme, kms, kme,   &
 934                        its, ite, jts, jte, kts, kte      )
 935      ENDIF
 936
 937      rk_step_1: IF( rk_step == 1 ) THEN
 938
 939        diff_opt1 : IF (config_flags%diff_opt .eq. 1) THEN
 940
 941        CALL horizontal_diffusion ( 'm', scalar(ims,kms,jms,im),            &
 942                                         scalar_tends(ims,kms,jms,im), mut, &
 943                                         config_flags,                      &
 944                                         msfux, msfuy, msfvx, msfvx_inv,    &
 945                                         msfvy, msftx, msfty,               &
 946                                         khdq , xkmhd, rdx, rdy,            &
 947                                         ids, ide, jds, jde, kds, kde,      &
 948                                         ims, ime, jms, jme, kms, kme,      &
 949                                         its, ite, jts, jte, kts, kte      )
 950
 951        pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
 952
 953          IF( (moist_step) .and. ( im == P_QV)) THEN
 954
 955             CALL vertical_diffusion_mp ( scalar(ims,kms,jms,im),       &
 956                                          scalar_tends(ims,kms,jms,im), &
 957                                          config_flags, base,           &
 958                                          alt, mut, rdn, rdnw, kvdq ,   &
 959                                          ids, ide, jds, jde, kds, kde, &
 960                                          ims, ime, jms, jme, kms, kme, &
 961                                          its, ite, jts, jte, kts, kte )
 962
 963          ELSE
 964
 965             CALL vertical_diffusion (  'm', scalar(ims,kms,jms,im),       &
 966                                             scalar_tends(ims,kms,jms,im), &
 967                                             config_flags,                 &
 968                                             alt, mut, rdn, rdnw, kvdq,    &
 969                                             ids, ide, jds, jde, kds, kde, &
 970                                             ims, ime, jms, jme, kms, kme, &
 971                                             its, ite, jts, jte, kts, kte )
 972
 973          END IF
 974
 975       ENDIF pbl_test
 976
 977     ENDIF diff_opt1
 978
 979     IF ( diff_6th_opt .NE. 0 )                                        &
 980       CALL sixth_order_diffusion( 'm', scalar(ims,kms,jms,im),        &
 981                                        scalar_tends(ims,kms,jms,im),  &
 982                                        mut, dt, config_flags,         &
 983                                        diff_6th_opt, diff_6th_factor, &
 984                                        ids, ide, jds, jde, kds, kde,  &
 985                                        ims, ime, jms, jme, kms, kme,  &
 986                                        its, ite, jts, jte, kts, kte )
 987
 988   ENDIF rk_step_1
 989
 990  END DO scalar_loop
 991
 992 END SUBROUTINE rk_scalar_tend
 993
 994 !-------------------------------------------------------------------------------
 995
 996 SUBROUTINE rk_update_scalar( scs, sce,                      &
 997                              scalar_1, scalar_2, sc_tend,   &
 998                              advect_tend, msftx, msfty,     &
 999                              mu_old, mu_new, mu_base,       &
1000                              rk_step, dt, spec_zone,        &
1001                              config_flags,                  &
1002                              ids, ide, jds, jde, kds, kde,  &
1003                              ims, ime, jms, jme, kms, kme,  &
1004                              its, ite, jts, jte, kts, kte  )
1005
1006    IMPLICIT NONE
1007
1008    !  Input data.
1009
1010    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
1011
1012    INTEGER ,                INTENT(IN   ) :: scs, sce, rk_step, spec_zone
1013    INTEGER ,                INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1014                                              ims, ime, jms, jme, kms, kme, &
1015                                              its, ite, jts, jte, kts, kte
1016
1017    REAL,                    INTENT(IN   ) :: dt
1018
1019    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
1020          INTENT(INOUT)                                  :: scalar_1,    &
1021                                                            scalar_2,  &
1022                                                            sc_tend
1023
1024    REAL, DIMENSION(ims:ime, kms:kme, jms:jme ),                &
1025          INTENT(IN)                                  :: advect_tend
1026
1027    REAL, DIMENSION(ims:ime, jms:jme  ), INTENT(IN   ) ::  mu_old,  &
1028                                                           mu_new,  &
1029                                                           mu_base, &
1030                                                           msftx,   &
1031                                                           msfty
1032
1033    INTEGER :: i,j,k,im
1034    REAL    :: sc_middle, msfsq
1035    REAL, DIMENSION(its:ite) :: muold, r_munew
1036
1037    REAL, DIMENSION(its:ite, kts:kte, jts:jte  ) :: tendency
1038
1039    INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
1040    INTEGER :: i_start_spc,i_end_spc,j_start_spc,j_end_spc,k_start_spc,k_end_spc
1041
1042 !<DESCRIPTION>
1043 !
1044 !  rk_scalar_update advances the scalar equation given the time t value
1045 !  of the scalar and the scalar tendency.
1046 !
1047 !</DESCRIPTION>
1048
1049
1050 !
1051 !  set loop limits.
1052
1053       i_start = its
1054       i_end   = ite
1055       j_start = jts
1056       j_end   = jte
1057       k_start = kts
1058       k_end   = kte-1
1059       IF(j_end == jde) j_end = j_end - 1
1060       IF(i_end == ide) i_end = i_end - 1
1061
1062       i_start_spc = i_start
1063       i_end_spc   = i_end
1064       j_start_spc = j_start
1065       j_end_spc   = j_end
1066       k_start_spc = k_start
1067       k_end_spc   = k_end
1068
1069     IF( config_flags%nested .or. config_flags%specified ) THEN
1070       IF( .NOT. config_flags%periodic_x)i_start = max( its,ids+spec_zone )
1071       IF( .NOT. config_flags%periodic_x)i_end   = min( ite,ide-spec_zone-1 )
1072       j_start = max( jts,jds+spec_zone )
1073       j_end   = min( jte,jde-spec_zone-1 )
1074       k_start = kts
1075       k_end   = min( kte, kde-1 )
1076     ENDIF
1077
1078     IF ( rk_step == 1 ) THEN
1079
1080       !  replace t-dt values (in scalar_1) with t values scalar_2,
1081       !  then compute new values by adding tendency to values at t
1082
1083       DO  im = scs,sce
1084
1085        DO  j = jts, min(jte,jde-1)
1086        DO  k = kts, min(kte,kde-1)
1087        DO  i = its, min(ite,ide-1)
1088            tendency(i,k,j) = 0.
1089        ENDDO
1090        ENDDO
1091        ENDDO
1092
1093        DO  j = j_start,j_end
1094        DO  k = k_start,k_end
1095        DO  i = i_start,i_end
1096           ! scalar was coupled with my
1097            tendency(i,k,j) = advect_tend(i,k,j) * msfty(i,j)
1098        ENDDO
1099        ENDDO
1100        ENDDO
1101
1102        DO  j = j_start_spc,j_end_spc
1103        DO  k = k_start_spc,k_end_spc
1104        DO  i = i_start_spc,i_end_spc
1105            tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1106        ENDDO
1107        ENDDO
1108        ENDDO
1109
1110       DO  j = jts, min(jte,jde-1)
1111
1112       DO  i = its, min(ite,ide-1)
1113         muold(i) = mu_old(i,j) + mu_base(i,j)
1114         r_munew(i) = 1./(mu_new(i,j) + mu_base(i,j))
1115       ENDDO
1116
1117       DO  k = kts, min(kte,kde-1)
1118       DO  i = its, min(ite,ide-1)
1119
1120         scalar_1(i,k,j,im) = scalar_2(i,k,j,im)
1121         scalar_2(i,k,j,im) = (muold(i)*scalar_1(i,k,j,im)   &
1122                              + dt*tendency(i,k,j))*r_munew(i)
1123
1124       ENDDO
1125       ENDDO
1126       ENDDO
1127
1128       ENDDO
1129
1130     ELSE
1131
1132       !  just compute new values, scalar_1 already at time t.
1133
1134       DO  im = scs, sce
1135
1136        DO  j = jts, min(jte,jde-1)
1137        DO  k = kts, min(kte,kde-1)
1138        DO  i = its, min(ite,ide-1)
1139            tendency(i,k,j) = 0.
1140        ENDDO
1141        ENDDO
1142        ENDDO
1143
1144        DO  j = j_start,j_end
1145        DO  k = k_start,k_end
1146        DO  i = i_start,i_end
1147            ! scalar was coupled with my
1148            tendency(i,k,j) = advect_tend(i,k,j) * msfty(i,j)
1149        ENDDO
1150        ENDDO
1151        ENDDO
1152
1153        DO  j = j_start_spc,j_end_spc
1154        DO  k = k_start_spc,k_end_spc
1155        DO  i = i_start_spc,i_end_spc
1156            tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1157        ENDDO
1158        ENDDO
1159        ENDDO
1160
1161       DO  j = jts, min(jte,jde-1)
1162
1163       DO  i = its, min(ite,ide-1)
1164         muold(i) = mu_old(i,j) + mu_base(i,j)
1165         r_munew(i) = 1./(mu_new(i,j) + mu_base(i,j))
1166       ENDDO
1167
1168       DO  k = kts, min(kte,kde-1)
1169       DO  i = its, min(ite,ide-1)
1170
1171         scalar_2(i,k,j,im) = (muold(i)*scalar_1(i,k,j,im)   &
1172                              + dt*tendency(i,k,j))*r_munew(i)
1173
1174       ENDDO
1175       ENDDO
1176       ENDDO
1177
1178       ENDDO
1179
1180     END IF
1181
1182 END SUBROUTINE rk_update_scalar
1183
1184 !-------------------------------------------------------------------------------
1185
1186 SUBROUTINE rk_update_scalar_pd( scs, sce,                      &
1187                                 scalar, sc_tend,               &
1188                                 mu_old, mu_new, mu_base,       &
1189                                 rk_step, dt, spec_zone,        &
1190                                 config_flags,                  &
1191                                 ids, ide, jds, jde, kds, kde,  &
1192                                 ims, ime, jms, jme, kms, kme,  &
1193                                 its, ite, jts, jte, kts, kte  )
1194
1195    IMPLICIT NONE
1196
1197    !  Input data.
1198
1199    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
1200
1201    INTEGER ,                INTENT(IN   ) :: scs, sce, rk_step, spec_zone
1202    INTEGER ,                INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1203                                              ims, ime, jms, jme, kms, kme, &
1204                                              its, ite, jts, jte, kts, kte
1205
1206    REAL,                    INTENT(IN   ) :: dt
1207
1208    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
1209          INTENT(INOUT)                                  :: scalar,      &
1210                                                            sc_tend
1211
1212    REAL, DIMENSION(ims:ime, jms:jme  ), INTENT(IN   ) ::  mu_old,  &
1213                                                           mu_new,  &
1214                                                           mu_base
1215
1216    INTEGER :: i,j,k,im
1217    REAL    :: sc_middle, msfsq
1218    REAL, DIMENSION(its:ite) :: muold, r_munew
1219
1220    REAL, DIMENSION(its:ite, kts:kte, jts:jte  ) :: tendency
1221
1222    INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
1223    INTEGER :: i_start_spc,i_end_spc,j_start_spc,j_end_spc,k_start_spc,k_end_spc
1224
1225 !<DESCRIPTION>
1226 !
1227 !  rk_scalar_update advances the scalar equation given the time t value
1228 !  of the scalar and the scalar tendency.
1229 !
1230 !</DESCRIPTION>
1231
1232
1233 !
1234 !  set loop limits.
1235
1236       i_start = its
1237       i_end   = ite
1238       j_start = jts
1239       j_end   = jte
1240       k_start = kts
1241       k_end   = kte-1
1242       IF(j_end == jde) j_end = j_end - 1
1243       IF(i_end == ide) i_end = i_end - 1
1244
1245       i_start_spc = i_start
1246       i_end_spc   = i_end
1247       j_start_spc = j_start
1248       j_end_spc   = j_end
1249       k_start_spc = k_start
1250       k_end_spc   = k_end
1251
1252     IF( config_flags%nested .or. config_flags%specified ) THEN
1253       IF( .NOT. config_flags%periodic_x)i_start = max( its,ids+spec_zone )
1254       IF( .NOT. config_flags%periodic_x)i_end   = min( ite,ide-spec_zone-1 )
1255       j_start = max( jts,jds+spec_zone )
1256       j_end   = min( jte,jde-spec_zone-1 )
1257       k_start = kts
1258       k_end   = min( kte, kde-1 )
1259     ENDIF
1260
1261       DO  im = scs, sce
1262
1263        DO  j = jts, min(jte,jde-1)
1264        DO  k = kts, min(kte,kde-1)
1265        DO  i = its, min(ite,ide-1)
1266            tendency(i,k,j) = 0.
1267        ENDDO
1268        ENDDO
1269        ENDDO
1270
1271        DO  j = j_start_spc,j_end_spc
1272        DO  k = k_start_spc,k_end_spc
1273        DO  i = i_start_spc,i_end_spc
1274            tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1275            sc_tend(i,k,j,im) = 0.
1276        ENDDO
1277        ENDDO
1278        ENDDO
1279
1280       DO  j = jts, min(jte,jde-1)
1281
1282       DO  i = its, min(ite,ide-1)
1283         muold(i) = mu_old(i,j) + mu_base(i,j)
1284         r_munew(i) = 1./(mu_new(i,j) + mu_base(i,j))
1285       ENDDO
1286
1287       DO  k = kts, min(kte,kde-1)
1288       DO  i = its, min(ite,ide-1)
1289
1290         scalar(i,k,j,im) = (muold(i)*scalar(i,k,j,im)   &
1291                              + dt*tendency(i,k,j))*r_munew(i)
1292       ENDDO
1293       ENDDO
1294       ENDDO
1295
1296       ENDDO
1297
1298 END SUBROUTINE rk_update_scalar_pd
1299
1300 !------------------------------------------------------------
1301
1302 SUBROUTINE init_zero_tendency(ru_tendf, rv_tendf, rw_tendf, ph_tendf,  &
1303                               t_tendf,  tke_tendf, mu_tendf,           &
1304                               moist_tendf,chem_tendf,scalar_tendf,     &
1305                               n_moist,n_chem,n_scalar,rk_step,         &
1306                               ids, ide, jds, jde, kds, kde,            &
1307                               ims, ime, jms, jme, kms, kme,            &
1308                               its, ite, jts, jte, kts, kte             )
1309 !-----------------------------------------------------------------------
1310    IMPLICIT NONE
1311 !-----------------------------------------------------------------------
1312
1313    INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1314                                     ims, ime, jms, jme, kms, kme, &
1315                                     its, ite, jts, jte, kts, kte
1316
1317    INTEGER ,       INTENT(IN   ) :: n_moist,n_chem,n_scalar,rk_step
1318
1319    REAL , DIMENSION( ims:ime , kms:kme, jms:jme  ) , INTENT(INOUT) ::  &
1320                                                              ru_tendf, &
1321                                                              rv_tendf, &
1322                                                              rw_tendf, &
1323                                                              ph_tendf, &
1324                                                               t_tendf, &
1325                                                             tke_tendf
1326
1327    REAL , DIMENSION( ims:ime , jms:jme  ) , INTENT(INOUT) ::  mu_tendf
1328
1329    REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_moist),INTENT(INOUT)::&
1330                                                           moist_tendf
1331
1332    REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_chem ),INTENT(INOUT)::&
1333                                                           chem_tendf
1334
1335    REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_scalar ),INTENT(INOUT)::&
1336                                                           scalar_tendf
1337
1338 ! LOCAL VARS
1339
1340    INTEGER :: im, ic, is
1341
1342 !<DESCRIPTION>
1343 !
1344 ! init_zero_tendency
1345 ! sets tendency arrays to zero for all prognostic variables.
1346 !
1347 !</DESCRIPTION>
1348
1349
1350    CALL zero_tend ( ru_tendf,                        &
1351                     ids, ide, jds, jde, kds, kde,    &
1352                     ims, ime, jms, jme, kms, kme,    &
1353                     its, ite, jts, jte, kts, kte     )
1354
1355    CALL zero_tend ( rv_tendf,                        &
1356                     ids, ide, jds, jde, kds, kde,    &
1357                     ims, ime, jms, jme, kms, kme,    &
1358                     its, ite, jts, jte, kts, kte     )
1359
1360    CALL zero_tend ( rw_tendf,                        &
1361                     ids, ide, jds, jde, kds, kde,    &
1362                     ims, ime, jms, jme, kms, kme,    &
1363                     its, ite, jts, jte, kts, kte     )
1364
1365    CALL zero_tend ( ph_tendf,                        &
1366                     ids, ide, jds, jde, kds, kde,    &
1367                     ims, ime, jms, jme, kms, kme,    &
1368                     its, ite, jts, jte, kts, kte     )
1369
1370    CALL zero_tend ( t_tendf,                         &
1371                     ids, ide, jds, jde, kds, kde,    &
1372                     ims, ime, jms, jme, kms, kme,    &
1373                     its, ite, jts, jte, kts, kte     )
1374
1375    CALL zero_tend ( tke_tendf,                       &
1376                     ids, ide, jds, jde, kds, kde,    &
1377                     ims, ime, jms, jme, kms, kme,    &
1378                     its, ite, jts, jte, kts, kte     )
1379
1380    CALL zero_tend ( mu_tendf,                        &
1381                     ids, ide, jds, jde, kds, kds,    &
1382                     ims, ime, jms, jme, kms, kms,    &
1383                     its, ite, jts, jte, kts, kts     )
1384
1385 !   DO im=PARAM_FIRST_SCALAR,n_moist
1386    DO im=1,n_moist                      ! make sure first one is zero too
1387       CALL zero_tend ( moist_tendf(ims,kms,jms,im),  &
1388                        ids, ide, jds, jde, kds, kde, &
1389                        ims, ime, jms, jme, kms, kme, &
1390                        its, ite, jts, jte, kts, kte  )
1391    ENDDO
1392
1393 !   DO ic=PARAM_FIRST_SCALAR,n_chem
1394    DO ic=1,n_chem                       ! make sure first one is zero too
1395       CALL zero_tend ( chem_tendf(ims,kms,jms,ic),   &
1396                        ids, ide, jds, jde, kds, kde, &
1397                        ims, ime, jms, jme, kms, kme, &
1398                        its, ite, jts, jte, kts, kte  )
1399    ENDDO
1400
1401 !   DO ic=PARAM_FIRST_SCALAR,n_scalar
1402    DO ic=1,n_scalar                       ! make sure first one is zero too
1403       CALL zero_tend ( scalar_tendf(ims,kms,jms,ic),   &
1404                        ids, ide, jds, jde, kds, kde, &
1405                        ims, ime, jms, jme, kms, kme, &
1406                        its, ite, jts, jte, kts, kte  )
1407    ENDDO
1408
1409 END SUBROUTINE init_zero_tendency
1410
1411 !===================================================================
1412
1413
1414 SUBROUTINE dump_data( a, field, io_unit,            &
1415                       ims, ime, jms, jme, kms, kme, &
1416                       ids, ide, jds, jde, kds, kde )
1417 implicit none
1418 integer ::  ims, ime, jms, jme, kms, kme, &
1419             ids, ide, jds, jde, kds, kde
1420 real, dimension(ims:ime, kms:kme, jds:jde) :: a
1421 character :: field
1422 integer :: io_unit
1423
1424 integer :: is,ie,js,je,ks,ke
1425
1426 !<DESCRIPTION
1427 !
1428 ! quick and dirty debug io utility
1429 !
1430 !</DESCRIPTION
1431
1432 is = ids
1433 ie = ide-1
1434 js = jds
1435 je = jde-1
1436 ks = kds
1437 ke = kde-1
1438
1439 if(field == 'u') ie = ide
1440 if(field == 'v') je = jde
1441 if(field == 'w') ke = kde
1442
1443 write(io_unit) is,ie,ks,ke,js,je
1444 write(io_unit) a(is:ie, ks:ke, js:je)
1445
1446 end subroutine dump_data
1447
1448 !-----------------------------------------------------------------------
1449
1450 SUBROUTINE calculate_phy_tend (config_flags,mu,muu,muv,pi3d,           &
1451                      RTHRATEN,                                         &
1452                      RUBLTEN,RVBLTEN,RTHBLTEN,                         &
1453                      RQVBLTEN,RQCBLTEN,RQIBLTEN,                       &
1454                      RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN,              &
1455                      RQICUTEN,RQSCUTEN,                                &
1456                      RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,        &
1457                      RMUNDGDTEN,                                       &
1458                      ids,ide, jds,jde, kds,kde,                        &
1459                      ims,ime, jms,jme, kms,kme,                        &
1460                      its,ite, jts,jte, kts,kte                         )
1461 !-----------------------------------------------------------------------
1462       IMPLICIT NONE
1463
1464       TYPE(grid_config_rec_type), INTENT(IN)     ::      config_flags
1465
1466       INTEGER,  INTENT(IN   )   ::          ids,ide, jds,jde, kds,kde, &
1467                                             ims,ime, jms,jme, kms,kme, &
1468                                             its,ite, jts,jte, kts,kte
1469
1470       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )               , &
1471                 INTENT(IN   )   ::                               pi3d
1472
1473       REAL,     DIMENSION( ims:ime, jms:jme )                        , &
1474                 INTENT(IN   )   ::                                 mu, &
1475                                                                   muu, &
1476                                                                   muv
1477
1478
1479 ! radiation
1480
1481       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ),                &
1482                 INTENT(INOUT)   ::                           RTHRATEN
1483
1484 ! cumulus
1485
1486       REAL,     DIMENSION( ims:ime , kms:kme , jms:jme ),              &
1487                 INTENT(INOUT)   ::                                     &
1488                                                              RTHCUTEN, &
1489                                                              RQVCUTEN, &
1490                                                              RQCCUTEN, &
1491                                                              RQRCUTEN, &
1492                                                              RQICUTEN, &
1493                                                              RQSCUTEN
1494 ! pbl
1495
1496       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )               , &
1497                 INTENT(INOUT)   ::                            RUBLTEN, &
1498                                                               RVBLTEN, &
1499                                                              RTHBLTEN, &
1500                                                              RQVBLTEN, &
1501                                                              RQCBLTEN, &
1502                                                              RQIBLTEN
1503
1504 ! fdda
1505
1506       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )               , &
1507                 INTENT(INOUT)   ::                            RUNDGDTEN, &
1508                                                               RVNDGDTEN, &
1509                                                              RTHNDGDTEN, &
1510                                                              RQVNDGDTEN
1511       REAL,     DIMENSION( ims:ime, jms:jme )               , &
1512                 INTENT(INOUT)   ::                           RMUNDGDTEN
1513
1514       INTEGER :: i,k,j
1515       INTEGER :: itf,ktf,jtf,itsu,jtsv
1516
1517 !-----------------------------------------------------------------------
1518
1519 !<DESCRIPTION>
1520 !
1521 !  calculate_phy_tend couples the physics tendencies to the column mass (mu),
1522 !  because prognostic equations are in flux form, but physics tendencies are
1523 !  computed for uncoupled variables.
1524 !
1525 !</DESCRIPTION>
1526
1527       itf=MIN(ite,ide-1)
1528       jtf=MIN(jte,jde-1)
1529       ktf=MIN(kte,kde-1)
1530       itsu=MAX(its,ids+1)
1531       jtsv=MAX(jts,jds+1)
1532
1533 ! radiation
1534
1535    IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
1536
1537       DO J=jts,jtf
1538       DO K=kts,ktf
1539       DO I=its,itf
1540          RTHRATEN(I,K,J)=mu(I,J)*RTHRATEN(I,K,J)
1541       ENDDO
1542       ENDDO
1543       ENDDO
1544
1545    ENDIF
1546
1547 ! cumulus
1548
1549    IF (config_flags%cu_physics .gt. 0) THEN
1550
1551       DO J=jts,jtf
1552       DO I=its,itf
1553       DO K=kts,ktf
1554          RTHCUTEN(I,K,J)=mu(I,J)*RTHCUTEN(I,K,J)
1555          RQVCUTEN(I,K,J)=mu(I,J)*RQVCUTEN(I,K,J)
1556       ENDDO
1557       ENDDO
1558       ENDDO
1559
1560       IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
1561          DO J=jts,jtf
1562          DO I=its,itf
1563          DO K=kts,ktf
1564             RQCCUTEN(I,K,J)=mu(I,J)*RQCCUTEN(I,K,J)
1565          ENDDO
1566          ENDDO
1567          ENDDO
1568       ENDIF
1569
1570       IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
1571          DO J=jts,jtf
1572          DO I=its,itf
1573          DO K=kts,ktf
1574             RQRCUTEN(I,K,J)=mu(I,J)*RQRCUTEN(I,K,J)
1575          ENDDO
1576          ENDDO
1577          ENDDO
1578       ENDIF
1579
1580       IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
1581          DO J=jts,jtf
1582          DO I=its,itf
1583          DO K=kts,ktf
1584             RQICUTEN(I,K,J)=mu(I,J)*RQICUTEN(I,K,J)
1585          ENDDO
1586          ENDDO
1587          ENDDO
1588       ENDIF
1589
1590       IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
1591          DO J=jts,jtf
1592          DO I=its,itf
1593          DO K=kts,ktf
1594             RQSCUTEN(I,K,J)=mu(I,J)*RQSCUTEN(I,K,J)
1595          ENDDO
1596          ENDDO
1597          ENDDO
1598       ENDIF
1599
1600    ENDIF
1601
1602 ! pbl
1603
1604    IF (config_flags%bl_pbl_physics .gt. 0) THEN
1605
1606       DO J=jts,jtf
1607       DO K=kts,ktf
1608       DO I=its,itf
1609          RUBLTEN(I,K,J) =mu(I,J)*RUBLTEN(I,K,J)
1610          RVBLTEN(I,K,J) =mu(I,J)*RVBLTEN(I,K,J)
1611          RTHBLTEN(I,K,J)=mu(I,J)*RTHBLTEN(I,K,J)
1612       ENDDO
1613       ENDDO
1614       ENDDO
1615
1616       IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
1617          DO J=jts,jtf
1618          DO K=kts,ktf
1619          DO I=its,itf
1620             RQVBLTEN(I,K,J)=mu(I,J)*RQVBLTEN(I,K,J)
1621          ENDDO
1622          ENDDO
1623          ENDDO
1624       ENDIF
1625
1626       IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
1627          DO J=jts,jtf
1628          DO K=kts,ktf
1629          DO I=its,itf
1630            RQCBLTEN(I,K,J)=mu(I,J)*RQCBLTEN(I,K,J)
1631          ENDDO
1632          ENDDO
1633          ENDDO
1634       ENDIF
1635
1636       IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
1637          DO J=jts,jtf
1638          DO K=kts,ktf
1639          DO I=its,itf
1640             RQIBLTEN(I,K,J)=mu(I,J)*RQIBLTEN(I,K,J)
1641          ENDDO
1642          ENDDO
1643          ENDDO
1644       ENDIF
1645
1646     ENDIF
1647
1648 ! fdda
1649 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
1650 !   so only couple those
1651
1652    IF (config_flags%grid_fdda .gt. 0) THEN
1653
1654       DO J=jts,jtf
1655       DO K=kts,ktf
1656       DO I=itsu,itf
1657 !     if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1658 !     write(*,'(a,3i6,e15.5)') 'u_ten before=',i,k,j, RUNDGDTEN(i,k,j)
1659          RUNDGDTEN(I,K,J) =muu(I,J)*RUNDGDTEN(I,K,J)
1660 !        if( i == itf/2 .AND. j == jtf/2 .AND. k==ktf/2 ) &
1661 !          write(*,'(a,2f15.5)') 'mu, muu=',mu(i,j), muu(i,j)
1662 !     if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1663 !     write(*,'(a,3i6,e15.5)') 'u_ten after=',i,k,j, RUNDGDTEN(i,k,j)
1664 !     if( RUNDGDTEN(i,k,j) > 30.0 ) write(*,*) 'IKJ=',i,k,j
1665       ENDDO
1666       ENDDO
1667       ENDDO
1668 !     write(*,'(a,e15.5)') 'u_ten MAXIMUM after=', maxval(RUNDGDTEN)
1669       DO J=jtsv,jtf
1670       DO K=kts,ktf
1671       DO I=its,itf
1672          RVNDGDTEN(I,K,J) =muv(I,J)*RVNDGDTEN(I,K,J)
1673       ENDDO
1674       ENDDO
1675       ENDDO
1676       DO J=jts,jtf
1677       DO K=kts,ktf
1678       DO I=its,itf
1679 !     if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1680 !     write(*,'(a,3i6,e15.5)') 'th before=',i,k,j, RTHNDGDTEN(I,K,J)
1681          RTHNDGDTEN(I,K,J)=mu(I,J)*RTHNDGDTEN(I,K,J)
1682 !        RMUNDGDTEN(I,J) - no coupling
1683 !     if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1684 !     write(*,'(a,3i6,e15.5)') 'th after=',i,k,j, RTHNDGDTEN(I,K,J)
1685       ENDDO
1686       ENDDO
1687       ENDDO
1688
1689       IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
1690          DO J=jts,jtf
1691          DO K=kts,ktf
1692          DO I=its,itf
1693             RQVNDGDTEN(I,K,J)=mu(I,J)*RQVNDGDTEN(I,K,J)
1694          ENDDO
1695          ENDDO
1696          ENDDO
1697       ENDIF
1698
1699     ENDIF
1700
1701 END SUBROUTINE calculate_phy_tend
1702
1703 !-----------------------------------------------------------------------
1704
1705 SUBROUTINE positive_definite_filter ( a,                          &
1706                                       ids,ide, jds,jde, kds,kde,  &
1707                                       ims,ime, jms,jme, kms,kme,  &
1708                                       its,ite, jts,jte, kts,kte  )
1709
1710   IMPLICIT NONE
1711
1712   INTEGER,  INTENT(IN   )   ::          ids,ide, jds,jde, kds,kde, &
1713                                         ims,ime, jms,jme, kms,kme, &
1714                                         its,ite, jts,jte, kts,kte
1715
1716   REAL, DIMENSION( ims:ime , kms:kme , jms:jme  ), INTENT(INOUT) :: a
1717
1718   INTEGER :: i,k,j
1719
1720 !<DESCRIPTION>
1721 !
1722 ! debug and testing code for bounding a variable
1723 !
1724 !</DESCRIPTION>
1725
1726   DO j=jts,min(jte,jde-1)
1727   DO k=kts,kte-1
1728   DO i=its,min(ite,ide-1)
1729 !    a(i,k,j) = max(a(i,k,j),0.)
1730     a(i,k,j) = min(1000.,max(a(i,k,j),0.))
1731   ENDDO
1732   ENDDO
1733   ENDDO
1734
1735   END SUBROUTINE positive_definite_filter
1736
1737 !-----------------------------------------------------------------------
1738
1739 SUBROUTINE bound_tke ( tke, tke_upper_bound,       &
1740                        ids,ide, jds,jde, kds,kde,  &
1741                        ims,ime, jms,jme, kms,kme,  &
1742                        its,ite, jts,jte, kts,kte  )
1743
1744   IMPLICIT NONE
1745
1746   INTEGER,  INTENT(IN   )   ::          ids,ide, jds,jde, kds,kde, &
1747                                         ims,ime, jms,jme, kms,kme, &
1748                                         its,ite, jts,jte, kts,kte
1749
1750   REAL, DIMENSION( ims:ime , kms:kme , jms:jme  ), INTENT(INOUT) :: tke
1751   REAL, INTENT(   IN) :: tke_upper_bound
1752
1753   INTEGER :: i,k,j
1754
1755 !<DESCRIPTION>
1756 !
1757 ! bounds tke between zero and tke_upper_bound.
1758 !
1759 !</DESCRIPTION>
1760
1761   DO j=jts,min(jte,jde-1)
1762   DO k=kts,kte-1
1763   DO i=its,min(ite,ide-1)
1764     tke(i,k,j) = min(tke_upper_bound,max(tke(i,k,j),0.))
1765   ENDDO
1766   ENDDO
1767   ENDDO
1768
1769   END SUBROUTINE bound_tke
1770
1771
1772
1773 END MODULE module_em