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) THEN
 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      ENDIF
 474
 475      CALL w_damp   ( rw_tend, max_vert_cfl,            &
 476                       max_horiz_cfl,                  &
 477                       u, v, ww, w, mut, rdnw,         &
 478                       rdx, rdy, msfux, msfuy, msfvx,  &
 479                       msfvy, dt, config_flags,        &
 480                       ids, ide, jds, jde, kds, kde,   &
 481                       ims, ime, jms, jme, kms, kme,   &
 482                       its, ite, jts, jte, kts, kte   )
 483
 484      IF(config_flags%pert_coriolis) THEN
 485
 486           CALL perturbation_coriolis ( ru, rv, rw,                   &
 487                                        ru_tend,  rv_tend,  rw_tend,  &
 488                                        config_flags,                 &
 489                                        u_base, v_base, z_base,       &
 490                                        muu, muv, phb, ph,            &
 491                                        msftx, msfty, msfux, msfuy,   &
 492                                        msfvx, msfvy,                 &
 493                                        f, e, sina, cosa, fnm, fnp,   &
 494                                        ids, ide, jds, jde, kds, kde, &
 495                                        ims, ime, jms, jme, kms, kme, &
 496                                        its, ite, jts, jte, kts, kte )
 497      ELSE
 498           CALL coriolis ( ru, rv, rw,                   &
 499                           ru_tend,  rv_tend,  rw_tend,  &
 500                           config_flags,                 &
 501                           msftx, msfty, msfux, msfuy,   &
 502                           msfvx, msfvy,                 &
 503                           f, e, sina, cosa, fnm, fnp,   &
 504                           ids, ide, jds, jde, kds, kde, &
 505                           ims, ime, jms, jme, kms, kme, &
 506                           its, ite, jts, jte, kts, kte )
 507
 508      END IF
 509
 510      CALL curvature ( ru, rv, rw, u, v, w,            &
 511                        ru_tend,  rv_tend,  rw_tend,    &
 512                        config_flags,                   &
 513                        msfux, msfuy, msfvx, msfvy,     &
 514                        msftx, msfty,                   &
 515                        xlat, fnm, fnp, rdx, rdy,       &
 516                        ids, ide, jds, jde, kds, kde,   &
 517                        ims, ime, jms, jme, kms, kme,   &
 518                        its, ite, jts, jte, kts, kte   )
 519
 520 ! Damping option added for Held-Suarez test (also uses lw option HELDSUAREZ)
 521       IF (config_flags%ra_lw_physics == HELDSUAREZ) THEN
 522          CALL held_suarez_damp ( ru_tend, rv_tend,               &
 523                                  ru,rv,p,pb,                     &
 524                                  ids, ide, jds, jde, kds, kde,   &
 525                                  ims, ime, jms, jme, kms, kme,   &
 526                                  its, ite, jts, jte, kts, kte   )
 527       END IF
 528
 529 !**************************************************************
 530 !
 531 !  Next, the terms that we integrate only with forward-in-time
 532 !  (evaluate with time t variables).
 533 !
 534 !**************************************************************
 535
 536   forward_step: IF( rk_step == 1 ) THEN
 537
 538     diff_opt1 : IF (config_flags%diff_opt .eq. 1) THEN
 539
 540         CALL horizontal_diffusion ('u', u, ru_tendf, mut, config_flags, &
 541                                         msfux, msfuy, msfvx, msfvx_inv, &
 542                                         msfvy,msftx, msfty,             &
 543                                         khdif, xkmhd, rdx, rdy,         &
 544                                         ids, ide, jds, jde, kds, kde,   &
 545                                         ims, ime, jms, jme, kms, kme,   &
 546                                         its, ite, jts, jte, kts, kte   )
 547
 548         CALL horizontal_diffusion ('v', v, rv_tendf, mut, config_flags, &
 549                                         msfux, msfuy, msfvx, msfvx_inv, &
 550                                         msfvy,msftx, msfty,             &
 551                                         khdif, xkmhd, rdx, rdy,         &
 552                                         ids, ide, jds, jde, kds, kde,   &
 553                                         ims, ime, jms, jme, kms, kme,   &
 554                                         its, ite, jts, jte, kts, kte   )
 555
 556         CALL horizontal_diffusion ('w', w, rw_tendf, mut, config_flags, &
 557                                         msfux, msfuy, msfvx, msfvx_inv, &
 558                                         msfvy,msftx, msfty,             &
 559                                         khdif, xkmhd, rdx, rdy,         &
 560                                         ids, ide, jds, jde, kds, kde,   &
 561                                         ims, ime, jms, jme, kms, kme,   &
 562                                         its, ite, jts, jte, kts, kte   )
 563
 564         khdq = 3.*khdif
 565         CALL horizontal_diffusion_3dmp ( 'm', t, t_tendf, mut,            &
 566                                          config_flags, t_init,            &
 567                                          msfux, msfuy, msfvx, msfvx_inv,  &
 568                                          msfvy, msftx, msfty,             &
 569                                          khdq , xkhh, rdx, rdy,           &
 570                                          ids, ide, jds, jde, kds, kde,    &
 571                                          ims, ime, jms, jme, kms, kme,    &
 572                                          its, ite, jts, jte, kts, kte    )
 573
 574         pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
 575
 576           CALL vertical_diffusion_u ( u, ru_tendf, config_flags,      &
 577                                       u_base,                         &
 578                                       alt, muu, rdn, rdnw, kvdif,     &
 579                                       ids, ide, jds, jde, kds, kde,   &
 580                                       ims, ime, jms, jme, kms, kme,   &
 581                                       its, ite, jts, jte, kts, kte   )
 582
 583           CALL vertical_diffusion_v ( v, rv_tendf, config_flags,      &
 584                                       v_base,                         &
 585                                       alt, muv, rdn, rdnw, kvdif,     &
 586                                       ids, ide, jds, jde, kds, kde,   &
 587                                       ims, ime, jms, jme, kms, kme,   &
 588                                       its, ite, jts, jte, kts, kte   )
 589
 590           IF (non_hydrostatic)                                           &
 591           CALL vertical_diffusion ( 'w', w, rw_tendf, config_flags,      &
 592                                     alt, mut, rdn, rdnw, kvdif,          &
 593                                     ids, ide, jds, jde, kds, kde,        &
 594                                     ims, ime, jms, jme, kms, kme,        &
 595                                     its, ite, jts, jte, kts, kte        )
 596
 597           kvdq = 3.*kvdif
 598           CALL vertical_diffusion_3dmp ( t, t_tendf, config_flags, t_init,     &
 599                                          alt, mut, rdn, rdnw, kvdq ,           &
 600                                          ids, ide, jds, jde, kds, kde,         &
 601                                          ims, ime, jms, jme, kms, kme,         &
 602                                          its, ite, jts, jte, kts, kte         )
 603
 604         ENDIF pbl_test
 605
 606    !  Theta tendency computations.
 607
 608     END IF diff_opt1
 609
 610     IF ( diff_6th_opt .NE. 0 ) THEN
 611
 612       CALL sixth_order_diffusion( 'u', u, ru_tendf, mut, dt,          &
 613                                        config_flags,                  &
 614                                        diff_6th_opt, diff_6th_factor, &
 615                                        ids, ide, jds, jde, kds, kde,  &
 616                                        ims, ime, jms, jme, kms, kme,  &
 617                                        its, ite, jts, jte, kts, kte )
 618
 619       CALL sixth_order_diffusion( 'v', v, rv_tendf, mut, dt,          &
 620                                        config_flags,                  &
 621                                        diff_6th_opt, diff_6th_factor, &
 622                                        ids, ide, jds, jde, kds, kde,  &
 623                                        ims, ime, jms, jme, kms, kme,  &
 624                                        its, ite, jts, jte, kts, kte )
 625
 626       IF (non_hydrostatic)                                            &
 627       CALL sixth_order_diffusion( 'w', w, rw_tendf, mut, dt,          &
 628                                        config_flags,                  &
 629                                        diff_6th_opt, diff_6th_factor, &
 630                                        ids, ide, jds, jde, kds, kde,  &
 631                                        ims, ime, jms, jme, kms, kme,  &
 632                                        its, ite, jts, jte, kts, kte )
 633
 634       CALL sixth_order_diffusion( 'm', t,  t_tendf, mut, dt,          &
 635                                        config_flags,                  &
 636                                        diff_6th_opt, diff_6th_factor, &
 637                                        ids, ide, jds, jde, kds, kde,  &
 638                                        ims, ime, jms, jme, kms, kme,  &
 639                                        its, ite, jts, jte, kts, kte )
 640
 641     ENDIF
 642
 643     IF( damp_opt .eq. 2 )                                      &
 644        CALL rk_rayleigh_damp( ru_tendf, rv_tendf,              &
 645                               rw_tendf, t_tendf,               &
 646                               u, v, w, t, t_init,              &
 647                               mut, muu, muv, ph, phb,          &
 648                               u_base, v_base, t_base, z_base,  &
 649                               dampcoef, zdamp,                 &
 650                               ids, ide, jds, jde, kds, kde,    &
 651                               ims, ime, jms, jme, kms, kme,    &
 652                               its, ite, jts, jte, kts, kte   )
 653
 654   END IF forward_step
 655
 656 END SUBROUTINE rk_tendency
 657
 658 !-------------------------------------------------------------------------------
 659
 660 SUBROUTINE rk_addtend_dry ( ru_tend, rv_tend, rw_tend, ph_tend, t_tend,      &
 661                             ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf, &
 662                             u_save, v_save, w_save, ph_save, t_save,         &
 663                             mu_tend, mu_tendf, rk_step,                      &
 664                             h_diabatic, mut, msftx, msfty, msfux, msfuy,     &
 665                             msfvx, msfvx_inv, msfvy,                         &
 666                             ids,ide, jds,jde, kds,kde,                       &
 667                             ims,ime, jms,jme, kms,kme,                       &
 668                             ips,ipe, jps,jpe, kps,kpe,                       &
 669                             its,ite, jts,jte, kts,kte                       )
 670
 671    IMPLICIT NONE
 672
 673    !  Input data.
 674
 675    INTEGER ,               INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 676                                             ims, ime, jms, jme, kms, kme, &
 677                                             ips, ipe, jps, jpe, kps, kpe, &
 678                                             its, ite, jts, jte, kts, kte
 679    INTEGER ,               INTENT(IN   ) :: rk_step
 680
 681    REAL , DIMENSION( ims:ime , kms:kme, jms:jme  ) , INTENT(INOUT) :: ru_tend, &
 682                                                                       rv_tend, &
 683                                                                       rw_tend, &
 684                                                                       ph_tend, &
 685                                                                       t_tend,  &
 686                                                                       ru_tendf, &
 687                                                                       rv_tendf, &
 688                                                                       rw_tendf, &
 689                                                                       ph_tendf, &
 690                                                                       t_tendf
 691
 692    REAL , DIMENSION( ims:ime , jms:jme  ) , INTENT(INOUT) :: mu_tend, &
 693                                                              mu_tendf
 694
 695    REAL , DIMENSION( ims:ime , kms:kme, jms:jme  ) , INTENT(IN   ) ::  u_save,  &
 696                                                                        v_save,  &
 697                                                                        w_save,  &
 698                                                                       ph_save,  &
 699                                                                        t_save,  &
 700                                                                       h_diabatic
 701
 702    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: mut,       &
 703                                                                     msftx,     &
 704                                                                     msfty,     &
 705                                                                     msfux,     &
 706                                                                     msfuy,     &
 707                                                                     msfvx,     &
 708                                                                     msfvx_inv, &
 709                                                                     msfvy
 710
 711
 712 ! Local
 713    INTEGER :: i, j, k
 714
 715 !<DESCRIPTION>
 716 !
 717 ! rk_addtend_dry constructs the full large-timestep tendency terms for
 718 ! momentum (u,v,w), theta and geopotential equations.   This is accomplished
 719 ! by combining the physics tendencies (in *tendf; these are computed
 720 ! the first RK substep, held fixed thereafter) with the RK tendencies
 721 ! (in *tend, these include advection, pressure gradient, etc;
 722 ! these change each rk substep).  Output is in *tend.
 723 !
 724 !</DESCRIPTION>
 725
 726 !  Finally, add the forward-step tendency to the rk_tendency
 727
 728 ! u/v/w/save contain bc tendency that needs to be multiplied by msf
 729 ! (u by msfuy, v by msfvx)
 730 !  before adding it to physics tendency (*tendf)
 731 ! For momentum we need the final tendency to include an inverse msf
 732 ! physics/bc tendency needs to be divided, advection tendency already has it
 733
 734 ! For scalars we need the final tendency to include an inverse msf (msfty)
 735 ! advection tendency is OK, physics/bc tendency needs to be divided by msf
 736
 737    DO j = jts,MIN(jte,jde-1)
 738    DO k = kts,kte-1
 739    DO i = its,ite
 740      ! multiply by my to uncouple u
 741      IF(rk_step == 1)ru_tendf(i,k,j) = ru_tendf(i,k,j) +  u_save(i,k,j)*msfuy(i,j)
 742      ! divide by my to couple u
 743      ru_tend(i,k,j) = ru_tend(i,k,j) + ru_tendf(i,k,j)/msfuy(i,j)
 744    ENDDO
 745    ENDDO
 746    ENDDO
 747
 748    DO j = jts,jte
 749    DO k = kts,kte-1
 750    DO i = its,MIN(ite,ide-1)
 751      ! multiply by mx to uncouple v
 752      IF(rk_step == 1)rv_tendf(i,k,j) = rv_tendf(i,k,j) +  v_save(i,k,j)*msfvx(i,j)
 753      ! divide by mx to couple v
 754      rv_tend(i,k,j) = rv_tend(i,k,j) + rv_tendf(i,k,j)*msfvx_inv(i,j)
 755    ENDDO
 756    ENDDO
 757    ENDDO
 758
 759    DO j = jts,MIN(jte,jde-1)
 760    DO k = kts,kte
 761    DO i = its,MIN(ite,ide-1)
 762      ! multiply by my to uncouple w
 763      IF(rk_step == 1)rw_tendf(i,k,j) = rw_tendf(i,k,j) +  w_save(i,k,j)*msfty(i,j)
 764      ! divide by my to couple w
 765      rw_tend(i,k,j) = rw_tend(i,k,j) + rw_tendf(i,k,j)/msfty(i,j)
 766      IF(rk_step == 1)ph_tendf(i,k,j) = ph_tendf(i,k,j) +  ph_save(i,k,j)
 767      ! divide by my to couple scalar
 768      ph_tend(i,k,j) = ph_tend(i,k,j) + ph_tendf(i,k,j)/msfty(i,j)
 769    ENDDO
 770    ENDDO
 771    ENDDO
 772
 773    DO j = jts,MIN(jte,jde-1)
 774    DO k = kts,kte-1
 775    DO i = its,MIN(ite,ide-1)
 776      IF(rk_step == 1)t_tendf(i,k,j) = t_tendf(i,k,j) +  t_save(i,k,j)
 777      ! divide by my to couple theta
 778       t_tend(i,k,j) =  t_tend(i,k,j) +  t_tendf(i,k,j)/msfty(i,j)  &
 779                                      +  mut(i,j)*h_diabatic(i,k,j)/msfty(i,j)
 780      ! divide by my to couple heating
 781    ENDDO
 782    ENDDO
 783    ENDDO
 784
 785    DO j = jts,MIN(jte,jde-1)
 786    DO i = its,MIN(ite,ide-1)
 787 ! mu tendencies not coupled with 1/msf
 788       mu_tend(i,j) =  mu_tend(i,j) +  mu_tendf(i,j)
 789    ENDDO
 790    ENDDO
 791
 792 END SUBROUTINE rk_addtend_dry
 793
 794 !-------------------------------------------------------------------------------
 795
 796 SUBROUTINE rk_scalar_tend ( scs, sce, config_flags,          &
 797                             rk_step, dt,                     &
 798                             ru, rv, ww, mut, mub, mu_old,    &
 799                             alt,                             &
 800                             scalar_old, scalar,              &
 801                             scalar_tends, advect_tend,       &
 802                             RQVFTEN,                         &
 803                             base, moist_step, fnm, fnp,      &
 804                             msfux, msfuy, msfvx, msfvx_inv,  &
 805                             msfvy, msftx, msfty,             &
 806                             rdx, rdy, rdn, rdnw,             &
 807                             khdif, kvdif, xkmhd,             &
 808                             diff_6th_opt, diff_6th_factor,   &
 809                             adv_opt,                         &
 810                             ids, ide, jds, jde, kds, kde,    &
 811                             ims, ime, jms, jme, kms, kme,    &
 812                             its, ite, jts, jte, kts, kte    )
 813
 814    IMPLICIT NONE
 815
 816    !  Input data.
 817
 818    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
 819
 820    INTEGER ,                INTENT(IN   ) :: rk_step, scs, sce
 821    INTEGER ,                INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
 822                                              ims, ime, jms, jme, kms, kme, &
 823                                              its, ite, jts, jte, kts, kte
 824
 825    LOGICAL , INTENT(IN   ) :: moist_step
 826
 827    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ),                &
 828                                          INTENT(IN   )  :: scalar, scalar_old
 829
 830    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ),                      &
 831                                          INTENT(INOUT)  :: scalar_tends
 832
 833    REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(INOUT) :: advect_tend
 834
 835    REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(OUT  ) :: RQVFTEN
 836
 837    REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(IN   ) ::     ru,  &
 838                                                                       rv,  &
 839                                                                       ww,  &
 840                                                                       xkmhd,  &
 841                                                                       alt
 842
 843
 844    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fnm,  &
 845                                                                   fnp,  &
 846                                                                   rdn,  &
 847                                                                   rdnw, &
 848                                                                   base
 849
 850    REAL , DIMENSION( ims:ime , jms:jme ) ,       INTENT(IN   ) :: msfux,    &
 851                                                                   msfuy,    &
 852                                                                   msfvx,    &
 853                                                                   msfvx_inv,    &
 854                                                                   msfvy,    &
 855                                                                   msftx,    &
 856                                                                   msfty,    &
 857                                                                   mub,     &
 858                                                                   mut,     &
 859                                                                   mu_old
 860
 861    REAL ,                                        INTENT(IN   ) :: rdx,     &
 862                                                                   rdy,     &
 863                                                                   khdif,   &
 864                                                                   kvdif
 865
 866    INTEGER, INTENT( IN ) :: diff_6th_opt
 867    REAL,    INTENT( IN ) :: diff_6th_factor
 868
 869    REAL ,                                        INTENT(IN   ) :: dt
 870
 871    INTEGER, INTENT(IN   ) :: adv_opt
 872
 873    ! Local data
 874
 875    INTEGER :: im, i,j,k
 876    INTEGER :: time_step
 877
 878    REAL    :: khdq, kvdq, tendency
 879
 880 !<DESCRIPTION>
 881 !
 882 ! rk_scalar_tend calls routines that computes scalar tendency from advection
 883 ! and 3D mixing (TKE or fixed eddy viscosities).
 884 !
 885 !</DESCRIPTION>
 886
 887
 888    khdq = khdif/prandtl
 889    kvdq = kvdif/prandtl
 890
 891    scalar_loop : DO im = scs, sce
 892
 893      CALL zero_tend ( advect_tend(ims,kms,jms),     &
 894                       ids, ide, jds, jde, kds, kde, &
 895                       ims, ime, jms, jme, kms, kme, &
 896                       its, ite, jts, jte, kts, kte )
 897
 898      CALL nl_get_time_step ( 1, time_step )
 899
 900       IF( (rk_step == 3) .and. (adv_opt == POSITIVEDEF) ) THEN
 901
 902         CALL advect_scalar_pd       ( scalar(ims,kms,jms,im),             &
 903                                       scalar_old(ims,kms,jms,im),         &
 904                                       advect_tend(ims,kms,jms),           &
 905                                       ru, rv, ww, mut, mub, mu_old,       &
 906                                       time_step, config_flags,            &
 907                                       msfux, msfuy, msfvx, msfvy,         &
 908                                       msftx, msfty, fnm, fnp,             &
 909                                       rdx, rdy, rdnw,dt,                  &
 910                                       ids, ide, jds, jde, kds, kde,       &
 911                                       ims, ime, jms, jme, kms, kme,       &
 912                                       its, ite, jts, jte, kts, kte     )
 913
 914       ELSE IF( (rk_step == 3) .and. (adv_opt == MONOTONIC) ) THEN
 915
 916         CALL advect_scalar_mono       ( scalar(ims,kms,jms,im),             &
 917                                         scalar_old(ims,kms,jms,im),         &
 918                                         advect_tend(ims,kms,jms),           &
 919                                         ru, rv, ww, mut, mub, mu_old,       &
 920                                         config_flags,                       &
 921                                         msfux, msfuy, msfvx, msfvy,         &
 922                                         msftx, msfty, fnm, fnp,             &
 923                                         rdx, rdy, rdnw,dt,                  &
 924                                         ids, ide, jds, jde, kds, kde,       &
 925                                         ims, ime, jms, jme, kms, kme,       &
 926                                         its, ite, jts, jte, kts, kte     )
 927
 928       ELSE
 929
 930         CALL advect_scalar     ( scalar(ims,kms,jms,im),        &
 931                                  scalar(ims,kms,jms,im),        &
 932                                  advect_tend(ims,kms,jms),      &
 933                                  ru, rv, ww, mut, time_step,    &
 934                                  config_flags,                  &
 935                                  msfux, msfuy, msfvx, msfvy,    &
 936                                  msftx, msfty, fnm, fnp,        &
 937                                  rdx, rdy, rdnw,                &
 938                                  ids, ide, jds, jde, kds, kde,  &
 939                                  ims, ime, jms, jme, kms, kme,  &
 940                                  its, ite, jts, jte, kts, kte  )
 941
 942       END IF
 943
 944      IF((config_flags%cu_physics == GDSCHEME .OR. config_flags%cu_physics == G3SCHEME) &
 945                      .and. moist_step .and. ( im == P_QV) ) THEN
 946
 947         CALL set_tend( RQVFTEN, advect_tend, msfty,    &
 948                        ids, ide, jds, jde, kds, kde,   &
 949                        ims, ime, jms, jme, kms, kme,   &
 950                        its, ite, jts, jte, kts, kte      )
 951      ENDIF
 952
 953      rk_step_1: IF( rk_step == 1 ) THEN
 954
 955        diff_opt1 : IF (config_flags%diff_opt .eq. 1) THEN
 956
 957        CALL horizontal_diffusion ( 'm', scalar(ims,kms,jms,im),            &
 958                                         scalar_tends(ims,kms,jms,im), mut, &
 959                                         config_flags,                      &
 960                                         msfux, msfuy, msfvx, msfvx_inv,    &
 961                                         msfvy, msftx, msfty,               &
 962                                         khdq , xkmhd, rdx, rdy,            &
 963                                         ids, ide, jds, jde, kds, kde,      &
 964                                         ims, ime, jms, jme, kms, kme,      &
 965                                         its, ite, jts, jte, kts, kte      )
 966
 967        pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
 968
 969          IF( (moist_step) .and. ( im == P_QV)) THEN
 970
 971             CALL vertical_diffusion_mp ( scalar(ims,kms,jms,im),       &
 972                                          scalar_tends(ims,kms,jms,im), &
 973                                          config_flags, base,           &
 974                                          alt, mut, rdn, rdnw, kvdq ,   &
 975                                          ids, ide, jds, jde, kds, kde, &
 976                                          ims, ime, jms, jme, kms, kme, &
 977                                          its, ite, jts, jte, kts, kte )
 978
 979          ELSE
 980
 981             CALL vertical_diffusion (  'm', scalar(ims,kms,jms,im),       &
 982                                             scalar_tends(ims,kms,jms,im), &
 983                                             config_flags,                 &
 984                                             alt, mut, rdn, rdnw, kvdq,    &
 985                                             ids, ide, jds, jde, kds, kde, &
 986                                             ims, ime, jms, jme, kms, kme, &
 987                                             its, ite, jts, jte, kts, kte )
 988
 989          END IF
 990
 991       ENDIF pbl_test
 992
 993     ENDIF diff_opt1
 994
 995     IF ( diff_6th_opt .NE. 0 )                                        &
 996       CALL sixth_order_diffusion( 'm', scalar(ims,kms,jms,im),        &
 997                                        scalar_tends(ims,kms,jms,im),  &
 998                                        mut, dt, config_flags,         &
 999                                        diff_6th_opt, diff_6th_factor, &
1000                                        ids, ide, jds, jde, kds, kde,  &
1001                                        ims, ime, jms, jme, kms, kme,  &
1002                                        its, ite, jts, jte, kts, kte )
1003
1004   ENDIF rk_step_1
1005
1006  END DO scalar_loop
1007
1008 END SUBROUTINE rk_scalar_tend
1009
1010 !-------------------------------------------------------------------------------
1011
1012 SUBROUTINE rk_update_scalar( scs, sce,                      &
1013                              scalar_1, scalar_2, sc_tend,   &
1014                              advect_tend, msftx, msfty,     &
1015                              mu_old, mu_new, mu_base,       &
1016                              rk_step, dt, spec_zone,        &
1017                              config_flags,                  &
1018                              ids, ide, jds, jde, kds, kde,  &
1019                              ims, ime, jms, jme, kms, kme,  &
1020                              its, ite, jts, jte, kts, kte  )
1021
1022    IMPLICIT NONE
1023
1024    !  Input data.
1025
1026    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
1027
1028    INTEGER ,                INTENT(IN   ) :: scs, sce, rk_step, spec_zone
1029    INTEGER ,                INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1030                                              ims, ime, jms, jme, kms, kme, &
1031                                              its, ite, jts, jte, kts, kte
1032
1033    REAL,                    INTENT(IN   ) :: dt
1034
1035    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
1036          INTENT(INOUT)                                  :: scalar_1,    &
1037                                                            scalar_2
1038
1039    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
1040          INTENT(IN)                                     :: sc_tend
1041
1042    REAL, DIMENSION(ims:ime, kms:kme, jms:jme ),                &
1043          INTENT(IN)                                  :: advect_tend
1044
1045    REAL, DIMENSION(ims:ime, jms:jme  ), INTENT(IN   ) ::  mu_old,  &
1046                                                           mu_new,  &
1047                                                           mu_base, &
1048                                                           msftx,   &
1049                                                           msfty
1050
1051    INTEGER :: i,j,k,im
1052    REAL    :: sc_middle, msfsq
1053    REAL, DIMENSION(its:ite) :: muold, r_munew
1054
1055    REAL, DIMENSION(its:ite, kts:kte, jts:jte  ) :: tendency
1056
1057    INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
1058    INTEGER :: i_start_spc,i_end_spc,j_start_spc,j_end_spc,k_start_spc,k_end_spc
1059
1060 !<DESCRIPTION>
1061 !
1062 !  rk_scalar_update advances the scalar equation given the time t value
1063 !  of the scalar and the scalar tendency.
1064 !
1065 !</DESCRIPTION>
1066
1067
1068 !
1069 !  set loop limits.
1070
1071       i_start = its
1072       i_end   = min(ite,ide-1)
1073       j_start = jts
1074       j_end   = min(jte,jde-1)
1075       k_start = kts
1076       k_end   = kte-1
1077
1078       i_start_spc = i_start
1079       i_end_spc   = i_end
1080       j_start_spc = j_start
1081       j_end_spc   = j_end
1082       k_start_spc = k_start
1083       k_end_spc   = k_end
1084
1085     IF( config_flags%nested .or. config_flags%specified ) THEN
1086       IF( .NOT. config_flags%periodic_x)i_start = max( its,ids+spec_zone )
1087       IF( .NOT. config_flags%periodic_x)i_end   = min( ite,ide-spec_zone-1 )
1088       j_start = max( jts,jds+spec_zone )
1089       j_end   = min( jte,jde-spec_zone-1 )
1090       k_start = kts
1091       k_end   = min( kte, kde-1 )
1092     ENDIF
1093
1094     IF ( rk_step == 1 ) THEN
1095
1096       !  replace t-dt values (in scalar_1) with t values scalar_2,
1097       !  then compute new values by adding tendency to values at t
1098
1099       DO  im = scs,sce
1100
1101        DO  j = jts, min(jte,jde-1)
1102        DO  k = kts, min(kte,kde-1)
1103        DO  i = its, min(ite,ide-1)
1104            tendency(i,k,j) = 0.
1105        ENDDO
1106        ENDDO
1107        ENDDO
1108
1109        DO  j = j_start,j_end
1110        DO  k = k_start,k_end
1111        DO  i = i_start,i_end
1112           ! scalar was coupled with my
1113            tendency(i,k,j) = advect_tend(i,k,j) * msfty(i,j)
1114        ENDDO
1115        ENDDO
1116        ENDDO
1117
1118        DO  j = j_start_spc,j_end_spc
1119        DO  k = k_start_spc,k_end_spc
1120        DO  i = i_start_spc,i_end_spc
1121            tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1122        ENDDO
1123        ENDDO
1124        ENDDO
1125
1126       DO  j = jts, min(jte,jde-1)
1127
1128       DO  i = its, min(ite,ide-1)
1129         muold(i) = mu_old(i,j) + mu_base(i,j)
1130         r_munew(i) = 1./(mu_new(i,j) + mu_base(i,j))
1131       ENDDO
1132
1133       DO  k = kts, min(kte,kde-1)
1134       DO  i = its, min(ite,ide-1)
1135
1136         scalar_1(i,k,j,im) = scalar_2(i,k,j,im)
1137         scalar_2(i,k,j,im) = (muold(i)*scalar_1(i,k,j,im)   &
1138                              + dt*tendency(i,k,j))*r_munew(i)
1139
1140       ENDDO
1141       ENDDO
1142       ENDDO
1143
1144       ENDDO
1145
1146     ELSE
1147
1148       !  just compute new values, scalar_1 already at time t.
1149
1150       DO  im = scs, sce
1151
1152        DO  j = jts, min(jte,jde-1)
1153        DO  k = kts, min(kte,kde-1)
1154        DO  i = its, min(ite,ide-1)
1155            tendency(i,k,j) = 0.
1156        ENDDO
1157        ENDDO
1158        ENDDO
1159
1160        DO  j = j_start,j_end
1161        DO  k = k_start,k_end
1162        DO  i = i_start,i_end
1163            ! scalar was coupled with my
1164            tendency(i,k,j) = advect_tend(i,k,j) * msfty(i,j)
1165        ENDDO
1166        ENDDO
1167        ENDDO
1168
1169        DO  j = j_start_spc,j_end_spc
1170        DO  k = k_start_spc,k_end_spc
1171        DO  i = i_start_spc,i_end_spc
1172            tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1173        ENDDO
1174        ENDDO
1175        ENDDO
1176
1177       DO  j = jts, min(jte,jde-1)
1178
1179       DO  i = its, min(ite,ide-1)
1180         muold(i) = mu_old(i,j) + mu_base(i,j)
1181         r_munew(i) = 1./(mu_new(i,j) + mu_base(i,j))
1182       ENDDO
1183
1184       DO  k = kts, min(kte,kde-1)
1185       DO  i = its, min(ite,ide-1)
1186
1187         scalar_2(i,k,j,im) = (muold(i)*scalar_1(i,k,j,im)   &
1188                              + dt*tendency(i,k,j))*r_munew(i)
1189
1190       ENDDO
1191       ENDDO
1192       ENDDO
1193
1194       ENDDO
1195
1196     END IF
1197
1198 END SUBROUTINE rk_update_scalar
1199
1200 !-------------------------------------------------------------------------------
1201
1202 SUBROUTINE rk_update_scalar_pd( scs, sce,                      &
1203                                 scalar, sc_tend,               &
1204                                 mu_old, mu_new, mu_base,       &
1205                                 rk_step, dt, spec_zone,        &
1206                                 config_flags,                  &
1207                                 ids, ide, jds, jde, kds, kde,  &
1208                                 ims, ime, jms, jme, kms, kme,  &
1209                                 its, ite, jts, jte, kts, kte  )
1210
1211    IMPLICIT NONE
1212
1213    !  Input data.
1214
1215    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
1216
1217    INTEGER ,                INTENT(IN   ) :: scs, sce, rk_step, spec_zone
1218    INTEGER ,                INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1219                                              ims, ime, jms, jme, kms, kme, &
1220                                              its, ite, jts, jte, kts, kte
1221
1222    REAL,                    INTENT(IN   ) :: dt
1223
1224    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
1225          INTENT(INOUT)                                  :: scalar,      &
1226                                                            sc_tend
1227
1228    REAL, DIMENSION(ims:ime, jms:jme  ), INTENT(IN   ) ::  mu_old,  &
1229                                                           mu_new,  &
1230                                                           mu_base
1231
1232    INTEGER :: i,j,k,im
1233    REAL    :: sc_middle, msfsq
1234    REAL, DIMENSION(its:ite) :: muold, r_munew
1235
1236    REAL, DIMENSION(its:ite, kts:kte, jts:jte  ) :: tendency
1237
1238    INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
1239    INTEGER :: i_start_spc,i_end_spc,j_start_spc,j_end_spc,k_start_spc,k_end_spc
1240
1241 !<DESCRIPTION>
1242 !
1243 !  rk_scalar_update advances the scalar equation given the time t value
1244 !  of the scalar and the scalar tendency.
1245 !
1246 !</DESCRIPTION>
1247
1248
1249 !
1250 !  set loop limits.
1251
1252       i_start = its
1253       i_end   = min(ite,ide-1)
1254       j_start = jts
1255       j_end   = min(jte,jde-1)
1256       k_start = kts
1257       k_end   = kte-1
1258
1259       i_start_spc = i_start
1260       i_end_spc   = i_end
1261       j_start_spc = j_start
1262       j_end_spc   = j_end
1263       k_start_spc = k_start
1264       k_end_spc   = k_end
1265
1266     IF( config_flags%nested .or. config_flags%specified ) THEN
1267       IF( .NOT. config_flags%periodic_x)i_start = max( its,ids+spec_zone )
1268       IF( .NOT. config_flags%periodic_x)i_end   = min( ite,ide-spec_zone-1 )
1269       j_start = max( jts,jds+spec_zone )
1270       j_end   = min( jte,jde-spec_zone-1 )
1271       k_start = kts
1272       k_end   = min( kte, kde-1 )
1273     ENDIF
1274
1275       DO  im = scs, sce
1276
1277        DO  j = jts, min(jte,jde-1)
1278        DO  k = kts, min(kte,kde-1)
1279        DO  i = its, min(ite,ide-1)
1280            tendency(i,k,j) = 0.
1281        ENDDO
1282        ENDDO
1283        ENDDO
1284
1285        DO  j = j_start_spc,j_end_spc
1286        DO  k = k_start_spc,k_end_spc
1287        DO  i = i_start_spc,i_end_spc
1288            tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1289            sc_tend(i,k,j,im) = 0.
1290        ENDDO
1291        ENDDO
1292        ENDDO
1293
1294       DO  j = jts, min(jte,jde-1)
1295
1296       DO  i = its, min(ite,ide-1)
1297         muold(i) = mu_old(i,j) + mu_base(i,j)
1298         r_munew(i) = 1./(mu_new(i,j) + mu_base(i,j))
1299       ENDDO
1300
1301       DO  k = kts, min(kte,kde-1)
1302       DO  i = its, min(ite,ide-1)
1303
1304         scalar(i,k,j,im) = (muold(i)*scalar(i,k,j,im)   &
1305                              + dt*tendency(i,k,j))*r_munew(i)
1306       ENDDO
1307       ENDDO
1308       ENDDO
1309
1310       ENDDO
1311
1312 END SUBROUTINE rk_update_scalar_pd
1313
1314 !------------------------------------------------------------
1315
1316 SUBROUTINE init_zero_tendency(ru_tendf, rv_tendf, rw_tendf, ph_tendf,  &
1317                               t_tendf,  tke_tendf, mu_tendf,           &
1318                               moist_tendf,chem_tendf,scalar_tendf,     &
1319                               tracer_tendf,n_tracer,                   &
1320                               n_moist,n_chem,n_scalar,rk_step,         &
1321                               ids, ide, jds, jde, kds, kde,            &
1322                               ims, ime, jms, jme, kms, kme,            &
1323                               its, ite, jts, jte, kts, kte             )
1324 !-----------------------------------------------------------------------
1325    IMPLICIT NONE
1326 !-----------------------------------------------------------------------
1327
1328    INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1329                                     ims, ime, jms, jme, kms, kme, &
1330                                     its, ite, jts, jte, kts, kte
1331
1332    INTEGER ,       INTENT(IN   ) :: n_moist,n_chem,n_scalar,n_tracer,rk_step
1333
1334    REAL , DIMENSION( ims:ime , kms:kme, jms:jme  ) , INTENT(INOUT) ::  &
1335                                                              ru_tendf, &
1336                                                              rv_tendf, &
1337                                                              rw_tendf, &
1338                                                              ph_tendf, &
1339                                                               t_tendf, &
1340                                                             tke_tendf
1341
1342    REAL , DIMENSION( ims:ime , jms:jme  ) , INTENT(INOUT) ::  mu_tendf
1343
1344    REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_moist),INTENT(INOUT)::&
1345                                                           moist_tendf
1346
1347    REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_chem ),INTENT(INOUT)::&
1348                                                           chem_tendf
1349    REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_tracer ),INTENT(INOUT)::&
1350                                                           tracer_tendf
1351
1352    REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_scalar ),INTENT(INOUT)::&
1353                                                           scalar_tendf
1354
1355 ! LOCAL VARS
1356
1357    INTEGER :: im, ic, is
1358
1359 !<DESCRIPTION>
1360 !
1361 ! init_zero_tendency
1362 ! sets tendency arrays to zero for all prognostic variables.
1363 !
1364 !</DESCRIPTION>
1365
1366
1367    CALL zero_tend ( ru_tendf,                        &
1368                     ids, ide, jds, jde, kds, kde,    &
1369                     ims, ime, jms, jme, kms, kme,    &
1370                     its, ite, jts, jte, kts, kte     )
1371
1372    CALL zero_tend ( rv_tendf,                        &
1373                     ids, ide, jds, jde, kds, kde,    &
1374                     ims, ime, jms, jme, kms, kme,    &
1375                     its, ite, jts, jte, kts, kte     )
1376
1377    CALL zero_tend ( rw_tendf,                        &
1378                     ids, ide, jds, jde, kds, kde,    &
1379                     ims, ime, jms, jme, kms, kme,    &
1380                     its, ite, jts, jte, kts, kte     )
1381
1382    CALL zero_tend ( ph_tendf,                        &
1383                     ids, ide, jds, jde, kds, kde,    &
1384                     ims, ime, jms, jme, kms, kme,    &
1385                     its, ite, jts, jte, kts, kte     )
1386
1387    CALL zero_tend ( t_tendf,                         &
1388                     ids, ide, jds, jde, kds, kde,    &
1389                     ims, ime, jms, jme, kms, kme,    &
1390                     its, ite, jts, jte, kts, kte     )
1391
1392    CALL zero_tend ( tke_tendf,                       &
1393                     ids, ide, jds, jde, kds, kde,    &
1394                     ims, ime, jms, jme, kms, kme,    &
1395                     its, ite, jts, jte, kts, kte     )
1396
1397    CALL zero_tend ( mu_tendf,                        &
1398                     ids, ide, jds, jde, kds, kds,    &
1399                     ims, ime, jms, jme, kms, kms,    &
1400                     its, ite, jts, jte, kts, kts     )
1401
1402 !   DO im=PARAM_FIRST_SCALAR,n_moist
1403    DO im=1,n_moist                      ! make sure first one is zero too
1404       CALL zero_tend ( moist_tendf(ims,kms,jms,im),  &
1405                        ids, ide, jds, jde, kds, kde, &
1406                        ims, ime, jms, jme, kms, kme, &
1407                        its, ite, jts, jte, kts, kte  )
1408    ENDDO
1409
1410 !   DO ic=PARAM_FIRST_SCALAR,n_chem
1411    DO ic=1,n_chem                       ! make sure first one is zero too
1412       CALL zero_tend ( chem_tendf(ims,kms,jms,ic),   &
1413                        ids, ide, jds, jde, kds, kde, &
1414                        ims, ime, jms, jme, kms, kme, &
1415                        its, ite, jts, jte, kts, kte  )
1416    ENDDO
1417
1418 !   DO ic=PARAM_FIRST_SCALAR,n_tracer
1419    DO ic=1,n_tracer                     ! make sure first one is zero too
1420       CALL zero_tend ( tracer_tendf(ims,kms,jms,ic), &
1421                        ids, ide, jds, jde, kds, kde, &
1422                        ims, ime, jms, jme, kms, kme, &
1423                        its, ite, jts, jte, kts, kte  )
1424    ENDDO
1425
1426 !   DO ic=PARAM_FIRST_SCALAR,n_scalar
1427    DO ic=1,n_scalar                       ! make sure first one is zero too
1428       CALL zero_tend ( scalar_tendf(ims,kms,jms,ic),   &
1429                        ids, ide, jds, jde, kds, kde, &
1430                        ims, ime, jms, jme, kms, kme, &
1431                        its, ite, jts, jte, kts, kte  )
1432    ENDDO
1433
1434 END SUBROUTINE init_zero_tendency
1435
1436 !===================================================================
1437
1438
1439 SUBROUTINE dump_data( a, field, io_unit,            &
1440                       ims, ime, jms, jme, kms, kme, &
1441                       ids, ide, jds, jde, kds, kde )
1442 implicit none
1443 integer ::  ims, ime, jms, jme, kms, kme, &
1444             ids, ide, jds, jde, kds, kde
1445 real, dimension(ims:ime, kms:kme, jds:jde) :: a
1446 character :: field
1447 integer :: io_unit
1448
1449 integer :: is,ie,js,je,ks,ke
1450
1451 !<DESCRIPTION
1452 !
1453 ! quick and dirty debug io utility
1454 !
1455 !</DESCRIPTION
1456
1457 is = ids
1458 ie = ide-1
1459 js = jds
1460 je = jde-1
1461 ks = kds
1462 ke = kde-1
1463
1464 if(field == 'u') ie = ide
1465 if(field == 'v') je = jde
1466 if(field == 'w') ke = kde
1467
1468 write(io_unit) is,ie,ks,ke,js,je
1469 write(io_unit) a(is:ie, ks:ke, js:je)
1470
1471 end subroutine dump_data
1472
1473 !-----------------------------------------------------------------------
1474
1475 SUBROUTINE calculate_phy_tend (config_flags,mu,muu,muv,pi3d,           &
1476                      RTHRATEN,                                         &
1477                      RUBLTEN,RVBLTEN,RTHBLTEN,                         &
1478                      RQVBLTEN,RQCBLTEN,RQIBLTEN,                       &
1479                      RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN,              &
1480                      RQICUTEN,RQSCUTEN,                                &
1481                      RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,        &
1482                      RMUNDGDTEN,                                       &
1483                      ids,ide, jds,jde, kds,kde,                        &
1484                      ims,ime, jms,jme, kms,kme,                        &
1485                      its,ite, jts,jte, kts,kte                         )
1486 !-----------------------------------------------------------------------
1487       IMPLICIT NONE
1488
1489       TYPE(grid_config_rec_type), INTENT(IN)     ::      config_flags
1490
1491       INTEGER,  INTENT(IN   )   ::          ids,ide, jds,jde, kds,kde, &
1492                                             ims,ime, jms,jme, kms,kme, &
1493                                             its,ite, jts,jte, kts,kte
1494
1495       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )               , &
1496                 INTENT(IN   )   ::                               pi3d
1497
1498       REAL,     DIMENSION( ims:ime, jms:jme )                        , &
1499                 INTENT(IN   )   ::                                 mu, &
1500                                                                   muu, &
1501                                                                   muv
1502
1503
1504 ! radiation
1505
1506       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ),                &
1507                 INTENT(INOUT)   ::                           RTHRATEN
1508
1509 ! cumulus
1510
1511       REAL,     DIMENSION( ims:ime , kms:kme , jms:jme ),              &
1512                 INTENT(INOUT)   ::                                     &
1513                                                              RTHCUTEN, &
1514                                                              RQVCUTEN, &
1515                                                              RQCCUTEN, &
1516                                                              RQRCUTEN, &
1517                                                              RQICUTEN, &
1518                                                              RQSCUTEN
1519 ! pbl
1520
1521       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )               , &
1522                 INTENT(INOUT)   ::                            RUBLTEN, &
1523                                                               RVBLTEN, &
1524                                                              RTHBLTEN, &
1525                                                              RQVBLTEN, &
1526                                                              RQCBLTEN, &
1527                                                              RQIBLTEN
1528
1529 ! fdda
1530
1531       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )               , &
1532                 INTENT(INOUT)   ::                            RUNDGDTEN, &
1533                                                               RVNDGDTEN, &
1534                                                              RTHNDGDTEN, &
1535                                                              RQVNDGDTEN
1536       REAL,     DIMENSION( ims:ime, jms:jme )               , &
1537                 INTENT(INOUT)   ::                           RMUNDGDTEN
1538
1539       INTEGER :: i,k,j
1540       INTEGER :: itf,ktf,jtf,itsu,jtsv
1541
1542 !-----------------------------------------------------------------------
1543
1544 !<DESCRIPTION>
1545 !
1546 !  calculate_phy_tend couples the physics tendencies to the column mass (mu),
1547 !  because prognostic equations are in flux form, but physics tendencies are
1548 !  computed for uncoupled variables.
1549 !
1550 !</DESCRIPTION>
1551
1552       itf=MIN(ite,ide-1)
1553       jtf=MIN(jte,jde-1)
1554       ktf=MIN(kte,kde-1)
1555       itsu=MAX(its,ids+1)
1556       jtsv=MAX(jts,jds+1)
1557
1558 ! radiation
1559
1560    IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
1561
1562       DO J=jts,jtf
1563       DO K=kts,ktf
1564       DO I=its,itf
1565          RTHRATEN(I,K,J)=mu(I,J)*RTHRATEN(I,K,J)
1566       ENDDO
1567       ENDDO
1568       ENDDO
1569
1570    ENDIF
1571
1572 ! cumulus
1573
1574    IF (config_flags%cu_physics .gt. 0) THEN
1575
1576       DO J=jts,jtf
1577       DO I=its,itf
1578       DO K=kts,ktf
1579          RTHCUTEN(I,K,J)=mu(I,J)*RTHCUTEN(I,K,J)
1580          RQVCUTEN(I,K,J)=mu(I,J)*RQVCUTEN(I,K,J)
1581       ENDDO
1582       ENDDO
1583       ENDDO
1584
1585       IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
1586          DO J=jts,jtf
1587          DO I=its,itf
1588          DO K=kts,ktf
1589             RQCCUTEN(I,K,J)=mu(I,J)*RQCCUTEN(I,K,J)
1590          ENDDO
1591          ENDDO
1592          ENDDO
1593       ENDIF
1594
1595       IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
1596          DO J=jts,jtf
1597          DO I=its,itf
1598          DO K=kts,ktf
1599             RQRCUTEN(I,K,J)=mu(I,J)*RQRCUTEN(I,K,J)
1600          ENDDO
1601          ENDDO
1602          ENDDO
1603       ENDIF
1604
1605       IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
1606          DO J=jts,jtf
1607          DO I=its,itf
1608          DO K=kts,ktf
1609             RQICUTEN(I,K,J)=mu(I,J)*RQICUTEN(I,K,J)
1610          ENDDO
1611          ENDDO
1612          ENDDO
1613       ENDIF
1614
1615       IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
1616          DO J=jts,jtf
1617          DO I=its,itf
1618          DO K=kts,ktf
1619             RQSCUTEN(I,K,J)=mu(I,J)*RQSCUTEN(I,K,J)
1620          ENDDO
1621          ENDDO
1622          ENDDO
1623       ENDIF
1624
1625    ENDIF
1626
1627 ! pbl
1628
1629    IF (config_flags%bl_pbl_physics .gt. 0) THEN
1630
1631       DO J=jts,jtf
1632       DO K=kts,ktf
1633       DO I=its,itf
1634          RUBLTEN(I,K,J) =mu(I,J)*RUBLTEN(I,K,J)
1635          RVBLTEN(I,K,J) =mu(I,J)*RVBLTEN(I,K,J)
1636          RTHBLTEN(I,K,J)=mu(I,J)*RTHBLTEN(I,K,J)
1637       ENDDO
1638       ENDDO
1639       ENDDO
1640
1641       IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
1642          DO J=jts,jtf
1643          DO K=kts,ktf
1644          DO I=its,itf
1645             RQVBLTEN(I,K,J)=mu(I,J)*RQVBLTEN(I,K,J)
1646          ENDDO
1647          ENDDO
1648          ENDDO
1649       ENDIF
1650
1651       IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
1652          DO J=jts,jtf
1653          DO K=kts,ktf
1654          DO I=its,itf
1655            RQCBLTEN(I,K,J)=mu(I,J)*RQCBLTEN(I,K,J)
1656          ENDDO
1657          ENDDO
1658          ENDDO
1659       ENDIF
1660
1661       IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
1662          DO J=jts,jtf
1663          DO K=kts,ktf
1664          DO I=its,itf
1665             RQIBLTEN(I,K,J)=mu(I,J)*RQIBLTEN(I,K,J)
1666          ENDDO
1667          ENDDO
1668          ENDDO
1669       ENDIF
1670
1671     ENDIF
1672
1673 ! fdda
1674 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
1675 !   so only couple those
1676
1677    IF (config_flags%grid_fdda .gt. 0) THEN
1678
1679       DO J=jts,jtf
1680       DO K=kts,ktf
1681       DO I=itsu,itf
1682 !     if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1683 !     write(*,'(a,3i6,e15.5)') 'u_ten before=',i,k,j, RUNDGDTEN(i,k,j)
1684          RUNDGDTEN(I,K,J) =muu(I,J)*RUNDGDTEN(I,K,J)
1685 !        if( i == itf/2 .AND. j == jtf/2 .AND. k==ktf/2 ) &
1686 !          write(*,'(a,2f15.5)') 'mu, muu=',mu(i,j), muu(i,j)
1687 !     if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1688 !     write(*,'(a,3i6,e15.5)') 'u_ten after=',i,k,j, RUNDGDTEN(i,k,j)
1689 !     if( RUNDGDTEN(i,k,j) > 30.0 ) write(*,*) 'IKJ=',i,k,j
1690       ENDDO
1691       ENDDO
1692       ENDDO
1693 !     write(*,'(a,e15.5)') 'u_ten MAXIMUM after=', maxval(RUNDGDTEN)
1694       DO J=jtsv,jtf
1695       DO K=kts,ktf
1696       DO I=its,itf
1697          RVNDGDTEN(I,K,J) =muv(I,J)*RVNDGDTEN(I,K,J)
1698       ENDDO
1699       ENDDO
1700       ENDDO
1701       DO J=jts,jtf
1702       DO K=kts,ktf
1703       DO I=its,itf
1704 !     if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1705 !     write(*,'(a,3i6,e15.5)') 'th before=',i,k,j, RTHNDGDTEN(I,K,J)
1706          RTHNDGDTEN(I,K,J)=mu(I,J)*RTHNDGDTEN(I,K,J)
1707 !        RMUNDGDTEN(I,J) - no coupling
1708 !     if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1709 !     write(*,'(a,3i6,e15.5)') 'th after=',i,k,j, RTHNDGDTEN(I,K,J)
1710       ENDDO
1711       ENDDO
1712       ENDDO
1713
1714       IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
1715          DO J=jts,jtf
1716          DO K=kts,ktf
1717          DO I=its,itf
1718             RQVNDGDTEN(I,K,J)=mu(I,J)*RQVNDGDTEN(I,K,J)
1719          ENDDO
1720          ENDDO
1721          ENDDO
1722       ENDIF
1723
1724     ENDIF
1725
1726 END SUBROUTINE calculate_phy_tend
1727
1728 !-----------------------------------------------------------------------
1729
1730 SUBROUTINE positive_definite_filter ( a,                          &
1731                                       ids,ide, jds,jde, kds,kde,  &
1732                                       ims,ime, jms,jme, kms,kme,  &
1733                                       its,ite, jts,jte, kts,kte  )
1734
1735   IMPLICIT NONE
1736
1737   INTEGER,  INTENT(IN   )   ::          ids,ide, jds,jde, kds,kde, &
1738                                         ims,ime, jms,jme, kms,kme, &
1739                                         its,ite, jts,jte, kts,kte
1740
1741   REAL, DIMENSION( ims:ime , kms:kme , jms:jme  ), INTENT(INOUT) :: a
1742
1743   INTEGER :: i,k,j
1744
1745 !<DESCRIPTION>
1746 !
1747 ! debug and testing code for bounding a variable
1748 !
1749 !</DESCRIPTION>
1750
1751   DO j=jts,min(jte,jde-1)
1752   DO k=kts,kte-1
1753   DO i=its,min(ite,ide-1)
1754 !    a(i,k,j) = max(a(i,k,j),0.)
1755     a(i,k,j) = min(1000.,max(a(i,k,j),0.))
1756   ENDDO
1757   ENDDO
1758   ENDDO
1759
1760   END SUBROUTINE positive_definite_filter
1761
1762 !-----------------------------------------------------------------------
1763
1764 SUBROUTINE bound_tke ( tke, tke_upper_bound,       &
1765                        ids,ide, jds,jde, kds,kde,  &
1766                        ims,ime, jms,jme, kms,kme,  &
1767                        its,ite, jts,jte, kts,kte  )
1768
1769   IMPLICIT NONE
1770
1771   INTEGER,  INTENT(IN   )   ::          ids,ide, jds,jde, kds,kde, &
1772                                         ims,ime, jms,jme, kms,kme, &
1773                                         its,ite, jts,jte, kts,kte
1774
1775   REAL, DIMENSION( ims:ime , kms:kme , jms:jme  ), INTENT(INOUT) :: tke
1776   REAL, INTENT(   IN) :: tke_upper_bound
1777
1778   INTEGER :: i,k,j
1779
1780 !<DESCRIPTION>
1781 !
1782 ! bounds tke between zero and tke_upper_bound.
1783 !
1784 !</DESCRIPTION>
1785
1786   DO j=jts,min(jte,jde-1)
1787   DO k=kts,kte-1
1788   DO i=its,min(ite,ide-1)
1789     tke(i,k,j) = min(tke_upper_bound,max(tke(i,k,j),0.))
1790   ENDDO
1791   ENDDO
1792   ENDDO
1793
1794   END SUBROUTINE bound_tke
1795
1796
1797
1798 END MODULE module_em