wrf svn trunk commit r4103
[wrffire.git] / wrfv2_fire / dyn_em / module_em.F
blobf4785fe874b45eff96e012c49c2de8bb0ef96d1c
1 !WRF:MODEL_LAYER:DYNAMICS
4 MODULE module_em
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
12 CONTAINS
14 !------------------------------------------------------------------------
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    )
32    IMPLICIT NONE
35    !  Input data.
37    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
39    INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
40                                     ims, ime, jms, jme, kms, kme, &
41                                     its, ite, jts, jte, kts, kte
43    INTEGER ,       INTENT(IN   ) :: n_moist, rk_step
45    REAL ,          INTENT(IN   ) :: rdx, rdy
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
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
69    REAL , DIMENSION(  ims:ime , kms:kme, jms:jme ) ,                      &
70                                                INTENT(IN   ) ::  p
71                                                                  
75    REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist ), INTENT(   IN) :: &
76                                                            moist
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
88    REAL , DIMENSION( ims:ime , jms:jme ) ,    INTENT(  OUT) :: muu,    &
89                                                                muv,    &
90                                                                mut
92    REAL , DIMENSION( kms:kme ) ,    INTENT(IN   ) :: fnm, fnp, dnw
94    integer :: k
97 !<DESCRIPTION>
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
103 !  (1) total column dry air mass (mut, call to calculate_full)
105 !  (2) total column dry air mass at u and v points 
106 !      (muu, muv, call to calculate_mu_uv)
108 !  (3) mass-coupled velocities for advection
109 !      (ru, rv, and rw, call to couple_momentum)
111 !  (4) omega (call to calc_ww_cp)
113 !  (5) moisture coefficients (cqu, cqv, cqw, call to calc_cq)
115 !  (6) inverse density (alt, call to calc_alt)
117 !  (7) geopotential at pressure points (php, call to calc_php)
119 !</DESCRIPTION>
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 )
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  )
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  )
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    )
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  )
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 )
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 )
163 END SUBROUTINE rk_step_prep
165 !-------------------------------------------------------------------------------
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)
194    IMPLICIT NONE
196    !  Input data.
198    TYPE(grid_config_rec_type)    ,           INTENT(IN   ) :: config_flags
200    INTEGER ,               INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
201                                             ims, ime, jms, jme, kms, kme, &
202                                             its, ite, jts, jte, kts, kte
204    LOGICAL ,               INTENT(IN   ) :: non_hydrostatic, top_lid
206    INTEGER ,               INTENT(IN   ) :: n_moist, rk_step
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
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
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
257    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(  OUT) :: mu_tend, &
258                                                                     mu_save
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
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
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
298    INTEGER, INTENT( IN ) :: damp_opt
300    REAL, INTENT( IN ) :: zdamp, dampcoef
302    REAL, INTENT( OUT ) :: max_horiz_cfl
303    REAL, INTENT( OUT ) :: max_vert_cfl
305    REAL    :: kdift, khdq, kvdq, cfn, cfn1, cf1, cf2, cf3
306    INTEGER :: i,j,k
307    INTEGER :: time_step
309 !<DESCRIPTION>
311 !  rk_tendency computes the large-timestep tendency terms in the 
312 !  momentum, thermodynamic (theta), and geopotential equations.  
313 !  These terms include:
315 !  (1) advection (for u, v, w, theta - calls to advect_u, advect_v,
316 !                 advect_w, and advact_scalar).
318 !  (2) geopotential equation terms (advection and "gw" - call to rhs_ph).
320 !  (3) buoyancy term in vertical momentum equation (call to pg_buoy_w).
322 !  (4) Coriolis and curvature terms in u,v,w momentum equations
323 !      (calls to subroutines coriolis, curvature)
325 !  (5) 3D diffusion on coordinate surfaces.
327 !</DESCRIPTION>
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 )
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 )
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 )
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 )
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 )
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 )
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 )
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 )
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 )
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 )
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 )
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 )
389    !  advection tendencies
390    CALL nl_get_time_step ( 1, time_step )
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 )
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 )
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 )
420 !  theta flux divergence
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 ) 
431      IF ( config_flags%cu_physics == GDSCHEME  .OR.     &
432           config_flags%cu_physics == G3SCHEME ) THEN
434      ! theta advection only:
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     )
441      END IF
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      )
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   )
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
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   )
484      IF(config_flags%pert_coriolis) THEN
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 )
508      END IF
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   )
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
529 !**************************************************************
531 !  Next, the terms that we integrate only with forward-in-time
532 !  (evaluate with time t variables).
534 !**************************************************************
536   forward_step: IF( rk_step == 1 ) THEN
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   )
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   )
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   )
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    )
574         pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
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   )
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   )
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        )
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         )
604         ENDIF pbl_test
606    !  Theta tendency computations.
608     END IF diff_opt1
610     IF ( diff_6th_opt .NE. 0 ) THEN
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 )
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 )
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 )
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 )
641     ENDIF
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   )
654   END IF forward_step
656 END SUBROUTINE rk_tendency
658 !-------------------------------------------------------------------------------
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                       )
671    IMPLICIT NONE
673    !  Input data.
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
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
692    REAL , DIMENSION( ims:ime , jms:jme  ) , INTENT(INOUT) :: mu_tend, &
693                                                              mu_tendf
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
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
712 ! Local
713    INTEGER :: i, j, k
715 !<DESCRIPTION>
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.
724 !</DESCRIPTION>
726 !  Finally, add the forward-step tendency to the rk_tendency
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
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
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
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
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
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
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
792 END SUBROUTINE rk_addtend_dry
794 !-------------------------------------------------------------------------------
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    )
814    IMPLICIT NONE
816    !  Input data.
818    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
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
825    LOGICAL , INTENT(IN   ) :: moist_step
827    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ),                &
828                                          INTENT(IN   )  :: scalar, scalar_old
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
835    REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(OUT  ) :: RQVFTEN
837    REAL, DIMENSION(ims:ime, kms:kme, jms:jme  ), INTENT(IN   ) ::     ru,  &
838                                                                       rv,  &
839                                                                       ww,  &
840                                                                       xkmhd,  &
841                                                                       alt
844    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fnm,  &
845                                                                   fnp,  &
846                                                                   rdn,  &
847                                                                   rdnw, &
848                                                                   base
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
861    REAL ,                                        INTENT(IN   ) :: rdx,     &
862                                                                   rdy,     &
863                                                                   khdif,   &
864                                                                   kvdif
866    INTEGER, INTENT( IN ) :: diff_6th_opt
867    REAL,    INTENT( IN ) :: diff_6th_factor
869    REAL ,                                        INTENT(IN   ) :: dt
871    INTEGER, INTENT(IN   ) :: adv_opt          
873    ! Local data
874   
875    INTEGER :: im, i,j,k
876    INTEGER :: time_step
878    REAL    :: khdq, kvdq, tendency
880 !<DESCRIPTION>
882 ! rk_scalar_tend calls routines that computes scalar tendency from advection 
883 ! and 3D mixing (TKE or fixed eddy viscosities).
885 !</DESCRIPTION>
888    khdq = khdif/prandtl
889    kvdq = kvdif/prandtl
891    scalar_loop : DO im = scs, sce
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 )
898      CALL nl_get_time_step ( 1, time_step )
900       IF( (rk_step == 3) .and. (adv_opt == POSITIVEDEF) ) THEN
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     )
914       ELSE IF( (rk_step == 3) .and. (adv_opt == MONOTONIC) ) THEN
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     )
928       ELSE
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  )
942       END IF
944      IF((config_flags%cu_physics == GDSCHEME .OR. config_flags%cu_physics == G3SCHEME) & 
945                      .and. moist_step .and. ( im == P_QV) ) THEN
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
953      rk_step_1: IF( rk_step == 1 ) THEN
955        diff_opt1 : IF (config_flags%diff_opt .eq. 1) THEN
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      )
967        pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
969          IF( (moist_step) .and. ( im == P_QV)) THEN
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 )
979          ELSE 
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 )
989          END IF
991       ENDIF pbl_test
993     ENDIF diff_opt1
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 )
1004   ENDIF rk_step_1
1006  END DO scalar_loop
1008 END SUBROUTINE rk_scalar_tend
1010 !-------------------------------------------------------------------------------
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  )
1022    IMPLICIT NONE
1024    !  Input data.
1026    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
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
1033    REAL,                    INTENT(IN   ) :: dt
1035    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
1036          INTENT(INOUT)                                  :: scalar_1,    &
1037                                                            scalar_2
1039    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
1040          INTENT(IN)                                     :: sc_tend
1042    REAL, DIMENSION(ims:ime, kms:kme, jms:jme ),                &
1043          INTENT(IN)                                  :: advect_tend
1045    REAL, DIMENSION(ims:ime, jms:jme  ), INTENT(IN   ) ::  mu_old,  &
1046                                                           mu_new,  &
1047                                                           mu_base, &
1048                                                           msftx,   &
1049                                                           msfty
1051    INTEGER :: i,j,k,im
1052    REAL    :: sc_middle, msfsq
1053    REAL, DIMENSION(its:ite) :: muold, r_munew
1055    REAL, DIMENSION(its:ite, kts:kte, jts:jte  ) :: tendency
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
1060 !<DESCRIPTION>
1062 !  rk_scalar_update advances the scalar equation given the time t value
1063 !  of the scalar and the scalar tendency.  
1065 !</DESCRIPTION>
1069 !  set loop limits.
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
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
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
1094     IF ( rk_step == 1 ) THEN
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
1099       DO  im = scs,sce
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)
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
1133       DO  k = kts, min(kte,kde-1)
1134       DO  i = its, min(ite,ide-1)
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)
1140       ENDDO
1141       ENDDO
1142       ENDDO
1144       ENDDO
1146     ELSE
1148       !  just compute new values, scalar_1 already at time t.
1150       DO  im = scs, sce
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
1177       DO  j = jts, min(jte,jde-1)
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
1184       DO  k = kts, min(kte,kde-1)
1185       DO  i = its, min(ite,ide-1)
1187         scalar_2(i,k,j,im) = (muold(i)*scalar_1(i,k,j,im)   &
1188                              + dt*tendency(i,k,j))*r_munew(i)
1190       ENDDO
1191       ENDDO
1192       ENDDO
1194       ENDDO
1196     END IF
1198 END SUBROUTINE rk_update_scalar
1200 !-------------------------------------------------------------------------------
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  )
1211    IMPLICIT NONE
1213    !  Input data.
1215    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
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
1222    REAL,                    INTENT(IN   ) :: dt
1224    REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce),                &
1225          INTENT(INOUT)                                  :: scalar,      &
1226                                                            sc_tend
1228    REAL, DIMENSION(ims:ime, jms:jme  ), INTENT(IN   ) ::  mu_old,  &
1229                                                           mu_new,  &
1230                                                           mu_base
1232    INTEGER :: i,j,k,im
1233    REAL    :: sc_middle, msfsq
1234    REAL, DIMENSION(its:ite) :: muold, r_munew
1236    REAL, DIMENSION(its:ite, kts:kte, jts:jte  ) :: tendency
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
1241 !<DESCRIPTION>
1243 !  rk_scalar_update advances the scalar equation given the time t value
1244 !  of the scalar and the scalar tendency.  
1246 !</DESCRIPTION>
1250 !  set loop limits.
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
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
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
1275       DO  im = scs, sce
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
1294       DO  j = jts, min(jte,jde-1)
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
1301       DO  k = kts, min(kte,kde-1)
1302       DO  i = its, min(ite,ide-1)
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
1310       ENDDO
1312 END SUBROUTINE rk_update_scalar_pd
1314 !------------------------------------------------------------
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 !-----------------------------------------------------------------------
1328    INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1329                                     ims, ime, jms, jme, kms, kme, &
1330                                     its, ite, jts, jte, kts, kte
1332    INTEGER ,       INTENT(IN   ) :: n_moist,n_chem,n_scalar,n_tracer,rk_step
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
1342    REAL , DIMENSION( ims:ime , jms:jme  ) , INTENT(INOUT) ::  mu_tendf
1344    REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_moist),INTENT(INOUT)::&
1345                                                           moist_tendf
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
1352    REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_scalar ),INTENT(INOUT)::&
1353                                                           scalar_tendf
1355 ! LOCAL VARS
1357    INTEGER :: im, ic, is
1359 !<DESCRIPTION>
1361 ! init_zero_tendency 
1362 ! sets tendency arrays to zero for all prognostic variables.
1364 !</DESCRIPTION>
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     )
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     )
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     )
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     )
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     )
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     )
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     )
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
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
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
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
1434 END SUBROUTINE init_zero_tendency
1436 !===================================================================
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
1449 integer :: is,ie,js,je,ks,ke
1451 !<DESCRIPTION
1453 ! quick and dirty debug io utility
1455 !</DESCRIPTION
1457 is = ids
1458 ie = ide-1
1459 js = jds
1460 je = jde-1
1461 ks = kds
1462 ke = kde-1
1464 if(field == 'u') ie = ide
1465 if(field == 'v') je = jde
1466 if(field == 'w') ke = kde
1468 write(io_unit) is,ie,ks,ke,js,je
1469 write(io_unit) a(is:ie, ks:ke, js:je)
1471 end subroutine dump_data
1473 !-----------------------------------------------------------------------
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
1489       TYPE(grid_config_rec_type), INTENT(IN)     ::      config_flags
1491       INTEGER,  INTENT(IN   )   ::          ids,ide, jds,jde, kds,kde, &
1492                                             ims,ime, jms,jme, kms,kme, &
1493                                             its,ite, jts,jte, kts,kte
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
1506       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme ),                &
1507                 INTENT(INOUT)   ::                           RTHRATEN
1509 ! cumulus
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
1521       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )               , &
1522                 INTENT(INOUT)   ::                            RUBLTEN, &
1523                                                               RVBLTEN, &
1524                                                              RTHBLTEN, &
1525                                                              RQVBLTEN, &
1526                                                              RQCBLTEN, &
1527                                                              RQIBLTEN
1529 ! fdda
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
1539       INTEGER :: i,k,j
1540       INTEGER :: itf,ktf,jtf,itsu,jtsv
1542 !-----------------------------------------------------------------------
1544 !<DESCRIPTION>
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.
1550 !</DESCRIPTION>
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)
1558 ! radiation
1560    IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
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
1570    ENDIF
1572 ! cumulus
1574    IF (config_flags%cu_physics .gt. 0) THEN
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
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
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
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
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
1625    ENDIF
1627 ! pbl
1629    IF (config_flags%bl_pbl_physics .gt. 0) THEN
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
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
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
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
1671     ENDIF
1673 ! fdda
1674 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
1675 !   so only couple those
1677    IF (config_flags%grid_fdda .gt. 0) THEN
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
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
1724     ENDIF
1726 END SUBROUTINE calculate_phy_tend
1728 !-----------------------------------------------------------------------
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  )
1735   IMPLICIT NONE
1737   INTEGER,  INTENT(IN   )   ::          ids,ide, jds,jde, kds,kde, &
1738                                         ims,ime, jms,jme, kms,kme, &
1739                                         its,ite, jts,jte, kts,kte
1741   REAL, DIMENSION( ims:ime , kms:kme , jms:jme  ), INTENT(INOUT) :: a
1743   INTEGER :: i,k,j
1745 !<DESCRIPTION>
1747 ! debug and testing code for bounding a variable
1749 !</DESCRIPTION>
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
1760   END SUBROUTINE positive_definite_filter
1762 !-----------------------------------------------------------------------
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  )
1769   IMPLICIT NONE
1771   INTEGER,  INTENT(IN   )   ::          ids,ide, jds,jde, kds,kde, &
1772                                         ims,ime, jms,jme, kms,kme, &
1773                                         its,ite, jts,jte, kts,kte
1775   REAL, DIMENSION( ims:ime , kms:kme , jms:jme  ), INTENT(INOUT) :: tke
1776   REAL, INTENT(   IN) :: tke_upper_bound
1778   INTEGER :: i,k,j
1780 !<DESCRIPTION>
1782 ! bounds tke between zero and tke_upper_bound.
1784 !</DESCRIPTION>
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
1794   END SUBROUTINE bound_tke
1798 END MODULE module_em