source: trunk/WRF.COMMON/WRFV3/dyn_em/module_em.F @ 2761

Last change on this file since 2761 was 2761, checked in by aslmd, 2 years ago

Applied planetary adaptation changes to WRFV3. job done previously by LMD_LES_MARS_install. Moved Registry.EM.newphys to put it simply in Registry.EM

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