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