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